基于AMReX 框架的重疊網(wǎng)格方法研究

向 輝

（南京航空航天大學(xué)，南京 210016）

在計(jì)算流體力學(xué)（Computational Fluid Dynamics，CFD）中，復(fù)雜外形與相對(duì)運(yùn)動(dòng)部件的流動(dòng)模擬問(wèn)題十分常見(jiàn)，這類(lèi)問(wèn)題的難點(diǎn)在于，傳統(tǒng)的網(wǎng)格技術(shù)無(wú)法生成高質(zhì)量的計(jì)算網(wǎng)格。重疊網(wǎng)格方法（overset grid method）是解決該類(lèi)問(wèn)題的重要研究方向。

重疊網(wǎng)格方法對(duì)計(jì)算域各區(qū)域分別生成網(wǎng)格，這大大降低了網(wǎng)格生成的難度。然后將各部分子網(wǎng)格裝配為計(jì)算網(wǎng)格。重疊網(wǎng)格方法最早由Steger 等[1]提出，發(fā)展了結(jié)構(gòu)網(wǎng)格的重疊網(wǎng)格方法。Nakahashi 等[2]將重疊網(wǎng)格方法應(yīng)用到非結(jié)構(gòu)網(wǎng)格上，并采用自動(dòng)裝配方法進(jìn)行網(wǎng)格裝配，提高了計(jì)算效率。這之后，Togashi 等[3]將重疊網(wǎng)格方法應(yīng)用于存在接觸的多個(gè)物體的流動(dòng)模擬問(wèn)題中，并著重對(duì)重疊區(qū)域的處理進(jìn)行了研究。張玉東等[4]對(duì)子母彈分離進(jìn)行了模擬，田書(shū)玲等[5]使用動(dòng)態(tài)重疊網(wǎng)格對(duì)三維多體運(yùn)動(dòng)進(jìn)行了模擬。

重疊網(wǎng)格方法中主要有2 類(lèi)網(wǎng)格：子網(wǎng)格與背景網(wǎng)格。子網(wǎng)格是對(duì)不同部件分別生成的網(wǎng)格，背景網(wǎng)格在整個(gè)計(jì)算域生成，用來(lái)將子網(wǎng)格裝配為整體。背景網(wǎng)格不包含物面，因此可以使用直角坐標(biāo)網(wǎng)格。直角坐標(biāo)網(wǎng)格是一種與直角坐標(biāo)軸嚴(yán)格對(duì)齊的網(wǎng)格。直角坐標(biāo)網(wǎng)格的優(yōu)點(diǎn)在于生成簡(jiǎn)單、守恒性好、計(jì)算量小。為了提高計(jì)算的分辨率，常常對(duì)直角坐標(biāo)網(wǎng)格進(jìn)行加密，文獻(xiàn)[6]指出，對(duì)于層流邊界層，直角坐標(biāo)網(wǎng)格單元數(shù)與貼體網(wǎng)格單元數(shù)之比與Re0.5n成正比。為了減少單元數(shù)量，可以使用自適應(yīng)直角坐標(biāo)網(wǎng)格。Epstein 等[7]提出了局部網(wǎng)格加密的方法優(yōu)化解的重構(gòu)，Aftosmis[8]在自適應(yīng)直角坐標(biāo)網(wǎng)格上使用了迎風(fēng)差分格式，這之后，Da rren 等[9]提出了直角坐標(biāo)網(wǎng)格自適應(yīng)的四叉樹(shù)數(shù)據(jù)結(jié)構(gòu)，在物體物面和流場(chǎng)中梯度較大的區(qū)域都進(jìn)行了網(wǎng)格自適應(yīng)。對(duì)二維定常流動(dòng)進(jìn)行了模擬，結(jié)果表明該自適應(yīng)方法兼顧了靈活性與計(jì)算精度。

本文基于AMReX 框架[10]開(kāi)展重疊網(wǎng)格方法的研究，發(fā)展了一套基于自適應(yīng)直角網(wǎng)格的重疊網(wǎng)格算法。

1 流動(dòng)數(shù)值解法

1.1 控制方程

黏性流體的流動(dòng)控制方程的積分形式[11]為

式中：W，F(xiàn)c（W），F(xiàn)v（W）分別為守恒變量、對(duì)流通量和黏性通量；ρ 為密度；u，v，w 為速度分量；E，H 分別為單位質(zhì)量流體的總能與總焓；nx，ny，nz為控制體單元表面的外法向量的分量；V 為逆變速度，V=nxu+nyv+nzw。為了使方程封閉，式（1）要和理想氣體狀態(tài)方程聯(lián)立

式中：R 為氣體常數(shù)，T 為溫度。為了模擬湍流，使用Spalart-Allmars（S-A）一方程模型[12]作為湍流模型計(jì)算黏性通量。

1.2 空間離散方法

本文采用有限體積方法，使用格心格式，在控制體i 上對(duì)式（1）進(jìn)行空間離散

式中：Ri表示控制體i 的殘差，Sij表示控制體i 的第j個(gè)面的面矢量，NF表示控制體表面的數(shù)量。對(duì)流通量使用HLLC 格式[13]計(jì)算。

1.3 時(shí)間離散方法

1.3.1 定常流動(dòng)

直角坐標(biāo)網(wǎng)格上使用Runge-Kutta 方法[14]進(jìn)行時(shí)間離散，對(duì)式（3）進(jìn)行時(shí)間離散

1.3.2 非定常流動(dòng)

非定常流動(dòng)是CFD 中經(jīng)常遇到的一類(lèi)問(wèn)題，一般使用雙時(shí)間步長(zhǎng)方法，該方法最早由Jameson 等[15]提出。對(duì)式（3）的時(shí)間導(dǎo)數(shù)項(xiàng)使用向后三點(diǎn)差分

式中：n-1，n，n+1 分別表示上一時(shí)間步、當(dāng)前時(shí)間步和下一時(shí)間步。Δt 表示全局時(shí)間步長(zhǎng)。設(shè)置偽時(shí)間步，在偽時(shí)間步上，式（4）可以寫(xiě)成

式中：t*表示偽時(shí)間，W*是的估計(jì)值。當(dāng)?W*/?t*=0 時(shí)，式（6）變?yōu)槭剑?），這說(shuō)明式（6）的穩(wěn)態(tài)解就是式（5）的解。本文使用四步Runge-Kutta 方法進(jìn)行偽時(shí)間推進(jìn)。

2 直角坐標(biāo)網(wǎng)格自適應(yīng)方法

本文基于AMReX 框架進(jìn)行了直角坐標(biāo)網(wǎng)格自適應(yīng)方法的研究，采用多層嵌套的矩形網(wǎng)格，不同層的網(wǎng)格之間是層層加密的。

在網(wǎng)格初始化階段，設(shè)置網(wǎng)格加密層數(shù)，并對(duì)所有層進(jìn)行加密，并在每一層選擇需要的計(jì)算區(qū)域。其中，最底層網(wǎng)格覆蓋全部計(jì)算域，其余各層網(wǎng)格的計(jì)算區(qū)域應(yīng)在上一層計(jì)算區(qū)域之內(nèi)。網(wǎng)格的加密主要分為2部分：初始網(wǎng)格加密與網(wǎng)格自適應(yīng)。初始網(wǎng)格加密是指在計(jì)算之前，人為對(duì)物面附近的網(wǎng)格進(jìn)行加密，可以減少迭代步數(shù)。網(wǎng)格自適應(yīng)是指在計(jì)算中，對(duì)流場(chǎng)劇烈變化的區(qū)域或者渦等流動(dòng)特征出現(xiàn)的區(qū)域網(wǎng)格進(jìn)行加密，對(duì)流動(dòng)平緩的區(qū)域網(wǎng)格進(jìn)行稀疏，自適應(yīng)操作只發(fā)生在物理時(shí)間步。

在進(jìn)行時(shí)間推進(jìn)時(shí)，首先在粗網(wǎng)格上計(jì)算，然后下一層網(wǎng)格上計(jì)算，一般使用兩步推進(jìn)，即在粗網(wǎng)格上推進(jìn)Δt，細(xì)網(wǎng)格上先推進(jìn)0.5Δt，再推進(jìn)0.5Δt。

3 重疊網(wǎng)格方法

3.1 網(wǎng)格單元分類(lèi)

本文發(fā)展的重疊網(wǎng)格方法，其子網(wǎng)格使用混合非結(jié)構(gòu)網(wǎng)格，背景網(wǎng)格使用自適應(yīng)直角坐標(biāo)網(wǎng)格。網(wǎng)格裝配方法使用基于物面距的自動(dòng)裝配方法[16]。網(wǎng)格之間的重疊關(guān)系主要有2 種：子網(wǎng)格與子網(wǎng)格、子網(wǎng)格與背景網(wǎng)格。

子網(wǎng)格與子網(wǎng)格的重疊區(qū)域如圖1 所示，淺色單元i 與深色單元j 都屬于網(wǎng)格grid1，Si，Sdi分別表示單元i 格心到grid1 和grid2 的物面距，Sj，Sdj分別表示單元j 格心到grid1 和grid2 的物面距。其中，Si＜Sdi，淺色單元是活動(dòng)單元，參與流場(chǎng)計(jì)算；Sj＜Sdj，深色單元是非活動(dòng)單元，不參與流場(chǎng)計(jì)算，需要去除。

圖1 子網(wǎng)格重疊情況

背景網(wǎng)格與子網(wǎng)格之間的重疊與上一種情況有所不同。背景網(wǎng)格不包含物面，需要人為設(shè)置參考距離ΔL。分別對(duì)子網(wǎng)格與背景網(wǎng)格單元分類(lèi)方法進(jìn)行說(shuō)明。

所有網(wǎng)格（包括子網(wǎng)格與背景網(wǎng)格）單元到背景網(wǎng)格的物面距均設(shè)置為ΔL。記子網(wǎng)格單元到自身的物面距為d0，如果d0≤ΔL，表明該子網(wǎng)格單元離物面較近，是活動(dòng)單元，否則為非活動(dòng)單元；記背景網(wǎng)格單元到其他子網(wǎng)格的物面距為dw，如果ΔL＜dw，表明該背景單元離物面較遠(yuǎn)，是活動(dòng)單元，否則為非活動(dòng)單元。如圖2所示，P、A、B、C 都是背景網(wǎng)格單元中心，實(shí)線圓表示物體的物面，矩形區(qū)域（包括被遮擋的區(qū)域）表示背景網(wǎng)格。按照前面的分類(lèi)方法，P、C 所在單元為活動(dòng)單元，A、B 所在單元為非活動(dòng)單元。其中點(diǎn)C 所在單元位于物體內(nèi)部，實(shí)際上不參與流場(chǎng)計(jì)算，對(duì)于此類(lèi)點(diǎn)應(yīng)該予以去除。

另外，為了保證一定可以搜索到貢獻(xiàn)單元，可以將與活動(dòng)單元相鄰的非活動(dòng)單元也設(shè)為活動(dòng)單元。

3.2 插值單元與貢獻(xiàn)單元

不同網(wǎng)格上的流場(chǎng)信息在重疊區(qū)域進(jìn)行交流，插值單元的作用就是通過(guò)插值得到其他網(wǎng)格的流場(chǎng)信息。本文將與活動(dòng)單元相鄰的非活動(dòng)單元定義為插值單元，為了實(shí)現(xiàn)高階數(shù)值格式，可以設(shè)置多層插值單元。圖3 展示了定義2 層插值單元的情形，1 表示活動(dòng)單元，0 表示非活動(dòng)單元，-1 表示第一層插值單元，-2表示第二層插值單元，洞邊界是活動(dòng)單元與非活動(dòng)單元的分界面。

圖3 插值單元示意圖

插值單元中流場(chǎng)變量的計(jì)算需要相應(yīng)的貢獻(xiàn)單元的流場(chǎng)值，貢獻(xiàn)單元與插值單元屬于不同網(wǎng)格。定義插值單元的貢獻(xiàn)單元，首先要找到插值單元的宿主網(wǎng)格單元，宿主單元指包含插值單元中心，且與插值單元不屬于同一個(gè)部件網(wǎng)格的網(wǎng)格單元。宿主單元的搜索使用交替數(shù)字二叉樹(shù)（alternating digital tree，ADT）方法[17]。在混合網(wǎng)格中，貢獻(xiàn)單元就是宿主單元。在直角坐標(biāo)網(wǎng)格中，貢獻(xiàn)單元是由距離插值單元中心最近的4 個(gè)或8 個(gè)單元中心組成的矩形或長(zhǎng)方體。

4 算例驗(yàn)證

本節(jié)對(duì)混合網(wǎng)格上的計(jì)算方法與基于自適應(yīng)直角坐標(biāo)網(wǎng)格的重疊網(wǎng)格算法進(jìn)行數(shù)值驗(yàn)證。

4.1 ONERA M6

使用標(biāo)準(zhǔn)的ONERA M6 機(jī)翼算例[18]對(duì)混合網(wǎng)格流場(chǎng)計(jì)算方法進(jìn)行數(shù)值驗(yàn)證。計(jì)算條件為來(lái)流馬赫數(shù)Ma∞=0.839 5，機(jī)翼攻角a=3.06°，雷諾數(shù)Re=1.172×107。將計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比，分別繪制機(jī)翼展向0.2、0.65、0.95 處壓強(qiáng)系數(shù)分布曲線。如圖4 所示。

圖4 壓強(qiáng)系數(shù)分布曲線

從結(jié)果可以看出，數(shù)值模擬的結(jié)果與實(shí)驗(yàn)數(shù)據(jù)符合得很好，這表明混合網(wǎng)格上的計(jì)算方法具有較高的精度。

4.2 NACA0015 低速流動(dòng)

為了驗(yàn)證基于自適應(yīng)直角坐標(biāo)網(wǎng)格的重疊網(wǎng)格算法的有效性，使用由NACA0015 翼型拉伸而成的直機(jī)翼，展弦比為6.6。網(wǎng)格裝配情況如圖5（a）所示，其中深色區(qū)域是混合網(wǎng)格，淺色區(qū)域?yàn)橹丿B區(qū)域，其余區(qū)域?yàn)橹苯亲鴺?biāo)網(wǎng)格。計(jì)算條件如下：來(lái)流馬赫數(shù)Ma∞=0.21，雷諾數(shù)Re=1.5×106，機(jī)翼攻角a=12°，側(cè)滑角β=0°。

圖5 重疊網(wǎng)格裝配情況及模擬結(jié)果云圖

該算例屬于非定常流動(dòng)，計(jì)算總時(shí)間t=0.58 s。在進(jìn)行時(shí)間推進(jìn)時(shí)，需要設(shè)置全局時(shí)間步長(zhǎng)。在混合網(wǎng)格上采用雙時(shí)間步長(zhǎng)方法，全局時(shí)間步長(zhǎng)△t=0.000 05 s，將混合網(wǎng)格的時(shí)間步定為物理時(shí)間步；在自適應(yīng)直角坐標(biāo)網(wǎng)格上使用顯式Runge-Kutta 方法。

計(jì)算結(jié)果如圖5（b）所示，計(jì)算時(shí)，翼尖渦附近進(jìn)行了網(wǎng)格自適應(yīng)，算法對(duì)渦這類(lèi)流動(dòng)特征實(shí)現(xiàn)了準(zhǔn)確捕捉，這表明本文發(fā)展的基于自適應(yīng)直角坐標(biāo)網(wǎng)格的重疊網(wǎng)格方法是有效的，具有較高的實(shí)用價(jià)值。

5 結(jié)束語(yǔ)

本文發(fā)展了基于自適應(yīng)直角坐標(biāo)網(wǎng)格的重疊網(wǎng)格方法，背景網(wǎng)格使用自適應(yīng)直角坐標(biāo)網(wǎng)格，減少了網(wǎng)格生成的難度，同時(shí)增強(qiáng)了對(duì)流動(dòng)特征的捕捉能力，子網(wǎng)格使用混合網(wǎng)格，減少了網(wǎng)格總數(shù)，提高了計(jì)算精度，具有較高的工程應(yīng)用價(jià)值。