唐志共，陳浩，畢林，袁先旭

中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

隨著計算流體動力學(xué)在工程實際中的深入應(yīng)用，所面臨的幾何外形和流場越來越復(fù)雜，高質(zhì)量的網(wǎng)格生成也變得越來越困難，往往需要占用大量人力資源。在針對復(fù)雜外形的網(wǎng)格生成方法中，相對于傳統(tǒng)的分區(qū)結(jié)構(gòu)網(wǎng)格方法和非結(jié)構(gòu)網(wǎng)格生成方法[1-3]，自適應(yīng)笛卡爾網(wǎng)格[4-6]在自動化生成高質(zhì)量網(wǎng)格方面具有很大優(yōu)勢。此外，自適應(yīng)笛卡爾網(wǎng)格還具備天然的多重網(wǎng)格特性，可以在計算中加速收斂；無需存儲矩陣、運算操作少，有利于高精度算法的應(yīng)用。自20世紀80年代以來，計算機技術(shù)和性能的飛速發(fā)展促使了自適應(yīng)笛卡爾網(wǎng)格技術(shù)的應(yīng)用和推廣，鑒于上述優(yōu)勢，越來越多的CFD工作者投入到該方法技術(shù)的研究和使用之中。

自適應(yīng)笛卡爾網(wǎng)格生成技術(shù)的主要不足之處在于黏性物面邊界處理較復(fù)雜，這是由笛卡爾網(wǎng)格的非貼體特性導(dǎo)致的。笛卡爾網(wǎng)格一般與物面相交，產(chǎn)生不同形狀的切割單元。這些切割單元有的尺寸很小，容易引起在物面附近解的非物理振蕩，如果數(shù)值模擬過程中采用的差分格式是顯式的，還會限制時間步長以及影響計算的穩(wěn)定性。目前主要的處理方法有：混合網(wǎng)格方法(Hybrid Cell Method)[7-8]、切割單元方法(Cut Cell Method)[9-10]、嵌入邊界方法(Embedded Boundary Approach)[11]以及浸入邊界方法(Immersed Boundary Method)[12-13]。

目前，基于自適應(yīng)笛卡爾網(wǎng)格的超聲速黏性流動數(shù)值模擬尚未工程實用化，有進一步發(fā)展的需求。為了實現(xiàn)對二維含強間斷的高超聲速黏性流動問題自動、高效、高精度的模擬，本文發(fā)展了自適應(yīng)笛卡爾網(wǎng)格生成技術(shù)，以及笛卡爾網(wǎng)格下可有效模擬黏性物面邊界條件的方法，并構(gòu)建了相應(yīng)的黏性數(shù)值求解器。

1 網(wǎng)格自適應(yīng)技術(shù)

1.1 四叉樹數(shù)據(jù)結(jié)構(gòu)

在四叉樹數(shù)據(jù)結(jié)構(gòu)中，父親網(wǎng)格單元加密得到4個子網(wǎng)格單元，父單元和子單元通過指針建立聯(lián)系，如圖1所示。這種數(shù)據(jù)結(jié)構(gòu)存儲量較小，自適應(yīng)加密和粗化時方便高效：粗化時，只需將4個子指針置空，設(shè)為無指向即可；加密時，則只要對子指針分別分配內(nèi)存，并指向子單元。

圖1 四叉樹數(shù)據(jù)結(jié)構(gòu)Fig.1 Quadtree data structure

圖2 網(wǎng)格信息存儲Fig.2 Grid information storage

1.2 幾何自適應(yīng)

幾何自適應(yīng)加密是根據(jù)物體幾何特征進行自適應(yīng)加密，要求能夠準確反映模型的外形結(jié)構(gòu)特征。主要從兩個方面進行：首先，基于物體幾何外形，將與物面相交的網(wǎng)格單元進行一定次數(shù)的加密，得到初步的加密網(wǎng)格單元，如圖3(a)所示；接著，基于物體表面曲率，在物體表面曲率變化較大的位置進行進一步加密，如圖3(b)所示。此外，對于相交單元鄰近的單元也需要進行一定次數(shù)的加密，以確保網(wǎng)格過渡均勻。對于表面曲率，通過相鄰單元內(nèi)物面曲率變化來判斷，判據(jù)為:|kcell-kcell_neighbor|>Δk，kcell和kcell_neighbor分別代表目標(biāo)單元和相應(yīng)的鄰近單元內(nèi)的物面曲率，Δk一般不大于0.1，具體數(shù)值根據(jù)網(wǎng)格加密需求確定。

圖3 幾何自適應(yīng)加密網(wǎng)格Fig.3 Grid adaptation based on geometric feature

1.3 流場自適應(yīng)

在流場計算到一定程度后，還需要根據(jù)流場解的特征進行自適應(yīng)。為了較好地捕捉流場中的激波、剪切層等特征結(jié)構(gòu)，本文將速度的散度和旋度判據(jù)相結(jié)合來進行流場解自適應(yīng)[14]，其中速度的散度主要用來捕捉激波，速度的旋度則主要用于剪切層的捕捉。網(wǎng)格單元的速度散度Dcell和旋度Rcell可表示為

(1)

(2)

式中：l為單元cell的邊長；系數(shù)a為給定值，二維時取2，三維時取3。

流場解自適應(yīng)過程中需要對每個單元進行速度的散度和旋度計算，令計算域內(nèi)所有網(wǎng)格單元的總數(shù)記作N，在得到N個單元的Dcell以及Rcell之后，可得到速度散度的標(biāo)準差SD以及速度旋度的標(biāo)準差SR，即

(3)

(4)

如果Dcell>δ1SD或者Rcell>δ2SR,則該網(wǎng)格單元需要加密；如果Dcell<0.1α1SD且Rcell<0.1α2SR,則該網(wǎng)格單元需要粗化。對于系數(shù)δ1、δ2、α1、α2的選取要結(jié)合流場特性給出[4]，這是由于不同情況下流場中占主導(dǎo)地位的流動現(xiàn)象也是不同的。如果流場中沒有明顯主導(dǎo)的流動特征，則上述系數(shù)可以都取1；若流場中剪切層的作用較大，激波作用較小，則相應(yīng)的系數(shù)可以取α1=δ1=2，α2=δ2=0.5；反之，若激波的作用較大，剪切層的作用較小，則可以取α1=δ1=0.5，α2=δ2=2。這些系數(shù)可能并不是最佳值，在具體流場計算中還需要根據(jù)流場特征進行調(diào)整，以使自適應(yīng)后的網(wǎng)格分布更加合理。

2 Navier-Stokes方程數(shù)值求解器

2.1 方程離散

針對二維Navier-Stokes方程：

(5)

式中：E、F為無黏流項；Ev、Fv為黏性流項。

為了避免激波間斷附近的非物理振蕩，對于對流項的離散采用張涵信院士構(gòu)造的無波動、無自由參數(shù)的耗散差分(NND)格式，即[15]

(6)

式中：minmod為限制器；本文E±采用Van Leer通量分裂[16]計算。

黏性項的離散采用中心差分格式，時間上則采用二階Runge-Kutta方法[17-18]：

(7)

式中：R(U)為網(wǎng)格單元的殘差。該格式要具備總變差衰減(Total Variation Diminishing，TVD)性質(zhì)還需要滿足穩(wěn)定性約束條件：

(8)

式中：C為柯朗(Courant)數(shù)；Δt為當(dāng)?shù)貢r間步長；Δx為網(wǎng)格單元步長。式(8)為柯朗-弗里德里奇-列維(CFL)條件，在空間、時間均為二階精度的情況下，常量const取值為1。

通過上述工作，本文針對Navier-Stokes方程構(gòu)造的全離散格式具有二階精度，保持TVD性質(zhì)，可以避免激波間斷附近的非物理振蕩。

2.2 懸掛網(wǎng)格

在網(wǎng)格自適應(yīng)之后，對于粗細網(wǎng)格的過渡區(qū)域，計算模版點流場值的獲取需要特殊處理，即處理懸掛網(wǎng)格問題。懸掛網(wǎng)格問題分為兩類：

第1類，計算模板點具有子單元。以圖4為例，對于目標(biāo)單元Ω，它的計算模板點X具有葉子單元，并非葉子結(jié)點，X的流場值的獲取需要通過子單元插值得到，即

(9)

式中：qa、qb、qc、qd和qX分別代表網(wǎng)格單元a、b、c、d和X處的流場值，r1、r2、r3、r4分別為網(wǎng)格單元a、b、c、d與網(wǎng)格單元X中心點之間的距離。

第2類，計算模板點并非網(wǎng)格單元。以圖5為例，對于目標(biāo)單元Ω，它的最右側(cè)計算模板點并非完整的網(wǎng)格單元，而是葉子單元L的一部分。在這種情況下，就需要對L進行剖分，以確定計算模板點L1的位置信息，再通過鄰近單元插值方法獲取L1處的流場信息，即

圖4 第1類懸掛網(wǎng)格問題Fig.4 First type of hang-grid problem

(10)

圖5 第2類懸掛網(wǎng)格問題Fig.5 Second type of hang-grid problem

3 物面邊界處理

3.1 虛擬鏡像對稱方法

對于自適應(yīng)笛卡爾網(wǎng)格，在物面邊界附近進行插值計算時，可能需要用到物體內(nèi)部的單元，即虛擬單元(Ghost Cell)。對于這些單元流場值的獲取，需要采用虛擬鏡像對稱方法。以圖6為例，為了獲得虛擬單元A的流場值，首先要確定A在流場中的鏡像點MA，通過鄰近單元插值的方法獲得MA的流場值，之后就是根據(jù)物面的邊界條件(壁面速度的無穿透、無滑移以及絕熱壁或等溫壁等條件)，確定虛擬點A和鏡像點MA流場值之間的關(guān)系，從而得到虛擬點A處的流場信息。

圖6 虛擬鏡像對稱方法Fig.6 Ghost mirror-symmetry method

鏡像點流場信息的獲得方法是通過鄰近流場單元插值，對于虛擬點A的鏡像點MA，其鄰近的4個單元結(jié)點分別記作C、D、E、F，則MA點處的流場值為

(11)

式中：r1、r2、r3、r4分別為點MA到網(wǎng)格點C、D、E、F的距離。

3.2 曲率修正技術(shù)

在曲率修正對稱方法(Curvature Corrected Symmetry Technique，CCST)[19]中，插值虛擬點的流場值是通過在物面附近虛構(gòu)的渦流場得到的，渦流場的熵和總焓在法向上對稱分布。Dadone和Grossman在針對貼體網(wǎng)格所發(fā)展的曲率修正對稱方法的基礎(chǔ)上，將其推廣到非貼體的笛卡爾網(wǎng)格，即虛擬體單元方法(Ghost Body-Cell Method，GBCM)[20]。

如圖7所示，虛擬點A和鏡像點B的流場值滿足：

圖7 虛擬體單元方法Fig.7 Ghost body-cell method

(12)

值得一提的是，鏡像點處的流場值是通過臨近的流場單元插值得到的，若是4個臨近單元中有2個或者2個以上在物體內(nèi)部，可能會影響插值點處流場值的精度，因此要盡可能地使鏡像點所有的插值點都在流場中。本文借鑒Forrer等提出的虛擬單元方法(Forrer’s Ghost Cell Method，F(xiàn)GCM)[21]對于鏡像點的處理思路，針對該問題提出了新的鏡像點位置取法，如圖8所示，鏡像點的位置不是取在和虛擬點A關(guān)于物面對稱的位置A′，而是將線段AA′再延伸一個網(wǎng)格長度Δl，從而確保鏡像點周圍的插值點都在流場之中。Forrer和Berger等[21]的研究表明，這種位置取法可以提高計算的穩(wěn)定性，而且不會降低邊界處理的精度，其不足之處是局部的質(zhì)量和動量守恒性難以保證。

Dadone等通過研究表明，在CCST方法基礎(chǔ)上推廣得到的GBCM在保證熵的增加和總焓的減少方面，相對于傳統(tǒng)的物面邊界處理方法有較大優(yōu)勢，有利于計算的穩(wěn)定性和收斂性。

圖8 鏡像點位置調(diào)整Fig.8 Adjustment of position of mirroring point

3.3 多值點問題

使用浸入邊界方法模擬研究具有尖點或薄體結(jié)構(gòu)的外形時，會遇到多值點問題[20]。對于多值點問題的有效處理是自適應(yīng)笛卡爾網(wǎng)格方法自動化處理任意復(fù)雜外形的關(guān)鍵技術(shù)之一，對于解算器的可靠度和精度有重要影響。

多值點單元的類型主要分為兩種：流場中的多值點單元和物面內(nèi)的多值點單元。在圖9(a)中，網(wǎng)格點X作為目標(biāo)單元Ω的計算模板點被調(diào)用，但是X與Ω在流場中異側(cè)，不能夠直接調(diào)用網(wǎng)格點X存儲的流場值，X在此時即為流場多值點單元；在圖9(b)中，上下兩側(cè)的目標(biāo)單元Ω1、Ω2的計算模板點均用到虛擬單元G，該單元就存在兩組流場值中，稱之為物體內(nèi)的多值點單元。

對于流場中的多值點單元，以圖10中的網(wǎng)格點X為例，由于其與目標(biāo)單元Ω異側(cè)，不能直接使用其本身的流場值，而是將其看作上側(cè)流場的虛擬單元，將網(wǎng)格點X關(guān)于上邊界l1取對稱點MX，通過虛擬鏡像對稱方法確定模版點X的流場值。

對于物體內(nèi)的多值點單元，以圖11中的單元G為例，上下兩側(cè)目標(biāo)單元Ω1、Ω2的計算模板點均用到它，需要對G分別關(guān)于l1、l2取鏡像點M1、M2，通過鏡像虛擬對稱方法確定其對應(yīng)的兩組流場值，并分別存儲，根據(jù)目標(biāo)單元的位置進行調(diào)用。

圖9 多值點單元類型Fig.9 Types of multi-valued grids

對于多值點單元的信息的存儲，需要在網(wǎng)格數(shù)據(jù)結(jié)構(gòu)中多分配一塊內(nèi)存，以存儲網(wǎng)格單元不同情況下的流場信息，根據(jù)具體情況選擇要調(diào)用的數(shù)據(jù)，如圖12和圖13所示。這種方式雖然增加了存儲量，但是相對于每次調(diào)用都重新計算的方式，計算量大大減少，調(diào)用方便高效。

圖10 流場多值點單元處理方法Fig.10 Processing method for multi-valued grid in flow field

圖11 物體多值點單元處理方法Fig.11 Processing method for multi-valued grid inside the body

圖12 流場多值點單元信息存儲Fig.12 Information storage of multi-valued grids in flow field

圖13 物體多值點單元信息存儲Fig.13 Information storage of multi-valued grids inside the body

4 算例考核

下面針對圓柱繞流、前臺階流動和三角尖劈繞流以及壓縮拐角流動等典型算例進行數(shù)值模擬，考核本文所發(fā)展的自適應(yīng)笛卡爾網(wǎng)格生成技術(shù)和構(gòu)建的黏性數(shù)值求解器的準確可靠性。

4.1 圓柱繞流

參照實驗[22]進行參數(shù)設(shè)置，來流馬赫數(shù)Ma=1.7，圓柱模型直徑D=2.0 m，基于直徑的雷諾數(shù)Re=2.0×105。初始計算網(wǎng)格邊長為0.4 m，最大加密層數(shù)為8層，計算過程中基于速度的散度和旋度判據(jù)進行了4次自適應(yīng)加密(Adaptive Mesh Refinement,AMR)，即AMR=4，相應(yīng)的自適應(yīng)網(wǎng)格如圖14(a)所示，圖14(b)為馬赫數(shù)云圖，圖15為密度(ρ)等值線和流場流線。

結(jié)合流場自適應(yīng)網(wǎng)格和馬赫數(shù)云圖、密度等值線以及流線圖可以看出，網(wǎng)格分布合理，層次清晰，曲線過渡光滑，外形信息保真度高，對激波間斷、剪切層以及尾流渦等結(jié)構(gòu)捕捉能力較強，網(wǎng)格自適應(yīng)效果較好。

圖16中，本文得到的表面壓力系數(shù)Cp分布與實驗及文獻[23]結(jié)果吻合得較好。通過表1中分離點位置以及阻力系數(shù)等數(shù)據(jù)的對比，可以看出，本文所得到的結(jié)果與實驗值更為接近，證明本文所發(fā)展的網(wǎng)格技術(shù)和構(gòu)建的求解器對于超聲速圓柱繞流問題的求解具有較高的精度和可靠性。

圖14 流場自適應(yīng)Cartesian網(wǎng)格和馬赫數(shù)云圖Fig.14 Adaptive Cartesian grid and Mach number contour of flow field

圖15 密度等值線和圓柱尾流區(qū)流線Fig.15 Density contour lines and streamlines in wake region of cylinder

圖16 表面壓力系數(shù)分布對比Fig.16 Comparison of surface pressure coefficients distribution

表1本文數(shù)據(jù)與文獻及實驗結(jié)果對比

Table1Comparisonoftheproposeddataandreferenceandexperimentresults

參數(shù)文獻[23]實驗本文網(wǎng)格數(shù)目21600不存在32542分離點θ/(°)133．0112．0120．0阻力系數(shù)CD1．501．431．47

4.2 前臺階流動

前臺階激波反射流動問題包含激波的多次壁面反射，激波與激波相互作用，激波與中心膨脹波、接觸間斷的碰撞以及波系的相互干擾，是檢驗計算格式性能和網(wǎng)格方法的經(jīng)典算例。

臺階的尺寸為3 m(長)×1 m(寬)，臺階位于離風(fēng)洞入口0.6 m處，高度為0.2 m。遠場水平來流馬赫數(shù)Ma=3.0，進、出口邊界分別設(shè)為超聲速入、出流邊界，上、下壁面均設(shè)為滑移絕熱固壁邊界。相應(yīng)的模型構(gòu)建和初始網(wǎng)格生成如圖17所示。其中初始單元長0.06 m，最大加密層數(shù)為5層。

計算在進行到t=4.0 s時得到的自適應(yīng)網(wǎng)格和流場結(jié)果分別如圖18和圖19所示。結(jié)合兩組圖中的自適應(yīng)網(wǎng)格和流場結(jié)果，自適應(yīng)網(wǎng)格的分布清晰地捕捉了經(jīng)壁面反射形成的λ型激波和馬赫桿(Mach Stem)、三叉點處的滑移線、臺階前緣奇點處的中心膨脹波以及與壁面相撞形成的反射激波等結(jié)構(gòu)，網(wǎng)格分布均勻有序，激波線過渡光滑，自適應(yīng)效果較好。

圖17 前臺階及初始計算網(wǎng)格Fig.17 Front step and initial computation grid

圖18 t=4.0 s自適應(yīng)網(wǎng)格Fig.18 Adaptive grid at t=4.0 s

圖19 t=4.0 s流場計算結(jié)果Fig.19 Calculated results of flow field at t=4.0 s

圖20 激波線位置與文獻[24]結(jié)果對比Fig.20 Comparison of the shock wave line in this paper and Ref.[24]

在圖20中，流場等值線分布與文獻[24]結(jié)果符合較好。該文獻基于非結(jié)構(gòu)網(wǎng)格，采用的是時空三階精度有限元格式(Third-order Taylor Galerkin Convection，TTGC)，本文與其h=0.01(網(wǎng)格數(shù)為56 892)的結(jié)果相對比，在前緣弓形激波靠近水平壁面處，本文的激波線與壁面垂直性更好，更貼近理論情況，并且對于前緣弓形激波和馬赫桿的捕捉更加精細；在網(wǎng)格數(shù)目方面，本文網(wǎng)格數(shù)目為26 238，說明在流場波系較復(fù)雜時，本文發(fā)展的網(wǎng)格方法網(wǎng)格數(shù)目較少。綜上，本文所發(fā)展的網(wǎng)格生成技術(shù)和構(gòu)建的數(shù)值求解器在計算前臺階流動問題的可靠性和優(yōu)勢方面得到驗證。

4.3 尖劈繞流

為了考察本文所發(fā)展的方法技術(shù)對于多值點問題的處理能力，選取三角尖劈模型，在超聲速來流情況下進行流場的模擬研究。

尖劈模型頂角θ=30°，高度H=1.0 m，來流馬赫數(shù)Ma=2.2，基于尖劈高度的雷諾數(shù)Re=1.0×105。初始計算網(wǎng)格長0.2 m，最大加密次數(shù)為7，AMR=5，相應(yīng)的流場自適應(yīng)網(wǎng)格和流場結(jié)果如圖21所示。

圖21 流場自適應(yīng)網(wǎng)格及流場結(jié)果Fig.21 Adaptive grid and results of flow field

根據(jù)流場解的速度散度和旋度判據(jù)生成的自適應(yīng)網(wǎng)格，對于附體斜激波、膨脹波、剪切層以及分離區(qū)等流場結(jié)構(gòu)和特征捕捉清晰，網(wǎng)格分布均勻，層次分明，冗余網(wǎng)格和網(wǎng)格空洞很少，進一步證明了本文網(wǎng)格生成方法具有很好的自適應(yīng)能力和較高的質(zhì)量，并且能夠有效處理多值點問題。

下面將不同馬赫數(shù)下的附體斜激波的激波角與理論值進行對比，其理論公式[25]為

式中：η為斜面偏轉(zhuǎn)角；β為斜激波角，即附體斜激波的激波線與水平方向的夾角。

在表2中，對于附體斜激波的激波角，本文得到的計算值和理論值相差較小，均在4%以下，證明了本文方法能夠較準確地捕捉激波位置，具有較高的精度和可靠性。

表2斜激波角計算值與理論值對比

Table2Comparisonofcalculatedvalueandtheoreticalvalueofobliqueshockwaveangles

馬赫數(shù)2．25．07．0斜激波角理論值/(°)41．2524．2921．54斜激波角計算值/(°)42．3725．0722．11誤差/%2．713．212．65

4.4 壓縮拐角流動

壓縮拐角流動問題的幾何邊界雖然簡單，但包含了剪切層、流動分離和再附等多種復(fù)雜流動結(jié)構(gòu)和現(xiàn)象。本文針對該問題進一步考核本文的網(wǎng)格技術(shù)和構(gòu)建的數(shù)值求解器對于復(fù)雜流動問題求解的準確可靠性。

選取Ma=3.0，基于壓縮拐角模型的雷諾數(shù)Re=7.5×104，壓縮角度為24°，總溫T0=300 K，壁面溫度Tw=297 K。網(wǎng)格最大加密層數(shù)為9層，AMR=5。相應(yīng)的自適應(yīng)網(wǎng)格和流場結(jié)果如圖22所示。

從圖22可以看出，本文所發(fā)展的自適應(yīng)笛卡爾網(wǎng)格技術(shù)對于壓縮拐角的波系結(jié)構(gòu)可以清晰的捕捉，自適應(yīng)效果較好。由壓縮坡面引起的逆壓梯度影響了層流邊界層，對于流動產(chǎn)生了阻礙，導(dǎo)致入口邊界附近產(chǎn)生了一道誘導(dǎo)激波。在拐點附近的上游區(qū)域，流動產(chǎn)生分離，形成了分離激波和分離區(qū)，之后在壓縮坡面上的某一位置發(fā)生再附現(xiàn)象，邊界層重新發(fā)展，并形成再附激波。分離激波和再附激波相互作用，沿壓縮坡面向上傳播。

改變模型壓縮面角度，在壓縮角分別為15°、18°和24°的情況下進行數(shù)值模擬，并將計算得到的表面壓力系數(shù)Cp分別與實驗值進行對比，如圖23所示。

通過不同壓縮角情況下表面壓力系數(shù)的計算結(jié)果和實驗結(jié)果的對比發(fā)現(xiàn)，本文計算結(jié)果和實驗符合較好，進一步證明了本文方法技術(shù)和求解器對于壓縮拐角流動問題求解的可行性和可靠性。另外，在壓縮角為24°時，兩者相差相對較大，Rudy等[26]在數(shù)值模擬過程中也遇到過類似問題，認為是由于實驗中存在三維效應(yīng)，對于實驗的測量值有所影響。

圖22 流場自適應(yīng)網(wǎng)格及流場結(jié)果Fig.22 Adaptive grid and results of flow field

圖23 不同壓縮角下計算與實驗表面壓力系數(shù)對比Fig.23 Comparison of surface pressure coefficients at different compression corners in calculation and experiment

5 結(jié) 論

本文針對二維情況下含強間斷的超聲速黏性流動問題，發(fā)展了自適應(yīng)笛卡爾網(wǎng)格生成技術(shù)和黏性物面邊界的處理方法，并構(gòu)建了相應(yīng)的Navier-Stokes方程數(shù)值求解器，在此基礎(chǔ)上進行了典型算例的考核驗證。

1) 在笛卡爾網(wǎng)格自適應(yīng)技術(shù)方面，本文選用叉樹數(shù)據(jù)結(jié)構(gòu)進行網(wǎng)格信息的存儲，分別從物體外形和流場特征出發(fā)進行網(wǎng)格的自適應(yīng)加密和粗化，發(fā)展了一種自動高效的自適應(yīng)笛卡爾網(wǎng)格生成方法，具有流場特征捕捉效果好、自動化程度高等特點。

2) 對于黏性物面邊界的處理，本文基于浸入邊界方法，結(jié)合虛擬鏡像對稱方法和曲率修正技術(shù)進行黏性物面邊界的處理，同時建立了多值點問題的處理技術(shù)，發(fā)展了一種在笛卡爾網(wǎng)格下可有效模擬黏性物面邊界條件的方法，經(jīng)算例考核，證明該方法具有較高的精度和可靠性。

3) 在數(shù)值求解器方面，針對Navier-Stokes方程，建立了具有二階精度和保持TVD性質(zhì)的全離散格式，避免了解在激波間斷附近的非物理振蕩；構(gòu)建了粗細網(wǎng)格間懸掛網(wǎng)格的插值方法，可以有效應(yīng)用于離散方程的求解。

目前，基于自適應(yīng)笛卡爾網(wǎng)格的超聲速黏性流動數(shù)值模擬尚未工程實用化，有進一步發(fā)展的需求。后續(xù)工作的首要目標(biāo)是將本文方法推廣到三維任意復(fù)雜外形的流動問題模擬中，重點研究三維復(fù)雜幾何外形物體的自動化處理方法，建立相應(yīng)的自適應(yīng)笛卡爾網(wǎng)格生成技術(shù)，并針對Navier-Stokes方程構(gòu)建黏性數(shù)值求解器。同時，通過自適應(yīng)笛卡爾網(wǎng)格生成技術(shù)與人工智能技術(shù)的結(jié)合，發(fā)展智能化、自動化生成高質(zhì)量網(wǎng)格的方法；發(fā)展各向異性的笛卡爾網(wǎng)格生成技術(shù)，減少網(wǎng)格數(shù)量和數(shù)據(jù)存儲量；將生成笛卡爾網(wǎng)格的軟件直接與模型生成和處理軟件對接，盡可能地形成一體化、軟件化和商品化。

參 考 文 獻

[1] MAVRIPLIS D J. An advancing front Delaunay triangulation algorithm designed for robustness[J]. Journal of Computational Physics, 1995, 117(1): 90-101.

[2] ITO Y, NAKAHASHI K. Unstructured hybrid grid generation baseds on isotropic tetrahedral grids: AIAA-2002-0861[R]. Reston, VA: AIAA, 2002.

[3] KALLINDERIS Y, WARD S. Prismatic grid generation with an efficient algebraic method for aircraft configuration: AIAA-1992-2721[R]. Reston, VA: AIAA, 1992.

[4] DE ZEEUW D L. A quadtree-based adaptively-refined Cartesian-grid algorithm for solution of the Euler equations[D]. Michigan: University of Michigan, 1993.

[5] MURAT B. Quadtree-based adaptively-refined Cartesian-grid algorithm for solution of the Euler equations[D]. Michigan: University of Michigan, 2005.

[6] ROBERT R A, VLADIMIR I K. Direct simulation Monte Carlo with octree Cartesian mesh[J]. AIAA Journal, 2012, 20(2): 25-28.

[7] FUJIMOTO K, FUJII K, WANG Z J. Improvements in the reliability and efficiency of body-fitted Cartesian grid method: AIAA-2009-1173[R]. Reston, VA: AIAA, 2009.

[8] UDAYKUMAR H S. Multiphase dynamics in arbitrary geometries on fixed Cartesian grids[J]. Journal of Computational Physics, 1997, 137(2): 366-405.

[9] UDAYKUMAR H S. A sharp interface Cartesian grid method for simulating flow with complex moving boundaries[J]. Journal of Computational Physics, 2001, 174(1): 345-380.

[10] MARSHALL D D. Extending the functionalities of Cartesian grid solvers: Viscous effects modeling and MPI parallelization[D]. Georgia: Georgia Institute of Technology, 2002.

[11] FORRER H, JELTSCH R. A higher order boundary treatment for Cartesian-grid method[J]. Journal of Computational Physics, 1998, 140(2): 259-277.

[12] JAE-DOO L, RUFFIN S M. Development of a turbulent wall-function based viscous Cartesian-grid methodology: AIAA-2007-1326[R]. Reston, VA: AIAA, 2007.

[13] JAE-DOO 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.

[14] DE ZEEUW D L. A wave-model-based refinement criterion for adaptive-grid computation of compressible flow: AIAA-1992-0332[R]. Reston, VA: AIAA, 1992.

[15] 張涵信. 無波動、無自由參數(shù)的耗散差分格式[J]. 空氣動力學(xué)學(xué)報, 1988, 6(2): 143-165.

ZHANG H X. A non-oscillatory, containing no free parameters and dissipative scheme[J]. Acta Aerodynamica Sinica, 1988, 6(2): 143-165 (in Chinese).

[16] VAN LEER B. Flux vector splitting for Euler equations[J]. Lecture Notes in Physics, 1982: 507-512.

[17] ENGQUIST B, MAJDA A. Absorbing boundary conditions for the numerical simulation of waves[J]. Mathematics of Computation, 1977, 31: 629-651.

[18] GUSIAFSSON B, SANDSTROM A. Incompletely parabolic problems in fluid dynamics[J]. SIAM Journal on Applied Mathematics, 1978, 35: 343-357.

[19] DADONE A, GROSSMAN B. An immersed body methodology for inviscid flows on Cartesian grids: AIAA-2002-1059[R]. Reston, VA: AIAA, 2002.

[20] DADONE A, GROSSMAN B. Surface boundary conditions for the numerical solution of the Euler equations[J]. AIAA Journal, 1994, 32(2): 285-293.

[21] FORRER H, BERGER M. Flow simulations on Cartesian grids involving complex moving geometries[C]∥Proceedings of 6th International Conference on Hyperbolic Problems, 1998.

[22] BASHKIN V A, VAGANOV A V, EGOROV I V, et al. Comparison of calculated and experimental data on supersonic flow past a circular cylinder[J]. Fluid Dynamics, 2002, 37(3): 473-483.

[23] MUNIKRISHNA N, LIOU M S. A Cartesian based body-fitted adaptive grid method for compressible viscous flows[R]. Reston, VA: AIAA, 2009.

[24] 魯陽, 鄒建鋒, 鄭耀. 基于非結(jié)構(gòu)網(wǎng)格的TTGC有限元格式的實現(xiàn)及在超聲速流動中的應(yīng)用[J]. 計算力學(xué)學(xué)報, 2013, 30(5): 712-722.

LU Y, ZOU J F, ZHENG Y. Unstructured grids based on TTGC finite element scheme and its application to supersonic flows[J]. Chinese Journal of Computational Mechanics, 2013, 30(5): 712-722 (in Chinese).

[25] BOIRON O, CHIAVASSA G, DONAT R. A high-resolution penalization method for large Mach number flows in the presence of obstacles[J]. Computers & Fluids, 2009, 38(3): 703-714.

[26] RUDY D H, THOMAS J L, KUMAR A. Computation of laminar viscous-inviscid interactions in high-speed internal flows[J]. AIAA Journal, 1991, 29(7): 1108-1113.