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

中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

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

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

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

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

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

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

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

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

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

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

2 翹曲現(xiàn)象

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

2.1 測(cè)試外形

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

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

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

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

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

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

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

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

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

3 翹曲形成機(jī)制

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

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

3.1 分析模型

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

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

3.2 位移擾動(dòng)信號(hào)

推進(jìn)位移沿著各個(gè)方向的分量隨著物面曲率或推進(jìn)法矢發(fā)生變化，推進(jìn)方向變化越急，推進(jìn)位移分量改變?cè)娇?，推進(jìn)方向變化越慢，推進(jìn)位移分量改變?cè)铰?。因此，不同曲率或法矢改變引起的推進(jìn)位移分量變化采用不同波長(zhǎng)的正弦擾動(dòng)信號(hào)來(lái)進(jìn)行模擬，長(zhǎng)波信號(hào)用來(lái)模擬曲率變化緩慢的情況，短波信號(hào)用于模擬曲率變化劇烈的情況。

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

圖6 位移擾動(dòng)信號(hào)

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

3.3 擾動(dòng)測(cè)試

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

圖7 長(zhǎng)波擾動(dòng)信號(hào)下的平板棱柱網(wǎng)格生成

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

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

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

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

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

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

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

圖10 不同的短波擾動(dòng)信號(hào)

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

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

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

4 翹曲消除算法

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

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

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

拉普拉斯光順?lè)椒▽?duì)不同波長(zhǎng)信號(hào)的耗散特性存在明顯的差異，波長(zhǎng)越短，信號(hào)衰減越快，波長(zhǎng)越長(zhǎng)，信號(hào)衰減越慢[28]。因此，本文采用基于拉普拉斯的思路發(fā)展面向邊界棱線空間推進(jìn)位移分布的光順?biāo)惴ā?/p>

圖13給出了Laplace光順?biāo)惴ㄓ糜谄矫媲€的示意圖，圖中Pi為平面曲線上的第i個(gè)點(diǎn)，li-1, i表示點(diǎn)Pi-1與點(diǎn)Pi之間的距離，li, i+1表示點(diǎn)Pi與點(diǎn)Pi+1之間的距離。

圖13 Laplace光順?biāo)惴ㄊ疽鈭D

Laplace算法的迭代過(guò)程可以表示為

(1)

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

(2)

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

(3)

本文中待光順的Pi分別為(s, dx)、(s, dy)、(s, dz)，其中s為邊界棱線網(wǎng)格點(diǎn)相對(duì)起始點(diǎn)的曲線距離，通過(guò)分布線段線性疊加得到，dx、dy、dz分別為邊界棱線推進(jìn)位移在3個(gè)坐標(biāo)方向的分量。

5 結(jié)果與討論

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

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

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

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

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

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

6 結(jié) 論

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

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

2) 基于拉普拉斯的邊界棱線推進(jìn)位移分布光順?lè)椒軌蛳庵W(wǎng)格生成中的翹曲，但不同區(qū)域的光順次數(shù)具有各自的最優(yōu)值。

3) 添加光順?lè)椒ê蟮睦庵W(wǎng)格生成算法具有在工程復(fù)雜外形網(wǎng)格生成中應(yīng)用的潛力。