蒲賽虎，陳紅全

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

0 引 言

早期數(shù)值方法采用的直角網(wǎng)格具有簡(jiǎn)單、規(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)格不易控制的正交性對(duì)解有較大的影響。目前，許多學(xué)者都在對(duì)直角網(wǎng)格方法［1］和貼體網(wǎng)格方法［2］進(jìn)行研究和改進(jìn)。

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

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

1 控制方程

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

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

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

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

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

本文采用的混合算法由無(wú)網(wǎng)格方法和直角網(wǎng)格方法構(gòu)成，其中無(wú)網(wǎng)格方法的基本計(jì)算單元稱為點(diǎn)云。采用無(wú)網(wǎng)格方法，需要對(duì)計(jì)算區(qū)域進(jìn)行布點(diǎn)離散，再對(duì)每一點(diǎn)構(gòu)成合適的點(diǎn)云。圖1（a）給出了一個(gè)七點(diǎn)無(wú)網(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］的方法，對(duì)二次極小曲面逼近進(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)的距離。需要指出的是，對(duì)于本文粘性流動(dòng)計(jì)算所涉及到的各向異性點(diǎn)云（圖1b中所示），切向點(diǎn)距與法向點(diǎn)距之比（拉伸比）可能非常大，甚至達(dá)到100倍以上，若不采用加權(quán)（權(quán)系數(shù)wik＝1），則在采用二次極小曲面逼近式（5）時(shí)系數(shù)矩陣會(huì)出現(xiàn)病態(tài)（比如當(dāng)拉伸比為100時(shí)，系數(shù)矩陣的條件數(shù)為10000），而本文在最小二乘方法中加入上述權(quán)系數(shù)則避免了系數(shù)矩陣病態(tài)的問(wèn)題，比如采用權(quán)系數(shù)后，當(dāng)拉伸比為100時(shí)，系數(shù)矩陣的條件數(shù)約為4，表明系數(shù)矩陣是健康的。

圖1 無(wú)網(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)格歩長(zhǎng)（本文生成的直角網(wǎng)格每一點(diǎn)上x(chóng)方向和y方向的網(wǎng)格歩長(zhǎng)相同），fiE、fiW、fiS、fiN分別是i點(diǎn)和其東、西、南、北方向上鄰近點(diǎn)之間中點(diǎn)處的函數(shù)值。

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

采用無(wú)網(wǎng)格／直角網(wǎng)格混合算法，計(jì)算區(qū)域整體采用直角網(wǎng)格，只在物面附近采用局部無(wú)網(wǎng)格處理。求解無(wú)粘流動(dòng)問(wèn)題時(shí)，物面附近的局部無(wú)網(wǎng)格通常采用圖1（a）所示的各向同性點(diǎn)云［3，9］，但對(duì)本文涉及到的粘性流動(dòng)，特別是高雷諾數(shù)粘性流動(dòng)而言，為了有效模擬邊界層內(nèi)較大的法向速度梯度，沿物面法線方向布點(diǎn)需較密，此時(shí)若仍采用各向同性點(diǎn)云，勢(shì)必造成整體布點(diǎn)量太大，不利于減少計(jì)算時(shí)間?；诖耍疚脑谖锩娓浇捎昧藞D1（b）所示的各向異性點(diǎn)云。采用各向異性點(diǎn)云，物面法線方向布點(diǎn)較密，從而能夠模擬邊界層內(nèi)的法向速度梯度，而物面切線方向點(diǎn)距相對(duì)較大，有利于減少布點(diǎn)總量，節(jié)約計(jì)算時(shí)間。對(duì)單個(gè)物體，本文的布點(diǎn)步驟是：首先將物面邊界劃分成離散點(diǎn)的形式，然后沿物面法線方向按層推進(jìn)的方法進(jìn)行布點(diǎn)［9］。在層推進(jìn)過(guò)程中，我們引入了指數(shù)函數(shù)來(lái)控制每一層點(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é)距離增長(zhǎng)速度的參數(shù)，本文取為0.1。對(duì)于流場(chǎng)中存在多個(gè)物體相互干擾的情況，本文在對(duì)每個(gè)物體單獨(dú)采用上述方法進(jìn)行布點(diǎn)的同時(shí)，引入了無(wú)網(wǎng)格區(qū)之間的布點(diǎn)干涉機(jī)制，以解決各個(gè)無(wú)網(wǎng)格區(qū)之間可能相互重疊的問(wèn)題，即在對(duì)每個(gè)物體同時(shí)進(jìn)行層推進(jìn)布點(diǎn)過(guò)程中，若發(fā)現(xiàn)本無(wú)網(wǎng)格區(qū)某點(diǎn)周?chē)木嚯x范圍內(nèi)存在其它無(wú)網(wǎng)格區(qū)的點(diǎn)，則在該點(diǎn)處，層推進(jìn)布點(diǎn)停止，而其它點(diǎn)上的層推進(jìn)布點(diǎn)操作繼續(xù)進(jìn)行。δ是一個(gè)距離閥值，可根據(jù)具體計(jì)算對(duì)象確定。

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

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

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

圖2 GA（W）－1翼型周?chē)臒o(wú)網(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)格生成過(guò)程確定。k為縮放系數(shù)，本文研究發(fā)現(xiàn)取1.5是一個(gè)合適的值。由于本文直角網(wǎng)格生成采用了四叉樹(shù)數(shù)據(jù)結(jié)構(gòu)，因此上述候選點(diǎn)可以根據(jù)既有的四叉樹(shù)數(shù)據(jù)結(jié)構(gòu)快速確定；

2）對(duì)于某一候選點(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)過(guò)上述判斷，圖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)過(guò)步驟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)則確定出來(lái)的k點(diǎn)的點(diǎn)云結(jié)構(gòu)。可見(jiàn)采用上述選點(diǎn)準(zhǔn)則，選出的衛(wèi)星點(diǎn)圍繞中心點(diǎn)的分布比較均勻，這有助于保證點(diǎn)云上空間導(dǎo)數(shù)的離散精度。

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

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

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

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

而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)格方法離散過(guò)程：

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

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

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

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

式中σ為CFL數(shù)，C為常數(shù)，本文取值為4。對(duì)于無(wú)網(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í)間歩長(zhǎng)也可按式（14）計(jì)算得到。數(shù)值求解時(shí)還需補(bǔ)充邊界條件，本文物面邊界采用無(wú)滑移邊界條件，遠(yuǎn)場(chǎng)則采用無(wú)反射邊界條件［11］。

4 算例與分析

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

4.1 NACA0012翼型亞音速粘性繞流

由于有實(shí)驗(yàn)數(shù)據(jù)可供對(duì)照，NACA0012翼型亞音速粘性繞流計(jì)算經(jīng)常用來(lái)考核算法。此算例的來(lái)流條件是Ma＝0.5，α＝3.51°，Re＝2.93×106。本文混合算法采用的無(wú)網(wǎng)格點(diǎn)和直角網(wǎng)格如圖4中所示，其中翼型表面布置了298個(gè)點(diǎn)，直角網(wǎng)格點(diǎn)和無(wú)網(wǎng)格點(diǎn)總數(shù)為20235個(gè)，第一層無(wú)網(wǎng)格點(diǎn)距物面的距離為翼型弦長(zhǎng)的5.0×10－5倍。圖5中給出了計(jì)算得到的流場(chǎng)等密度線分布，可以看到等密度線在網(wǎng)格區(qū)和無(wú)網(wǎng)格區(qū)之間過(guò)渡光滑（圖中翼型周?chē)暮诰€是無(wú)網(wǎng)格區(qū)外邊界），說(shuō)明本文采用的兩個(gè)區(qū)之間的流場(chǎng)信息交換方法是可行的。圖6給出了翼型上表面30%弦長(zhǎng)附近的速度型，該速度型反映出高雷諾數(shù)粘性流動(dòng)邊界層內(nèi)速度梯度較大的特點(diǎn)。圖7給出了計(jì)算得到的翼型表面壓力系數(shù)分布，可以看到計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)［12］吻合較好。另外我們還對(duì)迎角較大的情況（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翼型周?chē)臒o(wú)網(wǎng)格點(diǎn)和直角網(wǎng)格Fig.4 Points and Cartesian grid around NACA0012airfoil

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

圖9 RAE2822翼型周?chē)臒o(wú)網(wǎng)格點(diǎn)和直角網(wǎng)格Fig.9 Points and Cartesian grid around RAE2822airfoil

圖10 流場(chǎng)等密度線Fig.10 Density contours of the flow field

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

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

在采用前兩個(gè)算例對(duì)算法進(jìn)行驗(yàn)證的基礎(chǔ)上，本文還對(duì)GA（W）－1兩段翼型的粘性流動(dòng)進(jìn)行了數(shù)值模擬。GA（W）－1兩段翼型由主翼和襟翼構(gòu)成，兩者之間有一個(gè)狹縫。文獻(xiàn)［2］在塊對(duì)接結(jié)構(gòu)網(wǎng)格上求解了該翼型的粘性流動(dòng)問(wèn)題，從該文塊結(jié)構(gòu)網(wǎng)格生成過(guò)程可以看出，對(duì)此兩段翼型所在流場(chǎng)生成塊結(jié)構(gòu)網(wǎng)格時(shí)，需要首先根據(jù)主翼、襟翼之間的相對(duì)位置對(duì)計(jì)算區(qū)域進(jìn)行分塊，再在每一個(gè)塊中生成結(jié)構(gòu)網(wǎng)格，而分塊操作目前還不能完全實(shí)現(xiàn)自動(dòng)化，需要技巧。本文采用無(wú)網(wǎng)格／直角網(wǎng)格混合方法處理，不需要預(yù)先對(duì)計(jì)算區(qū)域進(jìn)行分塊，而是采用第2.1節(jié)中描述的方法進(jìn)行無(wú)網(wǎng)格布點(diǎn)及直角網(wǎng)格劃分，在布點(diǎn)過(guò)程和直角網(wǎng)格劃分過(guò)程中不需要人工干預(yù)。圖2已給出了翼型周?chē)臒o(wú)網(wǎng)格點(diǎn)和直角網(wǎng)格，其中主翼上布置了250個(gè)點(diǎn)，襟翼上布置了190個(gè)點(diǎn)，無(wú)網(wǎng)格點(diǎn)和直角網(wǎng)格點(diǎn)共計(jì)29768個(gè)，第一層無(wú)網(wǎng)格點(diǎn)距物面的距離為主翼弦長(zhǎng)的5.0×10－5倍。圖12給出了來(lái)流條件為Ma＝0.21，α＝10°，Re＝2.3×106時(shí)計(jì)算得到的流場(chǎng)等密度線，可以看到等密度線在直角網(wǎng)格區(qū)和無(wú)網(wǎng)格區(qū)以及兩個(gè)無(wú)網(wǎng)格區(qū)之間過(guò)渡都比較光滑。圖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 流場(chǎng)等密度線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ú)網(wǎng)格的無(wú)網(wǎng)格／直角網(wǎng)格混合算法發(fā)展用于求解粘性流動(dòng)問(wèn)題。采用局部四叉樹(shù)加密方法生成整體非結(jié)構(gòu)直角網(wǎng)格，實(shí)現(xiàn)了對(duì)流場(chǎng)大部分區(qū)域的快速網(wǎng)格劃分，簡(jiǎn)化了流場(chǎng)計(jì)算的前置處理。算例表明：算法涉及的各向異性點(diǎn)云處理，能有效捕捉邊界層的流動(dòng)特征，捕捉的激波及邊界層清晰。所發(fā)展的混合算法不僅適合處理單個(gè)翼型的粘性流動(dòng)問(wèn)題，也適合處理多段翼型等物體間存在相互干擾的復(fù)雜流動(dòng)問(wèn)題。

［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ú)網(wǎng)格及網(wǎng)格和無(wú)網(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ú)網(wǎng)格算法及其應(yīng)用研究［J］.空氣動(dòng)力學(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.