重疊網(wǎng)格中割補(bǔ)法的改進(jìn)及應(yīng)用

陳小毛, 宣傳偉

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

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

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

1 割補(bǔ)法

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

1.1 切割過程

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

1.2 填補(bǔ)過程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表1 不同方法效率對比

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

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

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

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

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

圖7 挖洞效果

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

3 數(shù)值模擬

3.1 計(jì)算方法

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

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個(gè)部分組成:縫翼網(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ù)值與試驗(yàn)結(jié)果對比如圖11所示。圖11中：Cp為壓力系數(shù);x/c為x坐標(biāo)與弦長c的比值。計(jì)算結(jié)果與試驗(yàn)值[19]吻合良好,表明本文洞面優(yōu)化方法準(zhǔn)確可靠。

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

圖9 30P30翼面壓力云圖

圖10 30P30翼面等馬赫線

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

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

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

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

圖13 Titan Ⅳ挖洞效果

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

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

圖15 芯級中心線壓力分布

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

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

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