網(wǎng)格生長動力學(xué)模型的建立和求解及其混合黏性網(wǎng)格生成中的應(yīng)用1)

任一鵬 劉 楓 林 其 胡華偉

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

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

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

真實飛行器外形復(fù)雜，所涉及的流動結(jié)構(gòu)和現(xiàn)象對數(shù)值模擬中所需的計算網(wǎng)格提出了較高的要求。一方面，飛行器本身構(gòu)型多樣，各種部件無論從尺度還是位置都存在著多種變化和組合，將傳統(tǒng)的參數(shù)化的網(wǎng)格生成方法應(yīng)用其中存在重大障礙；另一方面，真實飛行條件多變，所產(chǎn)生的流動結(jié)構(gòu)性質(zhì)各異，對計算網(wǎng)格的要求也不盡相同。在流動的數(shù)值模擬中，一個熟知的原則就是網(wǎng)格分布應(yīng)該盡可能符合當(dāng)?shù)氐牧鲃訝顟B(tài)。因此，建立能夠適應(yīng)復(fù)雜外形的網(wǎng)格生成方法是必要的。

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

1 網(wǎng)格的類型和性質(zhì)

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

圖1 網(wǎng)格類型圖

2 網(wǎng)格生長動力學(xué)模型

2.1 網(wǎng)格生長受力分析

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

圖2 網(wǎng)格生長示意圖

假設(shè)網(wǎng)格點之間存在力的相互作用，并進一步假設(shè)作用力滿足胡克定律，即

式中，k為彈性系數(shù)，為網(wǎng)格點之間的矢徑。

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

圖3 凸殼角點受力分析

圖4 非凸殼角點受力分析

2.2 網(wǎng)格生長動力學(xué)模型

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

圖5 網(wǎng)格坐標系

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

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

式中，Kv為阻尼函數(shù)。初始條件為

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

對整個網(wǎng)格生長陣面而言，則有

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

2.3 模型求解

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

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

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

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

不難找出

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

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

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

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

對比式(14)和式(15)，不難看出兩種不同的生長模式對網(wǎng)格生長高度和方向的影響。下面分別考慮凸殼和非凸殼中網(wǎng)格生長隨時間的變化模式。在凸殼(如圖3)中，令則圖6和圖7分別給出了凸殼中模式一(阻尼模式)和模式二(非阻尼模式)中網(wǎng)格生長高度隨時間的變化。在非凸殼(如圖4)中，令則圖8和圖9分別給出了非凸殼中模式一(阻尼模式)和模式二(非阻尼模式)中網(wǎng)格生長高度隨時間變化?？梢钥闯觯枘崮Ｊ綍种凭W(wǎng)格生長高度的過快增長，從而增強網(wǎng)格生成的魯棒性。

圖6 凸殼角點網(wǎng)格生長高度r x隨時間變化

圖7 凸殼角點網(wǎng)格生長高度r y隨時間變化

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

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

2.4 法向量求解與阻尼函數(shù)選取

法向量存在于網(wǎng)格推進陣面的面心處，而網(wǎng)格格點的法向量有

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

式(16)等價于線性系統(tǒng)

其中

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

根據(jù)式(4)的定義，阻尼函數(shù)的選取應(yīng)該與網(wǎng)格生長距離l的模成反比，即

3 各向異性網(wǎng)格生成算法

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

步驟1：計算每層網(wǎng)格的預(yù)計生長高度。每層網(wǎng)格預(yù)計生長高度按指數(shù)增加計算，即

式中，α為層增長率，β為預(yù)設(shè)網(wǎng)格層數(shù)，h0為第一層網(wǎng)格高度。

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

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

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

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

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

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

其主要步驟包括：

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

步驟2：然后通過運用改進的最差面細化算法將網(wǎng)格細化，其間運用Delaunay面交換算法優(yōu)化網(wǎng)格質(zhì)量；

步驟3：求解基于頂點彈簧的線性系統(tǒng)進一步提高網(wǎng)格質(zhì)量。

5 算例驗證

5.1 雙橢球算例

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

圖10 表面網(wǎng)格

圖11 頭部網(wǎng)格

圖13 底部轉(zhuǎn)角網(wǎng)格

5.2 空天飛機算例

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

圖12 拐角網(wǎng)格

圖14 表面網(wǎng)格

圖15 截面x/L=0.25網(wǎng)格

圖17 截面y/L=0.2網(wǎng)格

圖16 截面x/L=0.9網(wǎng)格

6 結(jié)論

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

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