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

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

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

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

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

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

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

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

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

1 交互式棱柱網(wǎng)格生成方法

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

圖1 交互式棱柱網(wǎng)格生成方法

2 翹曲現(xiàn)象

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

2.1 測試外形

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

圖2 DLR-F6翼身組合體構(gòu)型[25]

2.2 棱柱網(wǎng)格生成

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

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

圖3 DLR-F6模型機翼表面網(wǎng)格

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

圖4 DLR-F6機翼棱柱網(wǎng)格生成

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

3 翹曲形成機制

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

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

3.1 分析模型

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

圖5 平板棱柱網(wǎng)格生成測試算例

3.2 位移擾動信號

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

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

圖6 位移擾動信號

位移擾動信號添加到y(tǒng)=1.0位置邊界棱線的空間推進面網(wǎng)格的最外層網(wǎng)格點上，空間推進面網(wǎng)格的內(nèi)部網(wǎng)格點的位移擾動根據(jù)推進方向上的網(wǎng)格點分布線性插值得到。

3.3 擾動測試

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

圖7 長波擾動信號下的平板棱柱網(wǎng)格生成

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

圖8 短波擾動信號下的平板棱柱網(wǎng)格生成

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

圖9 混合擾動信號下的平板棱柱網(wǎng)格生成

3.4 翹曲產(chǎn)生機制分析

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

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

圖10 不同的短波擾動信號

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

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

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

4 翹曲消除算法

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

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

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

拉普拉斯光順方法對不同波長信號的耗散特性存在明顯的差異，波長越短，信號衰減越快，波長越長，信號衰減越慢[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為邊界棱線網(wǎng)格點相對起始點的曲線距離，通過分布線段線性疊加得到，dx、dy、dz分別為邊界棱線推進位移在3個坐標方向的分量。

5 結(jié)果與討論

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

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

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

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

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

圖16 DLR-F6翼身組合體構(gòu)型棱柱網(wǎng)格生成

6 結(jié) 論

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

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

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

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