陳小毛, 宣傳偉

(1.南京航空航天大學(xué) 無人機(jī)研究院,江蘇 南京 210016; 2.上海宇航系統(tǒng)工程研究所,上海 201100)

由于精度高、邏輯關(guān)系簡單、流場計算效率高等特點,結(jié)構(gòu)網(wǎng)格技術(shù)是計算流體動力學(xué)(computational fluid dynamics,CFD)的首選技術(shù)。但是,復(fù)雜外形流場結(jié)構(gòu)網(wǎng)格的生成是一項費(fèi)時并需要經(jīng)驗積累的工作。重疊網(wǎng)格技術(shù)[1]將計算域劃分為若干個簡單的子域,針對各個子域單獨劃分網(wǎng)格,從而大大降低了結(jié)構(gòu)網(wǎng)格的生成難度,同時保證了網(wǎng)格生成質(zhì)量,也提高了網(wǎng)格生成效率。因此,結(jié)構(gòu)重疊網(wǎng)格技術(shù)得以廣泛應(yīng)用于復(fù)雜外形流場[2]及多體相對運(yùn)動[3]數(shù)值模擬中。結(jié)構(gòu)重疊網(wǎng)格技術(shù)的核心是“挖洞”,即確定各子域網(wǎng)格之間的插值關(guān)系并剔除不需要計算的網(wǎng)格。洞面優(yōu)化[4]則是針對“挖洞”的一種優(yōu)化技術(shù),其目的是改善重疊區(qū)域,使得網(wǎng)格重疊區(qū)域最小、插值邊界遠(yuǎn)離物面,提高流場求解精度。

重疊網(wǎng)格技術(shù)是求解復(fù)雜外形流場及多體相對運(yùn)動問題的有效途徑,而準(zhǔn)確高效的洞面優(yōu)化方法是其關(guān)鍵技術(shù)之一。常用的洞面優(yōu)化方法包括陣面推進(jìn)技術(shù)[5]、隱式挖洞技術(shù)[6]以及割補(bǔ)法[7]。陣面推進(jìn)技術(shù)是利用2個陣面反復(fù)迭代推進(jìn)挖洞面來進(jìn)行重疊網(wǎng)格挖洞的過程。陣面的推進(jìn)距離依賴于陣面和相鄰網(wǎng)格點之間的相對距離,且一定程度上需要人工操作,自動化程度稍低。隱式挖洞技術(shù)是一種基于網(wǎng)格單元品質(zhì)的插值單元選擇方法[8]。該方法沒有單獨的初始挖洞過程,僅僅需要尋找最佳品質(zhì)的貢獻(xiàn)單元即可[9],因此網(wǎng)格裝配效率和自動化程度較高。但正是由于挖洞過程的缺失,插值單元可能出現(xiàn)在物體內(nèi)部,從而影響計算結(jié)果[10]。作為陣面推進(jìn)技術(shù)的改進(jìn)方法,割補(bǔ)法以洞邊界點為起點,沿著網(wǎng)格線離散地做推進(jìn)運(yùn)動,每段網(wǎng)格線的長度即為該次推進(jìn)的距離。與陣面推進(jìn)方法相比,割補(bǔ)法的迭代推進(jìn)次數(shù)少,自動化程度高,因此得到了較好的發(fā)展。但是,割補(bǔ)法仍然存在如下2個問題:一是當(dāng)壁面距離很近時,割補(bǔ)法的填補(bǔ)過程容易因少數(shù)洞邊界點找點過程的失敗而導(dǎo)致網(wǎng)格進(jìn)入壁面內(nèi);二是割補(bǔ)完后洞邊界的位置不可預(yù)測,洞邊界參差不齊。針對第1類問題,文獻(xiàn)[11]提出了Volume Method和Vertex Method 2種方法。其實質(zhì)是去掉填補(bǔ)過程,將割補(bǔ)法的2個過程變?yōu)?個過程。但由于網(wǎng)格體積的非單調(diào)性,Volume Method可能會導(dǎo)致單元屬性轉(zhuǎn)換混亂,出現(xiàn)網(wǎng)格孤島。文獻(xiàn)[12]通過引入2類洞內(nèi)點概念,使得位于另一物體內(nèi)部的網(wǎng)格點的屬性在填補(bǔ)過程中不被修改,同樣解決了狹縫重疊中洞邊界易挖進(jìn)物面的問題。文獻(xiàn)[13-14]提出將割補(bǔ)法與物面距準(zhǔn)則[15]相結(jié)合的方法解決了第2類問題,該方法在網(wǎng)格割補(bǔ)之前增加了一個以物面距離為準(zhǔn)則的網(wǎng)格分類過程。

針對割補(bǔ)過程中產(chǎn)生的2類問題,本文在網(wǎng)格切割過程中將物面距離作為判別標(biāo)準(zhǔn)。與文獻(xiàn)[14]不同,本文方法僅在切割過程中搜索洞邊界點的貢獻(xiàn)單元,避免了對全局所有網(wǎng)格點進(jìn)行貢獻(xiàn)單元的搜索,提高了網(wǎng)格裝配效率。當(dāng)洞邊界點的物面距小于或等于其貢獻(xiàn)單元的物面距時停止切割,從而使得洞邊界處于壁面中間位置;在網(wǎng)格填補(bǔ)過程中,以切割結(jié)束時的洞邊界作為初始陣面進(jìn)行網(wǎng)格回填,直至產(chǎn)生兩層插值邊界。本文首先對傳統(tǒng)割補(bǔ)法進(jìn)行簡單介紹,然后詳細(xì)闡述改進(jìn)后的割補(bǔ)法過程,最后通過2個典型流動算例驗證本文所提方法的可靠性。

1 割補(bǔ)法

挖洞和貢獻(xiàn)單元搜索是結(jié)構(gòu)重疊網(wǎng)格技術(shù)的2個關(guān)鍵過程,而其中的挖洞過程本質(zhì)上是一個將網(wǎng)格節(jié)點分類的過程,即直接參與流場計算的點為正常點,不參與計算的點為洞內(nèi)點,參與插值計算的點為洞邊界點。挖洞又可以細(xì)分為初始挖洞和洞邊界優(yōu)化2個步驟。初始挖洞旨在標(biāo)記出所有落入物體內(nèi)部的網(wǎng)格點,從而不參與流場計算。割補(bǔ)法在初始挖洞的基礎(chǔ)上,經(jīng)過“切割”和“填補(bǔ)”2個階段，使得各子網(wǎng)格間的洞邊界最優(yōu)化。貢獻(xiàn)單元的搜索也在該過程中完成。

1.1 切割過程

經(jīng)過初始挖洞之后的網(wǎng)格重疊區(qū)域通常十分巨大,洞邊界可能位于流場梯度較大的邊界層內(nèi),流場計算精度低。在該階段,所有洞邊界點依次尋找貢獻(xiàn)單元。若貢獻(xiàn)單元的所有節(jié)點都為正常點,則該洞邊界點轉(zhuǎn)變?yōu)槎磧?nèi)點,與其相鄰的正常點轉(zhuǎn)變?yōu)槎催吔琰c,由此將洞邊界推離物面;若貢獻(xiàn)單元的節(jié)點不全為正常點,則停止該過程。

1.2 填補(bǔ)過程

經(jīng)切割之后,沒有網(wǎng)格單元重疊,而填補(bǔ)階段本質(zhì)上是切割階段的逆過程,其目的是將洞邊界逐步收縮從而實現(xiàn)網(wǎng)格單元的重疊。在該階段,所有洞邊界點依次尋找貢獻(xiàn)單元。若貢獻(xiàn)單元的節(jié)點不全為正常點,則該洞邊界點轉(zhuǎn)變?yōu)檎｜c,與其相鄰的洞內(nèi)點轉(zhuǎn)變?yōu)槎催吔琰c,由此收縮洞邊界,直至再無洞內(nèi)點轉(zhuǎn)變?yōu)槎催吔琰c。

割補(bǔ)法過程示意圖如圖1所示,圖1中:N表示正常網(wǎng)格點(normal point),參與流場計算;H表示挖去的洞內(nèi)點(hole point)不參與計算;F表示洞邊界點(fringe point),用于插值計算。

圖1 割補(bǔ)法過程示意圖

2 割補(bǔ)法的改進(jìn)

在網(wǎng)格切割過程中,洞邊界的推進(jìn)是按照網(wǎng)格線來進(jìn)行的,因此最終洞邊界的位置嚴(yán)重依賴于網(wǎng)格單元尺度在各方向上的分布特征。對于同一個問題,經(jīng)過割補(bǔ)后,不同的網(wǎng)格可能會產(chǎn)生不同的洞邊界;而且當(dāng)網(wǎng)格匹配性較差時,洞邊界甚至可能深入物面邊界層,從而降低插值精度。另外,洞邊界在空間上的分布也會參差不齊[14],如圖2所示,機(jī)翼和油箱分離模型經(jīng)過割補(bǔ)后,洞邊界靠近物面而且參差不齊。

圖2 割補(bǔ)法洞邊界參差不齊現(xiàn)象

在網(wǎng)格填補(bǔ)過程中,當(dāng)物面距離很近時,填補(bǔ)過程容易由于少數(shù)洞邊界點的尋點過程失敗而導(dǎo)致網(wǎng)格進(jìn)入物面內(nèi)部,如圖3所示[12]。

圖3 割補(bǔ)法的縫隙問題

針對上述2種問題,本文對傳統(tǒng)割補(bǔ)法進(jìn)行改進(jìn),具體步驟如下:

(1) 物面距計算。將所有網(wǎng)格點標(biāo)記為正常點,并計算所有網(wǎng)格點的物面距(若該點屬于背景網(wǎng)格,則將其物面距離設(shè)置為一個大值)。

(2) 初始挖洞。使用洞映射法[16]對所有網(wǎng)格進(jìn)行挖洞(洞邊界位于物面處),將位于物面內(nèi)部的網(wǎng)格點標(biāo)記為洞內(nèi)點,與洞內(nèi)點相鄰的正常點標(biāo)記為洞邊界點。

(3) 切割。針對所有洞邊界點,進(jìn)行ADT[17]尋點操作,搜索并保存其貢獻(xiàn)單元信息,通過三線性插值法得到貢獻(xiàn)單元的物面距。比較該洞邊界點與其貢獻(xiàn)單元的物面距大小。若前者大于后者,則將該洞邊界點標(biāo)記為洞內(nèi)點,而與其相鄰的正常點標(biāo)記為洞邊界點。

(4) 循環(huán)步驟(3),直至所有洞邊界點的物面距小于或等于其貢獻(xiàn)單元的物面距。

(5) 填補(bǔ)。將與洞邊界點相鄰的2個洞內(nèi)點標(biāo)記為插值點,然后將該洞邊界點標(biāo)記為正常點。

本文方法的過程示意圖如圖4所示,該套網(wǎng)格系統(tǒng)由4塊子網(wǎng)格和1塊背景網(wǎng)格組成。每塊子網(wǎng)格包含3 960個網(wǎng)格點,背景網(wǎng)格包含42 411個網(wǎng)格點。以物面邊界為初始挖洞邊界的結(jié)果如圖4a所示,從圖4a可以看出,位于物面內(nèi)部的網(wǎng)格均已被挖去,將已被挖去的網(wǎng)格點標(biāo)記為洞內(nèi)點,即上述步驟(2)。將與洞內(nèi)點相鄰的正常點標(biāo)記為洞邊界點對應(yīng)步驟(3)、步驟(4),使用ADT方法循環(huán)搜索洞邊界點的貢獻(xiàn)單元,比較洞邊界點與其貢獻(xiàn)單元的物面距大小,更新洞邊界,直至所有洞邊界點的物面距小于或等于其貢獻(xiàn)單元的物面距,結(jié)果如圖4b所示。此時,洞邊界位于物面中間位置。為保證二階計算精度,對應(yīng)上述步驟(5),將與洞邊界點相鄰的2個洞內(nèi)點標(biāo)記為插值點,如圖4c所示。

圖4 改進(jìn)型割補(bǔ)法過程

針對上述網(wǎng)格系統(tǒng),本文方法與文獻(xiàn)[14]方法的效率對比情況見表1所列。文獻(xiàn)[14]方法需要搜索所有網(wǎng)格點的貢獻(xiàn)單元,因此需要進(jìn)行尋點操作的網(wǎng)格點數(shù)為所有子網(wǎng)格點數(shù)目之和,即58 251個。而本文方法的尋點操作是在網(wǎng)格切割的同時進(jìn)行的,且只需對洞邊界點進(jìn)行尋點操作,累計點數(shù)為6 582個,明顯少于文獻(xiàn)[14]方法,耗時不足其1/2,具有更高的效率(本文所有計算均基于Intel Xeon CPU E5-2650平臺,主頻為2.00 GHz)。

表1 不同方法效率對比

如前所述,針對三段翼網(wǎng)格系統(tǒng),當(dāng)縫翼與主翼網(wǎng)格很近時,在填補(bǔ)過程中傳統(tǒng)割補(bǔ)法往往會因少數(shù)洞邊界點的尋點過程失敗,而導(dǎo)致洞邊界進(jìn)入物面內(nèi)部(見圖3)。使用本文方法對其進(jìn)行重疊網(wǎng)格裝配的效果如圖5所示。由于切割過程中保存了所有洞邊界點的貢獻(xiàn)單元信息,因此填補(bǔ)過程可不用搜索貢獻(xiàn)單元,從而避免了壁面距離很近時因少數(shù)點的貢獻(xiàn)單元搜索失敗而導(dǎo)致整個網(wǎng)格重疊過程失敗的可能性。從圖5可以看出,洞邊界位于縫翼和主翼面中間,而且主翼面網(wǎng)格也沒有進(jìn)入縫翼內(nèi)部。

圖5 縫翼距離主翼很近時重疊網(wǎng)格效果

S300C型旋翼結(jié)合Robin機(jī)身重疊網(wǎng)格系統(tǒng)如圖6所示。

圖6 旋翼機(jī)身重疊網(wǎng)格系統(tǒng)

該套網(wǎng)格系統(tǒng)由如下5個部分組成:背景網(wǎng)格、Robin機(jī)身網(wǎng)格、3個S300C旋翼網(wǎng)格,共3 588 672個單元。使用本文方法對其進(jìn)行網(wǎng)格裝配后的效果如圖7所示,不同子網(wǎng)格間插值邊界分布正確、合理。本文方法與文獻(xiàn)[14]方法效率對比情況見表2所列。從表2可以看出,本文方法具有更高的效率。

圖7 挖洞效果

表2 旋翼機(jī)身網(wǎng)格不同方法效率對比

3 數(shù)值模擬

3.1 計算方法

本文采用有限體積法對Navier-Stokes控制方程進(jìn)行離散求解,空間離散采用Roe-FDS(flux difference splitting)格式,具備二階空間離散精度;湍流模型采用SST模型;時間離散采用LUSGS(Lower-Upper Symmetric Gauss-Seidel)隱式計算方法。

3.2 30P30N三段式機(jī)翼

對經(jīng)典的二維三段式機(jī)翼30P30N[18]進(jìn)行數(shù)值模擬。來流馬赫數(shù)Ma=0.2,雷諾數(shù)Re=8.7×106,攻角α=8°。流場網(wǎng)格由以下4個部分組成:縫翼網(wǎng)格為300×101、主翼面網(wǎng)格為682×151、襟翼網(wǎng)格為300×141、背景網(wǎng)格為609×600。經(jīng)洞面優(yōu)化后重疊效果如圖8所示,重疊網(wǎng)格的洞邊界清晰整潔,位于不同翼面的中間位置。翼面壓力云圖和等馬赫線如圖9、圖10所示。從圖9、圖10可以看出流場分布合理,重疊邊界處等馬赫線和壓力線過渡光滑。翼面壓力系數(shù)的數(shù)值與試驗結(jié)果對比如圖11所示。圖11中：Cp為壓力系數(shù);x/c為x坐標(biāo)與弦長c的比值。計算結(jié)果與試驗值[19]吻合良好,表明本文洞面優(yōu)化方法準(zhǔn)確可靠。

圖8 30P30N機(jī)翼重疊網(wǎng)格效果圖

圖9 30P30翼面壓力云圖

圖10 30P30翼面等馬赫線

圖11 30P30N翼面壓力系數(shù)計算值與試驗值對比

3.3 Titan Ⅳ運(yùn)載火箭

為進(jìn)一步驗證本文方法在三維復(fù)雜外形流場中的適用性,對美國大力神Titan Ⅳ運(yùn)載火箭[20]超音速外流場進(jìn)行了數(shù)值仿真。TitanⅣ運(yùn)載火箭采用并聯(lián)式布局,芯級和助推級距離較近,是典型的狹縫重疊問題。Titan Ⅳ火箭重疊網(wǎng)格系統(tǒng)如圖12所示,由以下3個部分組成:芯級網(wǎng)格、2個助推器網(wǎng)格,共 4 905 743個六面體單元。采用本文方法進(jìn)行挖洞后的效果如圖13所示。從圖13可以看出,洞邊界整齊合理,位于物面中間,且不存在狹縫問題。

圖12 Titan Ⅳ重疊網(wǎng)格系統(tǒng)

圖13 Titan Ⅳ挖洞效果

本算例計算條件為:馬赫數(shù)Ma=1.6, 攻角α=0°,雷諾數(shù)Re=115×107。對稱面流場馬赫數(shù)分布如圖14所示。芯級中心線壓力數(shù)值結(jié)果與試驗值的比較結(jié)果如圖15所示。

圖14 Titan Ⅳ對稱面馬赫數(shù)云圖

圖15 芯級中心線壓力分布

圖15中:p為芯級中心線壓強(qiáng);pref為來流遠(yuǎn)場壓強(qiáng);x為沿芯級中心線位置;r為芯級半徑。從圖14、圖15可以看出,網(wǎng)格重疊區(qū)域等壓線分布合理、光滑, 計算結(jié)果與試驗值吻合得很好。由于芯級與助推級之間存在連接機(jī)構(gòu), 導(dǎo)致試驗的峰值點略高,且在連接機(jī)構(gòu)附近數(shù)值結(jié)果與實驗值之間的偏差稍大(圖15中,x/r值在25～35之間的區(qū)域)。本文所涉及的4種模型分別采用本文方法與文獻(xiàn)[14]方法進(jìn)行重疊網(wǎng)格裝配效率對比,結(jié)果如圖16所示。圖16中:橫坐標(biāo)表示4種模型編號;縱坐標(biāo)表示網(wǎng)格裝配所用時間。從圖16可以看出,本文方法具有更高的效率,且隨著網(wǎng)格數(shù)目的增多,本文方法具有更明顯的效率優(yōu)勢。

圖16 4種模型采用不同方法的效率對比

4 結(jié) 論

本文針對傳統(tǒng)割補(bǔ)法的狹縫問題和洞邊界位置不可預(yù)測性問題,發(fā)展了一種改進(jìn)型洞邊界優(yōu)化方法。

(1) 通過將壁面距離作為網(wǎng)格切割過程的判別標(biāo)準(zhǔn),使得洞邊界位于物面中間,解決了傳統(tǒng)割補(bǔ)法洞邊界不可預(yù)測且參差不齊的問題。

(2) 相比文獻(xiàn)[14]方法,本文方法更簡潔,在網(wǎng)格切割過程中只需搜索當(dāng)前洞邊界點的貢獻(xiàn)單元,避免了搜索所有網(wǎng)格點的貢獻(xiàn)單元,具備更高的網(wǎng)格裝配效率。

(3) 在網(wǎng)格填補(bǔ)過程中,由于不需要再次搜索貢獻(xiàn)單元,因此解決了因少數(shù)洞邊界點找點過程的失敗而導(dǎo)致網(wǎng)格進(jìn)入物面內(nèi)的問題。與文獻(xiàn)[11]的Volume Method方法相比,由于本文方法未涉及到網(wǎng)格體積大小,因此不存在“孤點”問題,具有更好的魯棒性。

(4) 通過經(jīng)典算例驗證了本文方法的準(zhǔn)確性與可靠性。