蒲賽虎，陳紅全

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

0 引 言

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

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

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

1 控制方程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

在形成的無網(wǎng)格區(qū)中，除外邊界點(diǎn)外，其它點(diǎn)的點(diǎn)云結(jié)構(gòu)可以在布點(diǎn)過程中根據(jù)層推進(jìn)結(jié)構(gòu)同時(shí)確定，邊界點(diǎn)點(diǎn)云結(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ū)交界處的點(diǎn)云結(jié)構(gòu)打下了基礎(chǔ)。圖2中給出了GA（W）－1兩段翼型周圍的無網(wǎng)格點(diǎn)和非結(jié)構(gòu)直角網(wǎng)格。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中，uik、vik、ρik、μL、μT、cik分別表示中心點(diǎn)i與其第k個(gè)衛(wèi)星點(diǎn)之間中點(diǎn)處的速度分量、密度、粘性系數(shù)以及音速。αik和βik為式（4）所確定的系數(shù)。在直角網(wǎng)格點(diǎn)上，當(dāng)?shù)貢r(shí)間歩長也可按式（14）計(jì)算得到。數(shù)值求解時(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í)驗(yàn)數(shù)據(jù)可供對照，NACA0012翼型亞音速粘性繞流計(jì)算經(jīng)常用來考核算法。此算例的來流條件是Ma＝0.5，α＝3.51°，Re＝2.93×106。本文混合算法采用的無網(wǎng)格點(diǎn)和直角網(wǎng)格如圖4中所示，其中翼型表面布置了298個(gè)點(diǎn)，直角網(wǎng)格點(diǎn)和無網(wǎng)格點(diǎn)總數(shù)為20235個(gè)，第一層無網(wǎng)格點(diǎn)距物面的距離為翼型弦長的5.0×10－5倍。圖5中給出了計(jì)算得到的流場等密度線分布，可以看到等密度線在網(wǎng)格區(qū)和無網(wǎng)格區(qū)之間過渡光滑（圖中翼型周圍的黑線是無網(wǎng)格區(qū)外邊界），說明本文采用的兩個(gè)區(qū)之間的流場信息交換方法是可行的。圖6給出了翼型上表面30%弦長附近的速度型，該速度型反映出高雷諾數(shù)粘性流動邊界層內(nèi)速度梯度較大的特點(diǎn)。圖7給出了計(jì)算得到的翼型表面壓力系數(shù)分布，可以看到計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)［12］吻合較好。另外我們還對迎角較大的情況（Ma＝0.301，α＝9.79°，Re＝2.6×106）進(jìn)行了計(jì)算，圖8給出了計(jì)算得到的翼型表面壓力系數(shù)分布，可以看到計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)［12］也比較吻合。

圖4 NACA0012翼型周圍的無網(wǎng)格點(diǎn)和直角網(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í)翼型表面壓力系數(shù)分布Fig.7 Cpdistribution around NACA0012 with the angle of attack 3.51°

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

4.2 RAE2822翼型跨音速粘性繞流

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

圖9 RAE2822翼型周圍的無網(wǎng)格點(diǎn)和直角網(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兩段翼型粘性繞流

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

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

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

5 結(jié) 論

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