新人的第一个物理网格NACA0012

rockandtcl

从大二开始看《计算流体力学入门》（清华版的）。结果老是看公式头都疼死了。这次总算有时间静下心来做个。呵呵。有点成就感。
这一周来我一直在学习做网格。首先是物理网格，最简单的就是几何法生成网格。
在翼型的生成上使用的是代数方法。通过一个函数来生成NACA0012。
呵呵，有点成就感。不过开始的时候出了点问题，因为是老书《计算流体动力学》。符号写错了。害我调试程序花了2天时间。今天总算改了过来，得到了师兄师姐姐的帮助（主要是TECPLOT的一些格式说明上）勉强把网格做出来了。
问题又来了
前缘总是往里凹。还望师兄师姐们告诉我改怎么解决。
主要思路是：
1：首先用函数法生成NACA0012
2、然后从后缘开始，第一个theta取0，后一个点的坐标与前一个点的坐标做运算求出两条直线的夹角，把夹角与前一个theta相加。同时算出直线长度R。
3、A=SQRT（R）*COS（THETA/2）；
S=SQRT（R）*SIN（THETA/2）
4、X=A*A-（S+F*）*（S+F*）；Y=2*A*（S+F*）；
F1〈F2〈F3〈F4.....
5、输出X，Y；
问题就出在了A，在接近前缘的时候R的数值很小，而THETA的变化量很大，但由于接近前缘的时候THETA是接近3.141592（弧度），所以微小的误差对于COS（THETA）来说是致命的，而且R由于很小，即使绝对误差也很小，但相对误差就大了。于是就出现了前缘的向里凹的现象。不知道如何处理才好。希望各位师兄师姐能提供个解决的办法。[我是分开计算了，所以如果合上的话可能就没有误差了]。
接下来要把物理转成计算域。呵呵
新手需要的话请留言，可以提供代码
C语言版本、FORTRAN版本（给我个可以用的序列号。呵呵）

无心向学

要一个fortran版本的
先谢过了
哈哈
china95189@yahoo.com.cn