網格生長動力學模型的建立和求解及其混合黏性網格生成中的應用1)

任一鵬 劉 楓 林 其 胡華偉

?(北京宇航系統工程研究所，北京100076)

?(中國空氣動力研究與發展中心，四川綿陽621000)

??(成都銘峰新源科技有限公司，成都610042)

真實飛行器外形復雜，所涉及的流動結構和現象對數值模擬中所需的計算網格提出了較高的要求。一方面，飛行器本身構型多樣，各種部件無論從尺度還是位置都存在著多種變化和組合，將傳統的參數化的網格生成方法應用其中存在重大障礙；另一方面，真實飛行條件多變，所產生的流動結構性質各異，對計算網格的要求也不盡相同。在流動的數值模擬中，一個熟知的原則就是網格分布應該盡可能符合當地的流動狀態。因此，建立能夠適應復雜外形的網格生成方法是必要的。

國內外一些學者在復雜外形的網格生成方法[1-2]上進行了諸多探索。在過去的二十年中，許多重要的網格生成技術得到建立和發展，比如多塊對接/拼接結構網格、重疊/嵌套網格、非結構(四面體/金字塔)網格等。近年來混合網格(包含有多種形態的網格)越來越受到關注，得到了較快的發展。Zhang等[3-4]建立了基于各向異性四面體凝聚的混合網格生成方法，并成功運用于大雷諾數的復雜飛行器數值模擬中；Esquieu[5]和Kallinderis等[6]建立了三棱柱/四面體混合網格生成算法，并運用于噴氣式客機繞流流場的計算；Timothy[7]建立了三棱柱/金字塔/四面體/六面體混合網格生成方法；Kannan等[8]和Wang[9]建立了自適應笛卡爾網格生成算法。眾所周知，邊界層網格的質量對計算結果有重大影響。三棱柱網格作為一種各向異性網格可以任意沿其法向拉伸或者壓縮而不影響其形狀，能夠滿足邊界層計算的需要。然而，傳統的三棱柱網格生成算法(層推進算法及其改進[10-11])在凹凸相鄰或相交區域都遇到了較大的困難，因此難以運用于復雜三維外形的網格生成中。本文試圖建立網格生長動力學模型，通過迭代求解動態優化調整網格點生長法向和長度，形成邊界層網格生成算法，然后結合各向同性四面體網格的陣面推進Delaunay網格細化算法，建立一種新型的混合黏性網格生成方法。

1 網格的類型和性質

根據網格的幾何特征，可以將網格類型分為各向同性網格和各向異性網格，如圖1所示。眾所周知，不同流動結構對網格尺度的要求是不同的，在流動變化劇烈的方向需要加密網格。因此，各向同性網格適合布置在流動梯度較小、物理量分布較光滑的區域；各向異性網格由于網格幾何特征決定了網格方向性的差異，不同幾何方向可以存在較大的尺寸比，可以用于流動梯度較大、物理量分布存在間斷的區域。

圖1 網格類型圖

2 網格生長動力學模型

2.1 網格生長受力分析

如圖2所示，根據生長邊界的幾何特性的不同，一般可以分為凸殼生長、非凸殼生長以及二者的組合。這三種生長類型都試圖將網格的法向“光順化”，由法向變化劇烈的最底層網格逐漸過渡到法向光滑的頂層網格。這種“光順化”在二維表現為圓弧，在三維情況表現為球面。

圖2 網格生長示意圖

假設網格點之間存在力的相互作用，并進一步假設作用力滿足胡克定律，即

式中，k為彈性系數，為網格點之間的矢徑。

不失一般性，分別考慮凸殼和非凸殼中法向變化最劇烈的角點受力情況，如圖3和圖4所示。顯然，在最底層網格中，除角點外的其他網格點所受合力為0。凸殼中角點受力沿著網格生長的方向，合力使得網格生長加速；非凸殼中角點受力沿著網格生長相反的方向，合力會抑制網格生長。這就使得凸殼中的角點在光順化條件下，生長總距離較其他網格點長，網格生長速度較快；在非凸殼中反之。一般來說，網格生成成敗的關鍵在于“角點”的處理，網格交叉、網格扭曲最嚴重的區域也在于此。

圖3 凸殼角點受力分析

圖4 非凸殼角點受力分析

2.2 網格生長動力學模型

如圖5所示，建立O?xyz笛卡爾網格坐標系。由上節分析可知，對于無論凸殼或是非凸殼情況，網格點上的合力是網格生長的動力，這種合力作用隨著網格生長而改變。由牛頓第二定律，即有

圖5 網格坐標系

不妨假設在生長過程中網格點作用合力不變，對式(2)積分，顯然有

其中，c0，c1為常數。隨著時間推進，|r|呈拋物線快速增長或者減小。這種生長模式簡單直接，然而由于網格尺度增加或縮減速度過快，對于復雜外形極易出現網格交叉和負體積，導致網格生成失敗。根據對圖1和圖2的分析，如果引入某種阻尼，對網格生長速度進行動態抑制，則有

式中，Kv為阻尼函數。初始條件為

式中，r0為初始位置向量，n0為初始法向量。

對整個網格生長陣面而言，則有

式(6)是一個N維二階常微分非線性系統。

2.3 模型求解

動力學系統公式(6)存在強烈的非線性且相互耦合，難以直接求得解析解。然而，我們知道在每一層網格陣面進行推進，特別是生成邊界層網格時，由于邊界層尺度η遠小于飛行器的特征尺度L，不妨假設在每一層網格推進過程中，各網格點生長過程相互獨立，則有

式(7)是一個二階非齊次常微分方程。不失一般性，去掉方程中的下標i，式(7)對應的齊次方程為

式(8)對應的特征方程是

方程(9)有兩個不相等的實數解M1,2=0,?Kv，則齊次方程(8)的通解形式為

不難找出

是方程(7)的一個特解。根據常微分方程解析理論，非齊次方程的通解為對應齊次方程的通解與非齊次方程的一個特解之和，則方程(7)的通解為

代入初始條件(5)，則有

將方程(13)代入方程(12)中，則方程(7)的解為

將初始條件式(5)代入式(3)中，則有

對比式(14)和式(15)，不難看出兩種不同的生長模式對網格生長高度和方向的影響。下面分別考慮凸殼和非凸殼中網格生長隨時間的變化模式。在凸殼(如圖3)中，令則圖6和圖7分別給出了凸殼中模式一(阻尼模式)和模式二(非阻尼模式)中網格生長高度隨時間的變化。在非凸殼(如圖4)中，令則圖8和圖9分別給出了非凸殼中模式一(阻尼模式)和模式二(非阻尼模式)中網格生長高度隨時間變化。可以看出，阻尼模式會抑制網格生長高度的過快增長，從而增強網格生成的魯棒性。

圖6 凸殼角點網格生長高度r x隨時間變化

圖7 凸殼角點網格生長高度r y隨時間變化

圖8 非凸殼角點網格生長高度r x隨時間變化

圖9 非凸殼角點網格生長高度r y隨時間變化

2.4 法向量求解與阻尼函數選取

法向量存在于網格推進陣面的面心處，而網格格點的法向量有

S i為格點p相鄰網格陣面的面積，n i F為相鄰網格陣面的法向量。然而，網格格點法向量與周圍格點法向量密切相關且相互作用，本文取周圍格點法向量的平均值，則有

式(16)等價于線性系統

其中

線性系統(17)對應的放大因子顯然小于1，因此該線性系統可以通過Jacobi迭代或者Gauss–Seidel迭代求解。迭代求解初值按照式(16)求得。

根據式(4)的定義，阻尼函數的選取應該與網格生長距離l的模成反比，即

3 各向異性網格生成算法

上節建立和求解了網格生長動力學模型，運用該模型可以得到各向異性網格生成算法，歸結起來可以分為以下幾步：

步驟1：計算每層網格的預計生長高度。每層網格預計生長高度按指數增加計算，即

式中，α為層增長率，β為預設網格層數，h0為第一層網格高度。

步驟2：計算網格推進陣面上各個網格點的法向矢量。根據式(16)計算法向矢量初值，求解式(17)得到各個網格點的法向矢量。

步驟3：計算當前時間步中網格推進陣面上各個網格點的阻尼函數Kv和所受合力F。分別根據式(1)和式(19)計算合力F和阻尼函數Kv。

步驟4：根據式(14)計算每個時間步中生長矢量r。

步驟5：判斷當前推進陣面中任意網格點的生長高度是否達到該網格點的預計生長高度。如果是，則轉入下一層網格生長；如果否，則返回步驟2繼續時間推進，直到當前推進陣面中任意網格點的生長高度達到該網格點的預計生長高度。

4 各向同性四面體生成算法

各向同性四面體生成算法，國內外已經有眾多學者進行了廣泛而深入的研究，已經相對成熟。各類算法的主要區別在于網格生成質量和效率的差異。陣面推進法的特點是搜索量大，每次搜索只能完成一個網格的生成，網格生成效率較低，生成網格不一定能滿足外接球準則即Delaunay準則；而Delaunay方法因其生成網格的效率較高，而且網格質量高，天然滿足Delaunay準則，然而需要高質量的邊界恢復算法才能保形。本文基于計算幾何的基本原理，建立了基于可視面的快速四面體初始化算法，然后建立了陣面推進Delaunay網格細化算法，并結合網格優化算法，發展了一套任意多面體的快速非結構網格生成算法。該算法一方面克服了陣面推進法網格生成效率低的缺點，另一方面克服了Delaunay方法在邊界恢復中不保形的缺陷。詳細算法描述及定義參考文獻[12-13]。

其主要步驟包括：

步驟1：通過建立的基于可視面的初始化四面體算法，將網格剖分區域進行初始化，形成初始四面體陣列；

步驟2：然后通過運用改進的最差面細化算法將網格細化，其間運用Delaunay面交換算法優化網格質量；

步驟3：求解基于頂點彈簧的線性系統進一步提高網格質量。

5 算例驗證

5.1 雙橢球算例

雙橢球算例是高超聲速風洞試驗的標準模型，外形由兩個不同形狀的橢球疊加而成。圖10給出了雙橢球外形的表面網格，共由14102個三角形構成。三棱柱網格共10層，網格量為141020。圖11～圖13分別給出了頭部、拐角和底部轉角處的局部網格，可以看出網格無交叉且光順性較好。

圖10 表面網格

圖11 頭部網格

圖13 底部轉角網格

5.2 空天飛機算例

空天飛機算例是翼身融合體模型，外形中有數個凹凸部件相連的區域，對網格生成算法的魯棒性有很高的要求。圖14給出了空天飛機外形的表面網格，共由80606個三角形構成。三棱柱網格共2層，網格量為1612120。圖15～圖17分別給出了不同截面的局部網格，可以看出網格光順性較好。

圖12 拐角網格

圖14 表面網格

圖15 截面x/L=0.25網格

圖17 截面y/L=0.2網格

圖16 截面x/L=0.9網格

6 結論

(1)本文提出并求解的網格生長動力學模型，可以用于黏性網格的邊界層網格生成中，具有較高的可靠性；

(2)基于各向異性網格生成算法和各向同性四面體網格生成算法的網格生成策略，可以有效地運用于飛行器的網格生成中。