孫巖，江盟，孟德虹，龐宇飛

中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

近幾十年，伴隨著計算機硬件與計算方法的快速發(fā)展，計算流體力學(Computational Fluid Dynamics, CFD)模擬能力迅速提升，已經能夠在部分階段代替試驗手段，成為飛行器快速精細化設計的重要工具，從根本上改變航空航天飛行器的設計模式[1-2]。

目前主流CFD方法基于網格對計算區(qū)域進行空間離散，然后迭代求解離散后的代數方程得到流場的數值解，并通過表面積分獲取飛行器的宏觀氣動載荷特性。因此，網格是開展CFD分析的前提和基礎，直接影響流場數值模擬的精準度與效率[3]。

CFD計算網格按照相鄰節(jié)點之間的拓撲關系可以分為結構網格和非結構網格兩類[4]。結構網格內部節(jié)點的相鄰節(jié)點數量是固定的，可以直接通過數組的索引序號來定義網格節(jié)點間的連接關系，而非結構網格內部節(jié)點的相鄰節(jié)點數量是變化的，需要通過單獨的數據結構對節(jié)點間的連接關系進行定義。相比于非結構網格，結構網格在效率、精度、分辨率和程序實現上均有明顯的優(yōu)勢。但針對復雜的工程外形，結構網格生成十分困難，需要投入大量的人力開展結構網格的拓撲設計[5]。而隨著計算機能力的快速提升，非結構網格相比結構網格的計算效率短板逐漸縮小，且在外形靈活性、生成速度方面的優(yōu)勢逐漸凸顯，已經成為主流商業(yè)軟件和部分in-house代碼的首選網格類型[6]。

非結構網格單元類型主要有二維的三角形和三維的四面體，均可以實現任意離散區(qū)域的自動化填充，但在模擬梯度占主導的邊界層流動時表現糟糕[7]。針對邊界層流動問題，棱柱/四面體混合網格是目前CFD模擬中主要采用的解決方案，其中半結構化的棱柱網格能夠彌補非結構三角形/四面體單元模擬黏性流動的不足，具有更好的流動分辨率。但復雜外形的棱柱網格生成仍然是具有挑戰(zhàn)性的難題，是混合網格生成的核心和難點，已經成為網格生成領域的重要研究方向[8]。

目前，邊界層棱柱網格主要沿著自動生成的方向發(fā)展，有不同的算法實現思路，例如各項異性四面體網格聚合[9-10]、法向推進[11-16]、Level Set界面追蹤[17-19]以及求解PDE方程[20-21]，但不同算法思路的核心均是通過確定前沿推進單元的法矢和步長，構造棱柱網格單元。自動化棱柱網格生成方法具有較高的生成效率，人工參與少，但針對復雜外形，很難確定普適的控制參數，最終生成的棱柱網格難以避免較差質量單元的存在，且難以對局部網格單元質量進行調整和優(yōu)化，影響整體CFD模擬的精度和效率[22]。

此外，還有研究者通過構造的思路人工交互生成棱柱網格，例如本文作者等發(fā)展的基于背景結構網格約束框架[23]、基于網格變形插值[24]的交互式棱柱網格生成技術，能夠對局部棱柱網格質量進行有效調整，從而生成高質量的邊界層棱柱網格。尤其是基于變形插值的棱柱網格生成方法，可以直接嵌入目前的主流網格生成軟件中，對航空航天飛行器流動模擬具有潛在的工程應用價值。

但在進一步的算例測試中，基于變形思路的棱柱網格生成會在部分情況下產生翹曲現象，影響網格的質量，嚴重情況下會產生負體積單元，導致棱柱網格生成失敗。因此，本文通過特定算例對棱柱網格生成中翹曲現象的形成機制進行研究，找出翹曲現象產生的原因，并探索相應的消除算法，改善交互式棱柱網格生成技術的魯棒性。

1 交互式棱柱網格生成方法

交互式棱柱網格生成借鑒了網格變形的思路，將棱柱網格生成看成是物面網格連續(xù)變形的過程。圖1展示了交互式棱柱網格的生成思路，首先生成物面網格邊界棱線的空間推進面網格，然后利用三維空間插值方法將邊界棱線網格點上的推進位移插值到內部網格點上，得到棱柱網格每一層網格點的坐標，最后組裝獲得整個棱柱網格。交互式棱柱網格生成中采用徑向基函數實現邊界棱線網格點位移到內部點位移的插值，基函數使用網格質量控制效果較好的Wendland’s C2型緊支函數，關于交互式棱柱網格生成的詳細算法過程可以參考文獻[24]。

圖1 交互式棱柱網格生成方法

2 翹曲現象

針對物面形狀比較復雜的情況，例如曲率或外法矢變化劇烈的區(qū)域，基于徑向基函數插值的交互式棱柱網格生成方法可能在局部位置形成翹曲現象，降低棱柱網格的質量，嚴重時可產生負體積網格單元，導致棱柱網格生成失敗。本節(jié)通過典型航空外形的棱柱網格生成算例對這種局部區(qū)域產生的翹曲現象進行介紹和描述。

2.1 測試外形

選擇德國宇航中心(DLR)設計的雙引擎現代寬體運輸機標模DLR-F6作為測試算例。DLR-F6模型被選為兩屆AIAA阻力預測會議DPWII、DPWIII的標準計算外形，用于評估世界范圍內CFD軟件的發(fā)展水平，具有十分廣泛的影響力。圖2給出了DLR-F6模型翼身組合體構型下的外形，采用下單翼的設計，機翼從機身的中下部穿出，形成一些復雜的交叉位置，給棱柱網格的生成帶來一定挑戰(zhàn)[25]。網格生成測試中采用DLR-F6風洞縮比模型的尺寸，翼展長度為1 173 mm，機身長度為1 192 mm，詳細幾何信息可見文獻[26-27]。

圖2 DLR-F6翼身組合體構型[25]

2.2 棱柱網格生成

棱柱網格生成測試中發(fā)現翹曲現象主要發(fā)生在DLR-F6翼身組合體模型的機翼區(qū)域，機身位置翹曲不明顯。為了更好地觀察和分析翹曲現象，減少多余因素的干擾，本節(jié)采用DLR-F6模型的機翼開展棱柱網格生成的測試。

圖3給出了DLR-F6模型機翼的表面網格，整個機翼表面依據前后緣區(qū)域分割成17個網格面，網格面由三角形和四邊形單元混合而成，后緣為了保證推進后的棱柱網格質量，采用全四邊形的單元。整個表面網格包含12 826個網格點，374個三角形單元，12 565個四邊形單元。

圖3 DLR-F6模型機翼表面網格

棱柱網格生成測試通過中國空氣動力研究與發(fā)展中心自主開發(fā)的棱柱網格生成程序PG2(Prismatic Grid Generation)開展。圖4給出了生成后的DLR-F6機翼棱柱網格，其中第1層網格高度Δh設定為0.001 mm，網格高度法向增長比例為1.2，沿法向分布有39個網格點?？梢钥闯?，PG2能夠成功生成機翼的棱柱網格，且外網格面比較規(guī)則，但在靠近前緣和后緣的位置產生了明顯的條狀翹曲，如圖4(a)所示。選取翹曲明顯的兩個剖面z=217 mm、z=245 mm，可以發(fā)現，翹曲現象在前后緣附近的上下表面均有發(fā)生，后緣附近的翹曲更加明顯，存在顯著的凸起，而在遠離前后緣的機翼中部位置，沒有發(fā)生翹曲問題，如圖4(b)和圖4(c)所示。

圖4 DLR-F6機翼棱柱網格生成

圖4的棱柱網格生成結果表明，在曲率或法矢變化明顯的區(qū)域(機翼前后緣)容易發(fā)生翹曲現象，而在曲率或推進法矢變化緩慢的區(qū)域(遠離前后緣的機翼中部)發(fā)生翹曲現象的可能性較小。

3 翹曲形成機制

基于徑向基函數插值的交互式棱柱網格生成方法是使用物面邊界網格點推進位移插值得到內部網格點的位移，然后組裝得出整體的棱柱網格。當物面曲率或推進法矢急劇變化時，推進位移在某個坐標方向的分量可能會隨著法矢的變化而劇烈改變，并通過徑向基函數插值傳遞給內部網格點。因此猜測棱柱網格生成中的翹曲現象是由法矢變化導致的邊界點推進位移分布不光滑引起的。

為了驗證這一猜測，本節(jié)以平板外形棱柱網格生成為例，通過添加位移擾動的方式模擬曲率或法矢改變引起的位移變化，進一步分析棱柱網格翹曲現象的形成機制。

3.1 分析模型

圖5給出了棱柱網格翹曲現象形成機制的分析模型，采用單位長度的方形平板。該平板模型左下角位于坐標原點，有4條邊界棱線，每條棱線均勻分布101個網格點，平面網格為了便于顯示，在平板中心位置進行了粗化處理，整個平面網格由4 932個 三角形單元組成。每條邊界棱線沿著平板法向(z軸)推進了20層，第1層高度為0.001，整個推進高度H為0.1。模擬曲率或法矢變化效應的位移擾動添加在y=1.0的邊界棱線上。

圖5 平板棱柱網格生成測試算例

3.2 位移擾動信號

推進位移沿著各個方向的分量隨著物面曲率或推進法矢發(fā)生變化，推進方向變化越急，推進位移分量改變越快，推進方向變化越慢，推進位移分量改變越慢。因此，不同曲率或法矢改變引起的推進位移分量變化采用不同波長的正弦擾動信號來進行模擬，長波信號用來模擬曲率變化緩慢的情況，短波信號用于模擬曲率變化劇烈的情況。

圖6給出了添加的位移擾動信號，這里使用了3種擾動信號，長波信號Signal 1、短波信號Signal 2和混合信號Signal 3。長波信號和短波信號的擾動幅值均為0.01，波長分別為0.5和0.04， 混合信號為長波信號與短波信號的線性疊加。

圖6 位移擾動信號

位移擾動信號添加到y=1.0位置邊界棱線的空間推進面網格的最外層網格點上，空間推進面網格的內部網格點的位移擾動根據推進方向上的網格點分布線性插值得到。

3.3 擾動測試

圖7給出了長波擾動信號下的平板棱柱網格生成?？梢钥闯?，長波擾動信號添加后，y=1.0位置邊界棱線推進位移緩慢地起伏變化，但沿x方向的整體分布依然保持著光滑的特性，如圖7(a)所示；擾動位移添加后生成的棱柱網格在靠近y=1.0位置存在微小的起伏變化，但擾動效應迅速衰減，并沒有在整體棱柱網格中誘導出明顯的翹曲，如圖7(b)所示。

圖7 長波擾動信號下的平板棱柱網格生成

圖8給出了短波擾動信號下的平板棱柱網格生成?？梢钥闯?，短波擾動信號添加后，y=1.0位置邊界棱線推進位移快速劇烈地變化，沿x方向的整體分布不再保持光滑，呈現鋸齒形的波動，如圖8(a)所示；擾動位移添加后生成的棱柱網格在靠近y=1.0位置并沒有明顯的變化，但擾動效應卻在整體棱柱網格中被放大，誘導出顯著的翹曲現象，靠近x=0的區(qū)域產生了明顯的凸起，在靠近x=1.0區(qū)域產生了明顯的凹陷，如圖8(b)所示。

圖8 短波擾動信號下的平板棱柱網格生成

圖9給出了混合擾動信號下的平板棱柱網格生成?？梢钥闯?，混合擾動信號和短波擾動信號的情況相似，y=1.0位置邊界棱線推進位移出現快速劇烈的變化，并伴隨緩慢的起伏，位移分布光滑性不再存在，呈現鋸齒形的波動，如圖9(a)所示；擾動位移添加后生成的棱柱網格也出現了相似的翹曲現象，靠近x=0的區(qū)域凸起，靠近x=1.0區(qū)域凹陷，如圖9(b)所示。

圖9 混合擾動信號下的平板棱柱網格生成

3.4 翹曲產生機制分析

從3.3節(jié)的平板棱柱網格生成測試中可以得出，邊界棱線推進位移分布不光滑將誘導棱柱網格產生明顯的翹曲現象。但平板與DLR-F6機翼棱柱網格生成中的翹曲現象有所差異，平板棱柱網格的翹曲不僅有凸起的鼓包，還存在凹陷，而DLR-F6機翼棱柱網格的翹曲體現為鼓包。

進一步測試發(fā)現，棱柱網格翹曲的表現形式受邊界擾動信號的影響。圖10給出了一組不同的短波擾動信號，其中Case I采用3.2節(jié)中的短波擾動信號Signal 2，Case II在Case I信號的基礎上通過減小波長使得整個邊界棱線上增加半個周期的擾動，Case III在Case II信號基礎上相位提前180°。

圖10 不同的短波擾動信號

圖11給出了不同短波擾動信號下的平板棱柱網格生成?？梢钥闯?，在Case I擾動信號作用下，平板棱柱網格翹曲表現為凸起和凹陷；在Case II擾動信號作用下，平板棱柱網格翹曲表現為凸起；在Case III擾動信號作用下，平板棱柱網格翹曲表現為凹陷。圖11表明棱柱網格的翹曲存在多種表現形式，受邊界棱線推進位移分布的影響。DLR-F6機翼棱柱網格的翹曲與Case II短波擾動信號作用下的平板棱柱網格的翹曲相似，是多種翹曲表現形式中的一種。

圖11 不同短波擾動信號下的平板棱柱網格生成

上述結果間接驗證了本節(jié)開始的猜測，從而闡釋清楚交互式棱柱網格生成中翹曲現象的形成機制：曲率或推進法矢的劇烈變化引起推進位移的快速改變，造成邊界棱線推進位移分布的不光滑，誘導棱柱網格生成出現翹曲現象。

4 翹曲消除算法

圖12給出了不同擾動信號下平板棱柱網格y=0.5位置的剖面形狀，圖中擾動信號仍然使用3.2節(jié)中圖6所示的3種信號?？梢钥闯?，疊加了Signal 1和Signal 2的Signal 3誘導產生的棱柱網格翹曲與Signal 2的結果十分接近，說明棱柱網格翹曲主要由短波信號分量引起，長波信號分量的貢獻十分小。

圖12 平板棱柱網格剖面形狀(y=0.5)

不同波長信號的貢獻差異性為解決棱柱網格生成中的翹曲問題提供了一種可行的思路：通過對推進位移分布中的短波分量進行過濾，使推進位移分布光滑，從而消除網格翹曲。

拉普拉斯光順方法對不同波長信號的耗散特性存在明顯的差異，波長越短，信號衰減越快，波長越長，信號衰減越慢[28]。因此，本文采用基于拉普拉斯的思路發(fā)展面向邊界棱線空間推進位移分布的光順算法。

圖13給出了Laplace光順算法用于平面曲線的示意圖，圖中Pi為平面曲線上的第i個點，li-1, i表示點Pi-1與點Pi之間的距離，li, i+1表示點Pi與點Pi+1之間的距離。

圖13 Laplace光順算法示意圖

Laplace算法的迭代過程可以表示為

(1)

式中：n表示當前迭代步；n+1表示下一時刻迭代步；λ為松弛因子，取值范圍為0～1；L(·)為Laplace算子，可以表示為

(2)

式(2)中的Laplace算子假定曲線上點均勻分布的情況，對于非均勻分布的點，可以表示為

(3)

本文中待光順的Pi分別為(s, dx)、(s, dy)、(s, dz)，其中s為邊界棱線網格點相對起始點的曲線距離，通過分布線段線性疊加得到，dx、dy、dz分別為邊界棱線推進位移在3個坐標方向的分量。

5 結果與討論

圖14給出了邊界棱線推進位移拉普拉斯光順后生成的DLR-F6機翼棱柱網格，每個網格塊邊界棱線的光順次數統一設為20，松弛因子λ選擇0.5。可以看出，前后緣附近區(qū)域的翹曲現象被消除，棱柱網格外推進面保持光滑規(guī)則的形狀，如圖14(a)所示；仍然選取z=217，245 mm兩個剖面，并和光順前的剖面形狀比較，可以發(fā)現，機翼前后緣附近位置棱柱網格的鼓包現象被抑制，光順后的棱柱剖面形狀更加光滑，且光順過程對未發(fā)生鼓包現象的中部區(qū)域的作用很小，如圖14(b)和圖14(c)所示。圖14的結果表明發(fā)展的基于拉普拉斯的光順算法能夠消除棱柱網格生成中的翹曲現象。

圖14 DLR-F6機翼棱柱網格生成(光順后)

圖15給出了不同光順次數下DLR-F6機翼棱柱網格z=217 mm位置的剖面形狀。可以看出，隨著光順次數的增加，光順后的外輪廓線更加趨向直線。對于后緣位置，曲率變化平緩，光順次數越多，外輪廓線愈光滑，但對于前緣位置，曲率變化明顯，光順次數太大，會引起部分前緣區(qū)域的棱柱網格被壓縮，導致網格正交性降低。針對本文DLR-F6機翼算例，前緣位置的光順次數在5附近，后緣位置的光順次數大于20，能夠獲得較好的棱柱網格生成。圖15結果表明，針對不同區(qū)域，棱柱網格生成中的光順次數具有不同的最優(yōu)值。

圖15 不同光順次數下DLR-F6機翼棱柱網格剖面形狀(z=217 mm)

圖16給出了DLR-F6翼身組合體構型的棱柱網格生成，該構型表面網格有18 371個網格點，16 077個四邊形單元，642個三角形單元。棱柱網格生成的參數，如第1層網格高度、網格增長速率、推進層數，均與前文DLR-F6機翼的相同，光順次數仍然統一選擇為20?？梢钥闯?，整個棱柱網格的外推進面保持規(guī)則的形狀，棱柱網格單元的連續(xù)性和正交性很好，整體網格具有較高的質量。圖16的算例表明，完善后的棱柱網格生成方法具有在工程復雜外形中廣泛應用的潛力。

圖16 DLR-F6翼身組合體構型棱柱網格生成

6 結 論

針對交互式棱柱網格生成中出現的翹曲現象，利用平板外形，通過擾動信號模擬邊界棱線推進位移的變化，弄清了翹曲現象的形成機制?；诶绽构忭樂椒òl(fā)展了消除翹曲現象的光順技術，并通過DLR-F6外形棱柱網格生成進行了測試?；跍y試和分析結果，可以得到以下結論：

1) 交互式棱柱網格生成中的翹曲現象是由曲率或推進法矢劇烈變化導致推進位移沿邊界棱線的分布不光滑引起的。

2) 基于拉普拉斯的邊界棱線推進位移分布光順方法能夠消除棱柱網格生成中的翹曲，但不同區(qū)域的光順次數具有各自的最優(yōu)值。

3) 添加光順方法后的棱柱網格生成算法具有在工程復雜外形網格生成中應用的潛力。