采用橢圓型方程數(shù)值解法生成 NACA0015 翼型 O 型（C 型）網(wǎng)格，并采用速度勢方程求解繞 NACA0015 翼型迎角為零時的低速不可壓無黏流場。
1.1.1 求解思路
首先通過 Laplace 方程生成翼型O型網(wǎng)格，然后通過速度勢方程求解流場。
1.2.1 控制方程
對稱翼型 NACA0015 翼型的表面(包括上下表面點)坐標計算公式：
網(wǎng)格生成的控制方程為 Laplace 方程，只考慮源項為0時的情況，其表達式為：
可列出計算平面變量 與實際物理平面x,y之間的關系：
其中 系數(shù)均由反變換系數(shù)組成：
1.2.2 求解原理
（1） 設定遠場邊界條件
對于 O 型網(wǎng)格，遠場一般設定為半徑是一定倍數(shù)弦長的圓，例如可以設定以翼型中心點為圓心，以 15 倍弦長為半徑的圓，作為遠場邊界。根據(jù)設置的離散點數(shù)，確定遠場網(wǎng)格點的(x，y)坐標值。
（2） 設定物面邊界條件
確定翼型表面點分布的原則是保證流動變化劇烈的區(qū)域分布密集，變化平緩的區(qū)域分布稀疏。翼型表面點的x坐標通過下列方程給定：
（3） 計算初值
根據(jù)給定的邊界值插值出每一個網(wǎng)格點的初始值，
（4） 迭代實現(xiàn)
根據(jù)差分方程選擇合適的迭代格式進行迭代求解，直至達到收斂標準，即所有網(wǎng)格點上的相鄰兩次迭代的結果差異均小于一個極小量 。遍歷所有x,y，直至滿足以下條件：
1.2.3 程序實現(xiàn)
（1） 參數(shù)輸入
采用將程序的輸入?yún)?shù)寫入文件的方式，將輸入文件命名為input.txt。輸入?yún)?shù)包括翼型表面與遠場的離散點個數(shù)N、遠場邊界條件的半徑r、*大迭代次數(shù)iterations以及收斂標準。取切向點數(shù)101，法向點數(shù)也是101，遠場半徑30倍弦長（弦長為1），*大迭代次數(shù)為次，收斂標準為 1e-10。程序開始后，先讀取輸入文件中的參數(shù)。
（2） 網(wǎng)格初始化
在翼型周圍生成O型網(wǎng)格，割線取為沿x軸正方向的直線，原點為翼型前緣點。翼型弦長為1，中心與遠場的圓心重合（x=0）。通過邊界條件給定翼型與遠場的表面坐標，初始化后遠場邊界上的網(wǎng)格點均勻分布，物面邊界翼型表面上網(wǎng)格點在前緣以及后緣處進行了適當加密，符合CFD在物體表面加密網(wǎng)格的思想。
在割縫上并沒有給出固定的坐標，這是因為在割縫上使用了周期性邊界條件，并且迭代時割縫上的網(wǎng)格點也參與了迭代，其初始值并不影響迭代的結果以及收斂性。查找資料可知，割縫上的點參與迭代后可以使網(wǎng)格的正交性更好。
（3） 網(wǎng)格迭代生成
松弛技術是一種適用于橢圓型差分方程的求解技術。通過分析題目，已知控制方程為橢圓型偏微分方程，因此可以運用松弛迭代方法來進行流場求解。首先寫出中心差分后的差分方程：
采用簡單迭代：
（4） 網(wǎng)格結果導出
將*終迭代的網(wǎng)格（x,y）坐標輸出到.dat文件中，便于用Tecplot軟件進行后處理，網(wǎng)格結果文件命名為O-grid.dat。
需要注意的是Tecplot軟件的文件輸入格式，導入的數(shù)據(jù)要嚴格按照點形式（POINT）的格式輸入，即要在.dat文件的開始加上標頭：variables=x,y；zone i=101,j=101,F=POINT其中i和j是網(wǎng)格離散點的坐標，i,j的數(shù)量代表x,y方向上各有101個點，這步操作在其他的程序輸出文件時也要用到。
1.2.4 網(wǎng)格生成結果與分析
（1） 結果
文件成功導出后，用Tecplot打開.dat文件進行后處理，將網(wǎng)格顯示出來。圖1是翼型的全局網(wǎng)格：
圖1翼型的全局網(wǎng)格
圖2是翼型表面的網(wǎng)格，圖3、圖4是進一步放大后翼型前緣與后緣的網(wǎng)格：
（2） 分析
放大到翼型表面后觀察，可以看到翼型前、后緣的網(wǎng)格滿足了網(wǎng)格局部加密的要求。*終翼型表面與前緣、后緣的網(wǎng)格正交性不錯，實現(xiàn)了理想的效果。
1.3.1 控制方程
流場求解的控制方程仍然是為 Laplace 方程，源項設定為0，具體形式為：
代表的是空氣動力學中的速度勢函數(shù)，可求出計算平面內(nèi)的變換形式：
其中系數(shù) 的形式如下：
1.3.2 求解原理
（1） 設定遠場邊界條件
遠場邊界條件是自由來流條件，即
（2） 設定物面邊界條件
物面邊界條件是無穿透條件，即
進一步可推導出如下形式：
計算平面變換后的控制方程（計算流體力學基礎及其應用中式（5-17））如下式所示：
（3） 速度勢迭代求解
由于速度勢求解的控制方程也是拉普拉斯方程，中心差分后的差分方程與簡單迭代的格式與1.2.3節(jié)第（3）小節(jié)的格式一致，這里就不贅述了。
（4） 流場參數(shù)求解
求解出速度勢函數(shù)后，根據(jù)以下公式便可得到流場中的速度：
求解速度分布需要用雅可比行列式進行逆變換，下面給出雅可比行列式的表示形式（書中式（5-22b））：
速度求解完成后，*后求解流場的壓強系數(shù)分布，計算式如下所示：
1.3.3 程序實現(xiàn)
（1） 程序輸入
在輸入文件中添加新的參數(shù)：自由來流速度，流場求解時其他參數(shù)與生成網(wǎng)格的參數(shù)保持一致。設定攻角為0度，無窮遠處來流速度為20m/s。在網(wǎng)格的迭代生成函數(shù)迭代完后，依次調用流場計算的子函數(shù)，先將網(wǎng)格（x,y）坐標、收斂標準等參數(shù)導入迭代求解速度勢的子函數(shù)。
（2） 迭代求解速度勢
調用迭代求解速度勢的子函數(shù)，運用逐次超松弛迭代方法迭代求解速度勢函數(shù)，*終達到收斂的標準：
仍取1e-10。
（3） 求解流場參數(shù)
調用求解流場參數(shù)的函數(shù)，輸入為速度勢函數(shù)，用*終收斂的速度勢分布求出流場的速度分布與壓力系數(shù)分布。
（4） 流場結果導出
先在.dat文件的前幾行加上標頭，將流場參數(shù)：u，v，Cp與（x,y）坐標輸出到.dat文件中，用Tecplot進行后處理，得到速度分布云圖與壓力系數(shù)的云圖。
將離散點的x坐標與對應的壓力系數(shù)輸出到.dat文件中，后處理后得到x方向上的壓力系數(shù)分布曲線。
1.3.4 流場結果與分析
（1） 結果
下面是流場的輸出結果（翼型的弦長為1），圖5是翼型的壓力云圖，圖6是翼型的速度云圖：
圖6是x方向上的翼型的壓力系數(shù)分布曲線：
用Python實現(xiàn)，主函數(shù)如下：
如需代碼可私
