面向粘性流動問題的無網(wǎng)格／直角網(wǎng)格混合算法研究

蒲賽虎，陳紅全

（南京航空航天大學(xué) 航空宇航學(xué)院，江蘇 南京 210016）

0 引 言

早期數(shù)值方法采用的直角網(wǎng)格具有簡單、規(guī)則的特點，良好的正交性避免了單元扭曲變形，因此被人們廣泛使用。但處理復(fù)雜外形時，由于網(wǎng)格線與物面通常并不一致，在處理邊界條件上存在較大困難。因此人們發(fā)展了貼體網(wǎng)格方法：將物理域內(nèi)的曲線貼體網(wǎng)格映射成計算域內(nèi)規(guī)則的網(wǎng)格，再進(jìn)行求解。但是計算域內(nèi)的微分方程往往比物理域內(nèi)原形式更加復(fù)雜，而且曲線網(wǎng)格不易控制的正交性對解有較大的影響。目前，許多學(xué)者都在對直角網(wǎng)格方法［1］和貼體網(wǎng)格方法［2］進(jìn)行研究和改進(jìn)。

無網(wǎng)格方法是近年來興起的一種新型數(shù)值計算方法，在剖分計算區(qū)域時，該方法只需要進(jìn)行布點填充，而無需生成網(wǎng)格單元，因此具有靈活性，易于處理復(fù)雜外形。但在計算效率方面，無網(wǎng)格方法與網(wǎng)格方法相比，目前還不具有競爭性［3］。基于此，馬志華等人［4］提出了一種計算區(qū)域整體采用直角網(wǎng)格，只在物面附近嵌入局部無網(wǎng)格處理的無網(wǎng)格／直角網(wǎng)格混合算法，力求無網(wǎng)格方法和直角網(wǎng)格方法的優(yōu)勢互補(bǔ)。目前該混合算法已被應(yīng)用于求解Euler方程，成功模擬了鼓包、圓柱及翼型等的無粘流動［4－5］，同時文獻(xiàn)［4］展示了混合算法相對于傳統(tǒng)網(wǎng)格方法在計算效率上有競爭性。

本文將無網(wǎng)格／直角網(wǎng)格混合算法發(fā)展用于求解NS方程，從而對粘性流動進(jìn)行數(shù)值模擬。為了處理粘性流動問題，首先必須解決計算區(qū)域的無網(wǎng)格／直角網(wǎng)格的離散問題。對于整體直角網(wǎng)格，本文改用直角坐標(biāo)與四叉樹加密的方法生成非結(jié)構(gòu)直角網(wǎng)格，力求通過四叉樹加密等措施使網(wǎng)格區(qū)與無網(wǎng)格區(qū)過渡可控。在無網(wǎng)格區(qū)，本文引入了點距控制函數(shù)及無網(wǎng)格區(qū)之間的布點干涉機(jī)制等措施（見第2節(jié)），使得生成的點云結(jié)構(gòu)能反映粘性流動邊界層特點，同時適合處理多個物體相互干擾的復(fù)雜流動問題。基于計算區(qū)域的無網(wǎng)格直角網(wǎng)格離散，給出了求解NS方程耦合SA湍流模型的無網(wǎng)格／直角網(wǎng)格混合算法的具體實施方法（見第3節(jié)），編程計算了典型翼型的粘性繞流，同時成功模擬了兩段翼型粘性繞流，展示了所發(fā)展的算法處理存在多個物體的粘性流動問題的可行性。

1 控制方程

本文采用的控制方程是直角坐標(biāo)系下的可壓縮雷諾平均NS方程，其微分形式可以寫為：

式中W為守恒矢變量，E和F為無粘矢通量，Ev和Fv為粘性矢通量，表達(dá)式如下：

Ev和Fv中的粘性應(yīng)力項表達(dá)式為：

其中ρ為密度，u、v分別是沿x、y方向的速度分量，p為壓強(qiáng)，E是單位質(zhì)量內(nèi)的總能。其它符號的含義可以參見文獻(xiàn)［6］。本文層流粘性系數(shù)μL利用Surtherland公式計算得到，湍流粘性系數(shù)μT通過SA湍流模型［7］得到。

2 計算區(qū)域的無網(wǎng)格／直角網(wǎng)格離散問題

本文采用的混合算法由無網(wǎng)格方法和直角網(wǎng)格方法構(gòu)成，其中無網(wǎng)格方法的基本計算單元稱為點云。采用無網(wǎng)格方法，需要對計算區(qū)域進(jìn)行布點離散，再對每一點構(gòu)成合適的點云。圖1（a）給出了一個七點無網(wǎng)格點云。其中i點是中心點，點1～點6是其衛(wèi)星點。以函數(shù)f為例，在i點處，其空間一階偏導(dǎo)數(shù)可近似為：

其中m表示點云中的衛(wèi)星點數(shù)。fik表示i點和其第k個衛(wèi)星點之間中點處的函數(shù)值。αik、βik是只與點的坐標(biāo)有關(guān)的系數(shù)，在本文中，這些系數(shù)是在點云上用加權(quán)的二次極小曲面逼近如下線性方程得到：

本文沿用文獻(xiàn)［8］的方法，對二次極小曲面逼近進(jìn)行了加權(quán)，權(quán)系數(shù)形式為wik＝（rref／rik）2。其中rref為i點點云的參考半徑，定義為i點到其最遠(yuǎn)衛(wèi)星點的距離；rik為i點與其第k個衛(wèi)星點的距離。需要指出的是，對于本文粘性流動計算所涉及到的各向異性點云（圖1b中所示），切向點距與法向點距之比（拉伸比）可能非常大，甚至達(dá)到100倍以上，若不采用加權(quán)（權(quán)系數(shù)wik＝1），則在采用二次極小曲面逼近式（5）時系數(shù)矩陣會出現(xiàn)病態(tài)（比如當(dāng)拉伸比為100時，系數(shù)矩陣的條件數(shù)為10000），而本文在最小二乘方法中加入上述權(quán)系數(shù)則避免了系數(shù)矩陣病態(tài)的問題，比如采用權(quán)系數(shù)后，當(dāng)拉伸比為100時，系數(shù)矩陣的條件數(shù)約為4，表明系數(shù)矩陣是健康的。

圖1 無網(wǎng)格點云Fig.1 Clouds of points for gridless method

在直角網(wǎng)格上，同樣以函數(shù)f為例，在i點處，其空間一階偏導(dǎo)數(shù)可以采用如下的差分格式逼近：

其中，hi是直角網(wǎng)格在i點處的網(wǎng)格歩長（本文生成的直角網(wǎng)格每一點上x方向和y方向的網(wǎng)格歩長相同），fiE、fiW、fiS、fiN分別是i點和其東、西、南、北方向上鄰近點之間中點處的函數(shù)值。

2.1 局部無網(wǎng)格布點和整體直角網(wǎng)格生成

采用無網(wǎng)格／直角網(wǎng)格混合算法，計算區(qū)域整體采用直角網(wǎng)格，只在物面附近采用局部無網(wǎng)格處理。求解無粘流動問題時，物面附近的局部無網(wǎng)格通常采用圖1（a）所示的各向同性點云［3，9］，但對本文涉及到的粘性流動，特別是高雷諾數(shù)粘性流動而言，為了有效模擬邊界層內(nèi)較大的法向速度梯度，沿物面法線方向布點需較密，此時若仍采用各向同性點云，勢必造成整體布點量太大，不利于減少計算時間?；诖?，本文在物面附近采用了圖1（b）所示的各向異性點云。采用各向異性點云，物面法線方向布點較密，從而能夠模擬邊界層內(nèi)的法向速度梯度，而物面切線方向點距相對較大，有利于減少布點總量，節(jié)約計算時間。對單個物體，本文的布點步驟是：首先將物面邊界劃分成離散點的形式，然后沿物面法線方向按層推進(jìn)的方法進(jìn)行布點［9］。在層推進(jìn)過程中，我們引入了指數(shù)函數(shù)來控制每一層點的推進(jìn)距離，該函數(shù)為：

式中，di為第i層點與第i－1層點之間的距離；d1是第一層點與物面的距離，可根據(jù)計算的雷諾數(shù)由經(jīng)驗確定；r是一個調(diào)節(jié)距離增長速度的參數(shù)，本文取為0.1。對于流場中存在多個物體相互干擾的情況，本文在對每個物體單獨采用上述方法進(jìn)行布點的同時，引入了無網(wǎng)格區(qū)之間的布點干涉機(jī)制，以解決各個無網(wǎng)格區(qū)之間可能相互重疊的問題，即在對每個物體同時進(jìn)行層推進(jìn)布點過程中，若發(fā)現(xiàn)本無網(wǎng)格區(qū)某點周圍δ距離范圍內(nèi)存在其它無網(wǎng)格區(qū)的點，則在該點處，層推進(jìn)布點停止，而其它點上的層推進(jìn)布點操作繼續(xù)進(jìn)行。δ是一個距離閥值，可根據(jù)具體計算對象確定。

在形成的無網(wǎng)格區(qū)中，除外邊界點外，其它點的點云結(jié)構(gòu)可以在布點過程中根據(jù)層推進(jìn)結(jié)構(gòu)同時確定，邊界點點云結(jié)構(gòu)的確定方法將在下節(jié)中介紹。對流場中除無網(wǎng)格區(qū)之外的區(qū)域，本文則采用非結(jié)構(gòu)直角網(wǎng)格進(jìn)行快速網(wǎng)格劃分。在此過程中，我們采用了直角坐標(biāo)與四叉樹加密的方法生成非結(jié)構(gòu)直角網(wǎng)格［1］，采用該方法，網(wǎng)格區(qū)與無網(wǎng)格區(qū)之間過渡可控，為下一步確定網(wǎng)格區(qū)和無網(wǎng)格區(qū)交界處的點云結(jié)構(gòu)打下了基礎(chǔ)。圖2中給出了GA（W）－1兩段翼型周圍的無網(wǎng)格點和非結(jié)構(gòu)直角網(wǎng)格。

2.2 無網(wǎng)格區(qū)與網(wǎng)格區(qū)邊界處點云結(jié)構(gòu)的確定

采用無網(wǎng)格／直角網(wǎng)格混合算法，計算區(qū)域被分成了無網(wǎng)格區(qū)和網(wǎng)格區(qū)，這兩個區(qū)之間需要進(jìn)行流場信息交換。本文繼承了原混合算法采用對偶節(jié)點方法實現(xiàn)無網(wǎng)格區(qū)與網(wǎng)格區(qū)之間流場信息交換的思想。所謂對偶節(jié)點方法，是將靠近無網(wǎng)格區(qū)的部分直角網(wǎng)格點（如圖3中的方形點）作為無網(wǎng)格點處理，從而實現(xiàn)兩個區(qū)之間的流場信息交換，關(guān)于該方法的細(xì)節(jié)描述可以參見文獻(xiàn)［5］。采用對偶節(jié)點方法，需要確定直角網(wǎng)格區(qū)和無網(wǎng)格區(qū)之間邊界上的點（如圖3中的i點，j點）的點云結(jié)構(gòu)。另外，對于多物體情況，不同無網(wǎng)格區(qū)之間的邊界點（如圖3中的k點）的點云結(jié)構(gòu)也需要確定。確定某一點的點云結(jié)構(gòu)，就是選擇合適的點作為該點的衛(wèi)星點。為此，本文提出了一種基于距離和夾角的點云選點準(zhǔn)則，下面以圖3中i點的衛(wèi)星點選取過程為例，介紹該選點準(zhǔn)則的實施過程如下：

圖2 GA（W）－1翼型周圍的無網(wǎng)格點和直角網(wǎng)格Fig.2 Points and Cartesian grid around GA（W）－1airfoil

圖3 交界處點云結(jié)構(gòu)Fig.3 Structure of the clouds of points at the interface

1）確定一個以i點為圓心，r為半徑的圓形區(qū)域為選點區(qū)域，將該區(qū)域內(nèi)除i點本身外的點都作為其衛(wèi)星點的候選點，如圖中的點1～點11。本文中r的取值為：r＝k·dlocal，dlocal為i點當(dāng)?shù)氐木W(wǎng)格尺度，由網(wǎng)格生成過程確定。k為縮放系數(shù)，本文研究發(fā)現(xiàn)取1.5是一個合適的值。由于本文直角網(wǎng)格生成采用了四叉樹數(shù)據(jù)結(jié)構(gòu)，因此上述候選點可以根據(jù)既有的四叉樹數(shù)據(jù)結(jié)構(gòu)快速確定；

2）對于某一候選點a，分別判斷由點a，中心點i和其它任意候選點b之間形成的角度α，若α＜αmin（本文取αmin＝30°），則比較點a和點b到中心點i的距離大小，若，則保留點a為候選點，而將點b排除出候選點；反之若，則保留點b為候選點，而將點a排除出候選點。經(jīng)過上述判斷，圖3中的點4、點6、點7、點8就被排除出候選點范圍。此步驟的作用一方面是盡量選擇離中心點較近的點作為該點衛(wèi)星點，另一方面避免距離太近的點以及在同一直線上的點同時作為衛(wèi)星點，因為這兩種情況容易造成在采用最小二乘方法確定式（1）中的系數(shù)時，系數(shù)矩陣出現(xiàn)病態(tài)，影響數(shù)值求解精度［10］。

經(jīng)過步驟1）、2），最終剩下的候選點即作為i點的衛(wèi)星點。圖3中標(biāo)出了采用上述方法確定的i點的點云結(jié)構(gòu)，其衛(wèi)星點為點1、2、3、5、9、10、11共7個點。圖3同時標(biāo)出了采用上述選點準(zhǔn)則確定出來的k點的點云結(jié)構(gòu)?？梢姴捎蒙鲜鲞x點準(zhǔn)則，選出的衛(wèi)星點圍繞中心點的分布比較均勻，這有助于保證點云上空間導(dǎo)數(shù)的離散精度。

3 NS方程的離散及求解格式

對NS方程中的空間導(dǎo)數(shù)，本文在無網(wǎng)格區(qū)中采用式（4）進(jìn)行離散，而在直角網(wǎng)格上采用式（6）進(jìn)行離散。由于NS方程的無粘通量項與Euler方程的無粘通量項完全一樣，因此本文在離散NS方程中無粘通量項時，采用了與原混合算法離散Euler方程中無粘通量項相同的方法，如在無網(wǎng)格區(qū)采用式（4）進(jìn)行離散過程如下：

其中αik和βik是式（4）中確定的系數(shù)，Gik是i點和其第k個衛(wèi)星點中點處的通量。本文與原混合算法一樣，都采用Roe的近似Riemann解來確定Gik，具體過程可參見文獻(xiàn)［4］。對于NS方程中出現(xiàn)的粘性通量項，以其中的項為例，本文采用無網(wǎng)格方法和直角網(wǎng)格方法進(jìn)行離散的過程如下：

1）采用無網(wǎng)格方法離散過程：

而i點和其第k個衛(wèi)星點中點處的一階導(dǎo)數(shù)按如下確定［8］：

其中，Δx＝xk－xi，Δy＝y(tǒng)k－yi，Δs2＝Δx2＋Δy2

2）采用直角網(wǎng)格方法離散過程：

其中而i點與其東、西、南、北四個方向上的鄰近點中點處的一階導(dǎo)數(shù)也采用式（10）確定。

采用無網(wǎng)格方法和直角網(wǎng)格方法離散NS方程的空間導(dǎo)數(shù)后，可以得到其半離散形式：

為了得到流動定常解，對于NS方程的半離散形式（12），本文采用五步Runge－Kutta顯式時間推進(jìn)進(jìn)行處理：

其中，上標(biāo)n表示當(dāng)前時間層，n＋1表示新的時間層，α為歩長因子，Δt為當(dāng)?shù)貢r間歩長，其表達(dá)式為：

式中σ為CFL數(shù)，C為常數(shù)，本文取值為4。對于無網(wǎng)格方法，式中Λc和Λv的表達(dá)式如下：

其中，uik、vik、ρik、μL、μT、cik分別表示中心點i與其第k個衛(wèi)星點之間中點處的速度分量、密度、粘性系數(shù)以及音速。αik和βik為式（4）所確定的系數(shù)。在直角網(wǎng)格點上，當(dāng)?shù)貢r間歩長也可按式（14）計算得到。數(shù)值求解時還需補(bǔ)充邊界條件，本文物面邊界采用無滑移邊界條件，遠(yuǎn)場則采用無反射邊界條件［11］。

4 算例與分析

本文已采用發(fā)展的無網(wǎng)格／直角網(wǎng)格混合算法進(jìn)行了編程，并對典型粘性流動問題進(jìn)行了數(shù)值模擬。這里首先給出NACA0012翼型和RAE2822翼型的流動考核算例，在此基礎(chǔ)上發(fā)展用于GA（W）－1兩段翼型的流場數(shù)值模擬。

4.1 NACA0012翼型亞音速粘性繞流

由于有實驗數(shù)據(jù)可供對照，NACA0012翼型亞音速粘性繞流計算經(jīng)常用來考核算法。此算例的來流條件是Ma＝0.5，α＝3.51°，Re＝2.93×106。本文混合算法采用的無網(wǎng)格點和直角網(wǎng)格如圖4中所示，其中翼型表面布置了298個點，直角網(wǎng)格點和無網(wǎng)格點總數(shù)為20235個，第一層無網(wǎng)格點距物面的距離為翼型弦長的5.0×10－5倍。圖5中給出了計算得到的流場等密度線分布，可以看到等密度線在網(wǎng)格區(qū)和無網(wǎng)格區(qū)之間過渡光滑（圖中翼型周圍的黑線是無網(wǎng)格區(qū)外邊界），說明本文采用的兩個區(qū)之間的流場信息交換方法是可行的。圖6給出了翼型上表面30%弦長附近的速度型，該速度型反映出高雷諾數(shù)粘性流動邊界層內(nèi)速度梯度較大的特點。圖7給出了計算得到的翼型表面壓力系數(shù)分布，可以看到計算結(jié)果與實驗數(shù)據(jù)［12］吻合較好。另外我們還對迎角較大的情況（Ma＝0.301，α＝9.79°，Re＝2.6×106）進(jìn)行了計算，圖8給出了計算得到的翼型表面壓力系數(shù)分布，可以看到計算結(jié)果與實驗數(shù)據(jù)［12］也比較吻合。

圖4 NACA0012翼型周圍的無網(wǎng)格點和直角網(wǎng)格Fig.4 Points and Cartesian grid around NACA0012airfoil

圖5 流場等密度線Fig.5 Density contours of the flow field

圖6 上翼面速度型Fig.6 Velocity profiles near the upper surface

圖7 迎角為3.51°時翼型表面壓力系數(shù)分布Fig.7 Cpdistribution around NACA0012 with the angle of attack 3.51°

圖8 迎角為9.79°時翼型表面壓力系數(shù)分布Fig.8 Cpdistribution around NACA0012 with the angle of attack 9.79°

4.2 RAE2822翼型跨音速粘性繞流

上一個算例是亞音速流動問題，本算例是跨音速流動問題，用以考核所發(fā)展的混合算法對激波的捕捉能力。來流條件是Ma＝0.73，α＝2.79°，Re＝6.5×106。圖9給出了計算采用的局部無網(wǎng)格點和非結(jié)構(gòu)直角網(wǎng)格，其中翼型表面布置了310個點，直角網(wǎng)格點和無網(wǎng)格點總數(shù)為23866個，第一層無網(wǎng)格點距物面的距離為翼型弦長的1.0×10－5倍。圖10展示了計算得到的流場等密度線分布，可以看到在翼型上翼面50%弦長附近存在較強(qiáng)的激波，同時還可以觀察到激波與邊界層相互作用的情況。圖11給出了計算得到的翼型表面壓力系數(shù)和上翼面摩擦力系數(shù)分布，并同時給出了實驗數(shù)據(jù)，可以看出計算出的激波位置和強(qiáng)度等主要特征與實驗數(shù)據(jù)［13］吻合較好。

圖9 RAE2822翼型周圍的無網(wǎng)格點和直角網(wǎng)格Fig.9 Points and Cartesian grid around RAE2822airfoil

圖10 流場等密度線Fig.10 Density contours of the flow field

圖11 翼型表面壓力系數(shù)和摩擦力系數(shù)分布Fig.11 Cpand Cfdistribution around RAE2822airfoil

4.3 GA（W）－1兩段翼型粘性繞流

在采用前兩個算例對算法進(jìn)行驗證的基礎(chǔ)上，本文還對GA（W）－1兩段翼型的粘性流動進(jìn)行了數(shù)值模擬。GA（W）－1兩段翼型由主翼和襟翼構(gòu)成，兩者之間有一個狹縫。文獻(xiàn)［2］在塊對接結(jié)構(gòu)網(wǎng)格上求解了該翼型的粘性流動問題，從該文塊結(jié)構(gòu)網(wǎng)格生成過程可以看出，對此兩段翼型所在流場生成塊結(jié)構(gòu)網(wǎng)格時，需要首先根據(jù)主翼、襟翼之間的相對位置對計算區(qū)域進(jìn)行分塊，再在每一個塊中生成結(jié)構(gòu)網(wǎng)格，而分塊操作目前還不能完全實現(xiàn)自動化，需要技巧。本文采用無網(wǎng)格／直角網(wǎng)格混合方法處理，不需要預(yù)先對計算區(qū)域進(jìn)行分塊，而是采用第2.1節(jié)中描述的方法進(jìn)行無網(wǎng)格布點及直角網(wǎng)格劃分，在布點過程和直角網(wǎng)格劃分過程中不需要人工干預(yù)。圖2已給出了翼型周圍的無網(wǎng)格點和直角網(wǎng)格，其中主翼上布置了250個點，襟翼上布置了190個點，無網(wǎng)格點和直角網(wǎng)格點共計29768個，第一層無網(wǎng)格點距物面的距離為主翼弦長的5.0×10－5倍。圖12給出了來流條件為Ma＝0.21，α＝10°，Re＝2.3×106時計算得到的流場等密度線，可以看到等密度線在直角網(wǎng)格區(qū)和無網(wǎng)格區(qū)以及兩個無網(wǎng)格區(qū)之間過渡都比較光滑。圖13給出了計算得到的翼型表面壓力系數(shù)分布，并與實驗數(shù)據(jù)［14］和文獻(xiàn)［2］采用塊結(jié)構(gòu)網(wǎng)格有限體積法得到的結(jié)果進(jìn)行了比較，可以看出，本文結(jié)果和實驗數(shù)據(jù)以及文獻(xiàn)計算結(jié)果符合得都比較好。

圖12 流場等密度線Fig.12 Density contours of the flow field

圖13 翼型表面壓力系數(shù)分布Fig.13 Cpdistribution around GA（W）－1airfoil

5 結(jié) 論

本文成功地將基于計算區(qū)域整體非結(jié)構(gòu)直角網(wǎng)格物面附近局部無網(wǎng)格的無網(wǎng)格／直角網(wǎng)格混合算法發(fā)展用于求解粘性流動問題。采用局部四叉樹加密方法生成整體非結(jié)構(gòu)直角網(wǎng)格，實現(xiàn)了對流場大部分區(qū)域的快速網(wǎng)格劃分，簡化了流場計算的前置處理。算例表明：算法涉及的各向異性點云處理，能有效捕捉邊界層的流動特征，捕捉的激波及邊界層清晰。所發(fā)展的混合算法不僅適合處理單個翼型的粘性流動問題，也適合處理多段翼型等物體間存在相互干擾的復(fù)雜流動問題。

［1］ZEEUW D D，POWELL K G.An adaptively refined Cartesian mesh solver for the Euler equations［J］.JournalofComputationalPhysics，1993，104（1）：56－68.

［2］GUO T Q.LU Z L.N－S equation calculations on multielement airfoils with zonal patched grids［J］.TransactionsofNanjingUniversityofAeronautics＆Astronautics，2003，20（2）：155－158.

［3］BATINA J T.A gridless Euler／Navier－Stokes solution algorithm for complex aircraft applications［R］.AIAA 93－0333，1993.

［4］馬志華.自適應(yīng)無網(wǎng)格及網(wǎng)格和無網(wǎng)格混合算法研究［D］.南京：南京航空航天大學(xué)，2008.

［5］MA Z H，EMAD N A，CHEN H Q.A local meshless method for solving comprissible Euler equations［J］.SpaceResearchJournal，2008，1（1）：1－16.

［6］BLAZEK J.Computational fluid dynamics：principles and applications［M］.Amsterdam：Elsevier Science Ltd publisher，2001：171－173.

［7］SPALART P R，ALLMARAS S R.A one－equation turbulence model for aerodynamic flows［R］.AIAA－92－0439.

［8］MORINISHI K.Gridless type solution for high Reynolds number multielement flow fields［R］.AIAA 95－1856，1995.

［9］陳紅全.隱式無網(wǎng)格算法及其應(yīng)用研究［J］.空氣動力學(xué)學(xué)報，2002，20（2）：133－140.

［10］ANG S J，YEP K S，CHEW C S，et al.A singular－value decomposition（SVD）－based generalized finite difference（GFD）method for close－interaction moving boundary flow problems［J］.InternationalJournalforNumericalMethodsinEngineering，2008，76：1892－1929.

［11］閻超.計算流體力學(xué)方法及應(yīng)用［M］.北京：北京航空航天大學(xué)出版社，2006：252－254.

［12］THIBERT J J，GRANJACQUES M，OHMAN L H.NACA0012airfoil［R］.AGARD AR 138，1979.

［13］COOK P H，McDONALD M A，F(xiàn)IRMIN M C P.Aerofoil RAE 2822－pressure distributions，and boundary layer and wake measurements［R］.AGARD AR 138，1979.

［14］WENTZ W H，SEETHARM H C Jr.Development of a fowler flap system for a high performance general aviation airfoil［R］.NASA CR－2443，1974.