張啟明，葉友達(dá)，*，蔣勤學(xué)，田 浩

(1. 北京流體動(dòng)力科學(xué)研究中心，北京 100120；2. 國(guó)家計(jì)算流體力學(xué)實(shí)驗(yàn)室，北京 100191)

0 引 言

計(jì)算流體動(dòng)力學(xué)技術(shù)的進(jìn)步已經(jīng)到了幾乎所有流動(dòng)問題都能在可接受時(shí)長(zhǎng)內(nèi)解決的階段。網(wǎng)格是數(shù)值計(jì)算的基礎(chǔ)。但高質(zhì)量的網(wǎng)格生成所占的時(shí)間漸漸比求解流動(dòng)方程更多。我國(guó)著名空氣動(dòng)力學(xué)家張涵信院士將CFD 領(lǐng)域概括為5M1A 模型，即網(wǎng)格（Mesh）、數(shù)值方法（Method）、計(jì)算機(jī)（Machine）、機(jī)理（Mechanics）、繪圖（Mapping）和應(yīng)用（Application）。經(jīng)驗(yàn)豐富的CFD 應(yīng)用者，在生成網(wǎng)格時(shí)，是要將具體的流動(dòng)應(yīng)用背景、流動(dòng)機(jī)理、數(shù)值方法以及后處理等多個(gè)方面通盤考慮，而CFD 開發(fā)人員，則需要額外的計(jì)算機(jī)軟硬件基礎(chǔ)。

棱柱網(wǎng)格由于在某幾個(gè)方面具有特殊的優(yōu)勢(shì)，因而廣泛地應(yīng)用在邊界層網(wǎng)格的生成中：從網(wǎng)格生成的角度，物面三角化網(wǎng)格靈活性強(qiáng)，適用性高，邊界層內(nèi)棱柱網(wǎng)格生成自動(dòng)化程度高，當(dāng)面對(duì)復(fù)雜的外形時(shí)，可以大大減少人工交互；從流動(dòng)機(jī)理上看，靠近物面的黏性區(qū)域在垂直物面方向上的網(wǎng)格具有方向性，在高雷諾數(shù)流動(dòng)或者含湍流的流動(dòng)中對(duì)黏性邊界層的分辨至關(guān)重要；在數(shù)值方法上，在遠(yuǎn)離物面的方向，網(wǎng)格具有結(jié)構(gòu)化編號(hào)，這種內(nèi)在的網(wǎng)格結(jié)構(gòu)可以應(yīng)用到通量格式中以提高計(jì)算精度，也可以應(yīng)用到隱式時(shí)間推進(jìn)中以加速收斂。

根據(jù)以上觀察，Meakin 等[1]于2007 年引入了串網(wǎng)格劃分方法。在這種方法中，通過直接從由三角形和四邊形組成的離散曲面向外推進(jìn)獲得適合黏性捕捉的棱柱網(wǎng)格，這種體網(wǎng)格生成過程完全自動(dòng)化。離散曲面的頂點(diǎn)向外生成一系列曲線，這些曲線稱為“串”，通常是由指向向量和長(zhǎng)度表示的直線（見圖1），可以僅用幾個(gè)參數(shù)來表示。邊界層內(nèi)捕捉流場(chǎng)細(xì)節(jié)的體網(wǎng)格是就是通過這些曲線上的點(diǎn)構(gòu)建而成。

圖1 “串”的幾何示意圖[1]Fig. 1 Geometry schematic of a strand

串網(wǎng)格的生成，早期采用非結(jié)構(gòu)網(wǎng)格邊界層內(nèi)棱柱網(wǎng)格生成的方法。棱柱網(wǎng)格的生成，可以分為基于幾何的推進(jìn)方法、基于各向異性四面體聚合的方法和基于場(chǎng)的推進(jìn)方法?；趲缀蔚姆椒?，根據(jù)當(dāng)?shù)貛缀翁卣?，確定推進(jìn)方向和推進(jìn)步長(zhǎng)，重點(diǎn)解決在模型凹特征區(qū)域和鄰近特征區(qū)域附近的邊界層前沿面上可能出現(xiàn)的相交問題?；诰W(wǎng)格聚合的方法[2]將自動(dòng)生成的各向異性四面體網(wǎng)格基于一定的判定準(zhǔn)則進(jìn)行聚合，得到物面附近的棱柱網(wǎng)格單元。

基于場(chǎng)的方法，根據(jù)某種模型建立物面網(wǎng)格對(duì)空間網(wǎng)格影響的代數(shù)方程或者偏微分方程，由該方程確定的空間網(wǎng)格、網(wǎng)格推進(jìn)方向或者網(wǎng)格推進(jìn)步長(zhǎng)自然滿足網(wǎng)格的光滑性、正交性等條件，能夠有效避免或者減輕基于幾何的層推進(jìn)方法中遇到的推進(jìn)方向不確定或者可能出現(xiàn)的相交問題。基于場(chǎng)的方法的核心在于建立有效的控制方程，該方程能模擬計(jì)算出從物面到外部空間某種物理量，這個(gè)物理量自身的性質(zhì)有助于確定網(wǎng)格推進(jìn)中所需的推進(jìn)方向、推進(jìn)步長(zhǎng)或者是直接的空間網(wǎng)格點(diǎn)。下邊對(duì)已有的基于場(chǎng)的網(wǎng)格生成方法進(jìn)行介紹。

Nakahashi[3]通過求解控制推進(jìn)面總面積最小的變分方程，來確定推進(jìn)方向上的網(wǎng)格步長(zhǎng)，進(jìn)而確定每一層的網(wǎng)格面。Sethian 等[4]設(shè)計(jì)的水平集方法是通過求解Hamilton-Jacobi 方程得到的一系列隨時(shí)間演化的前沿面。Dawes 等[5]將水平集方法應(yīng)用到空間邊界層網(wǎng)格的生成中。Park 等[6]結(jié)合經(jīng)典層進(jìn)法和水平集方法，前沿點(diǎn)處的層進(jìn)法向定義為水平集方法中隱式函數(shù)的梯度，并將前沿點(diǎn)沿著梯度方向投影到下一層前沿面上。在該方法中初始前沿面的幾何信息得以借助前沿點(diǎn)的推進(jìn)過程傳遞到各等值面上。Wang[7]將水平集方法的控制方程改為空間網(wǎng)格點(diǎn)距離物面最小歐氏距離所滿足的Eikonal 方程，通過求解該方程并計(jì)算其梯度確定推進(jìn)方向，進(jìn)而確定下一個(gè)前沿面。Sikara[8]通過類比靜電場(chǎng)，建立勢(shì)函數(shù)滿足的拉普拉斯方程，使用邊界積分方法求解拉普拉斯方程，勢(shì)函數(shù)的梯度即電場(chǎng)方向，就是推進(jìn)方向。Takanashi[9]在各個(gè)網(wǎng)格點(diǎn)放置不同的電荷，電荷量通過求解關(guān)于電場(chǎng)函數(shù)的最優(yōu)化問題得到，進(jìn)而通過電場(chǎng)方向確定推進(jìn)方向。Zheng[10]利用快速多極子方法加速的邊界積分方法求解三組拉普拉斯方程，確定推進(jìn)方向場(chǎng)的三個(gè)分量，其本質(zhì)為：勢(shì)函數(shù)滿足拉普拉斯方程，那么勢(shì)函數(shù)梯度的每一個(gè)分量也滿足拉普拉斯方程。孫巖[11]提出一種交互式棱柱網(wǎng)格生成方法，該方法交互地生成物面邊界點(diǎn)的空間推進(jìn)面網(wǎng)格，將邊界點(diǎn)空間推進(jìn)面網(wǎng)格的每層網(wǎng)格點(diǎn)看成邊界網(wǎng)格點(diǎn)運(yùn)動(dòng)得到, 計(jì)算對(duì)應(yīng)的邊界網(wǎng)格點(diǎn)位移，利用徑向基函數(shù)插值將邊界網(wǎng)格點(diǎn)位移光滑傳遞給內(nèi)部網(wǎng)格點(diǎn)，獲得內(nèi)部網(wǎng)格點(diǎn)的空間推進(jìn)網(wǎng)格。

雙曲網(wǎng)格生成方法歷史悠久。Steger 和Chan[12]提出并改進(jìn)了雙曲網(wǎng)格生成方法，采用中心差分加黏性耗散方式離散控制方程，成功應(yīng)用到結(jié)構(gòu)網(wǎng)格中。Tai[13]通過迎風(fēng)格式離散雙曲控制方程，并顯式加入自適應(yīng)的耗散。Steger[14]也提出了將雙曲網(wǎng)格生成方法應(yīng)用到非結(jié)構(gòu)網(wǎng)格里去的思路，但由于非結(jié)構(gòu)網(wǎng)格坐標(biāo)轉(zhuǎn)換所需的網(wǎng)格連接結(jié)構(gòu)生成的困難，Steger 的文章只有一個(gè)圍繞簡(jiǎn)單球體的網(wǎng)格結(jié)果，沒有證明能應(yīng)用于復(fù)雜實(shí)用幾何的外形。Matsuno[15]應(yīng)用二階迎風(fēng)TVD 的方法離散雙曲控制方程，取代Steger 和Chan 方法中的中心差分加耗散，避免了耗散系數(shù)的人工調(diào)節(jié)。Matsuno[16]進(jìn)一步將其應(yīng)用到三角化網(wǎng)格面的棱柱網(wǎng)格生成中，前沿點(diǎn)是控制單元的格心，局部坐標(biāo)通過與相鄰單元格心的坐標(biāo)數(shù)據(jù)計(jì)算得到，下一層格心點(diǎn)確定后，通過簡(jiǎn)單幾何平均得到格點(diǎn)處的網(wǎng)格坐標(biāo)。

Matsuno 的方法有三個(gè)局限性：第一，基于格心的局部坐標(biāo)系的指定具有隨機(jī)性；第二，所能生成的網(wǎng)格僅限于三角形的網(wǎng)格；第三，所采用的三角形網(wǎng)格必須滿足一定的條件，角度不能過分扭曲。此外，對(duì)于所有和球面同胚的多面體，根據(jù)歐拉公式[17]，F(xiàn)+V=E+2，F(xiàn)為面元個(gè)數(shù)，V為頂點(diǎn)個(gè)數(shù)，E為面元邊的個(gè)數(shù)，對(duì)于三角化的表面，2E= 3F，易得F= 2V-4，也就是說面元的個(gè)數(shù)是頂點(diǎn)個(gè)數(shù)的兩倍。因此Matsuno方法中需要求解的以控制單元的格心坐標(biāo)推進(jìn)量為未知量，求解所需的內(nèi)存和計(jì)算量都比較大。

本文基于Knupp[18]提出的定義在非結(jié)構(gòu)網(wǎng)格局部離散點(diǎn)上的邏輯空間來近似計(jì)算坐標(biāo)變換導(dǎo)數(shù)。進(jìn)而將Steger 和Chan 的雙曲網(wǎng)格生成方法[12]應(yīng)用到非結(jié)構(gòu)棱柱網(wǎng)格上。在非結(jié)構(gòu)網(wǎng)格面上建立雙曲型網(wǎng)格生成方程。通過隱式方法離散，得到線性方程組， 調(diào)用Petsc 庫(kù)[19]中的廣義最小殘差算法（generalized mminimum RESidual，GMRES）求解。本質(zhì)上講，本文的方法是將雙曲型網(wǎng)格生成方程在格點(diǎn)有限體積法的框架下離散并求解[20]。

1 雙曲網(wǎng)格生成方法

本文首先在結(jié)構(gòu)網(wǎng)格的框架下描述三維雙曲型網(wǎng)格生成原理，給出了雙曲網(wǎng)格生成方程及數(shù)值求解方法。然后針對(duì)非結(jié)構(gòu)網(wǎng)格面的格點(diǎn)，重新定義局部坐標(biāo)來近似計(jì)算坐標(biāo)變換導(dǎo)數(shù)，發(fā)展出非結(jié)構(gòu)網(wǎng)格的框架下的雙曲網(wǎng)格生成體系，利用格點(diǎn)有限體積法離散控制方程，并采用迎風(fēng)方法和隱式推進(jìn)求解。

1.1 雙曲型結(jié)構(gòu)網(wǎng)格生成方法

用r=(x,y,z)表 示網(wǎng)格點(diǎn)的坐標(biāo)行向量，而ξ、η、ζ 是貼體廣義坐標(biāo)系（分別使用j、k、l來索引）。三維雙曲型網(wǎng)格生成方程是基于以下正交性條件及體積約束：

應(yīng)用局部線化后，經(jīng)過簡(jiǎn)單的代數(shù)運(yùn)算，可以得到如下向量形式的網(wǎng)格生成方程：

由于物面單元面積均不為零，且單元格點(diǎn)排布方向使面法向朝向推進(jìn)空間一側(cè)，易知P-1存在，于是：

其中：Q?1=P-1Q1，Q?2=P-1Q2，可以證明是實(shí)對(duì)稱矩陣，因而該方程可以看作沿 ζ方向推進(jìn)的雙曲方程。

到目前為止，雙曲網(wǎng)格生成控制方程的推導(dǎo)都是精確的。下一步需要將該方程離散，可以采用中心差分加人工黏性的方法[12]，也可以使用迎風(fēng)格式[13]。這里參考Chan[12]使用中心差分加人工黏性的方法，rξ、rη通 過中心差分計(jì)算，而rζ采用兩點(diǎn)隱式后插的方法：

其中： εiξ是人工指定隱式耗散控制參數(shù)，其大小為顯式耗散控制參數(shù)的兩倍；De為顯式耗散項(xiàng)，具體討論見下文。此外，以 η方向?yàn)槔?/p>

該方程在結(jié)構(gòu)化網(wǎng)格上是五塊對(duì)角矩陣，可以通過近似分裂分別在兩個(gè)方向化為三塊對(duì)角矩陣分別求解。

接下來討論右端耗散項(xiàng):

標(biāo)量函數(shù)Sl是利用某種函數(shù)歸一化的距離壁面的累積距離。距離壁面近的區(qū)域耗散小，網(wǎng)格正交性良好；距離壁面遠(yuǎn)的區(qū)域，一些表面凹區(qū)域生長(zhǎng)出來網(wǎng)格線可能開始相交，需要增大耗散避免網(wǎng)格相交；距離壁面更遠(yuǎn)的地方，則耗散保持不變即可。

1.2 雙曲型非結(jié)構(gòu)網(wǎng)格生成方法

如引言所述，Steger[14]和Matsuno[16]都曾經(jīng)將雙曲型網(wǎng)格生成方法應(yīng)用到非結(jié)構(gòu)網(wǎng)格上。Steger 的方法比較直接，如圖2 所示，在非結(jié)構(gòu)的網(wǎng)格面上指定當(dāng)?shù)刈鴺?biāo)系，然后找出可能存在的網(wǎng)格序號(hào)結(jié)構(gòu)性，進(jìn)而采用結(jié)構(gòu)網(wǎng)格的方法進(jìn)行網(wǎng)格生成。具體來講，使用坐標(biāo)為(1，+1)和(1，-1)的點(diǎn)計(jì)算ξ 方向的坐標(biāo)導(dǎo)數(shù)，使用坐標(biāo)為(2，+1)和(2，-1)的點(diǎn)計(jì)算η方向的坐標(biāo)導(dǎo)數(shù)。但由于這種坐標(biāo)系指定的隨機(jī)性以及網(wǎng)格潛在的結(jié)構(gòu)性連接結(jié)構(gòu)生成困難，Steger 的文章只有一個(gè)圍繞簡(jiǎn)單球體的網(wǎng)格結(jié)果，沒有證明能應(yīng)用于復(fù)雜實(shí)用幾何的外形。

圖2 Steger[14]采用的局部坐標(biāo)系Fig. 2 Local coordinate system used by Steger[14]

Matsuno[16]則在格心建立局部坐標(biāo)系，通過單元和單元的連接關(guān)系指定局部坐標(biāo)系的方向，即使用相鄰單元格心的坐標(biāo)計(jì)算得到，如圖3 所示。這樣同樣將非結(jié)構(gòu)的問題轉(zhuǎn)化為局部結(jié)構(gòu)化問題進(jìn)行雙曲網(wǎng)格生成。下一層格心點(diǎn)確定后，通過簡(jiǎn)單幾何平均得到格點(diǎn)處的網(wǎng)格坐標(biāo)。他的方法，可以視為在非結(jié)構(gòu)網(wǎng)格上構(gòu)造通過將三角形剖分為四邊形，從而構(gòu)造出局部結(jié)構(gòu)網(wǎng)格，開展雙曲網(wǎng)格生成，也可以看為基于格心有限體積法的雙曲網(wǎng)格生成。Matsuno 方法的局限性已在引言中講明。論文中算例的表面網(wǎng)格均為三角形網(wǎng)格，且三角形內(nèi)角沒有扭曲，單個(gè)頂點(diǎn)周圍三角形的數(shù)量均為6 個(gè)。可見他的方法適用性上是有限制的。在某些流動(dòng)特征的數(shù)值模擬中，物面也會(huì)根據(jù)需要生成各向異性的網(wǎng)格，這個(gè)方法就會(huì)無能為力。

圖3 Matsuno[16]采用的局部坐標(biāo)系Fig. 3 Local coordinate system used by Matsuno[16]

本文提出的方法則是完全適應(yīng)于三角形、四邊形混合表面網(wǎng)格，且對(duì)網(wǎng)格分布沒有額外要求。因?yàn)楸疚牡姆椒◤母旧暇褪窃诨诟顸c(diǎn)有限體積法的框架下，對(duì)守恒型式的雙曲方程進(jìn)行推導(dǎo)離散。

將控制方程（6）改為守恒型：

圖4 格點(diǎn)有限體積法控制體幾何示意圖Fig. 4 Schematic of dual control volume for node-centered finitevolume schemes with unit normals associated with an edge with an edge {j,k}

圖5 局部計(jì)算坐標(biāo)變換導(dǎo)數(shù)所用模板示意圖Fig. 5 Grid stencils for coordinate local transformation computation

體積的指定在當(dāng)前層格點(diǎn)控制體面積已知的情況下，就是沿推進(jìn)方向的控制網(wǎng)格分布函數(shù)的指定，可以采用指數(shù)函數(shù)、雙指數(shù)函數(shù)或者雙曲正切函數(shù)，也可以從相似的簡(jiǎn)單參考外形體已有的網(wǎng)格上“拷貝”網(wǎng)格分布。作者所在課題組長(zhǎng)期從事高超聲速外形飛行器的設(shè)計(jì)及優(yōu)化，我們的經(jīng)驗(yàn)是采用球錐模型近似飛行器，根據(jù)近似公式得到激波位置并外擴(kuò)作為遠(yuǎn)場(chǎng)邊界，使用代數(shù)的方法生成參考網(wǎng)格。

為進(jìn)一步提高網(wǎng)格的光滑性要求，對(duì)推進(jìn)層的體積做拉普拉斯光順：

其中：L是拉普拉斯算子，在非結(jié)構(gòu)網(wǎng)格上，使用偽拉普拉斯權(quán)方法計(jì)算[21]； υ是人工指定的體積光順因子。

同樣的，也可以在方程（28）右端加入顯式的耗散項(xiàng)（14）。顯式的耗散項(xiàng)從結(jié)構(gòu)網(wǎng)格往非結(jié)構(gòu)網(wǎng)格上擴(kuò)展時(shí)，唯一需要注意的就是角度感受器中角度的定義。本文參考文獻(xiàn)[22]，如圖6 所示，將角度定義為

圖6 角度感受器中角度的定義示意圖Fig. 6 Illustration of angle definition in angle sensor function sensor function

2 網(wǎng)格生成舉例

采用本文算法對(duì)幾個(gè)典型外形生成了雙曲型網(wǎng)格，驗(yàn)證本文方法的可行性。首先是最基本的球面，圖7 是圓球的表面網(wǎng)格，圖8 是本文方法生成的圓球雙曲型非結(jié)構(gòu)網(wǎng)格。從圖8 可以看出，表面光滑的圓球，由于表面網(wǎng)格點(diǎn)分布的非均勻性，也會(huì)導(dǎo)致生成的網(wǎng)格前沿面并非絕對(duì)球面。隨著前沿面的繼續(xù)推進(jìn)，耗散將不再變化，最終會(huì)在足夠遠(yuǎn)處得到比較圓的外邊界，無論初始外形多么復(fù)雜。

圖7 圓球的表面網(wǎng)格Fig. 7 Surface mesh of sphere

圖8 圓球的雙曲型非結(jié)構(gòu)網(wǎng)格y = 0 切面圖Fig. 8 Hyperbolic mesh slice at y = 0 cutting plane for sphere

第二個(gè)網(wǎng)格生成例子是表面含有凹凸的旋成體，母線函數(shù)為：

其中：a、b是 控制凸凹高度的參數(shù)，k是控制凸凹個(gè)數(shù)的參數(shù)，本文取a=0.2，b=2.5，k=4。

圖9 是含有凹凸結(jié)構(gòu)的旋成體表面網(wǎng)格，圖10是含有凹凸結(jié)構(gòu)的旋成體雙曲型網(wǎng)格y= 0 切面圖，圖11 是在物面附近的放大圖，可以看出，無論是在外形的凸起處，還是在凹陷處，均可生成高質(zhì)量的棱柱網(wǎng)格。圖12 是旋成體雙曲型網(wǎng)格x切面圖，由于旋成體周向網(wǎng)格分布均勻，生成的空間棱柱網(wǎng)格具有良好的對(duì)稱性。

圖9 含有凹凸結(jié)構(gòu)的旋成體表面網(wǎng)格Fig. 9 Surface mesh of the spiral body with extreme concave-convex structure

圖10 含有凹凸結(jié)構(gòu)的旋成體雙曲型網(wǎng)格y = 0 切面圖Fig. 10 Hyperbolic mesh slice at y = 0 cutting plane of the spiral body with extreme concave-convex structure

圖11 旋成體凹凸結(jié)構(gòu)的處的雙曲型網(wǎng)格y = 0 切面放大圖Fig. 11 Close up view of the hyperbolic mesh slice at y = 0 cutting plane of the spiral body with extreme concave-convex structure

圖12 含有凹凸結(jié)構(gòu)的旋成體雙曲型網(wǎng)格x 切面圖Fig. 12 Hyperbolic mesh slice at x cutting plane of the spiral body with extreme concave-convex structure

接下來是測(cè)試本文發(fā)展的雙曲網(wǎng)格生成方法在含有尖角或尖銳邊緣外形的適用情況，選取正方體外形進(jìn)行測(cè)試。圖13 是立方體表面網(wǎng)格，在結(jié)構(gòu)網(wǎng)格的基礎(chǔ)上通過擾動(dòng)并將四邊形劃分為三角形得到。圖14、圖15 是本文方法生成的空間棱柱網(wǎng)格兩個(gè)視角下的切面圖，可以看出得到的網(wǎng)格結(jié)果令人滿意。

圖13 立方體表面網(wǎng)格Fig. 13 Surface mesh of the cube with sharp edges

圖14 立方體雙曲型網(wǎng)格y 切面圖Fig. 14 Hyperbolic mesh slice at y const cutting plane of the cube

圖15 立方體雙曲型網(wǎng)格斜切面圖Fig. 15 Oblique view of hyperbolic mesh slice at x+y = const cutting plane of the cube

最后是復(fù)雜飛行器外形X38。圖16 為X38 外形表面網(wǎng)格，根據(jù)流動(dòng)特征進(jìn)行了加密。鄰居點(diǎn)的最大數(shù)量達(dá)到12，即在局部表面網(wǎng)格比較扭曲，三角形內(nèi)角小于30°。表面一共有106 138 個(gè)點(diǎn)，212 272 個(gè)三角形網(wǎng)格。圖17、圖18 分別展示了x切面和z切面下內(nèi)部網(wǎng)格的正交性，圖19、圖20、圖21 分別展示了幾個(gè)l=const時(shí)的前沿面。如圖17 所示，本文發(fā)展的網(wǎng)格生成方法在壁面附近和遠(yuǎn)場(chǎng)生成的網(wǎng)格均保持了法向正交性。但在機(jī)身和兩翼間形成的凹區(qū)域，網(wǎng)格極度聚集，盡管網(wǎng)格沒有相交，沒有產(chǎn)生負(fù)體積的單元，但極薄的三棱柱大大降低了網(wǎng)格質(zhì)量。

圖16 X38 飛行器表面網(wǎng)格示意圖Fig. 16 Surface mesh of X38 Crew-Return-Vehicle

圖17 X38 雙曲型網(wǎng)格x 切面圖Fig. 17 Hyperbolic mesh slice at x const cutting plane of X38 CRV

圖18 X38 雙曲型網(wǎng)格z 切面圖Fig. 18 Hyperbolic mesh slice at z const cutting plane of X38 CRV

圖19 l =20時(shí)網(wǎng)格前沿面示意圖Fig. 19 Generated Hyperbolic frontier mesh at l=20

圖20 l =30時(shí)網(wǎng)格前沿面示意圖Fig. 20 Generated Hyperbolic frontier mesh atl=30

圖21 l =40時(shí)網(wǎng)格前沿面示意圖Fig. 21 Generated Hyperbolic frontier mesh at l=40

3 數(shù)值流場(chǎng)計(jì)算

本文基于完全氣體模型求解三維層流N-S 方程，使用格心有限體積法的隱式離散化，在通量計(jì)算上，利用本文網(wǎng)格的方向性，使用混合無黏通量格式。所采用的混合策略是，在壁面平行的面元上使用耗散足夠小、能精確捕捉邊界層的Godunov 格式，在其他面元上使用耗散大、能抑制激波不穩(wěn)定的旋轉(zhuǎn)迎風(fēng)格式。時(shí)間推進(jìn)上采用線隱LUSGS 求解。

3.1 X38 外形氣動(dòng)力數(shù)值模擬

X38 飛機(jī)是為國(guó)際空間站研制的升力體再入式乘員往返飛行器的原型機(jī)，作為宇航員緊急逃逸裝置使用，該模型主要用來考察本文發(fā)展的外形生成方法及數(shù)值仿真方法對(duì)型號(hào)外形基本氣動(dòng)力的模擬能力。

計(jì)算條件為：Ma=6，Re=2.275×105/m，來流靜溫T=216.7 K，攻角為20°～40°。壁面法向第一層網(wǎng)格雷諾數(shù)為10。

如前所述，表面非結(jié)構(gòu)網(wǎng)格在雙曲網(wǎng)格生成中，生成的體網(wǎng)格在局部凹陷的區(qū)域會(huì)極度聚集。需要發(fā)展更合適的耗散參數(shù)感受器來進(jìn)一步控制網(wǎng)格的質(zhì)量。因此本文對(duì)X38 外形進(jìn)行數(shù)值計(jì)算時(shí)，使用的是混合網(wǎng)格，即在網(wǎng)格表面生成棱柱網(wǎng)格，遠(yuǎn)場(chǎng)生成四面體網(wǎng)格。

為了說明本文方法生成的附面層網(wǎng)格的優(yōu)勢(shì)，將本文方法生成的附面層網(wǎng)格與使用軟件NNW-Gridstar生成的附面層網(wǎng)格進(jìn)行了對(duì)比。圖22 為本文方法生成的X38 飛行器附面層網(wǎng)格切面圖。圖23 為Gridstar生成的X38 飛行器附面層網(wǎng)格切面圖?？梢钥闯鲈诖蟛糠止饣锩妫瑑煞N方法生成附面層網(wǎng)格都具備良好的正交性。而在機(jī)身和機(jī)翼之間的凹區(qū)域，以及在凹凸關(guān)聯(lián)的區(qū)域，Gridstar 為了避免棱柱網(wǎng)格相交，在局部停止生成棱柱網(wǎng)格，而本文的方法在這兩個(gè)區(qū)域依然能夠生成良好的棱柱網(wǎng)格。換而言之，本文方法生成棱柱網(wǎng)格，從外形表面到附面層網(wǎng)格外緣，每一條“串”上具有相等數(shù)量的棱柱網(wǎng)格，這對(duì)線隱格式的加速求解具有重要意義。

圖22 本文方法生成的X38 飛行器附面層網(wǎng)格x 切面圖Fig. 22 Cut-out view of the prismatic mesh of X38 generated by proposed method

圖23 Gridstar 生成的X38 飛行器附面層網(wǎng)格x 切面圖Fig. 23 Cut-out view of the prismatic mesh of X38 generated by Gridstar software

圖24 為本文計(jì)算所采用的混合網(wǎng)格，壁面法向第一層網(wǎng)格雷諾數(shù)為10，法向增長(zhǎng)率為1.1，棱柱網(wǎng)格共30 層。遠(yuǎn)場(chǎng)為四面體網(wǎng)格。圖25 為X38 在攻角20°時(shí)有量綱壓力場(chǎng)云圖。圖26 為升力系數(shù)隨攻角變化計(jì)算結(jié)果與實(shí)驗(yàn)值的對(duì)比圖，兩者符合得比較好。

圖24 X38 飛行器混合網(wǎng)格z 切面圖Fig. 24 Cut-out view of the hybrid mesh of X38

圖25 X38 攻角20°有量綱壓力場(chǎng)云圖Fig. 25 Pressure contour of the X38 at α = 20°

圖26 X38 升力系數(shù)隨攻角的變化曲線Fig. 26 Comparison of experimental and CFD results on lift coefficient with angle of attack

3.2 雙橢球外形氣動(dòng)熱數(shù)值模擬

高超聲速飛行器常采用簡(jiǎn)單組合體外形，相貫雙橢球體是其中比較典型的一類，如美國(guó)航天飛機(jī)頭部就具有典型雙橢球特征。國(guó)內(nèi)外對(duì)此外形開展的風(fēng)洞實(shí)驗(yàn)研究較多，表面熱流密度分布及流場(chǎng)內(nèi)波系變化等實(shí)驗(yàn)數(shù)據(jù)較豐富，因此是研究高超聲速大氣再入問題的標(biāo)準(zhǔn)算例之一，以此可以考察本文方法生成的網(wǎng)格對(duì)相對(duì)復(fù)雜外形的氣動(dòng)熱預(yù)測(cè)能力。本文對(duì)0°攻角的雙橢球高超聲速繞流進(jìn)行了數(shù)值模擬。

計(jì)算條件為：Ma=10.02，Re=2.2×106/m，總壓P0=6.9 MPa， 總溫T0=1 457 K。壁面法向第一層網(wǎng)格雷諾數(shù)為1。圖27 是本文計(jì)算雙橢球的物面網(wǎng)格。圖28 是本文方法生成的雙橢球空間雙曲型網(wǎng)格。圖29 是對(duì)稱面和物面上的壓力分布云圖，具有清晰的流場(chǎng)和激波分辨率。圖30 為壁面熱流分布云圖，具有較好的光滑性。圖31 為上下兩條中心線的熱流分布，與實(shí)驗(yàn)值吻合較好，表明能夠模擬出二次分離對(duì)熱流分布的影響。

圖27 雙橢球表面網(wǎng)格Fig. 27 Surface mesh of the double ellipsoid

圖28 雙橢球雙曲型網(wǎng)格Fig. 28 3D mesh of the hyperbolic double ellipsoid

圖29 雙橢球模型無量綱壓力場(chǎng)云圖Fig. 29 Pressure contour of the double ellipsoid

圖30 雙橢球模型熱流場(chǎng)云圖Fig. 30 Heat flux contour of the double ellipsoid

圖31 雙橢球中心對(duì)稱線熱流分布Fig. 31 Heat flux distributions of upper and lower symmetric lines

4 結(jié) 論

本文提出了一種基于格點(diǎn)有限體積法的非結(jié)構(gòu)雙曲網(wǎng)格生成的方法。建立了基于非結(jié)構(gòu)網(wǎng)格面格點(diǎn)的局部坐標(biāo)系，以在非結(jié)構(gòu)網(wǎng)格下計(jì)算坐標(biāo)變換導(dǎo)數(shù)。通過給控制方程右端添加人工耗散，并利用偽拉普拉斯權(quán)進(jìn)行光順，很好地控制非結(jié)構(gòu)雙曲型網(wǎng)格的自動(dòng)生成。

本文發(fā)展的方法，在光滑的表面上、含有深凹高凸的外形上、含有尖銳邊緣外形上，以及真實(shí)復(fù)雜三維高超聲速飛行器外形上，均可自動(dòng)生成良好的棱柱網(wǎng)格。在復(fù)雜外形上，亦可生成附面層網(wǎng)格，遠(yuǎn)場(chǎng)使用四面體網(wǎng)格。通過對(duì)X38 外形的氣動(dòng)力數(shù)值模擬和對(duì)雙橢球外形的氣動(dòng)熱數(shù)值模擬，結(jié)果表明本文發(fā)展的雙曲型網(wǎng)格可以滿足高超聲速流場(chǎng)氣動(dòng)力熱數(shù)值仿真的需求。

下一步工作，可在以下方面展開：

1）在每一層推進(jìn)過程中，發(fā)展適用于非結(jié)構(gòu)網(wǎng)格面的橢圓型光順方法，特別是針對(duì)有尖角的三維外形。

2）將本文發(fā)展的方法，進(jìn)一步應(yīng)用于含三角形、四邊形物面網(wǎng)格的空間雙曲型網(wǎng)格生成。

3）進(jìn)一步發(fā)展適用于非結(jié)構(gòu)網(wǎng)格的耗散控制感應(yīng)器，提高雙曲網(wǎng)格生成方法在非結(jié)構(gòu)網(wǎng)格生成中的魯棒性。

4）將雙曲非結(jié)構(gòu)網(wǎng)格生成方法與交互的棱柱網(wǎng)格生成算法、空間四面體網(wǎng)格生成算法、低質(zhì)量網(wǎng)格重構(gòu)方法或重疊網(wǎng)格技術(shù)結(jié)合，發(fā)展適用于CFD 計(jì)算的混合網(wǎng)格自動(dòng)生成技術(shù)。