陳浩，華如豪，袁先旭，唐志共，畢林,*

1.中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000 2.中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

網(wǎng)格生成技術(shù)是計(jì)算流體動(dòng)力學(xué)(CFD)不可或缺的重要組成部分。計(jì)算網(wǎng)格的合理設(shè)計(jì)和高質(zhì)量生成是數(shù)值模擬的前提條件，直接影響數(shù)值模擬過(guò)程的穩(wěn)定性和數(shù)值計(jì)算結(jié)果的精準(zhǔn)度，關(guān)乎計(jì)算的成敗。在CFD計(jì)算中，網(wǎng)格生成需要頻繁的人機(jī)交互，是人工工作量最大的部分，且網(wǎng)格生成質(zhì)量嚴(yán)重依賴用戶經(jīng)驗(yàn)，是制約CFD效率的主要瓶頸。美國(guó)Sandia國(guó)家實(shí)驗(yàn)室的研究表明，前處理過(guò)程占據(jù)了數(shù)值模擬實(shí)踐中總用時(shí)的90%，而網(wǎng)格生成相關(guān)的步驟獨(dú)占總用時(shí)的76%。

目前，網(wǎng)格生成技術(shù)主要有3類(lèi)：貼體結(jié)構(gòu)網(wǎng)格、非結(jié)構(gòu)網(wǎng)格和非貼體的笛卡爾網(wǎng)格。相對(duì)前兩者而言，笛卡爾網(wǎng)格在自動(dòng)化程度方面有著天然的優(yōu)勢(shì)，且易于捕捉流場(chǎng)精細(xì)特征，能夠保證網(wǎng)格質(zhì)量。然而，對(duì)于自適應(yīng)笛卡爾網(wǎng)格方法而言，其普及性和推廣度不如傳統(tǒng)的結(jié)構(gòu)和非結(jié)構(gòu)網(wǎng)格，很大程度上是因?yàn)槠浞琴N體特性引起的高精度模擬時(shí)網(wǎng)格規(guī)模過(guò)大，導(dǎo)致計(jì)算量和存儲(chǔ)量的急劇增加，計(jì)算成本可能難以接受，這在當(dāng)前的計(jì)算機(jī)硬件條件下仍然是不容忽視的。例如，對(duì)于工程型號(hào)問(wèn)題中遇到的復(fù)雜構(gòu)型，往往存在精細(xì)和跨尺度結(jié)構(gòu)，此時(shí)使用自適應(yīng)笛卡爾網(wǎng)格時(shí)，必須通過(guò)足夠次數(shù)的自適應(yīng)加密，高保真地刻畫(huà)模型邊界信息，通常會(huì)比貼體類(lèi)網(wǎng)格花費(fèi)數(shù)倍的網(wǎng)格量，甚至有量級(jí)上的差距。如果復(fù)雜外形存在運(yùn)動(dòng)或變形，對(duì)于笛卡爾網(wǎng)格帶來(lái)的計(jì)算量過(guò)大的問(wèn)題就更加凸顯，計(jì)算周期會(huì)很長(zhǎng)甚至難以接受。

因此，對(duì)于自適應(yīng)笛卡爾網(wǎng)格方法而言，在實(shí)現(xiàn)工程實(shí)用化的前進(jìn)道路上，還需要聚焦到上述技術(shù)瓶頸，重點(diǎn)解決網(wǎng)格規(guī)模控制、網(wǎng)格生成高效算法設(shè)計(jì)、非貼體壁面高保真處理等關(guān)鍵問(wèn)題，以實(shí)現(xiàn)大規(guī)模自適應(yīng)笛卡爾網(wǎng)格完全自動(dòng)化高效生成和復(fù)雜外形/流場(chǎng)高效高精度地計(jì)算，為縮短CFD計(jì)算周期，實(shí)現(xiàn)工程型號(hào)快速迭代設(shè)計(jì)提供技術(shù)支撐。本文正是以此為出發(fā)點(diǎn)，對(duì)黏性自適應(yīng)笛卡爾網(wǎng)格方法進(jìn)行發(fā)展，并開(kāi)展復(fù)雜構(gòu)型流動(dòng)問(wèn)題的應(yīng)用研究。

本文所選取的應(yīng)用對(duì)象為小展弦比飛翼布局飛機(jī)，其采用翼身融合的全翼式設(shè)計(jì)，取消了平尾和垂尾，整體上具備較好的氣動(dòng)特性和隱身特性，是下一代超聲速高性能戰(zhàn)機(jī)的理想構(gòu)型。由于其氣動(dòng)布局和常規(guī)戰(zhàn)斗機(jī)存在較大差異，目前對(duì)其氣動(dòng)特性的認(rèn)知仍然不足，需要進(jìn)一步開(kāi)展研究。

由于小展弦比飛翼飛機(jī)構(gòu)型的特殊性，存在較多的尖點(diǎn)結(jié)構(gòu)，會(huì)給貼體網(wǎng)格的高質(zhì)量生成帶來(lái)困難。而采用自適應(yīng)笛卡爾網(wǎng)格方法可自動(dòng)化生成高質(zhì)量的空間網(wǎng)格，避免了上述問(wèn)題并能對(duì)空間流動(dòng)特征進(jìn)行精細(xì)化模擬。目前，從公開(kāi)發(fā)表的文獻(xiàn)來(lái)看，尚未見(jiàn)基于黏性自適應(yīng)笛卡爾網(wǎng)格的小展弦比飛翼飛機(jī)流動(dòng)數(shù)值模擬。因此，本文選取該構(gòu)型進(jìn)行模擬研究，考核本文的黏性自適應(yīng)笛卡爾網(wǎng)格方法，并深入探究其氣動(dòng)力特性。

1 自適應(yīng)笛卡爾網(wǎng)格生成技術(shù)

自適應(yīng)笛卡爾網(wǎng)格生成的基本流程如圖1所示，圖中：為自適應(yīng)加密時(shí)間長(zhǎng);Δ為時(shí)間步長(zhǎng)。

1) 讀取STL(STereoLithography)格式的幾何模型數(shù)據(jù)文件，并構(gòu)造物面離散單元的二叉樹(shù)，方便快速檢索。

2) 設(shè)定計(jì)算域，劃分初始網(wǎng)格。

3) 構(gòu)建笛卡爾網(wǎng)格數(shù)據(jù)結(jié)構(gòu)，建立網(wǎng)格信息存儲(chǔ)和鏈接關(guān)系(鄰居關(guān)系、父子關(guān)系、兄弟關(guān)系)。

4) 判斷網(wǎng)格單元類(lèi)型(物面內(nèi)網(wǎng)格單元、物面外網(wǎng)格單元以及物面相交網(wǎng)格單元)。

5) 對(duì)于物面相交單元加密一定的次數(shù)，直到達(dá)到模型保真度的需求。

6) 對(duì)于加密過(guò)程中新出現(xiàn)的網(wǎng)格單元，重復(fù)4)和5)。

7) 對(duì)于通過(guò)上述加密得到的網(wǎng)格進(jìn)行光順優(yōu)化。

8) 進(jìn)行初步流場(chǎng)計(jì)算，在進(jìn)行到一定時(shí)間步后，根據(jù)流場(chǎng)解特征進(jìn)行進(jìn)一步自適應(yīng)，以使空間網(wǎng)格分布更合理。

圖1 自適應(yīng)笛卡爾網(wǎng)格生成流程Fig.1 Flowchart of adaptive Cartesian mesh generation procedure

在上述步驟中，網(wǎng)格數(shù)據(jù)結(jié)構(gòu)尤其關(guān)鍵，直接決定了網(wǎng)格信息查詢效率；網(wǎng)格單元類(lèi)型判斷方法則對(duì)網(wǎng)格生成的準(zhǔn)確度有著重要影響；幾何自適應(yīng)和流場(chǎng)自適應(yīng)技術(shù)是本文網(wǎng)格自適應(yīng)加密和粗化的依據(jù)；此外，對(duì)于非貼體物面邊界的處理是構(gòu)建笛卡爾網(wǎng)格數(shù)值求解器的重要前提。下面，分別對(duì)上述技術(shù)方法展開(kāi)介紹。

1.1 網(wǎng)格數(shù)據(jù)結(jié)構(gòu)

叉樹(shù)結(jié)構(gòu)是自適應(yīng)笛卡爾網(wǎng)格使用最多的數(shù)據(jù)結(jié)構(gòu)形式，其易于自適應(yīng)的特點(diǎn)與笛卡爾網(wǎng)格天然契合。然而，由于笛卡爾網(wǎng)格自適應(yīng)后是分層式的組織框架，不同層的網(wǎng)格之間的數(shù)據(jù)交互往往需要遞歸查詢，甚至遍歷整個(gè)樹(shù)結(jié)構(gòu)，導(dǎo)致計(jì)算量較大，給網(wǎng)格生成和流場(chǎng)計(jì)算的高效性帶來(lái)阻礙。為此，發(fā)展了一種全線程樹(shù)數(shù)據(jù)結(jié)構(gòu)(Fully Threaded Tree，F(xiàn)TT)。

FTT最初是由Khokhlov提出，該方法的基本思想是：在傳統(tǒng)四叉樹(shù)或八叉樹(shù)的基礎(chǔ)上，將其閑置的葉子節(jié)點(diǎn)利用起來(lái)，存儲(chǔ)指向鄰居單元的父節(jié)點(diǎn)的指針，以建立它將叉樹(shù)數(shù)據(jù)結(jié)構(gòu)中葉子節(jié)點(diǎn)閑置的指針利用起來(lái)，指向鄰居單元的父單元，從而構(gòu)建鄰居單元之間快速交互的“線程”。這種方法不僅可以提高網(wǎng)格生成、信息調(diào)用的效率，還能夠提高網(wǎng)格存儲(chǔ)的利用率。

為了進(jìn)一步減少網(wǎng)格類(lèi)型判斷過(guò)程中相交信息的重復(fù)計(jì)算，提高網(wǎng)格生成的效率，發(fā)展了一種新型的FTT數(shù)據(jù)結(jié)構(gòu)。該方法的思路是：在傳統(tǒng)FTT數(shù)據(jù)結(jié)構(gòu)的基礎(chǔ)上，對(duì)于網(wǎng)格相交面元序號(hào)，進(jìn)行動(dòng)態(tài)存儲(chǔ)。事實(shí)上，父子單元相交面元間存在“繼承性”，即與子單元相交的面元肯定與父單元相交，基于這個(gè)原則，就可以把相交判斷過(guò)程中面元檢索范圍限制在很小的范圍，大幅提高網(wǎng)格相交判斷的效率，進(jìn)而提高整體的網(wǎng)格生成效率。

基于上述方式，所構(gòu)建的笛卡爾網(wǎng)格數(shù)據(jù)結(jié)構(gòu)如圖2所示，Oct和Cell之間的組織關(guān)系如圖3所示?？梢钥吹?，除了傳統(tǒng)FTT基本的Oct模塊(包含層次OctLv、中心坐標(biāo)(,,)、父指針OctPr、各方向鄰居的父單元OctNb(6)等信息)，還添加了動(dòng)態(tài)數(shù)組OctTC(：)，用于存儲(chǔ)與笛卡爾網(wǎng)格相交面元的序號(hào)。OctCh代表存儲(chǔ)的子單元；代表流動(dòng)矢量。當(dāng)存在相交面元時(shí)，該數(shù)組被分配內(nèi)存和賦值；在完成子單元的計(jì)算后，父單元的存儲(chǔ)數(shù)組就會(huì)被釋放內(nèi)存，以減少重復(fù)存儲(chǔ)。

圖2 FTT數(shù)據(jù)結(jié)構(gòu)Fig.2 FTT data structure

圖3 FTT數(shù)據(jù)結(jié)構(gòu)中Cell和Oct的關(guān)系[14]Fig.3 Relationship between Cell and Oct in FTT[14]

1.2 射線相交判斷

對(duì)于笛卡爾網(wǎng)格單元類(lèi)型的判斷，最常用的方法是射線相交方法。該方法的基本思路是：通過(guò)計(jì)算過(guò)目標(biāo)點(diǎn)的射線與封閉物面的相交點(diǎn)數(shù)目，來(lái)確定目標(biāo)點(diǎn)在物面的內(nèi)外(奇數(shù)時(shí)在內(nèi)部，偶數(shù)時(shí)在外部)。這種方法簡(jiǎn)單高效，但是存在特殊情況，如射線經(jīng)過(guò)幾何拐點(diǎn)時(shí)，若還按照上述判斷方法，就會(huì)出現(xiàn)錯(cuò)誤。因此，有必要準(zhǔn)確識(shí)別射線經(jīng)過(guò)拐點(diǎn)的情況，如圖4所示。圖中：、代表正常的射線；1、2代表經(jīng)過(guò)拐點(diǎn)的錯(cuò)誤射線。

圖4 基于射線相交方法判斷目標(biāo)點(diǎn)位置[14]Fig.4 Illustration of point inside or outside polygon using ray-casting method[14]

針對(duì)上述問(wèn)題，發(fā)展了一種奇點(diǎn)檢測(cè)算法，可以有效識(shí)別上述情形。將其與射線相交方法結(jié)合，可保證網(wǎng)格點(diǎn)類(lèi)型判斷的準(zhǔn)確性。奇異性檢測(cè)算法的思路是：通過(guò)射線方向與拐點(diǎn)位置幾何離散單元法向之間夾角的關(guān)系，來(lái)識(shí)別射線對(duì)應(yīng)的目標(biāo)點(diǎn)在幾何內(nèi)外。

以圖5為例，目標(biāo)點(diǎn)在封閉幾何外部所滿足的判據(jù)為

cos〈,〉·cos〈,〉>0

or cos〈,〉·cos〈,〉>0

〈,〉=〈,〉=>0

〈,〉=〈,〉=>0

目標(biāo)點(diǎn)在封閉幾何內(nèi)部所滿足的判據(jù)為

cos〈,〉·cos〈,〉>0

or cos〈,〉·cos〈,〉>0

〈,〉=〈,〉=>0

〈,〉=〈,〉=>0

上述公式的具體論述過(guò)程參照文獻(xiàn)[14]。本文針對(duì)射線經(jīng)過(guò)封閉幾何拐點(diǎn)的情況，以二維幾何問(wèn)題為例，具體處理方法如下：

1) 在射線與幾何離散單元存在交點(diǎn)時(shí)，判斷該交點(diǎn)是否為該幾何離散單元的端點(diǎn)或邊緣點(diǎn)，若是，即射線經(jīng)過(guò)封閉幾何拐點(diǎn)。

2) 計(jì)算該射線方向，幾何離散單元法向。

3) 通過(guò)上述判據(jù)公式，結(jié)合圖5，判斷目標(biāo)點(diǎn)在封閉幾何的內(nèi)部還是外部。

圖5 射線過(guò)拐點(diǎn)的情形[14]Fig.5 Illustration of ray passing inflection point[14]

通過(guò)上述步驟，可識(shí)別射線經(jīng)過(guò)封閉幾何拐點(diǎn)對(duì)應(yīng)的目標(biāo)點(diǎn)位置情況，避免傳統(tǒng)射線相交方法因相交次數(shù)計(jì)數(shù)不準(zhǔn)帶來(lái)的誤判。

1.3 網(wǎng)格數(shù)據(jù)結(jié)構(gòu)

在完成網(wǎng)格類(lèi)型的判斷后，需要根據(jù)幾何特征加密一定的次數(shù)，達(dá)到幾何信息保真度需求。具體來(lái)說(shuō)，主要包括相交單元自適應(yīng)和曲率自適應(yīng)，即分別對(duì)于與物面相交的網(wǎng)格單元和曲率變化較大位置的網(wǎng)格單元加密一定的次數(shù)，一般后者是在前者的基礎(chǔ)上進(jìn)行的。在針對(duì)具體流動(dòng)問(wèn)題進(jìn)行幾何自適應(yīng)加密時(shí)，還會(huì)考慮邊界層網(wǎng)格高度的要求，來(lái)確定具體的加密次數(shù)，以使得最大壁面網(wǎng)格高度滿足約束條件。

對(duì)于笛卡爾網(wǎng)格而言，在加密后，相鄰層次的粗網(wǎng)格與細(xì)網(wǎng)格的尺寸之比值恒定為2，即使在附面層內(nèi)。這種網(wǎng)格尺寸的大比例變化，會(huì)影響近壁面流動(dòng)模擬的精度和穩(wěn)定性。針對(duì)上述問(wèn)題，本文借鑒貼體結(jié)構(gòu)網(wǎng)格對(duì)于邊界層流動(dòng)特征捕捉的優(yōu)勢(shì)，引入虛擬層技術(shù)進(jìn)行近壁面網(wǎng)格的優(yōu)化，如圖6所示。首先根據(jù)物面外形數(shù)據(jù)點(diǎn)沿壁面法向射線方向按照一定規(guī)律在邊界層內(nèi)布置一系列點(diǎn)、、、…，由第一層數(shù)據(jù)點(diǎn)生成初始虛擬面，然后設(shè)定虛擬物面沿壁面法向距離層進(jìn)比例：

式中：為第層網(wǎng)格的高度。該比例的設(shè)定借鑒貼體結(jié)構(gòu)網(wǎng)格在邊界層內(nèi)的網(wǎng)格分布特點(diǎn)，得到邊界層內(nèi)多個(gè)距離變化均勻的虛擬面，其距離排布類(lèi)似貼體結(jié)構(gòu)網(wǎng)格。笛卡爾網(wǎng)格單元的加密次數(shù)和層級(jí)就結(jié)合這些虛擬物面確定。

圖6 虛擬層方法示意圖Fig.6 Illustration of virtual layer method

1.4 流場(chǎng)自適應(yīng)

為了使網(wǎng)格分布更合理，以節(jié)省網(wǎng)格量和精細(xì)捕捉流動(dòng)結(jié)構(gòu)，需要根據(jù)流場(chǎng)解的特征進(jìn)行進(jìn)一步的自適應(yīng)加密或粗化。

初始網(wǎng)格的分辨率往往不足以精細(xì)刻畫(huà)流動(dòng)特征，合理有效的流場(chǎng)解自適應(yīng)判據(jù)是開(kāi)展流動(dòng)自適應(yīng)的前提。根據(jù)不同的流動(dòng)特征，往往需要建立對(duì)應(yīng)的特征識(shí)別準(zhǔn)則。一般來(lái)說(shuō)，流動(dòng)的間斷區(qū)域，如激波，需要速度或壓力的梯度進(jìn)行識(shí)別；而渦結(jié)構(gòu)集中的區(qū)域，用渦量、速度的旋度等識(shí)別效果更好。針對(duì)本文的研究對(duì)象，即低速不可壓縮情況下的飛翼流動(dòng)模擬，為渦旋結(jié)構(gòu)主導(dǎo)，不存在激波等間斷特征，因此，本文選用的是以速度旋度作為流動(dòng)自適應(yīng)判據(jù)。

在笛卡爾網(wǎng)格框架下，單元cell的速度旋度可表示如下：

式中：為速度；為網(wǎng)格單元cell的邊長(zhǎng)，三維情況下=3。

相應(yīng)的速度的旋度的標(biāo)準(zhǔn)差為

式中：為計(jì)算域內(nèi)所有的網(wǎng)格單元總數(shù)；為單元cell的速度的旋度的標(biāo)準(zhǔn)差。具體判據(jù)如下：

式中：、為經(jīng)驗(yàn)系數(shù)，需要針對(duì)具體的流動(dòng)問(wèn)題給出。

2 物面邊界條件

基于浸入邊界方法處理非貼體笛卡爾網(wǎng)格物面邊界的基本思路是：

1) 首先，確定判斷物面附近計(jì)算模板點(diǎn)是否在物面內(nèi)部，并對(duì)內(nèi)部的計(jì)算模板點(diǎn)標(biāo)記為虛擬單元(Ghost Cell)。

2) 對(duì)于虛擬單元關(guān)于物面邊界取參考點(diǎn)。

3) 通過(guò)臨近單元插值方法獲取參考點(diǎn)的流場(chǎng)值。

4) 通過(guò)物面無(wú)滑移、無(wú)穿透條件建立虛擬點(diǎn)和參考點(diǎn)的對(duì)應(yīng)關(guān)系。

5) 考慮物面曲率的影響，進(jìn)行曲率修正。

以圖7為例，虛擬點(diǎn)在流場(chǎng)中的參考點(diǎn)為，為到壁面投影點(diǎn)的距離。對(duì)于參考點(diǎn)的取法，傳統(tǒng)方式一般是取關(guān)于物面的鏡像對(duì)稱點(diǎn)。這種方法雖然簡(jiǎn)單，但是在虛擬點(diǎn)很靠近物面的時(shí)候，參考點(diǎn)也會(huì)離物面很近，導(dǎo)致其插值用到的臨近單元會(huì)有部分在物體內(nèi)部，這會(huì)對(duì)其流場(chǎng)值得計(jì)算帶來(lái)誤差。本文的方式是借鑒Forrer和Jeltsch的取法，將參考點(diǎn)沿對(duì)稱線延長(zhǎng)一個(gè)當(dāng)?shù)鼐W(wǎng)的長(zhǎng)度，確保臨近插值單元均在

圖7 非貼體物面處理方法Fig.7 Processing method of non-body-fitted wall

流場(chǎng)之中。Forrer和Jeltsch的研究表明，這種位置取法可以提高邊界附近計(jì)算的穩(wěn)定性，而且并不會(huì)降低非貼體物面模擬的精度。

參考點(diǎn)流場(chǎng)值通過(guò)鄰近單元插值獲取，以參考點(diǎn)為例，其鄰近的4個(gè)單元結(jié)點(diǎn)分別記作、、、，則其流場(chǎng)值為

=

式中：、、、分別為點(diǎn)到網(wǎng)格點(diǎn)、、、的距離。

虛擬點(diǎn)和參考點(diǎn)的關(guān)系式如下：

式中：為氣體常數(shù)。

3 數(shù)值方法

采用三維守恒型Navier-Stokes方程作為可壓縮黏性流動(dòng)的控制方程，即

式中：和分別為來(lái)流馬赫數(shù)和雷諾數(shù)；為守恒變量；、、代表無(wú)黏通量；、、代表黏性通量。

對(duì)于Navier-Stokes方程的數(shù)值離散方法是通過(guò)張涵信院士建立的NND格式進(jìn)行無(wú)黏項(xiàng)的離散，以及通過(guò)中心差分方法進(jìn)行黏性項(xiàng)的離散，在避免出現(xiàn)非物理振蕩的同時(shí)，可實(shí)現(xiàn)較高精度的捕捉激波、剪切層、接觸間斷等流場(chǎng)結(jié)構(gòu)及現(xiàn)象。

以方向上的變量為例(其他方向的形式類(lèi)似)，采用NND格式的無(wú)黏通量空間離散如下：

()=+12--12

對(duì)于時(shí)間的推進(jìn)使用的是顯式推進(jìn)方法。顯式方法不需要進(jìn)行大量的矩陣運(yùn)算，存儲(chǔ)量也較小，同時(shí)方法的建立和程序?qū)崿F(xiàn)比較簡(jiǎn)單。為了保證所構(gòu)建的全離散格式的TVD性質(zhì)，本文選取的是具有TVD性質(zhì)的Runge-Kutta型二階時(shí)間離散格式，如下：

對(duì)于自適應(yīng)笛卡爾網(wǎng)格方法而言，在構(gòu)建數(shù)值離散框架時(shí)，需要考慮網(wǎng)格自適應(yīng)帶來(lái)的非均勻、粗細(xì)過(guò)渡的特征，即懸掛網(wǎng)格的處理。其思路是通過(guò)插值重構(gòu)的方式，構(gòu)造需要的計(jì)算模板點(diǎn)，具體可參考文獻(xiàn)[27]，這里不再贅述。

4 結(jié)果和討論

本節(jié)選取小展弦比飛翼布局飛行器這一復(fù)雜構(gòu)型飛行器，運(yùn)用本文發(fā)展的自適應(yīng)笛卡爾網(wǎng)格技術(shù)開(kāi)展模擬研究。該飛行器采用隱身布局，同時(shí)氣動(dòng)性能較好，在作戰(zhàn)飛機(jī)設(shè)計(jì)方面有著廣泛的應(yīng)用前景。該幾何標(biāo)模的外形如圖8所示，經(jīng)由國(guó)內(nèi)多個(gè)研究機(jī)構(gòu)聯(lián)合開(kāi)展了數(shù)值模擬和風(fēng)洞試驗(yàn)。本文采用該標(biāo)模來(lái)考核驗(yàn)證本文的笛卡網(wǎng)格生成技術(shù)和數(shù)值求解器，并開(kāi)展笛卡爾網(wǎng)格自適應(yīng)技術(shù)對(duì)于小展弦比飛翼布局飛行器氣動(dòng)力特性影響規(guī)律的研究。

圖8 小展弦比飛翼布局風(fēng)洞模型Fig.8 Wind tunnel model of low aspect ratio fly-wing aircraft

來(lái)流馬赫數(shù)為0.2，靜溫和靜壓分別為288.15 K 和101 325 Pa，參考實(shí)驗(yàn)設(shè)計(jì)，選取來(lái)流攻角分別為=20°、24°、28°、32°、36°、38°、40°、42°、44°、48°、52°，無(wú)側(cè)滑角。初始幾何自適應(yīng)網(wǎng)格如圖9所示，其中物面網(wǎng)格49 026， 為三角形離散單元，經(jīng)幾何自適應(yīng)加密后，空間笛卡爾網(wǎng)格數(shù)目為2 259 042。

圖9 幾何自適應(yīng)笛卡爾網(wǎng)格Fig.9 Geometry adaptive Cartesian mesh

圖10和圖11為流場(chǎng)計(jì)算結(jié)果的流動(dòng)自適應(yīng)網(wǎng)格和壓力分布云圖、流場(chǎng)流線?？梢钥闯?，本文的自適應(yīng)笛卡爾網(wǎng)格技術(shù)能夠?qū)τ诳臻g流場(chǎng)的特征結(jié)構(gòu)進(jìn)行很好的捕捉，并進(jìn)行自適應(yīng)加密，確保重要流動(dòng)結(jié)構(gòu)的高精度模擬。網(wǎng)格自適應(yīng)區(qū)域主要集中在翼尖和尾尖牽引的尾流區(qū)，同時(shí)對(duì)于大攻角分離區(qū)的流場(chǎng)梯度變化也能很敏銳的捕捉。通過(guò)截面圖，可以看到在背風(fēng)區(qū)分離渦也有很好的自適應(yīng)效果。相對(duì)28°攻角、40°攻角時(shí)在網(wǎng)格分布上展現(xiàn)出了更復(fù)雜的分布特性，考慮是由于大攻角流動(dòng)分離引起的渦結(jié)構(gòu)演化導(dǎo)致的。

圖10 來(lái)流攻角為28°和40°時(shí)的壓力和網(wǎng)格分布Fig.10 Pressure and mesh distribution at α=28°，40°

圖11 來(lái)流攻角為28°時(shí)流線和渦結(jié)構(gòu) 自適應(yīng)網(wǎng)格示意圖Fig.11 Stream lines and adaptive mesh around votex at α=28°

圖12為自適應(yīng)笛卡爾網(wǎng)格方法得到的不同攻角下的升力系數(shù)，并與不同風(fēng)洞試驗(yàn)數(shù)據(jù)的對(duì)比，呈現(xiàn)了較好的吻合度?？梢钥闯?，F(xiàn)D9、FL8、FL12這3組試驗(yàn)結(jié)果給出的失速攻角在40°附近，而本文計(jì)算結(jié)果在38°附近。在計(jì)算結(jié)果中，升力系數(shù)在40°時(shí)明顯下降，通過(guò)對(duì)比兩者的流場(chǎng)云圖可以看出，在40°攻角時(shí)，由于大攻角分離，分離再附區(qū)減小，渦結(jié)構(gòu)演化時(shí)橫向流動(dòng)變強(qiáng)，使得飛行器背風(fēng)區(qū)的壓力降低，從而導(dǎo)致升力系數(shù)減小。

圖12 升力系數(shù)的計(jì)算結(jié)果與風(fēng)洞試驗(yàn)結(jié)果對(duì)比Fig.12 Lift coefficient comparison of calculated and experimental results

對(duì)于自適應(yīng)前后流場(chǎng)特性和氣動(dòng)力特性的對(duì)比研究，從圖13和圖14中攻角為28°和40°上表面的壓力分布在自適應(yīng)前后的變化可以看出，在自適應(yīng)后，攻角為28°的狀態(tài)下，飛行器背風(fēng)區(qū)上表面的低壓區(qū)更大，這將導(dǎo)致更高的升力系數(shù)。而在40°攻角狀態(tài)下，背風(fēng)區(qū)壓力分布較復(fù)雜，上述規(guī)律并不明顯。自適應(yīng)對(duì)背風(fēng)區(qū)渦結(jié)構(gòu)有著更高精度的捕捉，降低了分離渦的數(shù)值耗散，呈現(xiàn)了更大的低壓背風(fēng)區(qū)，使得升力系數(shù)更高。

圖13 不同攻角下自適應(yīng)前后壓力系數(shù)對(duì)比Fig.13 Pressure coefficient comparison before and after adaptation at different angles of attack

圖14 自適應(yīng)前后渦結(jié)構(gòu)附近網(wǎng)格和流場(chǎng)云圖對(duì)比結(jié)果Fig.14 Mesh and pressure contour comparison before and after adaptation

下面結(jié)合渦量圖進(jìn)一步對(duì)比網(wǎng)格自適應(yīng)前后對(duì)于流動(dòng)結(jié)構(gòu)的捕捉能力。從圖15中可以清晰看出，在不添加網(wǎng)格自適應(yīng)時(shí)，背風(fēng)區(qū)的流動(dòng)分辨率很低，無(wú)法清晰辨識(shí)出渦結(jié)構(gòu)，只能看出大致輪廓。并且在尾部離物面較遠(yuǎn)處，由于笛卡爾網(wǎng)格粗細(xì)不均，導(dǎo)致出現(xiàn)了壓力云圖的間斷。而在引入網(wǎng)格自適應(yīng)后，對(duì)于背風(fēng)區(qū)分離流動(dòng)的精細(xì)渦結(jié)構(gòu)都可以較清晰的捕捉和呈現(xiàn)。本文中采用了NND格式和中心差分格式進(jìn)行空間離散，數(shù)值格式耗散較小，而自適應(yīng)技術(shù)的引入，又能夠充分保證流動(dòng)細(xì)節(jié)較多的位置的網(wǎng)格密度(見(jiàn)圖10，圖14～圖16)，從而能夠較高分辨率地捕捉流場(chǎng)渦結(jié)構(gòu)。此外，還可以看出不同來(lái)流攻角下的背風(fēng)區(qū)渦結(jié)構(gòu)演化特征。在來(lái)流攻角為28°時(shí)，能看到完整的渦條結(jié)構(gòu)，在靠近尾部時(shí)才破碎；在來(lái)流攻角為40°時(shí)，背風(fēng)區(qū)流場(chǎng)基本處于完全破碎的大分離狀態(tài)，僅在靠近頭部有較短、較細(xì)的渦條結(jié)構(gòu)。

圖15 自適應(yīng)前后渦量圖對(duì)比(α=28°)Fig.15 Comparison of vorticity contour before and after adaptation (α=28°)

圖16 自適應(yīng)前后渦量圖對(duì)比(α=40°)Fig.16 Comparison of vorticity contour before and after adaptation (α=40°)

圖17為自適應(yīng)前后的升力系數(shù)、阻力系數(shù)以及俯仰力矩系數(shù)的對(duì)比數(shù)據(jù)曲線?？梢钥闯觯啾戎?，自適應(yīng)加密后與試驗(yàn)結(jié)果吻合更好。整體上，自適應(yīng)加密會(huì)使計(jì)算結(jié)果值變大，這是因?yàn)閷?duì)于背風(fēng)區(qū)流動(dòng)結(jié)構(gòu)捕捉更精細(xì)，使得飛行器上表面的低壓區(qū)較大，結(jié)合前面的網(wǎng)格自適應(yīng)圖和壓力云圖可以清晰看出。在=38°以下和=40°以上，兩者差別相對(duì)較小，并且與試驗(yàn)結(jié)果也較為接近；在=38°和=40°附近，即接近失速攻角時(shí)，差別較為明顯。說(shuō)明在低速大攻角狀態(tài)下，空間流動(dòng)結(jié)構(gòu)模擬的精細(xì)化程度，對(duì)飛行器的氣動(dòng)力特性計(jì)算的精準(zhǔn)度影響相對(duì)較大。尤其是對(duì)于俯仰力矩系數(shù)，在自適應(yīng)后與試驗(yàn)結(jié)果吻合度很好，尤其是在攻角較小時(shí)，縱向靜穩(wěn)定度與試驗(yàn)基本相符。

圖17 不同來(lái)流攻角的氣動(dòng)力系數(shù)計(jì)算與試驗(yàn)結(jié)果對(duì)比Fig.17 Aerodynamic force coefficient comparison between calculated and experimental data at different angles of attach

圖18中選取攻角為40°狀態(tài)下，在小展弦比飛翼不同展向的截面位置上的壓力分布，并通過(guò)自適應(yīng)加密前后的數(shù)據(jù)進(jìn)行對(duì)比。由圖可以看出，對(duì)于飛翼模型的下表面，即迎風(fēng)面，表面壓力系數(shù)分布在自適應(yīng)前后幾乎無(wú)變化，這是因?yàn)樵谟L(fēng)的高壓面，流動(dòng)梯度變化并不大，自適應(yīng)起到的效果不明顯。而對(duì)于流動(dòng)梯度變化較大的背風(fēng)區(qū)，流場(chǎng)渦結(jié)構(gòu)較為復(fù)雜，是網(wǎng)格自適應(yīng)集中加密捕捉的區(qū)域。這使得自適應(yīng)后飛翼的上表面壓力相對(duì)更低，從而導(dǎo)致升力系數(shù)變大，俯仰力矩變小。這與前面的數(shù)據(jù)對(duì)比結(jié)果相一致，并與圖13中壓力分布云圖呈現(xiàn)出的規(guī)律相同。

圖18 攻角為40°時(shí)展向不同截面位置的壓力系數(shù)對(duì)比Fig.18 Pressure coefficient comparison of different spanwise locations at α=40°

5 結(jié) 論

在本文的工作中，發(fā)展了三維黏性自適應(yīng)笛卡爾網(wǎng)格方法，并開(kāi)展了小展弦比飛翼飛機(jī)低速流動(dòng)問(wèn)題的應(yīng)用研究。主要包括：

1) 通過(guò)網(wǎng)格數(shù)據(jù)的優(yōu)化和射線相交方法的改進(jìn)，以及近壁網(wǎng)格過(guò)渡區(qū)的優(yōu)化，保證了自適應(yīng)笛卡爾網(wǎng)格生成的高效性、魯棒性和優(yōu)質(zhì)性，提高了其對(duì)于復(fù)雜外形流動(dòng)問(wèn)題的適用性。

2) 發(fā)展了非貼體笛卡爾網(wǎng)格下黏性物面邊界處理方法，構(gòu)造了高保真的非貼體壁面邊界模擬方法。在此基礎(chǔ)上，構(gòu)建了自適應(yīng)笛卡爾網(wǎng)格下Navier-Stokes方程數(shù)值離散框架，建立了黏性流動(dòng)數(shù)值求解器。

3) 基于上述技術(shù)方法，開(kāi)展了自適應(yīng)笛卡爾網(wǎng)格方法在飛翼流動(dòng)問(wèn)題中的應(yīng)用研究。研究表明，本文的技術(shù)方法能夠?qū)崿F(xiàn)對(duì)于小展弦比飛翼大攻角流動(dòng)問(wèn)題自動(dòng)化、精細(xì)化模擬仿真，并且自適應(yīng)技術(shù)的引入，可有效改善其氣動(dòng)力預(yù)測(cè)精度。

黏性自適應(yīng)笛卡爾網(wǎng)格方法為飛行器自動(dòng)化設(shè)計(jì)和預(yù)測(cè)提供了一種可選手段，并且其精度可靠，具有廣闊的應(yīng)用前景。在后續(xù)工作中，將對(duì)網(wǎng)格技術(shù)進(jìn)一步優(yōu)化，并將構(gòu)建笛卡爾網(wǎng)格下的非定常數(shù)值模擬方法，為開(kāi)展飛行器飛行動(dòng)力學(xué)仿真、深入研究其氣動(dòng)和飛行特性奠定基礎(chǔ)。