沈志偉，趙寧，胡偶

（1．南京航空航天大學(xué) 航空宇航學(xué)院，江蘇 南京 210016；2．中航工業(yè)直升機(jī)設(shè)計(jì)研究所，江西 景德鎮(zhèn) 333000）

可壓縮粘性流動(dòng)笛卡爾網(wǎng)格虛擬單元方法研究

沈志偉1，趙寧1，胡偶2

（1．南京航空航天大學(xué) 航空宇航學(xué)院，江蘇 南京 210016；2．中航工業(yè)直升機(jī)設(shè)計(jì)研究所，江西 景德鎮(zhèn) 333000）

以二維高雷諾數(shù)可壓縮粘性流動(dòng)問題為背景，提出了一種全新的笛卡爾網(wǎng)格虛擬單元方法?；诒诿婧瘮?shù)基本假設(shè)，構(gòu)造了壁面函數(shù)-虛擬單元方法（WE-GCM），用于定義湍流壁面邊界條件。引入?yún)⒖键c(diǎn)的概念計(jì)算虛擬單元上的基本變量與湍流變量值，定義了“非貼體”笛卡爾網(wǎng)格下的湍流壁面邊界條件，并通過壁面函數(shù)模型修正近壁面單元與界面單元?；谧赃m應(yīng)笛卡爾網(wǎng)格體系，采用發(fā)展的具有二階精度的格心格式有限體積求解器，數(shù)值模擬了跨音速RAE2822翼型繞流問題與超音速圓柱繞流問題，計(jì)算結(jié)果與實(shí)驗(yàn)值吻合良好，顯示了WEGCM對(duì)高雷諾數(shù)可壓縮粘性問題是有效的。

虛擬單元方法；壁面函數(shù)；自適應(yīng)笛卡爾網(wǎng)格；湍流

0 引 言

網(wǎng)格生成是采用計(jì)算流體力學(xué)（Computational Eluid Dynamics，CED）進(jìn)行數(shù)值模擬的先行工作，目前較為成熟的網(wǎng)格形式有結(jié)構(gòu)網(wǎng)格、非結(jié)構(gòu)網(wǎng)格和笛卡爾網(wǎng)格。笛卡爾網(wǎng)格是一種區(qū)別于前兩者的非貼體網(wǎng)格類型。由于早期多采用結(jié)構(gòu)化的直角網(wǎng)格形式，在處理復(fù)雜外形時(shí)所需要的網(wǎng)格數(shù)目巨大，并未得到廣泛的應(yīng)用。隨著計(jì)算機(jī)水平和數(shù)值方法的發(fā)展，尤其是自適應(yīng)網(wǎng)格加密（Adaptive Mesh Refinement，AMR）技術(shù)的提出，使得笛卡爾網(wǎng)格在復(fù)雜流動(dòng)問題中凸顯出了其獨(dú)特的優(yōu)勢(shì)：1）簡化復(fù)雜幾何外形的網(wǎng)格生成過程，不需要過多的人為技巧和經(jīng)驗(yàn)；2）易于實(shí)現(xiàn)基于幾何外形和流場(chǎng)特征的網(wǎng)格自適應(yīng)技術(shù)，能夠減少計(jì)算網(wǎng)格的數(shù)量并且提高流場(chǎng)的分辨率；3）笛卡爾網(wǎng)格系統(tǒng)易于和高精度數(shù)值計(jì)算方法耦合。

采用純笛卡爾網(wǎng)格方法進(jìn)行數(shù)值模擬的最大難點(diǎn)在于物面條件的滿足。早期的學(xué)者們多采用切割單元［1］的方法，即切割與物體表面相交的笛卡爾網(wǎng)格單元，只保留浸沒在流場(chǎng)中的部分，使得笛卡爾網(wǎng)格具備了貼體的特性。但是，該方法具有自身的局限性：切割后的單元形式多樣，使原有的網(wǎng)格數(shù)據(jù)類型變得復(fù)雜，三維情況時(shí)會(huì)更加明顯；同時(shí)切割過程中可能形成微小的網(wǎng)格單元，造成方程系統(tǒng)的剛性問題，影響流場(chǎng)的收斂特性并在物面邊界處產(chǎn)生流場(chǎng)解的非物理振蕩。為此Udaykumar［2］提出了融合單元處理技術(shù)，Munikrishna等人［3］提出了網(wǎng)格縫合處理技術(shù)對(duì)切割單元方法進(jìn)行了改進(jìn)。另一方面，浸入類邊界處理方法也引起越來越多的關(guān)注，其與切割單元方法相比避免了復(fù)雜的幾何求交運(yùn)算以及小網(wǎng)格單元出現(xiàn)所帶來的諸多缺點(diǎn)。浸入邊界方法（Immersed Boundary Method，IBM）最初在低速不可壓流動(dòng)中廣泛應(yīng)用，成功地?cái)?shù)值模擬了流固耦合、大位移運(yùn)動(dòng)邊界等諸多問題。在可壓縮流動(dòng)方面，Eorrer和Jeltsch［4］首次基于笛卡爾網(wǎng)格引入虛擬單元方法處理無粘可壓縮流動(dòng)問題，該方法不考慮物面的存在，把物體作為邊界條件嵌入到流場(chǎng)中，而邊界條件通過某種重構(gòu)強(qiáng)加在物體內(nèi)部的虛擬單元上。隨后Dadone等人［5］對(duì)虛擬單元方法進(jìn)行了細(xì)致深入的研究，考慮物面曲率的影響引入了曲率修正和熵修正。目前，該類方法的研究熱點(diǎn)集中在可壓縮粘性流動(dòng)問題的數(shù)值模擬，尤其是面向高雷諾數(shù)的湍流流動(dòng)問題。相對(duì)于無粘和層流流動(dòng)問題，浸入類笛卡爾網(wǎng)格方法應(yīng)用到湍流流動(dòng)時(shí)面臨的難點(diǎn)有：1）高雷諾數(shù)下如何準(zhǔn)確定義湍流壁面邊界條件，模擬壁面湍流效應(yīng)；2）在不影響流場(chǎng)計(jì)算精度及邊界層分辨的情況下，如何降低計(jì)算網(wǎng)格數(shù)量，尤其是物面附近的網(wǎng)格數(shù)量，提高計(jì)算效率。一些學(xué)者基于浸入類笛卡爾網(wǎng)格方法使用十分細(xì)密的網(wǎng)格給出了湍流流動(dòng)的數(shù)值模擬結(jié)果，與貼體類網(wǎng)格方式相比，壁面附近網(wǎng)格需求量巨大導(dǎo)致了計(jì)算量很大，限制了方法在復(fù)雜幾何外形流動(dòng)和三維流動(dòng)問題中的使用。為此，從減少物面邊界附近網(wǎng)格數(shù)量的目的出發(fā)，使用壁面函數(shù)模型來定義湍流邊界條件成為一種可行的方式：Kalitzin和Iaccarino［6］提出了一種基于壁面律的浸入邊界方法處理高雷諾數(shù)湍流問題；Jae-doo Lee［7］基于虛擬單元浸入邊界方法，通過耦合壁面函數(shù)模型與k-ε湍流模型來定義壁面邊界條件；Capizzano［8］采用Two-layer wall modeling修正浸入邊界網(wǎng)格交接面上湍流變量值。

國內(nèi)有關(guān)浸入邊界笛卡爾網(wǎng)格方法在高雷諾數(shù)可壓縮流動(dòng)問題中的研究和應(yīng)用還并不多見，作者在前人的工作基礎(chǔ)上，基于浸入邊界虛擬單元方法，研究了貼體網(wǎng)格下壁面函數(shù)模型的作用和實(shí)現(xiàn)模式，從壁面函數(shù)的基本假設(shè)出發(fā)構(gòu)造了一種壁面函數(shù)-虛擬單元方法（Wall Eunction-Ghost Cell Method，WEGCM）。本文基于自適應(yīng)笛卡爾網(wǎng)格和有限體積方法，發(fā)展了一套可用于求解二維高雷諾數(shù)可壓縮粘性流動(dòng)的求解器，耦合SSTk-ω湍流模型與WE-GCM，成功定義了高雷諾數(shù)下的湍流壁面條件，并依據(jù)壁面函數(shù)模型修正近壁面單元和界面單元上的湍流變量值和湍流粘性系數(shù)。

1 數(shù)值方法

1.1 控制方程的求解

積分守恒形式的RANS方程表達(dá)式如下：

守恒變量W，對(duì)流通量Fc和黏性通量Fv的具體形式參見文獻(xiàn)［9］。

本文采用了基于非結(jié)構(gòu)網(wǎng)格數(shù)據(jù)形式的格心格式有限體積方法來求解上式。時(shí)間推進(jìn)采用LU-SGS（Lower-Upper Symmetric Gauss-Seidel）隱式格式［10］；空間離散上選取AUSM＋（Advection Upstream Splitting Method Plus）［11］近似Riemann求解器并加入Venkatakrishnan限制器［12］來抑制間斷區(qū)域的非物理振蕩，同時(shí)為獲得解的二階精度進(jìn)行了解的線性重構(gòu)。湍流模型選取SSTk-ω兩方程模型［13］。遠(yuǎn)場(chǎng)邊界采用特征邊界條件［9］；采用虛擬單元方法來定義湍流壁面邊界條件，具體內(nèi)容在第2節(jié)給出。

1.2 基于流場(chǎng)特征的網(wǎng)格自適應(yīng)方法

笛卡爾網(wǎng)格生成過程中，二維采用四叉樹（三維八叉數(shù)）的數(shù)據(jù)結(jié)構(gòu)來儲(chǔ)存網(wǎng)格信息，叉數(shù)狀數(shù)據(jù)結(jié)構(gòu)的最大優(yōu)勢(shì)就是能夠十分便利地進(jìn)行網(wǎng)格的加密和粗化，網(wǎng)格自適應(yīng)技術(shù)的產(chǎn)生和發(fā)展也直接推動(dòng)了笛卡爾網(wǎng)格在流動(dòng)模擬中的應(yīng)用。本文為了提高流場(chǎng)結(jié)構(gòu)的分辨率基于流場(chǎng)特征進(jìn)行了解自適應(yīng)操作。文獻(xiàn)［14］對(duì)比了不同的解自適應(yīng)判據(jù)，認(rèn)為結(jié)合速度散度和旋度的綜合判據(jù)是十分有效的，表達(dá)式如下：

計(jì)算過程中自適應(yīng)判斷標(biāo)準(zhǔn)如下：

1）如果網(wǎng)格單元滿足τci＞σc或者τdi＞σd，則網(wǎng)格需要加密；

2）如果網(wǎng)格單元滿足τci＜0．3σc或者τdi＜0．1σd，則網(wǎng)格需要粗化。

2 壁面函數(shù)-虛擬單元方法

2.1 壁面函數(shù)模型

壁面函數(shù)模型是一種由邊界層理論推導(dǎo)得到的經(jīng)驗(yàn)?zāi)Ｐ停浣⑿枰獫M足六條假設(shè)［15］：

1）邊界層底部區(qū)域的速度分布和溫度分布可以解析的表達(dá)；

2）近壁面第一層網(wǎng)格點(diǎn)上的湍流輸運(yùn)參數(shù)可以解析的表達(dá)；

3）邊界層底部區(qū)域壓力分布是常數(shù)；

4）邊界層底部區(qū)域剪切應(yīng)力分布是常數(shù)；

5）邊界層底部區(qū)域熱量傳遞是常數(shù)；

6）邊界層底部區(qū)域內(nèi)沒有化學(xué)反應(yīng)。

文獻(xiàn)［15］依據(jù)上述假設(shè)考慮流體可壓縮性以及物面熱傳導(dǎo)的影響，修正了Spalding的壁面函數(shù)模型［16］，給出了一種新的函數(shù)模型，定義如下：

式中y＋和u＋分別代表無量綱的邊界層厚度和邊界層內(nèi)的速度，定義如下：

式（4）式（5）中vt為切向速度大小，y為單元中心到壁面的垂直距離，ρw、μw、τw、uτ分別為壁面流體的密度、層流粘性系數(shù)、剪切應(yīng)力和摩擦速度，常數(shù)κ和B分別0．41和5．5。修正項(xiàng)表達(dá)式如下：

式中Г、β分別表示可壓縮性和熱傳導(dǎo)的影響。r＝Pr1／3（Pr＝0．72），Cp為空氣定壓比熱容，qw、Tw、kw分別為壁面熱流密度、壁面溫度和壁面熱傳導(dǎo)系數(shù)。本文考慮的粘性問題都為絕熱壁，因此β＝0，qw＝0。

2.2 虛擬單元方法

虛擬單元方法中，笛卡爾網(wǎng)格單元被分為三類：界面單元，固體單元和流體單元。將與界面單元相鄰的固體單元定義為虛擬單元（Ghost Cell，GC），將與界面單元相鄰的流體單元定義為近壁面單元（Near Wall Cell，NWC）。計(jì)算過程中流體單元和界面單元構(gòu)成計(jì)算單元，而固體單元不參與計(jì)算。無粘虛擬單元方法在文獻(xiàn)［17］中已有詳細(xì)的介紹，其基本思想是在壁面附近通過邊界條件的線性重構(gòu)獲得虛擬單元上的虛擬流場(chǎng)信息。但是高雷諾數(shù)可壓流動(dòng)問題中附面層結(jié)構(gòu)是非線性的，沿用原有的線性重構(gòu)會(huì)帶來如下的兩點(diǎn)問題：1）笛卡爾網(wǎng)格的非貼體性和虛擬單元的對(duì)稱點(diǎn)的不均勻分布可能引起附面層內(nèi)流動(dòng)的非物理振蕩。2）虛擬單元的對(duì)稱點(diǎn)必須落在附面層中的粘性底層區(qū)域，這就要求壁面附近的笛卡爾網(wǎng)格單元尺度很小，因此網(wǎng)格數(shù)量會(huì)劇烈增加。

粘性計(jì)算時(shí)，為解決上述兩點(diǎn)問題，引入?yún)⒖键c(diǎn)的概念來取代無粘情況下的對(duì)稱點(diǎn)。同時(shí)，基于壁面函數(shù)修正湍流壁面邊界條件從而改變物面上的通量計(jì)算。

2.2.1 參考點(diǎn)與壁面函數(shù)相關(guān)參數(shù)的確定

如圖1所示，點(diǎn)E，G分別為虛擬單元A，E的參考點(diǎn)，參考點(diǎn)的取法以A點(diǎn)為例：首先找到A點(diǎn)距離物面的最近點(diǎn)B，延長AB至E使得BE的長度為物面網(wǎng)格尺度h的倍。與無粘情況時(shí)取關(guān)于物面的對(duì)稱點(diǎn)不同，該種取法中所有參考點(diǎn)到物面的距離都一致，從而具有了結(jié)構(gòu)網(wǎng)格的一些特點(diǎn)，同時(shí)該種取法保證了插值獲得參考點(diǎn)流場(chǎng)基本變量值過程中用到的四個(gè)模板點(diǎn)S1、S2、S3、S4都為流體單元。根據(jù)壁面函數(shù)的基本假設(shè)（3）和Crocco-Busemann關(guān)系得到壁面上的壓力和溫度：

下標(biāo)ref表示參考點(diǎn)。由狀態(tài)方程和Sutherland公式可以計(jì)算得到壁面的密度ρw和層流粘性系數(shù)μw。參考點(diǎn)上法向和切向速度分量分別為：

將VTref、ρw、μw結(jié)合該點(diǎn)到物面的距離δr帶入式（4）式（5）求解壁面摩擦應(yīng)力τw與壁面摩擦速度uτ。

圖1 虛擬單元方法示意圖Fig.1 Example configuration of the ghost cell method

2.2.2 虛擬單元流場(chǎng)值的確定

壁面函數(shù)基本假設(shè)式（4）認(rèn)為壁面附近的剪切應(yīng)力為常數(shù)分布，根據(jù)牛頓摩擦阻力公式：

可以得到虛擬單元上的切向速度分量為：

同時(shí)虛擬單元上法向速度分量滿足無穿透條件：

虛擬單元上壓力和溫度的計(jì)算方法同上節(jié)中物面點(diǎn)類似：

密度ρg可由狀態(tài)方程解出。

結(jié)合式（4）式（5）將y＋對(duì)u＋求導(dǎo)得到近壁面處湍流粘性系數(shù)與層流粘性系數(shù)的關(guān)系式：

根據(jù)虛擬單元和參考點(diǎn)上總粘性系數(shù)相等，則可以求得虛擬單元A上的湍流粘性系數(shù)為：

SSTk-ω兩方程模型中k和ω可以采用下式得到：

其中常數(shù)Cμ和κ分別取值為0．09與0．41。

Nichols等人指出，采用壁面函數(shù)處理湍流問題時(shí)，除了上述提及的邊界條件的定義之外，靠近物面的網(wǎng)格單元還需要進(jìn)行湍流變量值和湍流粘性系數(shù)的修正，此處虛擬單元方法中對(duì)近壁面單元和界面單元都根據(jù)壁面函數(shù)模型做出了修正，具體修正方法參見文獻(xiàn)［15］。

2.2.3 笛卡爾網(wǎng)格WE-GCM的實(shí)施流程

綜上所述，可以將基于笛卡爾網(wǎng)格的WE-GCM流程總結(jié)如下：

1）劃分單元類型，確定流體單元、固體單元、界面單元，定義虛擬單元和近壁面單元，找到虛擬單元的參考點(diǎn)并確定各點(diǎn)的幾何信息；

2）插值計(jì)算參考點(diǎn)的流場(chǎng)信息，得到壁面壓力、溫度、速度、密度、層流粘性系數(shù)，并結(jié)合壁面函數(shù)模型公式（4）求解壁面剪切應(yīng)力和摩擦速度；

3）由式（10）～式（16）計(jì)算得到虛擬單元上的基本流場(chǎng)變量和湍流變量；

4）根據(jù)壁面函數(shù)模型修正近壁面單元和界面單元上的湍流變量值和湍流粘性系數(shù)。

3 數(shù)值算例與分析

3.1 平板絕熱流動(dòng)的數(shù)值模擬

該算例來源于文獻(xiàn)［7］，來流馬赫數(shù)為0．2，基于平板長度的雷諾數(shù)為1．0927×107。計(jì)算采用了貼體的結(jié)構(gòu)網(wǎng)格，通過疏密不同的四套網(wǎng)格以及開關(guān)壁面函數(shù)模型驗(yàn)證了壁面函數(shù)模型的準(zhǔn)確性，并分析了其與網(wǎng)格尺度之間的關(guān)系。平板展向布置120個(gè)網(wǎng)格，法向分別布置50、60、80、100個(gè)網(wǎng)格，對(duì)應(yīng)物面網(wǎng)格的y＋分別為1、5、20、40。從圖2給出的平板表面摩擦系數(shù)分布可以看出，當(dāng)y＋＝40時(shí)，關(guān)閉壁面函數(shù)模型，計(jì)算結(jié)果和實(shí)驗(yàn)值偏差很大，而采用壁面函數(shù)模型后，計(jì)算結(jié)果得到明顯的改善；同時(shí)可以看出采用壁面函數(shù)模型后不同網(wǎng)格尺度下的計(jì)算結(jié)果基本與實(shí)驗(yàn)值相近，對(duì)網(wǎng)格的依賴性減小。表1給出了四套網(wǎng)格下x＝0．685處的壁面摩擦系數(shù)Cf的數(shù)值。數(shù)據(jù)結(jié)果顯示y＋＝5、20、40時(shí)采用壁面函數(shù)模型的計(jì)算結(jié)果與y＋＝1時(shí)的計(jì)算結(jié)果以及實(shí)驗(yàn)值都比較吻合；而當(dāng)y＋＝40不采用壁面函數(shù)模型的情況下，Cf的相對(duì)誤差達(dá)到了26．5%。

3.2 RAE2822翼型跨聲速繞流數(shù)值模擬

Cook等人［18］對(duì)RAE2822翼型的繞流問題進(jìn)行了詳細(xì)的實(shí)驗(yàn)研究，故其常作為測(cè)試流場(chǎng)求解器準(zhǔn)確性的標(biāo)準(zhǔn)算例。本文選取的計(jì)算狀態(tài)為：來流馬赫數(shù)M∞＝0．73，攻角α＝2．79°，雷諾數(shù)Re∞＝6．2×106。

圖2 不同網(wǎng)格下平板表面的摩擦系數(shù)分布Fig.2 Frictiondrag coefficient on the sur face indifferent meshes

表1 不同網(wǎng)格下x＝0.685處的Cf比較Table 1 Comparison of Cfat x＝0.685 indifferent mehes

該狀態(tài)下，翼型上表面會(huì)產(chǎn)生一道明顯的激波，激波與附面層相互作用會(huì)使得附面層變厚并同時(shí)伴隨著激波誘導(dǎo)所產(chǎn)生的輕微流動(dòng)分離，這些都給準(zhǔn)確地?cái)?shù)值模擬增加了難度。計(jì)算域大小取弦長（c＝1．0）的20倍，將計(jì)算域劃分成80×80的網(wǎng)格單元，初始背景笛卡爾網(wǎng)格步長為Hmax＝0．25，依據(jù)翼型表面進(jìn)行了8次幾何自適應(yīng)加密后，最小網(wǎng)格尺度達(dá)到Hmin＝Hmax／28。計(jì)算過程中基于速度旋度和散度判據(jù)進(jìn)行了三次（AMR＝3）解自適應(yīng)加密，網(wǎng)格數(shù)量從初始的33054依次增加為81474、178760、379927。圖3給出了初始網(wǎng)格（圖3a）和三次解自適應(yīng)之后的計(jì)算網(wǎng)格（圖3b）；圖4給出了AMR＝3時(shí)計(jì)算所得的馬赫數(shù)云圖；圖5將采用WE-GCM得到的翼型表面壓力系數(shù)Cp分布與貼體網(wǎng)格的計(jì)算結(jié)果以及實(shí)驗(yàn)值進(jìn)行了比較。從圖5中可以看到，初始計(jì)算網(wǎng)格得到的前緣Cp峰值偏低，激波位置偏差較大，隨著解自適應(yīng)次數(shù)的增加，前緣Cp峰值和激波位置得到了明顯的改善，當(dāng)AMR＝3時(shí)的計(jì)算結(jié)果和貼體網(wǎng)格計(jì)算所得到的結(jié)果十分一致，但是過激波后的Cp分布與實(shí)驗(yàn)值都存在一定的誤差。表2給出了該狀態(tài)下計(jì)算所得的升阻力系數(shù)與實(shí)驗(yàn)值的對(duì)比，從結(jié)果可以看到，伴隨著解自適應(yīng)的進(jìn)行升阻力系數(shù)得到了明顯的改善，在AMR＝3時(shí)升力系數(shù)的誤差小于2%，阻力系數(shù)的誤差小于7%。

圖3 初始網(wǎng)格與解自適應(yīng)加密網(wǎng)格（AMR＝3）Fig.3 Initial mesh and feature-adapted mesh（AMR＝3）

圖4 馬赫數(shù)云圖的計(jì)算結(jié)果（AMR＝3）Fig.4 Computed result of Mach contours（AMR＝3）

圖5 壓力系數(shù)分布的比對(duì)Fig.5 Comparison of Cpdistribution

表2 RAE2822翼型粘性繞流氣動(dòng)力的比較Table 2 Comparison of CLand CDof transonic flow over RAE2822 airfoil

3.3 超聲速圓柱繞流數(shù)值模擬

該算例考察WE-GCM在處理超聲速流動(dòng)問題時(shí)的有效性和準(zhǔn)確性。取來流馬赫數(shù)為1．7，基于圓柱直徑的雷諾數(shù)為2×105，在該狀態(tài)下，圓柱上游會(huì)形成一道弓形激波，過激波后流動(dòng)加速會(huì)在圓柱后緣再次形成兩道激波，同時(shí)伴隨著圓柱后緣的流動(dòng)分離，上述復(fù)雜的流動(dòng)現(xiàn)象都增加了采用WE-GCM進(jìn)行數(shù)值模擬難度。計(jì)算區(qū)域大小取為［-8D，9D］×［-8D，8D］，其中D為圓柱直徑。計(jì)算過程中基于速度旋度和散度判據(jù)進(jìn)行了三次（AMR＝3）解自適應(yīng)加密，最終網(wǎng)格總數(shù)目達(dá)到184684，圖6給出了最終的計(jì)算網(wǎng)格（圖6a）與馬赫數(shù)云圖（圖6b）。圖7給出了圓柱尾緣流動(dòng)分離區(qū)域的流線圖（圖7a）和對(duì)應(yīng)位置的自適應(yīng)計(jì)算網(wǎng)格（圖7b），可以看出WE-GCM能夠捕捉圓柱流動(dòng)分離現(xiàn)象。圖8將WE-GCM計(jì)算所得的圓柱表面Cp分布與實(shí)驗(yàn)值［19］進(jìn)行了對(duì)比，從表3給出的數(shù)據(jù)可以看到，相對(duì)于文獻(xiàn)［20］中Munikrishna采用混合網(wǎng)格方法計(jì)算所得分離點(diǎn)的位置與阻力系數(shù)，WE-GCM計(jì)算所得的分離點(diǎn)位置與阻力系數(shù)與實(shí)驗(yàn)值都更為符合。

表3 圓柱阻力系數(shù)與分離點(diǎn)比較Table 3 Comparison ofθSand CD

圖6 超聲速圓柱繞流的計(jì)算結(jié)果Fig.6 Results of supersonic flow past cylinder

圖7 圓柱尾流附近流線圖及網(wǎng)格細(xì)節(jié)Fig.7 Streamlines and mesh in wake flow near cylinder

圖8 壓力系數(shù)分布的比對(duì)Fig.8 Comparison of Cpdistribution

4 結(jié) 論

針對(duì)高雷諾數(shù)可壓縮粘性流動(dòng)，從幾何角度與邊界重構(gòu)數(shù)值方法兩個(gè)方面出發(fā)，基于笛卡爾網(wǎng)格虛擬單元方法耦合壁面函數(shù)模型，提出了一種全新的湍流壁面邊界條件定義方法，用于定義虛擬單元上的虛擬流場(chǎng)信息，并給出了壁面函數(shù)模型在非貼體笛卡爾網(wǎng)格中的具體實(shí)現(xiàn)模式，在此基礎(chǔ)上結(jié)合有限體積方法結(jié)合發(fā)展了可用于高雷諾數(shù)可壓縮流動(dòng)問題的壁面函數(shù)-虛擬單元方法（WE-GCM）。

本文以貼體網(wǎng)格下平板絕熱流動(dòng)為例，驗(yàn)證了壁面函數(shù)模型的有效性與正確性，進(jìn)而使用WE-GCM數(shù)值模擬了跨聲速RAE2822翼型繞流與超聲速圓柱繞流，通過上述算例的研究可以得到如下結(jié)論：

1）壁面函數(shù)模型的使用能夠在確保數(shù)值模擬準(zhǔn)確性的前提下減少網(wǎng)格數(shù)量；

2）WE-GCM能夠有效準(zhǔn)確地?cái)?shù)值模擬高雷諾數(shù)可壓縮粘性繞流，耦合自適應(yīng)網(wǎng)格加密技術(shù)能夠提高流場(chǎng)模擬的分辨率與氣動(dòng)力參數(shù)計(jì)算的準(zhǔn)確性；

3）盡管壁面函數(shù)模型不能嚴(yán)格證明在流動(dòng)分離時(shí)依然有效，但從數(shù)值結(jié)果來看，WE-GCM處理湍流壁面邊界條件能夠較為準(zhǔn)確地捕捉流動(dòng)分離現(xiàn)象。

［1］COIRIER W J，POWELL K G．Solution-adaptive Cartesian cell approach for viscous and inviscid flows［J］．AIAA Journal，1996，34（5）：938-945．

［2］UDAYKUMAR H S，UDAYKUMAR H S，KAN H C，et al．Multiphasedynamics in arbitrary geometries on fixed cartesian grids［J］．Journal of Computational Physics，1997，137（2）：366-405．

［3］MUNIKRISHNA N，MONDAL P，BALAKRISHNAN N．Cartesian-like grids using a novel grid-stitching algorithm for viscous flow computations［J］．Journal of Aircraft，2007，44（5）：1598-1609．

［4］EORRER H，JELTSCH R．A higher order boundary treatment for Cartesian-grid method［J］．J．Comput．Phys．，1998，140：259-277．

［5］DADONE A，GROSSMAN B．Ghost-cell method for inviscid two-dimensional flows on cartesian grids［J］．AIAA Journal，2004，42（12）：2499-2507．

［6］KALITZIN G，IACCARINO G．Toward immersed boundary simulation of high Reynolds number flows［R］．In：Annual Research Briefs 2003．Center for Turbulence Research，Stanford University，2003．

［7］JAE-Too L．Development of an efficient viscous approach in a cartesian grid framework and application to rotor-fuselage interaction［D］．Georgia：Georgia Institute of Technology，2006．

［8］CAPIZZANO E．A turbulent wall model for immersed boundary methods［R］．AIAA paper 2010-712，2010．

［9］BLAZEK J．Computational fluiddynamics：principles and applications［M］．Amesterdam：ELSEVIER，2011．

［10］SH AROV D，NAKAHASHI K．Reordering of 3-D hybrid unstructured grids for vectorized LU-SGS Navier-Stokes calculations［R］．AIAA Paper 97-2102，1997．

［11］LIOU M S．A sequel to AUSM：AUSM＋［J］．Journal of Computational Physics，1996，129（2）：364-382．

［12］VENKATAKRISHNAN V．On the accuracy of limiters and convergence to Steady State Solutions［R］．AIAA Paper 93-0880，1993．

［13］MENTER E R．Two-equation eddy-viscosity turbulence models for engineering applications［J］．AIAA Journal，1994，32（8）：1598-1605．

［14］ZEEUW D L．A quadtree-based adaptively-refined cartesian-grid algorithm for solution of the Euler equations［D］．Michigan：U-niversity of Michigan，1993．

［15］NICHOLS R H，NELSON C C．Wall function boundary conditions including heat transfer and compressibility for transport turbulence models［R］．AIAA Paper 2004-582，2004．

［16］SPALDING D B．A single formula for the law of the wall［J］．J．App．Mech．，1961，28（3）：455-458．

［17］胡偶，趙寧，劉劍明，等．基于有限體積格式的自適應(yīng)笛卡爾網(wǎng)格虛擬單元方法及其應(yīng)用［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2011，29（4）：491-495．

［18］COOK P H，MCDONALD M A，EIRMIN M C．Aerofoil RAE 2822-pressuredistributions，and boundary layer and wake measurements［R］．In：Experimental Data Base for Computer Program Assessment．AGARD AR 138，1979．

［19］BASHKIN V A，VAGANOV A V，EGOROV I V，et al．Comparison of calculated and experimentaldata on supersonic flow past a circular cylinder［J］．Eluid Dynamic，2002，37（3）：473-483．

［20］MUNIKRISHNA N，LIOU M S．A Cartesian based body-fitted adaptive grid method for compressible viscous flows［R］．AIAA Paper 2009-1500，2009．

Numerical research of Cartesian based ghost cell method for compressible viscous flows

SHEN Zhiwei1，ZHAO Ning1，HU Ou2
（1．College of Aerospace Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China；2．AVICChina helicopter research anddevelopment institute，Jingdezheng 333000，China）

This work presents a new Cartesian-based ghost cell method for two-dimensional high Reynolds number compressible viscous flows．Based on the six fundamental assumptions used in the law-of-thewall，a wall function-ghost cell method（WE-GCM）isdeveloped to treat turbulent wall boundary conditions．Reference points are employed to compute primitive variables and turbulent properties at ghost cells．Meanwhile，the turbulent variables at the near wall cells and boundary cells are modified by using the wall function model．The turbulent boundary conditions are incorporated into a Reynolds average Navier-Stokes（RANS）finite volume solver that includes the SSTk-ωturbulence model．Einally，the transonic flow past a RAE2822 airfoil and supersonic flow past a circle cylinder are simulated with adaptive Cartesian grid．Good agreement with the experimentaldatas shows the accuracy and efficiency of the presented WE-GCM．

ghost cell method；wall function；adaptive Cartesian grid；turbulent flow

V211.3

Adoi：10.7638／kqdlxxb-2014.0105

0258-1825（2014）06-0748-07

2014-09-16；

2014-10-11

國家自然科學(xué)基金（11102179）；國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃（973計(jì)劃）（2014CB046201）

沈志偉（1988-），男，江蘇靖江，博士研究生，研究方向：計(jì)算流體力學(xué)．E-mail:shenzhiweitc＠163.com

沈志偉，趙寧，胡偶．可壓縮粘性流動(dòng)笛卡爾網(wǎng)格虛擬單元方法研究［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2014，32（6）：748-754．

10．7638／kqdlxxb-2014．0105 SHEN Z W，ZH AO N，HU O．Numerical research of Cartesian based ghost cell method for compressible viscous flows［J］．ACTA Aerodynamica Sinica，2014，32（6）：748-754．