采用椭圆型方程数值解法生成 NACA0015 翼型 O 型（C 型）网格，并采用速度势方程求解绕 NACA0015 翼型迎角为零时的低速不可压无黏流场。

　　1.1.1 求解思路

　　首先通过 Laplace 方程生成翼型O型网格，然后通过速度势方程求解流场。

　　1.2.1 控制方程

　　对称翼型 NACA0015 翼型的表面(包括上下表面点)坐标计算公式：

　　网格生成的控制方程为 Laplace 方程，只考虑源项为0时的情况，其表达式为：

　　可列出计算平面变量 与实际物理平面x,y之间的关系：

　　其中 系数均由反变换系数组成：

　　1.2.2 求解原理

　　（1） 设定远场边界条件

　　对于 O 型网格，远场一般设定为半径是一定倍数弦长的圆，例如可以设定以翼型中心点为圆心，以 15 倍弦长为半径的圆，作为远场边界。根据设置的离散点数，确定远场网格点的(x，y)坐标值。

　　（2） 设定物面边界条件

　　确定翼型表面点分布的原则是保证流动变化剧烈的区域分布密集，变化平缓的区域分布稀疏。翼型表面点的x坐标通过下列方程给定：

　　（3） 计算初值

　　根据给定的边界值插值出每一个网格点的初始值，

　　（4） 迭代实现

　　根据差分方程选择合适的迭代格式进行迭代求解，直至达到收敛标准，即所有网格点上的相邻两次迭代的结果差异均小于一个极小量 。遍历所有x,y，直至满足以下条件：

　　1.2.3 程序实现

　　（1） 参数输入

　　采用将程序的输入参数写入文件的方式，将输入文件命名为input.txt。输入参数包括翼型表面与远场的离散点个数N、远场边界条件的半径r、*大迭代次数iterations以及收敛标准。取切向点数101，法向点数也是101，远场半径30倍弦长（弦长为1），*大迭代次数为次，收敛标准为 1e-10。程序开始后，先读取输入文件中的参数。

　　（2） 网格初始化

　　在翼型周围生成O型网格，割线取为沿x轴正方向的直线，原点为翼型前缘点。翼型弦长为1，中心与远场的圆心重合（x=0）。通过边界条件给定翼型与远场的表面坐标，初始化后远场边界上的网格点均匀分布，物面边界翼型表面上网格点在前缘以及后缘处进行了适当加密，符合CFD在物体表面加密网格的思想。

　　在割缝上并没有给出固定的坐标，这是因为在割缝上使用了周期性边界条件，并且迭代时割缝上的网格点也参与了迭代，其初始值并不影响迭代的结果以及收敛性。查找资料可知，割缝上的点参与迭代后可以使网格的正交性更好。

　　（3） 网格迭代生成

　　松弛技术是一种适用于椭圆型差分方程的求解技术。通过分析题目，已知控制方程为椭圆型偏微分方程，因此可以运用松弛迭代方法来进行流场求解。首先写出中心差分后的差分方程：

　　采用简单迭代：

　　（4） 网格结果导出

　　将*终迭代的网格（x,y）坐标输出到.dat文件中，便于用Tecplot软件进行后处理，网格结果文件命名为O-grid.dat。

　　需要注意的是Tecplot软件的文件输入格式，导入的数据要严格按照点形式（POINT）的格式输入，即要在.dat文件的开始加上标头：variables=x,y；zone i=101,j=101,F=POINT其中i和j是网格离散点的坐标，i,j的数量代表x,y方向上各有101个点，这步操作在其他的程序输出文件时也要用到。

　　1.2.4 网格生成结果与分析

　　（1） 结果

　　文件成功导出后，用Tecplot打开.dat文件进行后处理，将网格显示出来。图1是翼型的全局网格：

　　图1翼型的全局网格

　　图2是翼型表面的网格，图3、图4是进一步放大后翼型前缘与后缘的网格：

　　（2） 分析

　　放大到翼型表面后观察，可以看到翼型前、后缘的网格满足了网格局部加密的要求。*终翼型表面与前缘、后缘的网格正交性不错，实现了理想的效果。

　　1.3.1 控制方程

　　流场求解的控制方程仍然是为 Laplace 方程，源项设定为0，具体形式为：

　　代表的是空气动力学中的速度势函数，可求出计算平面内的变换形式：

　　其中系数 的形式如下：

　　1.3.2 求解原理

　　（1） 设定远场边界条件

　　远场边界条件是自由来流条件，即

　　（2） 设定物面边界条件

　　物面边界条件是无穿透条件，即

　　进一步可推导出如下形式：

　　计算平面变换后的控制方程（计算流体力学基础及其应用中式（5-17））如下式所示：

　　（3） 速度势迭代求解

　　由于速度势求解的控制方程也是拉普拉斯方程，中心差分后的差分方程与简单迭代的格式与1.2.3节第（3）小节的格式一致，这里就不赘述了。

　　（4） 流场参数求解

　　求解出速度势函数后，根据以下公式便可得到流场中的速度：

　　求解速度分布需要用雅可比行列式进行逆变换，下面给出雅可比行列式的表示形式（书中式（5-22b））：

　　速度求解完成后，*后求解流场的压强系数分布，计算式如下所示：

　　1.3.3 程序实现

　　（1） 程序输入

　　在输入文件中添加新的参数：自由来流速度，流场求解时其他参数与生成网格的参数保持一致。设定攻角为0度，无穷远处来流速度为20m/s。在网格的迭代生成函数迭代完后，依次调用流场计算的子函数，先将网格（x,y）坐标、收敛标准等参数导入迭代求解速度势的子函数。

　　（2） 迭代求解速度势

　　调用迭代求解速度势的子函数，运用逐次超松弛迭代方法迭代求解速度势函数，*终达到收敛的标准：

　　仍取1e-10。

　　（3） 求解流场参数

　　调用求解流场参数的函数，输入为速度势函数，用*终收敛的速度势分布求出流场的速度分布与压力系数分布。

　　（4） 流场结果导出

　　先在.dat文件的前几行加上标头，将流场参数：u，v，Cp与（x,y）坐标输出到.dat文件中，用Tecplot进行后处理，得到速度分布云图与压力系数的云图。

　　将离散点的x坐标与对应的压力系数输出到.dat文件中，后处理后得到x方向上的压力系数分布曲线。

　　1.3.4 流场结果与分析

　　（1） 结果

　　下面是流场的输出结果（翼型的弦长为1），图5是翼型的压力云图，图6是翼型的速度云图：

　　图6是x方向上的翼型的压力系数分布曲线：

　　用Python实现，主函数如下：

　　如需代码可私