自適應(yīng)笛卡爾網(wǎng)格超聲速黏性流動(dòng)數(shù)值模擬

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

中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

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

自適應(yīng)笛卡爾網(wǎng)格生成技術(shù)的主要不足之處在于黏性物面邊界處理較復(fù)雜，這是由笛卡爾網(wǎng)格的非貼體特性導(dǎo)致的。笛卡爾網(wǎng)格一般與物面相交，產(chǎn)生不同形狀的切割單元。這些切割單元有的尺寸很小，容易引起在物面附近解的非物理振蕩，如果數(shù)值模擬過(guò)程中采用的差分格式是顯式的，還會(huì)限制時(shí)間步長(zhǎng)以及影響計(jì)算的穩(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)格的超聲速黏性流動(dòng)數(shù)值模擬尚未工程實(shí)用化，有進(jìn)一步發(fā)展的需求。為了實(shí)現(xiàn)對(duì)二維含強(qiáng)間斷的高超聲速黏性流動(dòng)問(wèn)題自動(dòng)、高效、高精度的模擬，本文發(fā)展了自適應(yīng)笛卡爾網(wǎng)格生成技術(shù)，以及笛卡爾網(wǎng)格下可有效模擬黏性物面邊界條件的方法，并構(gòu)建了相應(yīng)的黏性數(shù)值求解器。

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

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

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

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

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

1.2 幾何自適應(yīng)

幾何自適應(yīng)加密是根據(jù)物體幾何特征進(jìn)行自適應(yīng)加密，要求能夠準(zhǔn)確反映模型的外形結(jié)構(gòu)特征。主要從兩個(gè)方面進(jìn)行：首先，基于物體幾何外形，將與物面相交的網(wǎng)格單元進(jìn)行一定次數(shù)的加密，得到初步的加密網(wǎng)格單元，如圖3(a)所示；接著，基于物體表面曲率，在物體表面曲率變化較大的位置進(jìn)行進(jìn)一步加密，如圖3(b)所示。此外，對(duì)于相交單元鄰近的單元也需要進(jìn)行一定次數(shù)的加密，以確保網(wǎng)格過(guò)渡均勻。對(duì)于表面曲率，通過(guò)相鄰單元內(nèi)物面曲率變化來(lá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 流場(chǎng)自適應(yīng)

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

(1)

(2)

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

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

(3)

(4)

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

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

2.1 方程離散

針對(duì)二維Navier-Stokes方程：

(5)

式中：E、F為無(wú)黏流項(xiàng)；Ev、Fv為黏性流項(xiàng)。

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

(6)

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

黏性項(xiàng)的離散采用中心差分格式，時(shí)間上則采用二階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(shí)間步長(zhǎng)；Δx為網(wǎng)格單元步長(zhǎng)。式(8)為柯朗-弗里德里奇-列維(CFL)條件，在空間、時(shí)間均為二階精度的情況下，常量const取值為1。

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

2.2 懸掛網(wǎng)格

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

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

(9)

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

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

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

(10)

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

3 物面邊界處理

3.1 虛擬鏡像對(duì)稱方法

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

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

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

(11)

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

3.2 曲率修正技術(shù)

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

如圖7所示，虛擬點(diǎn)A和鏡像點(diǎn)B的流場(chǎng)值滿足：

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

(12)

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

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

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

3.3 多值點(diǎn)問(wèn)題

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

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

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

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

圖9 多值點(diǎn)單元類型Fig.9 Types of multi-valued grids

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

圖10 流場(chǎng)多值點(diǎn)單元處理方法Fig.10 Processing method for multi-valued grid in flow field

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

圖12 流場(chǎng)多值點(diǎn)單元信息存儲(chǔ)Fig.12 Information storage of multi-valued grids in flow field

圖13 物體多值點(diǎn)單元信息存儲(chǔ)Fig.13 Information storage of multi-valued grids inside the body

4 算例考核

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

4.1 圓柱繞流

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

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

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

圖14 流場(chǎng)自適應(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ù)分布對(duì)比Fig.16 Comparison of surface pressure coefficients distribution

表1本文數(shù)據(jù)與文獻(xiàn)及實(shí)驗(yàn)結(jié)果對(duì)比

Table1Comparisonoftheproposeddataandreferenceandexperimentresults

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

4.2 前臺(tái)階流動(dòng)

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

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

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

圖17 前臺(tái)階及初始計(jì)算網(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流場(chǎng)計(jì)算結(jié)果Fig.19 Calculated results of flow field at t=4.0 s

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

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

4.3 尖劈繞流

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

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

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

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

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

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

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

表2斜激波角計(jì)算值與理論值對(duì)比

Table2Comparisonofcalculatedvalueandtheoreticalvalueofobliqueshockwaveangles

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

4.4 壓縮拐角流動(dòng)

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

選取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)格和流場(chǎng)結(jié)果如圖22所示。

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

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

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

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

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

5 結(jié) 論

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

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

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

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

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

參 考 文 獻(xiàn)

[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] 張涵信. 無(wú)波動(dòng)、無(wú)自由參數(shù)的耗散差分格式[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 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] 魯陽(yáng), 鄒建鋒, 鄭耀. 基于非結(jié)構(gòu)網(wǎng)格的TTGC有限元格式的實(shí)現(xiàn)及在超聲速流動(dòng)中的應(yīng)用[J]. 計(jì)算力學(xué)學(xué)報(bào), 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.