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

1.中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000 2.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

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

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

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

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

由于小展弦比飛翼飛機構型的特殊性，存在較多的尖點結構，會給貼體網格的高質量生成帶來困難。而采用自適應笛卡爾網格方法可自動化生成高質量的空間網格，避免了上述問題并能對空間流動特征進行精細化模擬。目前，從公開發(fā)表的文獻來看，尚未見基于黏性自適應笛卡爾網格的小展弦比飛翼飛機流動數值模擬。因此，本文選取該構型進行模擬研究，考核本文的黏性自適應笛卡爾網格方法，并深入探究其氣動力特性。

1 自適應笛卡爾網格生成技術

自適應笛卡爾網格生成的基本流程如圖1所示，圖中：為自適應加密時間長;Δ為時間步長。

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

2) 設定計算域，劃分初始網格。

3) 構建笛卡爾網格數據結構，建立網格信息存儲和鏈接關系(鄰居關系、父子關系、兄弟關系)。

4) 判斷網格單元類型(物面內網格單元、物面外網格單元以及物面相交網格單元)。

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

6) 對于加密過程中新出現的網格單元，重復4)和5)。

7) 對于通過上述加密得到的網格進行光順優(yōu)化。

8) 進行初步流場計算，在進行到一定時間步后，根據流場解特征進行進一步自適應，以使空間網格分布更合理。

圖1 自適應笛卡爾網格生成流程Fig.1 Flowchart of adaptive Cartesian mesh generation procedure

在上述步驟中，網格數據結構尤其關鍵，直接決定了網格信息查詢效率；網格單元類型判斷方法則對網格生成的準確度有著重要影響；幾何自適應和流場自適應技術是本文網格自適應加密和粗化的依據；此外，對于非貼體物面邊界的處理是構建笛卡爾網格數值求解器的重要前提。下面，分別對上述技術方法展開介紹。

1.1 網格數據結構

叉樹結構是自適應笛卡爾網格使用最多的數據結構形式，其易于自適應的特點與笛卡爾網格天然契合。然而，由于笛卡爾網格自適應后是分層式的組織框架，不同層的網格之間的數據交互往往需要遞歸查詢，甚至遍歷整個樹結構，導致計算量較大，給網格生成和流場計算的高效性帶來阻礙。為此，發(fā)展了一種全線程樹數據結構(Fully Threaded Tree，FTT)。

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

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

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

圖2 FTT數據結構Fig.2 FTT data structure

圖3 FTT數據結構中Cell和Oct的關系[14]Fig.3 Relationship between Cell and Oct in FTT[14]

1.2 射線相交判斷

對于笛卡爾網格單元類型的判斷，最常用的方法是射線相交方法。該方法的基本思路是：通過計算過目標點的射線與封閉物面的相交點數目，來確定目標點在物面的內外(奇數時在內部，偶數時在外部)。這種方法簡單高效，但是存在特殊情況，如射線經過幾何拐點時，若還按照上述判斷方法，就會出現錯誤。因此，有必要準確識別射線經過拐點的情況，如圖4所示。圖中：、代表正常的射線；1、2代表經過拐點的錯誤射線。

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

針對上述問題，發(fā)展了一種奇點檢測算法，可以有效識別上述情形。將其與射線相交方法結合，可保證網格點類型判斷的準確性。奇異性檢測算法的思路是：通過射線方向與拐點位置幾何離散單元法向之間夾角的關系，來識別射線對應的目標點在幾何內外。

以圖5為例，目標點在封閉幾何外部所滿足的判據為

cos〈,〉·cos〈,〉>0

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

〈,〉=〈,〉=>0

〈,〉=〈,〉=>0

目標點在封閉幾何內部所滿足的判據為

cos〈,〉·cos〈,〉>0

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

〈,〉=〈,〉=>0

〈,〉=〈,〉=>0

上述公式的具體論述過程參照文獻[14]。本文針對射線經過封閉幾何拐點的情況，以二維幾何問題為例，具體處理方法如下：

1) 在射線與幾何離散單元存在交點時，判斷該交點是否為該幾何離散單元的端點或邊緣點，若是，即射線經過封閉幾何拐點。

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

3) 通過上述判據公式，結合圖5，判斷目標點在封閉幾何的內部還是外部。

圖5 射線過拐點的情形[14]Fig.5 Illustration of ray passing inflection point[14]

通過上述步驟，可識別射線經過封閉幾何拐點對應的目標點位置情況，避免傳統射線相交方法因相交次數計數不準帶來的誤判。

1.3 網格數據結構

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

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

式中：為第層網格的高度。該比例的設定借鑒貼體結構網格在邊界層內的網格分布特點，得到邊界層內多個距離變化均勻的虛擬面，其距離排布類似貼體結構網格。笛卡爾網格單元的加密次數和層級就結合這些虛擬物面確定。

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

1.4 流場自適應

為了使網格分布更合理，以節(jié)省網格量和精細捕捉流動結構，需要根據流場解的特征進行進一步的自適應加密或粗化。

初始網格的分辨率往往不足以精細刻畫流動特征，合理有效的流場解自適應判據是開展流動自適應的前提。根據不同的流動特征，往往需要建立對應的特征識別準則。一般來說，流動的間斷區(qū)域，如激波，需要速度或壓力的梯度進行識別；而渦結構集中的區(qū)域，用渦量、速度的旋度等識別效果更好。針對本文的研究對象，即低速不可壓縮情況下的飛翼流動模擬，為渦旋結構主導，不存在激波等間斷特征，因此，本文選用的是以速度旋度作為流動自適應判據。

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

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

相應的速度的旋度的標準差為

式中：為計算域內所有的網格單元總數；為單元cell的速度的旋度的標準差。具體判據如下：

式中：、為經驗系數，需要針對具體的流動問題給出。

2 物面邊界條件

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

1) 首先，確定判斷物面附近計算模板點是否在物面內部，并對內部的計算模板點標記為虛擬單元(Ghost Cell)。

2) 對于虛擬單元關于物面邊界取參考點。

3) 通過臨近單元插值方法獲取參考點的流場值。

4) 通過物面無滑移、無穿透條件建立虛擬點和參考點的對應關系。

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

以圖7為例，虛擬點在流場中的參考點為，為到壁面投影點的距離。對于參考點的取法，傳統方式一般是取關于物面的鏡像對稱點。這種方法雖然簡單，但是在虛擬點很靠近物面的時候，參考點也會離物面很近，導致其插值用到的臨近單元會有部分在物體內部，這會對其流場值得計算帶來誤差。本文的方式是借鑒Forrer和Jeltsch的取法，將參考點沿對稱線延長一個當地網的長度，確保臨近插值單元均在

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

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

參考點流場值通過鄰近單元插值獲取，以參考點為例，其鄰近的4個單元結點分別記作、、、，則其流場值為

=

式中：、、、分別為點到網格點、、、的距離。

虛擬點和參考點的關系式如下：

式中：為氣體常數。

3 數值方法

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

式中：和分別為來流馬赫數和雷諾數；為守恒變量；、、代表無黏通量；、、代表黏性通量。

對于Navier-Stokes方程的數值離散方法是通過張涵信院士建立的NND格式進行無黏項的離散，以及通過中心差分方法進行黏性項的離散，在避免出現非物理振蕩的同時，可實現較高精度的捕捉激波、剪切層、接觸間斷等流場結構及現象。

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

()=+12--12

對于時間的推進使用的是顯式推進方法。顯式方法不需要進行大量的矩陣運算，存儲量也較小，同時方法的建立和程序實現比較簡單。為了保證所構建的全離散格式的TVD性質，本文選取的是具有TVD性質的Runge-Kutta型二階時間離散格式，如下：

對于自適應笛卡爾網格方法而言，在構建數值離散框架時，需要考慮網格自適應帶來的非均勻、粗細過渡的特征，即懸掛網格的處理。其思路是通過插值重構的方式，構造需要的計算模板點，具體可參考文獻[27]，這里不再贅述。

4 結果和討論

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

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

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

圖9 幾何自適應笛卡爾網格Fig.9 Geometry adaptive Cartesian mesh

圖10和圖11為流場計算結果的流動自適應網格和壓力分布云圖、流場流線。可以看出，本文的自適應笛卡爾網格技術能夠對于空間流場的特征結構進行很好的捕捉，并進行自適應加密，確保重要流動結構的高精度模擬。網格自適應區(qū)域主要集中在翼尖和尾尖牽引的尾流區(qū)，同時對于大攻角分離區(qū)的流場梯度變化也能很敏銳的捕捉。通過截面圖，可以看到在背風區(qū)分離渦也有很好的自適應效果。相對28°攻角、40°攻角時在網格分布上展現出了更復雜的分布特性，考慮是由于大攻角流動分離引起的渦結構演化導致的。

圖10 來流攻角為28°和40°時的壓力和網格分布Fig.10 Pressure and mesh distribution at α=28°，40°

圖11 來流攻角為28°時流線和渦結構 自適應網格示意圖Fig.11 Stream lines and adaptive mesh around votex at α=28°

圖12為自適應笛卡爾網格方法得到的不同攻角下的升力系數，并與不同風洞試驗數據的對比，呈現了較好的吻合度。可以看出，FD9、FL8、FL12這3組試驗結果給出的失速攻角在40°附近，而本文計算結果在38°附近。在計算結果中，升力系數在40°時明顯下降，通過對比兩者的流場云圖可以看出，在40°攻角時，由于大攻角分離，分離再附區(qū)減小，渦結構演化時橫向流動變強，使得飛行器背風區(qū)的壓力降低，從而導致升力系數減小。

圖12 升力系數的計算結果與風洞試驗結果對比Fig.12 Lift coefficient comparison of calculated and experimental results

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

圖13 不同攻角下自適應前后壓力系數對比Fig.13 Pressure coefficient comparison before and after adaptation at different angles of attack

圖14 自適應前后渦結構附近網格和流場云圖對比結果Fig.14 Mesh and pressure contour comparison before and after adaptation

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

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

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

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

圖17 不同來流攻角的氣動力系數計算與試驗結果對比Fig.17 Aerodynamic force coefficient comparison between calculated and experimental data at different angles of attach

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

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

5 結 論

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

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

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

3) 基于上述技術方法，開展了自適應笛卡爾網格方法在飛翼流動問題中的應用研究。研究表明，本文的技術方法能夠實現對于小展弦比飛翼大攻角流動問題自動化、精細化模擬仿真，并且自適應技術的引入，可有效改善其氣動力預測精度。

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