向 輝

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

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

重疊網(wǎng)格方法對計算域各區(qū)域分別生成網(wǎng)格，這大大降低了網(wǎng)格生成的難度。然后將各部分子網(wǎng)格裝配為計算網(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)格上，并采用自動裝配方法進(jìn)行網(wǎng)格裝配，提高了計算效率。這之后，Togashi 等[3]將重疊網(wǎng)格方法應(yīng)用于存在接觸的多個物體的流動模擬問題中，并著重對重疊區(qū)域的處理進(jìn)行了研究。張玉東等[4]對子母彈分離進(jìn)行了模擬，田書玲等[5]使用動態(tài)重疊網(wǎng)格對三維多體運(yùn)動進(jìn)行了模擬。

重疊網(wǎng)格方法中主要有2 類網(wǎng)格：子網(wǎng)格與背景網(wǎng)格。子網(wǎng)格是對不同部件分別生成的網(wǎng)格，背景網(wǎng)格在整個計算域生成，用來將子網(wǎng)格裝配為整體。背景網(wǎng)格不包含物面，因此可以使用直角坐標(biāo)網(wǎng)格。直角坐標(biāo)網(wǎng)格是一種與直角坐標(biāo)軸嚴(yán)格對齊的網(wǎng)格。直角坐標(biāo)網(wǎng)格的優(yōu)點(diǎn)在于生成簡單、守恒性好、計算量小。為了提高計算的分辨率，常常對直角坐標(biāo)網(wǎng)格進(jìn)行加密，文獻(xiàn)[6]指出，對于層流邊界層，直角坐標(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ù)據(jù)結(jié)構(gòu)，在物體物面和流場中梯度較大的區(qū)域都進(jìn)行了網(wǎng)格自適應(yīng)。對二維定常流動進(jìn)行了模擬，結(jié)果表明該自適應(yīng)方法兼顧了靈活性與計算精度。

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

1 流動數(shù)值解法

1.1 控制方程

黏性流體的流動控制方程的積分形式[11]為

式中：W，F(xiàn)c（W），F(xiàn)v（W）分別為守恒變量、對流通量和黏性通量；ρ 為密度；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]作為湍流模型計算黏性通量。

1.2 空間離散方法

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

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

1.3 時間離散方法

1.3.1 定常流動

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

1.3.2 非定常流動

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

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

式中：t*表示偽時間，W*是的估計值。當(dāng)?W*/?t*=0 時，式（6）變?yōu)槭剑?），這說明式（6）的穩(wěn)態(tài)解就是式（5）的解。本文使用四步Runge-Kutta 方法進(jìn)行偽時間推進(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ù)，并對所有層進(jìn)行加密，并在每一層選擇需要的計算區(qū)域。其中，最底層網(wǎng)格覆蓋全部計算域，其余各層網(wǎng)格的計算區(qū)域應(yīng)在上一層計算區(qū)域之內(nèi)。網(wǎng)格的加密主要分為2部分：初始網(wǎng)格加密與網(wǎng)格自適應(yīng)。初始網(wǎng)格加密是指在計算之前，人為對物面附近的網(wǎng)格進(jìn)行加密，可以減少迭代步數(shù)。網(wǎng)格自適應(yīng)是指在計算中，對流場劇烈變化的區(qū)域或者渦等流動特征出現(xiàn)的區(qū)域網(wǎng)格進(jìn)行加密，對流動平緩的區(qū)域網(wǎng)格進(jìn)行稀疏，自適應(yīng)操作只發(fā)生在物理時間步。

在進(jìn)行時間推進(jìn)時，首先在粗網(wǎng)格上計算，然后下一層網(wǎng)格上計算，一般使用兩步推進(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)格單元分類

本文發(fā)展的重疊網(wǎng)格方法，其子網(wǎng)格使用混合非結(jié)構(gòu)網(wǎng)格，背景網(wǎng)格使用自適應(yīng)直角坐標(biāo)網(wǎng)格。網(wǎ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，淺色單元是活動單元，參與流場計算；Sj＜Sdj，深色單元是非活動單元，不參與流場計算，需要去除。

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

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

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

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

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

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

圖3 插值單元示意圖

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

4 算例驗證

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

4.1 ONERA M6

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

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

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

4.2 NACA0015 低速流動

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

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

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

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

5 結(jié)束語

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