康忠良 方媛媛

(1.中國建筑科學(xué)研究院建研科技股份有限公司，北京 100013；2.北京市燃?xì)饧瘓F(tuán)研究院，北京 100011)

混合網(wǎng)格重疊方法在多體相對(duì)運(yùn)動(dòng)中的應(yīng)用

康忠良 方媛媛

(1.中國建筑科學(xué)研究院建研科技股份有限公司，北京 100013；2.北京市燃?xì)饧瘓F(tuán)研究院，北京 100011)

本文發(fā)展了混合網(wǎng)格的動(dòng)態(tài)重疊方法?；诰€性重構(gòu)方法，給出了一種適用于混合網(wǎng)格間的插值策略。所建立的網(wǎng)格間邊界定義方法的性能優(yōu)秀可靠，網(wǎng)格裝配后的插值邊界光滑且網(wǎng)格匹配性較好。耦合求解非定常NS方程和剛體6DOF運(yùn)動(dòng)方程來模擬多體相對(duì)運(yùn)動(dòng)問題。外掛物分離問題的數(shù)值計(jì)算表明，所發(fā)展的動(dòng)網(wǎng)格方法對(duì)于處理多體相對(duì)運(yùn)動(dòng)模擬問題是快速準(zhǔn)確可靠的。

動(dòng)網(wǎng)格；重疊網(wǎng)格；混合網(wǎng)格；多體運(yùn)動(dòng)

1 引言

在計(jì)算流體力學(xué)(CFD)應(yīng)用中，采用混合網(wǎng)格方法不但對(duì)復(fù)雜構(gòu)型具有強(qiáng)大的幾何適應(yīng)能力，網(wǎng)格生成的人工工作量少，而且易于通過采用不同的單元類型及調(diào)整網(wǎng)格的疏密來適應(yīng)不同的流場特征，并容易生成整體網(wǎng)格、整體求解。但是，在實(shí)際應(yīng)用中，混合網(wǎng)格方法也遇到了一些困難，比如對(duì)于多體相對(duì)運(yùn)動(dòng)問題，就需要對(duì)網(wǎng)格進(jìn)行再生，處理非常復(fù)雜。

為了解決這一問題，Nakahashi[1]首先提出了非結(jié)構(gòu)網(wǎng)格重疊方法。該方法吸收了結(jié)構(gòu)網(wǎng)格重疊方法[2]的優(yōu)點(diǎn)，后來又被推廣到混合網(wǎng)格重疊及動(dòng)態(tài)重疊應(yīng)用中[3-5]，得到了較大的發(fā)展。混合網(wǎng)格動(dòng)態(tài)重疊方法的優(yōu)點(diǎn)就是在模擬多體間具有大幅相對(duì)運(yùn)動(dòng)問題時(shí)，在物體運(yùn)動(dòng)過程中不需要網(wǎng)格再生，算法實(shí)現(xiàn)簡單，并且可以獲得較高的動(dòng)網(wǎng)格處理效率。

盡管如此，混合網(wǎng)格重疊方法目前仍然面臨諸多難題，其中，網(wǎng)格間邊界的定義及其插值方法就是其中一個(gè)重要方面，已有研究均基于多層網(wǎng)格節(jié)點(diǎn)的流場變量進(jìn)行插值計(jì)算，該思路實(shí)現(xiàn)較為繁瑣，且容易在高速流動(dòng)強(qiáng)間斷附近產(chǎn)生較大的耗散。另一方面，在實(shí)際工程應(yīng)用中還應(yīng)該強(qiáng)調(diào)動(dòng)網(wǎng)格處理的自動(dòng)化程度、最大限度地減少人工干預(yù)，否則也不能體現(xiàn)混合網(wǎng)格重疊方法的優(yōu)勢。

本文在優(yōu)化網(wǎng)格間邊界定義的基礎(chǔ)上，給出一種適用于任意單元類型的插值策略，并給出一種完全自動(dòng)化的高效的網(wǎng)格動(dòng)態(tài)裝配策略。

2 數(shù)值方法

本文算法構(gòu)造基于任意多面體混合網(wǎng)格單元。算法不考慮具體的網(wǎng)格拓?fù)?，?duì)不同的單元類型統(tǒng)一處理?？臻g離散采用格心有限體積法，控制體取網(wǎng)格單元，控制面取單元表面，流場變量存儲(chǔ)在單元中心。

2.1 控制方程

舍去源項(xiàng)的三維非定?？蓧嚎sNS方程組在直角坐標(biāo)系下的守恒積分形式可表示為

(1)

式中：Ω為控制體；?Ω為控制面；Q為守恒變矢量；Fc為對(duì)流通量；Fv為粘性通量。各項(xiàng)的具體描述詳見文獻(xiàn)[6]。

2.2 數(shù)值離散方法

對(duì)流通量離散采用Roe[7]格式近似求解Riemann問題。假設(shè)控制面IJ的左右單元分別為I和J,則對(duì)流通量可表示為:

(2)

Barth[8]提出的線性重構(gòu)方法假定解在控制體內(nèi)呈線性分布，面左右兩側(cè)的值可表示為

(3)

式中：U為任意流場變量值；▽UI為單元I體心的梯度；Ψ為限制器函數(shù)；r為從體心到面心的矢量。▽UI采用最小二乘法[7]計(jì)算，Ψ采用收斂特性非常好的Venkatakrishnan[9]限制器。

另外，粘性通量采用中心格式離散，湍流模型采用Spalart-Allmaras[10]一方程模型，時(shí)間離散采用格式。非定常流動(dòng)控制方程采用雙時(shí)間步長法[11]求解，每一個(gè)物理時(shí)間步上的偽時(shí)間推進(jìn)采用LU-SGS隱式格式[6]。

3 混合網(wǎng)格動(dòng)態(tài)重疊方法

與現(xiàn)存的基于格點(diǎn)[3-5]的混合網(wǎng)格重疊方法不同，本文基于格心展開，通過網(wǎng)格間邊界定義后，初始網(wǎng)格單元被分為三類：活動(dòng)單元、插值單元和非活動(dòng)單元?；顒?dòng)單元在計(jì)算區(qū)域內(nèi)部參與流場計(jì)算，插值單元分布在網(wǎng)格邊界用于子網(wǎng)格間的信息交換，非活動(dòng)單元作為非計(jì)算單元被挖去。另外，本文定義包圍插值單元網(wǎng)格中心點(diǎn)的單元為其宿主單元。圖1給出了一個(gè)網(wǎng)格邊界定義示例。其中，子網(wǎng)格G2中單元B(1-2-3)是子網(wǎng)格G1中單元A(a-b-c)的宿主單元，這是因?yàn)閱卧狝的格心落在了單元B的內(nèi)部。

3.1 插值方法

在本文混合網(wǎng)格重疊系統(tǒng)的流場計(jì)算中，子網(wǎng)格間的信息交換需要插值單元通過其宿主單元來獲取其它子網(wǎng)格的相應(yīng)邊界信息，這一信息交換過程稱為插值。與傳統(tǒng)的完全依賴于格點(diǎn)流場變量進(jìn)行插值的方法不同，下面給出一種新型的插值策略。

該方法依賴于插值單元的格心流場變量值、單元梯度值和限制器函數(shù)值的計(jì)算來完成。比如圖1，插值單元A的狀態(tài)可以通過下式確定

▽UA=▽UB

ΨA=ΨB

(4)

通過上式計(jì)算，插值單元A可以充分利用宿主單元B的所有相鄰單元的信息。這是因?yàn)槭?4)與式(3)相類似，▽UB和ΨB都是通過單元B的所有相鄰單元構(gòu)造得來。顯然，這里要求宿主單元B在流場重構(gòu)時(shí)能夠構(gòu)造出足夠多的模板。

當(dāng)網(wǎng)格間邊界重疊區(qū)附近流動(dòng)梯度變化較小時(shí)，該插值方法與線性重構(gòu)一樣，可使空間離散精度達(dá)到二階。當(dāng)邊界附近流動(dòng)變化劇烈時(shí)，也可以通過調(diào)整重疊區(qū)網(wǎng)格的分辨率來提高插值精度，并且這時(shí)滿足插值單元與其宿主單元的尺寸匹配尤為重要。

因此，上述插值方法應(yīng)用簡潔方便并且能夠保證計(jì)算精度。它只需要一層插值單元，就可以充分利用其宿主及所有相鄰單元的信息。并且，因?yàn)榻柚松弦粫r(shí)間步已經(jīng)確定的梯度值，所以同梯度計(jì)算方法一樣，該插值方法適用于任意網(wǎng)格單元類型間的插值計(jì)算。

圖1 兩子網(wǎng)格間邊界定義示例

另外，應(yīng)用本文插值方法，插值邊界可以同其它邊界條件一樣，在復(fù)雜求解器中可以做到透明處理。即迭代計(jì)算中，插值單元與其它邊界虛網(wǎng)格單元功能類似。

3.2 網(wǎng)格間邊界定義方法

本文中網(wǎng)格間邊界定義過程可歸納為兩步：網(wǎng)格分類和邊界優(yōu)化。首先通過以物面距離作為網(wǎng)格分類參數(shù)，將網(wǎng)格單元分為活動(dòng)單元和非活動(dòng)單元，形成初始分類邊界。然后在初始分類的基礎(chǔ)上，對(duì)邊界進(jìn)行優(yōu)化，將網(wǎng)格單元分為活動(dòng)單元、插值單元和非活動(dòng)單元，最終形成插值邊界。邊界優(yōu)化過程主要考慮如下幾點(diǎn)因素：

(1)保證插值單元必須存在宿主單元。因?yàn)椴逯祮卧牧鲌鲂畔⑹菑乃谄渌W(wǎng)格的宿主單元獲得的，所以只有保證存在宿主單元，流場計(jì)算才能正確進(jìn)行。

(2)保證插值單元只存在一個(gè)宿主單元。因?yàn)閷?shí)際應(yīng)用中可能是任意多個(gè)子網(wǎng)格的任意單元類型重疊在一起，所以對(duì)一給定插值單元可能存在多個(gè)宿主單元。

(3)清除孤點(diǎn)。所謂孤點(diǎn)就是與周圍單元的類別不相同的局部小部分網(wǎng)格單元，它可能是活動(dòng)單元，也可能是非活動(dòng)單元。在流場計(jì)算中，孤點(diǎn)的存在不但會(huì)損害計(jì)算精度而且還會(huì)降低收斂效率。

(4)適當(dāng)把插值邊界外移，即保證適當(dāng)?shù)闹丿B區(qū)域，以保證流場計(jì)算中宿主單元重構(gòu)時(shí)足夠的模板需求，從而有利于提高流場重構(gòu)精度。

(5)光滑插值邊界，且盡量保持插值單元和宿主單元的網(wǎng)格尺寸一致，從而有利于提高插值精度。

圖2給出了六圓柱和五球體不同單元類型的邊界定義結(jié)果?？梢?，本文方法定義的網(wǎng)格間邊界光滑，網(wǎng)格尺寸大小匹配，網(wǎng)格重疊區(qū)域大小適當(dāng)，說明本文給出的網(wǎng)格間邊界定義方法是優(yōu)秀可靠的。

3.3 混合網(wǎng)格動(dòng)態(tài)重疊方法

本文混合網(wǎng)格動(dòng)態(tài)重疊方法具體算法如下：

(1)獲取所有子區(qū)域初始網(wǎng)格；

(2)采用2.2節(jié)的混合網(wǎng)格重疊方法裝配網(wǎng)格；

(3)對(duì)所有活動(dòng)區(qū)域，在當(dāng)前物理時(shí)間步上迭代求解非定常NS方程，并積分出運(yùn)動(dòng)物體受到的氣動(dòng)力和氣動(dòng)力矩；

(4)基于運(yùn)動(dòng)物體受到的氣動(dòng)力及氣動(dòng)力矩、其它外力及力矩、重力和其它相關(guān)已知量等，在當(dāng)前物理時(shí)間步求解剛體6DOF運(yùn)動(dòng)方程，得到下一時(shí)刻物體的位移和姿態(tài)角；

(5)根據(jù)步驟(4)的計(jì)算結(jié)果，移動(dòng)運(yùn)動(dòng)物體所在子網(wǎng)格；

(6)判斷非定常計(jì)算是否達(dá)到要求，如果達(dá)到，則結(jié)束計(jì)算；否則返回步驟(2)。

在上述計(jì)算流程中，步驟(2)到(6)不需要任何人工干預(yù)，完全是自動(dòng)完成的。并且實(shí)踐中發(fā)現(xiàn)，網(wǎng)格裝配消耗時(shí)間與網(wǎng)格量以及子網(wǎng)格分布結(jié)構(gòu)有關(guān)，但是一般情況下，網(wǎng)格裝配一次所用時(shí)間與偽時(shí)間步上迭代一次消耗的時(shí)間基本是同一量級(jí)。相對(duì)于整個(gè)非定常計(jì)算過程而言，網(wǎng)格裝配過程所消耗時(shí)間是一個(gè)小量。

圖2 重疊網(wǎng)格邊界定義示例

圖3 三維機(jī)翼與外掛物分離網(wǎng)格

圖4 三維外掛物分離不同時(shí)刻重疊網(wǎng)格及壓力分布云圖(x=5.3截面)

圖5 三維外掛物分離不同時(shí)刻壁面壓力分布云圖

因此，本文采用的基于混合網(wǎng)格動(dòng)態(tài)重疊方法的動(dòng)網(wǎng)格技術(shù)具有如下優(yōu)點(diǎn)：算法簡單可靠，網(wǎng)格裝配質(zhì)量好效率高，處理復(fù)雜外形及復(fù)雜運(yùn)動(dòng)方式能力強(qiáng)，且可以完全實(shí)現(xiàn)自動(dòng)化。

4 算例驗(yàn)證

本節(jié)通過模擬一個(gè)三維機(jī)翼與外掛物的分離過程對(duì)本文數(shù)值方法進(jìn)行驗(yàn)證。該算例是一個(gè)考核動(dòng)網(wǎng)格方法的標(biāo)準(zhǔn)算例，具有詳實(shí)可靠的風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)，其模型由機(jī)翼、掛架和外掛物組成，詳細(xì)幾何參數(shù)及實(shí)驗(yàn)數(shù)據(jù)見參考文獻(xiàn)[12]。

如圖3所示，采用的重疊網(wǎng)格包括三個(gè)子網(wǎng)格，分別是背景網(wǎng)格(共65萬)、繞機(jī)翼的子網(wǎng)格(共193萬)和繞外掛物的子網(wǎng)格(共287萬)。背景網(wǎng)格采用均勻的六面體單元覆蓋整個(gè)分離區(qū)，其它外圍區(qū)域采用四面體填充，繞機(jī)翼和繞外掛物的子網(wǎng)格在壁面附近均采用各向異性網(wǎng)格單元。

本算例模擬物理時(shí)間步長取為0.001s。圖4給出了不同時(shí)刻重疊后的網(wǎng)格及壓力分布云圖，圖5給出了不同時(shí)刻壁面壓力分布云圖?？梢姡S著時(shí)間的推進(jìn)，三個(gè)子網(wǎng)格在不同時(shí)刻自動(dòng)裝配，并且由于均勻背景子網(wǎng)格的存在，重疊區(qū)附近保證了網(wǎng)格單元尺寸大小匹配。計(jì)算的壓力在重疊區(qū)附近過度光滑，分布合理。

圖6給出沿外掛物軸向壓力取值的幾個(gè)不同位置，圖7給出了初始時(shí)刻不同位置的壁面壓力分布對(duì)比，可見數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果比較一致。

圖6 外掛物壓力取值的幾個(gè)不同位置

圖7 外掛物初始時(shí)刻壁面壓力分布

圖8 外掛物運(yùn)動(dòng)參數(shù)隨時(shí)間變化

圖8給出了外掛物運(yùn)動(dòng)參數(shù)隨時(shí)間的變化曲線。從圖中可以清楚地看到彈射力的消失對(duì)外掛物質(zhì)心速度和角速度的影響位置。計(jì)算的外掛物質(zhì)心位移、速度、角速度、滾轉(zhuǎn)角、偏航角及俯仰角均與實(shí)驗(yàn)結(jié)果吻合較好。

5 結(jié)論

(1)本文插值方法依賴于插值單元的格心流場變量值、單元梯度值和限制器函數(shù)值的計(jì)算來完成，只需要一層插值單元，就可以充分利用其宿主及所有相鄰單元的信息。

(2)本文給出的網(wǎng)格間邊界定義方法是優(yōu)秀可靠的，所定義的網(wǎng)格間邊界光滑，網(wǎng)格尺寸大小匹配，網(wǎng)格重疊區(qū)域大小適當(dāng)。

(3)本文采用的基于混合網(wǎng)格動(dòng)態(tài)重疊方法的動(dòng)網(wǎng)格技術(shù)具有如下優(yōu)點(diǎn)：算法簡單可靠，網(wǎng)格裝配質(zhì)量好效率高，處理復(fù)雜外形及復(fù)雜運(yùn)動(dòng)方式能力強(qiáng)，且可以完全實(shí)現(xiàn)自動(dòng)化。

(4)三維外掛物分離問題的數(shù)值計(jì)算表明，本文所發(fā)展的動(dòng)網(wǎng)格方法對(duì)于處理多體相對(duì)運(yùn)動(dòng)模擬問題是快速準(zhǔn)確可靠的。

[1]Kazuhiro Nakahash, Fumiya Togash, Dmitri Sharov. An intergrid-boundary definition method for overset unstructured grid approach. AIAA-1999-3304, 1999.

[2]Benek J A, Steger J L, Dougherty F.A flexible grid embedding technique with applications to the Euler equations.AIAA-1983-1944, 1983.

[3]Hong Luo, Dmitri Sharov,Joseph D Baum.Anoverlappingunstructured gridmethod for viscous flows.AIAA-2001-2603, 2001.

[4]Fumiya Togashi, Yasushi Ito,Kazuhiro Nakahashi. Extensions of overset unstructured grids to multiple bodies in contact.Journalof Aircraft 2006; 43(1): 52-57.

[5]Shun Takahashi, IchieMonjugawa, and Kazuhiro Nakahashi. Unsteady flow computation around moving multiple bodies using overset unstructured grids. AIAA-2006-2839, 2006.

[6]J Blazek.Computational fluid dynamics principles and application. Oxford: Elsevier Sci Ltd, 2001.

[7]Roe P L. Approximate Riemann solvers, parameter vectors, and difference schemes.Journal of Computational Physics 1981; 43: 357-372.

[8]Barth T J, Jespersen D C. The design and application of upwind schemes on unstructured meshes.AIAA-1989-0366, 1989.

[9]Venkatakrishnan V. On the accuracy of limiters and convergence to steady state solutions.AIAA-1993-0880, 1993.

[10]P R Spalart, S R Allmaras. A one-equation turbulence model for aerodynamic flows.AlAA-1992-0439, 1992.

[11]Jameson A..Time-dependent calculations using multigrid with application to unsteady flows past airfoils and wings[R]. AIAA-1991-1596.

[12]HallL. H.,ParthasarathV.Y..Validation of an automated chimera/6-dof methodology for multiple moving body problems[R]. AIAA-1998-0767.

A Dynamic Overset Mixed Grids Approach for Multiple Relative Motion

Kang ZhongLiang1, Fang Yuanyuan2

(1.CABRTECHNOLOGYCO.,LTD.，ChinaAcademyofBuildingResearch,Beijing100013,China；2.AcademyofBeijingGasGroupCo.,Ltd.,Beijing100011,China)

A dynamic overset mixed grids algorithm is proposed. An interpolation method based on the linear reconstruction is developed for the overset mixedgrids approach. The optimized intergrid-boundary definition method is shown to be excellent and reliable. The intergrid boundaries are very smooth and have similar cell sizes.The flows around multiple bodies in relative motion can be simulated by solving the unsteady Navier-Stokes equations and the 6DOF equations. The numerical results of store separation prove that the present dynamic overset grids method is fast, accurate and reliable.

Dynamic Grids; Overset Grids; Mixed Grids; Relative Motion

“十二五”國家科技支撐計(jì)劃(2012BAJ09B04)

康忠良(1981-)，工學(xué)博士。主要從事CFD技術(shù)和BIM技術(shù)研究。

V 211.3

A

1674-7461(2015)01-0069-06