三維復雜外形的雙曲型非結構網格生成方法

張啟明，葉友達，*，蔣勤學，田 浩

(1. 北京流體動力科學研究中心，北京 100120；2. 國家計算流體力學實驗室，北京 100191)

0 引 言

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

棱柱網格由于在某幾個方面具有特殊的優勢，因而廣泛地應用在邊界層網格的生成中：從網格生成的角度，物面三角化網格靈活性強，適用性高，邊界層內棱柱網格生成自動化程度高，當面對復雜的外形時，可以大大減少人工交互；從流動機理上看，靠近物面的黏性區域在垂直物面方向上的網格具有方向性，在高雷諾數流動或者含湍流的流動中對黏性邊界層的分辨至關重要；在數值方法上，在遠離物面的方向，網格具有結構化編號，這種內在的網格結構可以應用到通量格式中以提高計算精度，也可以應用到隱式時間推進中以加速收斂。

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

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

串網格的生成，早期采用非結構網格邊界層內棱柱網格生成的方法。棱柱網格的生成，可以分為基于幾何的推進方法、基于各向異性四面體聚合的方法和基于場的推進方法。基于幾何的方法，根據當地幾何特征，確定推進方向和推進步長，重點解決在模型凹特征區域和鄰近特征區域附近的邊界層前沿面上可能出現的相交問題。基于網格聚合的方法[2]將自動生成的各向異性四面體網格基于一定的判定準則進行聚合，得到物面附近的棱柱網格單元。

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

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

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

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

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

1 雙曲網格生成方法

本文首先在結構網格的框架下描述三維雙曲型網格生成原理，給出了雙曲網格生成方程及數值求解方法。然后針對非結構網格面的格點，重新定義局部坐標來近似計算坐標變換導數，發展出非結構網格的框架下的雙曲網格生成體系，利用格點有限體積法離散控制方程，并采用迎風方法和隱式推進求解。

1.1 雙曲型結構網格生成方法

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

應用局部線化后，經過簡單的代數運算，可以得到如下向量形式的網格生成方程：

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

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

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

其中： εiξ是人工指定隱式耗散控制參數，其大小為顯式耗散控制參數的兩倍；De為顯式耗散項，具體討論見下文。此外，以 η方向為例：

該方程在結構化網格上是五塊對角矩陣，可以通過近似分裂分別在兩個方向化為三塊對角矩陣分別求解。

接下來討論右端耗散項:

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

1.2 雙曲型非結構網格生成方法

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

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

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

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

本文提出的方法則是完全適應于三角形、四邊形混合表面網格，且對網格分布沒有額外要求。因為本文的方法從根本上就是在基于格點有限體積法的框架下，對守恒型式的雙曲方程進行推導離散。

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

圖4 格點有限體積法控制體幾何示意圖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 局部計算坐標變換導數所用模板示意圖Fig. 5 Grid stencils for coordinate local transformation computation

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

為進一步提高網格的光滑性要求，對推進層的體積做拉普拉斯光順：

其中：L是拉普拉斯算子，在非結構網格上，使用偽拉普拉斯權方法計算[21]； υ是人工指定的體積光順因子。

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

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

2 網格生成舉例

采用本文算法對幾個典型外形生成了雙曲型網格，驗證本文方法的可行性。首先是最基本的球面，圖7 是圓球的表面網格，圖8 是本文方法生成的圓球雙曲型非結構網格。從圖8 可以看出，表面光滑的圓球，由于表面網格點分布的非均勻性，也會導致生成的網格前沿面并非絕對球面。隨著前沿面的繼續推進，耗散將不再變化，最終會在足夠遠處得到比較圓的外邊界，無論初始外形多么復雜。

圖7 圓球的表面網格Fig. 7 Surface mesh of sphere

圖8 圓球的雙曲型非結構網格y = 0 切面圖Fig. 8 Hyperbolic mesh slice at y = 0 cutting plane for sphere

第二個網格生成例子是表面含有凹凸的旋成體，母線函數為：

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

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

圖9 含有凹凸結構的旋成體表面網格Fig. 9 Surface mesh of the spiral body with extreme concave-convex structure

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

圖11 旋成體凹凸結構的處的雙曲型網格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 含有凹凸結構的旋成體雙曲型網格x 切面圖Fig. 12 Hyperbolic mesh slice at x cutting plane of the spiral body with extreme concave-convex structure

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

圖13 立方體表面網格Fig. 13 Surface mesh of the cube with sharp edges

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

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

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

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

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

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

圖19 l =20時網格前沿面示意圖Fig. 19 Generated Hyperbolic frontier mesh at l=20

圖20 l =30時網格前沿面示意圖Fig. 20 Generated Hyperbolic frontier mesh atl=30

圖21 l =40時網格前沿面示意圖Fig. 21 Generated Hyperbolic frontier mesh at l=40

3 數值流場計算

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

3.1 X38 外形氣動力數值模擬

X38 飛機是為國際空間站研制的升力體再入式乘員往返飛行器的原型機，作為宇航員緊急逃逸裝置使用，該模型主要用來考察本文發展的外形生成方法及數值仿真方法對型號外形基本氣動力的模擬能力。

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

如前所述，表面非結構網格在雙曲網格生成中，生成的體網格在局部凹陷的區域會極度聚集。需要發展更合適的耗散參數感受器來進一步控制網格的質量。因此本文對X38 外形進行數值計算時，使用的是混合網格，即在網格表面生成棱柱網格，遠場生成四面體網格。

為了說明本文方法生成的附面層網格的優勢，將本文方法生成的附面層網格與使用軟件NNW-Gridstar生成的附面層網格進行了對比。圖22 為本文方法生成的X38 飛行器附面層網格切面圖。圖23 為Gridstar生成的X38 飛行器附面層網格切面圖。可以看出在大部分光滑物面，兩種方法生成附面層網格都具備良好的正交性。而在機身和機翼之間的凹區域，以及在凹凸關聯的區域，Gridstar 為了避免棱柱網格相交，在局部停止生成棱柱網格，而本文的方法在這兩個區域依然能夠生成良好的棱柱網格。換而言之，本文方法生成棱柱網格，從外形表面到附面層網格外緣，每一條“串”上具有相等數量的棱柱網格，這對線隱格式的加速求解具有重要意義。

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

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

圖24 為本文計算所采用的混合網格，壁面法向第一層網格雷諾數為10，法向增長率為1.1，棱柱網格共30 層。遠場為四面體網格。圖25 為X38 在攻角20°時有量綱壓力場云圖。圖26 為升力系數隨攻角變化計算結果與實驗值的對比圖，兩者符合得比較好。

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

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

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

3.2 雙橢球外形氣動熱數值模擬

高超聲速飛行器常采用簡單組合體外形，相貫雙橢球體是其中比較典型的一類，如美國航天飛機頭部就具有典型雙橢球特征。國內外對此外形開展的風洞實驗研究較多，表面熱流密度分布及流場內波系變化等實驗數據較豐富，因此是研究高超聲速大氣再入問題的標準算例之一，以此可以考察本文方法生成的網格對相對復雜外形的氣動熱預測能力。本文對0°攻角的雙橢球高超聲速繞流進行了數值模擬。

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

圖27 雙橢球表面網格Fig. 27 Surface mesh of the double ellipsoid

圖28 雙橢球雙曲型網格Fig. 28 3D mesh of the hyperbolic double ellipsoid

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

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

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

4 結 論

本文提出了一種基于格點有限體積法的非結構雙曲網格生成的方法。建立了基于非結構網格面格點的局部坐標系，以在非結構網格下計算坐標變換導數。通過給控制方程右端添加人工耗散，并利用偽拉普拉斯權進行光順，很好地控制非結構雙曲型網格的自動生成。

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

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

1）在每一層推進過程中，發展適用于非結構網格面的橢圓型光順方法，特別是針對有尖角的三維外形。

2）將本文發展的方法，進一步應用于含三角形、四邊形物面網格的空間雙曲型網格生成。

3）進一步發展適用于非結構網格的耗散控制感應器，提高雙曲網格生成方法在非結構網格生成中的魯棒性。

4）將雙曲非結構網格生成方法與交互的棱柱網格生成算法、空間四面體網格生成算法、低質量網格重構方法或重疊網格技術結合，發展適用于CFD 計算的混合網格自動生成技術。