基于AMReX 框架的重疊網格方法研究

向 輝

（南京航空航天大學，南京 210016）

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

重疊網格方法對計算域各區域分別生成網格，這大大降低了網格生成的難度。然后將各部分子網格裝配為計算網格。重疊網格方法最早由Steger 等[1]提出，發展了結構網格的重疊網格方法。Nakahashi 等[2]將重疊網格方法應用到非結構網格上，并采用自動裝配方法進行網格裝配，提高了計算效率。這之后，Togashi 等[3]將重疊網格方法應用于存在接觸的多個物體的流動模擬問題中，并著重對重疊區域的處理進行了研究。張玉東等[4]對子母彈分離進行了模擬，田書玲等[5]使用動態重疊網格對三維多體運動進行了模擬。

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

本文基于AMReX 框架[10]開展重疊網格方法的研究，發展了一套基于自適應直角網格的重疊網格算法。

1 流動數值解法

1.1 控制方程

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

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

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

1.2 空間離散方法

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

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

1.3 時間離散方法

1.3.1 定常流動

直角坐標網格上使用Runge-Kutta 方法[14]進行時間離散，對式（3）進行時間離散

1.3.2 非定常流動

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

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

式中：t*表示偽時間，W*是的估計值。當?W*/?t*=0 時，式（6）變為式（5），這說明式（6）的穩態解就是式（5）的解。本文使用四步Runge-Kutta 方法進行偽時間推進。

2 直角坐標網格自適應方法

本文基于AMReX 框架進行了直角坐標網格自適應方法的研究，采用多層嵌套的矩形網格，不同層的網格之間是層層加密的。

在網格初始化階段，設置網格加密層數，并對所有層進行加密，并在每一層選擇需要的計算區域。其中，最底層網格覆蓋全部計算域，其余各層網格的計算區域應在上一層計算區域之內。網格的加密主要分為2部分：初始網格加密與網格自適應。初始網格加密是指在計算之前，人為對物面附近的網格進行加密，可以減少迭代步數。網格自適應是指在計算中，對流場劇烈變化的區域或者渦等流動特征出現的區域網格進行加密，對流動平緩的區域網格進行稀疏，自適應操作只發生在物理時間步。

在進行時間推進時，首先在粗網格上計算，然后下一層網格上計算，一般使用兩步推進，即在粗網格上推進Δt，細網格上先推進0.5Δt，再推進0.5Δt。

3 重疊網格方法

3.1 網格單元分類

本文發展的重疊網格方法，其子網格使用混合非結構網格，背景網格使用自適應直角坐標網格。網格裝配方法使用基于物面距的自動裝配方法[16]。網格之間的重疊關系主要有2 種：子網格與子網格、子網格與背景網格。

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

圖1 子網格重疊情況

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

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

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

3.2 插值單元與貢獻單元

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

圖3 插值單元示意圖

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

4 算例驗證

本節對混合網格上的計算方法與基于自適應直角坐標網格的重疊網格算法進行數值驗證。

4.1 ONERA M6

使用標準的ONERA M6 機翼算例[18]對混合網格流場計算方法進行數值驗證。計算條件為來流馬赫數Ma∞=0.839 5，機翼攻角a=3.06°，雷諾數Re=1.172×107。將計算結果與實驗結果對比，分別繪制機翼展向0.2、0.65、0.95 處壓強系數分布曲線。如圖4 所示。

圖4 壓強系數分布曲線

從結果可以看出，數值模擬的結果與實驗數據符合得很好，這表明混合網格上的計算方法具有較高的精度。

4.2 NACA0015 低速流動

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

圖5 重疊網格裝配情況及模擬結果云圖

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

計算結果如圖5（b）所示，計算時，翼尖渦附近進行了網格自適應，算法對渦這類流動特征實現了準確捕捉，這表明本文發展的基于自適應直角坐標網格的重疊網格方法是有效的，具有較高的實用價值。

5 結束語

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