可壓縮多介質流動問題的高精度數值模擬方法

何志偉, 田保林,2,*, 李 理, 李海鋒, 張又升, 孟寶清

(1. 北京應用物理與計算數學研究所, 北京 100094; 2. 北京大學 應用物理與技術研究中心, 北京 100871)

0 引 言

可壓縮多介質流動問題，是指流場涉及多種物質，且物質之間存在相互耦合作用的流動現象。它廣泛存在于各種自然現象、工程應用中。如超新星爆發、超燃沖壓發動機中的燃料混合、慣性約束聚變和武器物理等[1]，有著非常重要的科學研究意義和工程應用價值。

在這種流動問題中，由于移動且可變形界面的存在，界面處流體性質不連續給物理建模與數值方法帶來了極大的困難。目前，文獻中存在多種具有不同復雜度和適用范圍的物理模型及數值計算方法。從界面處理方式來看，它們可以分為兩類[2]：第一類將界面處理為清晰的間斷，即清晰界面模型/方法；第二類將將界面處理為多物質相互擴散而成的混合區，即擴散界面模型/方法。

在擴散界面模型中，界面被處理為多物質相互擴散而成的混合區，就像數值捕獲空氣動力學中的接觸間斷一樣[2]。雖然這類擴散界面是一種人為做成的數值混合，但這類模型最重要的優點就是可以在固定網格上，所有計算單元都采用相同的數值方案來求解，而不用考慮流場中存在什么樣類型的波(沖擊波、稀疏波，界面等)以及它們之間的相互作用。這類模型/方法可進一步分為兩小類：基于歐拉/Navier-Stokes方程的擴散界面模型和基于多相流方程的擴散界面模型。

1) 基于單流體歐拉/Navier-Stokes方程的擴散界面模型。這類模型是由傳統的歐拉/Navier-Stokes方程，附加每種物質的(質量)濃度方程構成，每個計算網格即使含有多種流體組分仍具有相同的運動速度。這類方法雖然形式上簡單，但卻存在一些問題：①過界面非物理振蕩。Abgrall通過數值實驗表明[3]，對于兩種理想氣體的多介質問題，直接采用這類模型會在界面處產生壓力振蕩。為了消除這個問題，他提出了非守恒的四方程模型(擬守恒方法)[3]及等效四方程模型[4]，在恰當數值離散的基礎上，可以使用同種有限體積/差分算法來處理由界面分開的可壓縮多介質流動問題[5]。但是這樣的操作最終導致方程組喪失守恒性。② 等溫等壓假設限制此模型的應用范圍。這類模型通過兩種組分的狀態方程構造混合物的狀態方程。壓強均衡會導致兩相的聲速必須耦合，而且顯然無法滿足某些真實物理情況。

2) 基于多相流方程的擴散界面模型。多介質流動是一類含有可區分大尺度物質界面的多相流問題，因此很多學者也直接將多介質流動稱為多相流動。多相流方程的擴散界面模型是指將所有相視為連續介質，利用各種平均方法(時間平均、空間平均、系綜平均)對各相的場方程進行平均后得到的一類忽略了流動的一些微觀細節的宏觀均質或平均混合模型[6]。由于使用平均方法，最終的模型會出現非守恒項(即最終的模型無法寫成散度形式)和較多的附加項(物質之間的質量、動量和能量傳遞等過程的宏觀模化)。對這些附加項建立不同的唯象模式，就形成了不同的模型(對于兩種介質的流動，又稱雙流體模型或七方程模型)，如Baer-Nunziato模型(BN模型)[7]、Saurel-Abgrall模型(SA模型)[8]。這類模型具體如下一些特點：① 不同相之間是完全非平衡的，即每個相流體都有自己的密度、壓力、速度、溫度等，因此一維情況下共有7個方程。由于這類模型包含的方程個數多、波系復雜、存在非守恒項[10]，有時候并不能滿足嚴格雙曲特性，給數值求解帶來了很大挑戰。② 雙流體多相流模型一般只適用于稠密多相流情形，即可以離散可近似看做一種連續的流體(如大量顆粒、液滴等)。對于顆粒相稀疏或者處于稀疏到稠密過渡情形，雙流體多相流模型的適應性還存在問題(為克服這一困難，最近我們提出了一種適用于稀疏到稠密跨流態可壓縮兩相流的模擬方法(CMP-PIC)[9])。③ 對于一些特殊情況，這類模型還可以退化為一系列簡化模型[1]。

在實際的工程應用中，我們會根據具體問題的特殊性，針對性選擇一些簡化模型來開展相應的數值模擬工作。但不管選擇何種類型的模型，為建立一套高效的高精度數值模擬算法，我們遇到的問題都是基本類似的。因此，本文內容不局限于某種特殊模型，而將從數值框架的選擇、非守恒方程相容離散、高精度有界格式構造、界面壓縮、界面-單介質分區計算方法等5個維度出發，綜述近幾年我們在可壓縮多介質流動問題高精度數值模擬工作中的一些進展。

1 高精度數值模擬方法

1.1 基本要求

為了準確模擬可壓縮多介質流動問題，尤其是其中涉及的界面大變形運動，需要探索和建立一套高效的高精度數值算法。同面向其他可壓縮流動問題的數值算法一樣，這套算法同樣需要滿足一定的要求，如守恒性、高分辨率、健壯性等。但可壓縮多介質流動問題又有其特殊性，因此有一些額外的特殊要求，其中包括：1) 過界面不產生各種非物理的數值振蕩；2) 在長時間計算過程中，保持界面的銳利結構。界面不能隨著計算時間的增加而不斷彌散變寬，最終影響流場其他結構；3) 保界性。可壓縮多介質流動的物理模型中往往會涉及其他類型的物理量，如質量分數、體積分數。這類物理量有著明確的上界和下界。因此，在對可壓縮多介質流動問題進行模擬時，不僅僅需要考慮熱力學量的保正問題，還得同時考慮如質量分數這類物理量的保界問題。從而保證程序的穩定運行和界面混合封閉模式的精準執行。

為了實現構建滿足上述約束條件的數值方法，我們做了一系列探索和研究工作。下面詳細介紹研究的進展以及一些具體工作。

1.2 數值框架的選擇

流體力學數值方法通常可以分為拉格朗日方法和歐拉方法兩種類型。拉格朗日方法，計算網格隨流體一起運動，適宜計算一些變形不大的多介質界面問題。歐拉方法，網格固定不動，能夠模擬界面大變形問題。本文主要關注歐拉數值方法。目前歐拉數值離散方法大致分為有限差分方法、有限體積方法、有限元方法三種基本框架。這些數值框架能否適用于可壓縮多介質流動問題，需要做適應性分析。而對這類框架開展適用性分析的指導思想是Abgrall提出的界面無振蕩分析方法[3]——過界面速度-壓力需要保持為常數。

1.2.1 過界面非物理振蕩現象

考慮如下經典的滿足氣體狀態方程的界面運動問題(見圖1)。左邊是一種氣體，如空氣，其比熱比為γL=1.4。右邊是另一種氣體，如氦氣，其比熱比為γL=1.667。跨過界面時，兩種介質的壓力和速度均相同。

圖1 界面運動問題

對于這個問題，采用如下模型(一維形式)[3]：

其中，ρ、u、E、p、YL分別為密度、速度、總能、壓力和左邊物質的質量分數。其中質量分數滿足YL+YR=1，且混合物的內能和壓力之間滿足：

(2)

其中WL和WR分別是作用物質的摩爾質量。對于這個模型，可以采用經典的有限體積方法或有限差分方法進行數值離散。但是大量數值結果表明，過界面的壓力和速度無法繼續保持常數。此即所謂的界面壓力振蕩現象。

目前，此現象的數值機理已經很明確。即，傳統的數值方法具有數值黏性，會產生界面處的數值擴散，導致界面處的不同介質進行了偽混合。在這種偽混合區域，數值方法計算出的壓力和溫度是不準確的。此過程可以用圖2更清晰地進行闡述。如圖2所示，初始狀態1和狀態2在等壓線上，但是由于數值黏性的存在，界面處會產生新的混合狀態3，而這個新狀態如果不在等壓線上，那么一定會出現壓力振蕩現象。

圖2 數值混合的不兼容性示意圖

在文獻[11]中，作者證明，當介質具有p=Aρe+B(其中A和B是常數)這樣類型的狀態方程時，采用傳統的高精度差分方法是可以實現單介質接觸間斷處壓力無振蕩屬性。但是當介質具有非線性狀態方程，或者多介質混合問題(兩種氣體的混合造成混合物的狀態方程呈現非線性特點)，接觸間斷/界面處一定會出現壓力振蕩現象。

1.2.2 數值方案引起的非物理振蕩機制分析

除了上述根本原因外，不同的數值方案會額外導致過界面非物理振蕩。為此，我們以有限差分方法為例[12]，開展了數值方案引起的非物理數值振蕩的分析工作。

基于通量分裂的高精度差分方法[13]，由于算法簡單，容易實現高精度計算，因此得到了大量應用。我們在文獻[11-12,15]中詳細分析了此套方案應用于多介質流動問題時引起的非物理振蕩現象的數值機制。

考慮Abgrall提出的擬守恒方法[3]：

通過求解最后一個非守恒方程，避免了使用混合封閉模式(2)。此方法的本質在于將原有的等溫假設修正為等壓假設。這種修正對多介質界面運動問題是合理的。對于擬守恒方法的控制方程，我們可以將其寫為：

其中，

u=(ρ,ρu,ρE,1/(γ-1))T

f=(ρu,ρu2+p,ρEu+pu,0)T

h=diag(0,0,0,u)

(5)

采用傳統的基于通量分裂的高精度差分方法[13]，可以得到如下的半離散形式：

對于逐特征場有限差分方法，需要在通量分裂后，通過局部特征分解后在每個特征場上應用差分格式獲得相應的特征場數值通量，進而通過變換最終獲得物理空間的數值通量。

在文獻[12]，我們系統分析了逐方程有限差分方法在求解多介質界面問題時存在的問題。

首先，假定采用的差分格式為線性格式，通過推導發現，傳統的Steger-Warming通量分裂方法[14]無法保持界面處壓力不振蕩。只有Lax-Friedrichs類型通量分裂方法可以滿足質量-動量-能量方程的分裂得到的正通量之間/負通量之間保持相容性，其中global Lax-Friedrichs通量分裂方法[13]是最簡單有效的方法。具體推導參見文獻[12]。

其次，在固定使用global Lax-Friedrichs通量分裂方法的情況下，有限差分方法需要采用非線性WENO差分格式[13]來計算各個方程的數值通量。但是，由于每個方程的通量的量級大小不同，每個方程所計算出的WENO格式的非線性權系數不同，進而導致質量-動量-能量方程所獲得的數值通量之間不滿足相容性，進而造成界面處的壓力振蕩現象。具體推導參見文獻[12]。

進一步，為了使方程之間的非線性WENO離散具有兼容性，我們提出了一套“公共權”方案。即質量-動量-能量方程分別計算得到的WENO格式的非線性權修正為所有方程采用一套WENO權系數，這樣就可以保證質量-動量-能量方程所獲得的數值通量之間的兼容性。在文獻[12]中，我們給出了一種具體的設計方法。但是更為魯棒的公共權方法需要進一步研究。

最后，根據界面輸運問題中速度和壓力要保持常數這一原則，推導出了非守恒方程相容離散形式為：

其中，φ=1/(γ-1)，α為global Lax-Friedrichs通量分裂方法中所取的全局常數。

對于雙曲系統，逐特征場的有限差分方法更加匹配雙曲系統。當使用局部特征分解時，我們發現[11,15]：上述的WENO公共權方案的約束條件可以得到部分放松；質量-動量-能量方程之間的相容性會得到更好的保持。在文獻[11]中，詳細證明了只要非線性特征場和描述界面演化的線性場采用公共權技術，就可以實現過界面壓力無振蕩屬性。但是設計滿足這種條件的公共權技術仍然是一個懸而未決的問題。

基于前期的工作，我們在文獻[15]中最終給出了避免界面壓力振蕩的有限差分方法的一般性的處理方案[15]。具體如下：

首先，將擬守恒方程寫成“守恒+源項”的如下具體形式：

其中,

采用上述的方法結合局部特征分解技術可以保證界面輸運問題壓力-速度保持常數的屬性，具體證明見文獻[15]。

再其次，提出一個新的準則[15]“單介質保持準則”——多介質算法應能自動保持單介質流動屬性。利用這個性質，推導出了源項的兼容離散形式為：

至此，就完成了新算法的推導過程。此推導過程清楚地表明,新算法不需要對非保守項進行任何特殊處理。它成功地避免了設計WENO格式的公共權重技術。此外，它可以直接應用于任何差分格式，而不局限于WENO格式。更詳細的內容請參考文獻[15]。這里以如下的界面輸運問題為例給出具體效果。此問題的初始條件為：

此問題的具體設置見文獻[15]。

從圖3中可以看出，采用五階精度的WENO格式[13]的差分方法直接求解守恒型多組分方程時，速度和壓力都會產生非物理數值振蕩。相反，當使用新提出的算法求解擬守恒方法時，獲得了與精確解的吻合較好的數值結果。

圖3 界面輸運問題的數值結果[15]

相比于只能采用通量分裂形式的有限差分方法，有限體積框架則具有更多的靈活性。如在重構過程中，可以采用原始變量、守恒變量、特征量變等。在計算通量的過程中，可以采用各種黎曼求解器、矢通量分裂方法等[16]。Abgrall在最初的文章中就使用了基于原始變量重構的有限體積方法。Johnsen[17]進一步確認了采用原始變量重構的有限體積方法可以有效避免數值框架帶來的非物理數值振蕩。

這里需要特別指出，Deng等提出的基于半點值的緊致型有限差分方法[18-19]，由于可以借助有限體積中的重構和通量計算策略，因此也比較容易拓展應用于可壓縮多介質流動問題的數值模擬中。對此問題，Nonomura等做了相關分析[20]。

1.3 非守恒方程相容離散

在可壓縮多介質流動問題中，我們會遇到一些額外類型的方程。典型的如上節里提到的對流方程。這類方程與傳統的雙曲守恒方程不同，它們無法寫成守恒形式。從上一節的分析還可以看出：這類方法如果不恰當離散，同樣會導致過界面的非物理振蕩。因此對于一個包含守恒方程和非守恒方程的可壓縮多介質流動問題的控制方程，需要進一步解決的問題是：當守恒方程采用確定的雙曲守恒離散方法(如有限差分方法、有限體積方法)后，這些非守恒方程如何進行相容離散。

其實，在上一節數值框架的分析過程中，已經給出了一種有限差分方法框架下的離散方案。此外，對于有限體積方法，不同研究者也已經提出了一些不同方法[15,17,21]。最近，我們進一步拓展了此工作，給出了一套統一的離散方案。新方案自動兼容我們提出的有限差分離散方法[15]及其他部分研究者提出的有限體積離散方法[21]。此外，新方案還可以進一步導出一些黎曼求解器對應的新離散方案。關于這部分工作，詳見文獻[22]。

這里需要指出的是，針對一般性非守恒系統， Dal Maso、LeFloch和Murat的提出了DLM理論[49]，其中跨間斷的非守恒乘積項被定義為連接左右狀態的沿路徑積分。基于上述理論發展了所謂的路徑守恒離散方案[50-51]。但是，這類方法在路徑的選擇、收斂性等方面還沒有確切結論[52]，仍需進一步發展。

1.4 高精度有界數值格式構造

不同于單介質歐拉方程和N-S方程，可壓縮多介質流動的物理模型中會涉及其他類型的物理量，如質量分數、體積分數。這類物理量有著明確的上界和下界，直接影響界面混合封閉模式的精準執行。因此，在對可壓縮多介質流動問題進行模擬時，我們不僅需要考慮熱力學量的保正問題，還得同時考慮如質量分數這類物理量的保界問題。因此，除了要求具備高精度、基本無振蕩屬性等基本性質外，如何實現物理量的保界/保正要求是我們選擇和構造可壓縮多介質流動問題的高精度數值格式的額外指導原則。

目前，經數值實驗驗證比較成功的高精度數值格式為WENO類型格式[13]。這類格式在實現高精度計算的同時，也實現了對各種間斷結構的基本無振蕩捕捉。但是WENO格式同樣存在一些缺點：如WENO-JS格式耗散偏大，導致對一些多尺度問題(如激波-湍流邊界層干擾問題、聲學問題)的計算效率過低。為此，相關研究者進行了大量的探索。提出了諸如WENO-M[23]、WENO-Z[24]、TENO[25]、雜交方法[26-27]等各種方法。

除了上述經典方法外，其實還有其他一些類型格式。如：高精度保單調格式(monotonicty-preserving, MP)[28]。這種格式的基本思想為：首先采用任意類型高精度格式預估一個初始結果；其次通過構造局地的上界和下界，形成一個局地的區間；最終通過探測初始結果是否落入此區間內，決定是否對預估的初始結果進行進一步處理。這類格式繼承了TVD類型格式的保界思想，并進一步拓展至高精度計算[29-30]。

因為WENO格式不具備保界性，Balsara在2000年提出了MPWENO格式[31]。其基本思想是用WENO格式預估初始結果，然后用MP格式構造的局地上界和下界來對初始結果進行進一步限制，從而實現保界性。

在2016年，我們發現[32]，MP格式的局地上界和下界的設計方法會導致形成的區間過大，無法抑制一些數值振蕩。這些數值振蕩會持續傳播，最終導致MP格式的穩定性較差。上述發現也得到了進一步工作的證實[33]：“It was recently found that this limiting procedure is no longer accurate for long - term integrations. As indicated by He et al., the major problem arises mostly due to the use of considerable amounts of room for the MP constraint, resulting in either problem-dependent solutions or abundant room for the growth of spurious oscillations”。結合TVD有界理論，我們給出了一種改進版本的MP格式[32]，即MP-R格式。數值實驗表明改進的局地保界方案可以較為顯著地提高最終方案的健壯性。

由于保界/保正要求對于模擬多介質問題至關重要，我們自然希望計算法自動實現這個要求。因此，上述基于保界思想構造的MP格式是一類比較適用于可壓縮多介質流動問題的方法。但是，對于非線性問題而言，這種先驗性的保界/保正算法不具備普適性。因此，為了進一步提升算法的保界性和保正性，可以進一步采用最新發展的保界/保正處理方法，如Postivity-Preserving方法[34]、Multi-dimensional Optimal Order Detection(MOOD)方法[35]。

1.5 界面壓縮

擴散界面方法是將界面視為數值擴散的狹窄區域，比如氣體動力學中的數值模擬所捕捉的接觸間斷。因此造成的結果是——數值界面會一直彌散變寬。并且隨著計算時間的增加，這種彌散變寬現象會越來越嚴重，最終可能會導致影響流動的其他物理結構[36]。因此，如何避免或抑制界面的持續彌散變寬問題，一直是可壓縮多介質流動問題界面捕捉算法的一個研究重點。

一個自然的想法是發展高精度間斷捕捉數值格式。在同樣的網格上，高精度間斷捕捉格式可實現高階精度計算，對間斷的分辨能力也顯著提高。但是文獻[37]中發現，采用高精度間斷捕捉格式，雖然可以在一定程度上提高間斷/界面的分辨能力，但依舊無法抑制長時計算情況下的界面彌散變寬問題。此外，由于高精度間斷捕捉格式引入的數值耗散較小，因此結果具有更多的數值振蕩。為了抑制這些數值振蕩，雖然可以采用局部特征分解方法來獲得基本無振蕩結果，但如果物質由非常復雜的狀態方程(EOS)所描述，如何構造一個穩定且與物理模型兼容的局部特征分解方法仍需要進一步探索[37]。

總之，采用高精度格式依舊無法抑制長時計算情況下的界面彌散變寬問題。此外，還會帶來可能的數值振蕩[37]。

因此，即使采用了高精度間斷捕捉格式來進行計算，界面壓縮算法也是必須的[37]。目前，相關研究者提出了各種界面壓縮/反耗散算法[36-39,53]。這些方法可以大致分為四類：(1) 限制型下游格式方法[53]。這種方法的基本思想是在滿足TVD穩定性理論的基礎上盡可能使用下游格式。目前僅適用于Lagrange-Remap方案。(2) 反擴散方法[54]。這類方法的基本思想是直接求解反擴散方程從而實現界面壓縮。(3) 基于修改重構過程的代數型VOF方法，如基于雙曲正切函數發展的界面捕捉算法(THINC)[38-39]。(4) 界面重整化(又叫擴散-銳利化)方法[36]。這類方法的基本思想是添加恰當的擴散項和銳利化項，從而保證界面的數值解是收斂的(寬度不隨計算時間增加而變寬)。目前這幾類方法推廣應用到可壓縮多介質問題時，均會遇到一個共同且依舊沒有徹底解決的問題——當描述界面演化的方程(如體積分數方程)采用上述界面壓縮方法后，其他方程(質量-動量-能量方程)如何進行相容的壓縮。前三類方法為了實現相容壓縮，目前已有的方案僅有一階精度，我們的數值實驗證實過低的精度可能會抑制界面的大變形運動。第四類方法的相容壓縮方法已經從理論上得到解決(針對五方程模型[43])，但是目前這類方法添加的擴散項和銳利化項均來源于level set方法[55]或相場方法[56]，如何設計和離散這些項(穩定且保界)目前還沒有得到徹底解決[57]。

為了易于實現高精度計算且避免如何設計擴散項和銳利化項的問題，在文獻[37]中，針對五方程模型[40-41]，我們提出了一種基于界面壓縮方案(第三類和第四類方法的雜交方案)。

考慮如下兩種物質的五方程模型[40]：

(13)

其中，ρ、u、E、p分別為混合物密度、速度、總能和壓力，ρk、αk分為是第k種物質的相密度和體積分數。其中體積分數滿足α1+α2=1。上述方程可以寫成(以一維情況為例)：

其中V=(ρ1α1,ρ2α2,ρu,ρE)T。

采用標準的有限體積方案，上述方程可以離散為：

提出的界面壓縮方案主要包括兩部分[37]：

1) 界面壓縮物理量的選擇。針對五方程模型，采用了如下的原始變量進行(相密度、速度、壓力、體積分數)重構[37](可進一步投影至特征空間進行重構)，來獲得單元邊界的左右重構狀態：

(16)

對于這些物理量，可采用上述的各種高精度數值格式進行計算。具體內容詳見文獻[37]。

上述方案不具備界面壓縮效應。為了引入界面的壓縮效應(即采用THINC格式進行計算)，需要對上述的標準方案進行改進。對何種物理量采用界面壓縮是需要考慮的首要問題。在文獻[37]中提出了如下方案：只對體積分數采用THINC格式進行計算，從而獲得的單元邊界的左右重構狀態。

(17)

可以看到，在這里獲得的單元邊界的左右重構狀態時，除了體積分數采用了其他類型界面壓縮方法來計算，其他物理量均采用原高精度數值格式方案來進行計算。

提出上述方案的基本思想主要基于兩點考慮：① 我們知道，物質界面/接觸間斷是一類非常特殊的流動結構。各種物理量過界面的表現性質不一樣。我們可以把這些物理量分為快變量和慢變量。其中慢變量定義為過界面連續的物理量。與之相反，快變量定義為過界面存在間斷的變化劇烈的物理量。從這個分類可以看出，上述五方程模型中，相密度、速度、壓力均是慢變量，而只有體積分數是快變量。② 僅對快變量做界面壓縮方法。這個觀點是基于如下事實：我們發現，界面壓縮方法，尤其是THINC類型的具有反耗散性質的方法，存在一個問題：這類方法可以使光滑的物理結構被反物理地壓縮成鋸齒狀間斷結構[42]。如圖4所示，當采用THINC格式對Burgers方程進行計算時，反耗散方法對間斷結構的分辨率非常高。但如果在光滑區也采用THINC方法，光滑的物理結構亦被非物理地壓縮為間斷結構。 關于這個問題的討論，詳見最新的文獻[42]。基于這樣的實際情況，可以得出如下結論：界面壓縮只能應用于快變量。正是基于上述觀點，我們提出了僅對體積分數進行壓縮處理的界面壓縮方案。

圖4 采用THINC格式計算Burgers方程[42]

2) 相容性。大量數值實驗表明：上述方案在工程計算中可以取得較好的計算結果，這種方案也可以應用于其他模型。我們所需要做的僅僅是：① 選擇重構的原始變量；② 對這些原始變量進行區分：快變量還是慢變量；③ 對快變量進行界面壓縮處理。表1給出了四方程模型和五方程模型的一些典型快變量和慢變量。

(19)

有了上述的體積分數方程的源項，就可以直接利用Tiwari的理論結果：在質量方程、動量方程、能量方程上加入能嚴格保證相容的源項。具體表達式詳見文獻[37]。圖5給出了激波-氣泡干擾問題(具體設置參見文獻[37])的模擬結果(有/無界面壓縮情況)。

從圖5中可以看到，采用了界面壓縮方法后，界面的分辨率得到顯著提高。此外，數值結果亦表明界面的寬度基本不隨時間增加而變寬[37]。

圖5 Air-R22激波-氣泡干擾問題結果[37]

但是，上述方法也存在一定的缺陷——最終的方案不是嚴格守恒[43]。此外，文獻[37]中指出，在實際應用中，不同近似黎曼求解器帶來的數值耗散不一樣。而由這部分因素帶來的數值耗散不能僅通過修改重構方法而徹底回避。因此，如何構造守恒且相容的界面壓縮方法仍值得進一步深入研究。

1.6 界面-單介質分區計算方案

開展可壓縮多介質流動問題數值模擬時，由于我們采用的模型為擴散界面模型。其本質假設是每一空間點上同時存在多種介質，介質之間滿足一定的關系(如壓力平衡-溫度非平衡假設)。因此為了避免模型的退化，必須在純流體中添加一個小的體積分數[37]。

在文獻[44-45]中我們發現，這種人為處理方法會給計算帶來非物理的波。這種現象在考慮彈塑性-流體相互作用的問題時尤為明顯：① 在流體中添加小部分固體，會讓流體產生抵抗剪切的能力；② 固體和流體的聲速差異很大，即使是一小部分的固體，也會影響純流體的聲速。在此發現基礎上，我們提出了一個改進的壓力平衡-溫度非平衡擴散界面模型，該模型在等壓縮假設的基礎上，導出了和體積分數無關的相密度方程。采用該模型，即使體積分數為0，也不會導致方程退化，避免了非物理波的產生。圖6給出了新模型和舊模型模擬純剪切問題的對比，新模型可以保持住剪切間斷，不產生額外的剪應力。而采用原模型，必須在方程中添加小的體積分數，導致虛假剪切應力的產生，以及剪切波向兩側擴散。但上述方案是否能推廣到實際工程應用中，還需進一步探索。

圖6 純剪切流數值模擬結果[44-45]

另一種解決方案是采用區域分解這類典型的多尺度問題計算方案。即將計算區域分解為純物質區和界面區，不同區域采用不同的方案進行計算。其實，在上一節構造界面壓縮算法過程中，已經人為地區分了界面區域和純物質區域。相關數值實驗亦證實了此方案的可行性。但界面區域與純物質區域的邊界處是否會人為引入一些非物理波，這一問題還需要更深入的研究。

2 應 用

通過上述多個維度的工作，建立了一套適用于可壓縮多介質流動問題的低耗散、基本無數值振蕩的高精度歐拉方法。相關數值方法成果已集成至工程問題數值模擬軟件，為相關工程任務提供了重要技術支撐。在具體應用方面，不同的問題往往具有自身的特殊性，需要有針對性地選擇不同的具體離散方案進行應用。本小節給出幾個具體應用算例。

第一個算例是彈塑性介質的高速沖擊問題。金屬等固體材料在受到高速沖擊后，會發生彈性到塑性的轉變，甚至產生斷裂、破碎等復雜的物理現象，對數值模擬提出了極大的挑戰。我們基于上述五方程模型，采用有限差分方法、高精度(MP)WENO格式、非守恒離散方法以及結合界面壓縮技術，并進一步對彈塑性效應進行了建模和數值離散，成功地實現了對此類問題的模擬。圖7給出了如下問題的數值模擬結果：銅塊以800 m/s的初始速度沖擊一個靜止的銅靶，銅靶先后經歷彈性變形、塑性變形，最終發生斷裂和射流現象。相關結果和已有的文獻結果一致。

(a) t=0 ms (b) t=1.1 ms

第二個算例是含顆粒Richtmyer-Meshkov(RM)不穩定性問題[9]。激波驅動的多相流動(如氣固多相流)中經常伴隨著多相流RM不穩定性的的發生。在氣固兩相流動體系中，顆粒相在流動體系中顯著影響了界面的動力學行為。我們基于上述五方程模型，采用有限體積方法、高精度(MP)WENO格式、非守恒離散方法，并進一步對跨流態顆粒運動進行了數值建模(CMP-PIC 方法[9])，成功實現了對此類問題的數值模擬。圖8具體給出了模擬的含空氣/SF6的多相RM問題(其中顆粒均勻分布于空間)問題的具體情況：空氣以及SF6界面存在單模擾動，入射激波由空氣進入到SF6區域。顆粒均勻布置于計算域，其分布區域為：0.01 m

圖8 計算域示意圖

(a) t=0.0301 ms，界面演化過程

圖10 不同顆粒體積分數下混合區寬度隨時間演化

表2 不同體積分數下混合區寬度數值及理論預測比較

第三個算例為含激波二次沖擊的RM湍流混合問題[46]。含激波沖擊的湍流混合問題是工程應用及自然現象中的基本物理過程，其典型特點是波系與混合區的多次作用。針對該問題，開展了直接數值模擬研究，采用了四方程模型，并且由于考慮物理黏性，所以沒有使用界面壓縮技術。并且采用的網格規模達到11億。初始激波從輕流體向重流體傳播，透射激波在下游壁面反彈后二次作用于混合區，混合區快速湍流化(如圖11所示)，隨后混合區依次被稀疏波和壓縮波加速。稀疏波和壓縮波作用期間，混合區內形成了明顯的平均速度梯度，混合區被拉伸或者壓縮。采用混合寬度無量綱化的平均組分剖面在二次激波作用后基本保持不變。借助高精度直接數值模擬結果，對多介質湍流的生成、衰減過程進行了收支平衡分析。結果表明，在稀疏波和壓縮波作用階段，速度梯度生成項和壓力梯度生成項是主導項。壓力梯度生成項和速度梯度生成項的貢獻相反，前者的幅值更大。因此，稀疏波導致湍動能增強而壓縮波導致湍動能減弱。進一步統計了含能結構的尺度變化。結果表明，混合寬度的增長伴隨著大尺度結構的形成，稀疏波階段，結構被縱向拉伸，而壓縮波作用階段，結構被縱向壓縮。

圖11 三維質量組分等值面的時間演化[46]. (a) t=0.0005 s和(b) t=0.0002 s為二次激波作用前；(c) t=0.003 s 和 (d) t=0.0045 s 為二次激波作用后

最后一個算例是匯聚幾何不穩定及湍流混合問題。在實際工程應用中，RM湍流混合往往發生在不斷向中心匯聚(如柱面、球面等)的物質界面上。相比于平面界面上的RM湍流混合，匯聚界面上的RM湍流混合區在演化過程中會伴隨有其周向長度的改變，而這種改變將會深刻影響混合區的后續發展。針對該問題，基于集成上述方案所開發的模擬程序(具體方案同上一個算例)，對圓柱界面上的RM不穩定性開展了數值模擬研究。初始激波由外部的重流體向內部的輕流體傳播，并在隨后往復運動于幾何中心和混合區之間，使得激波多次作用于混合區，混合區快速湍流化(如圖12所示)。對于混合區寬度，其演化被分解為混合區整體的拉伸壓縮和兩種流體組分相互穿插兩部分的貢獻。通過與平面算例進行對比，我們發現向內匯聚的物質界面將整體拉伸混合區，并在同時促進兩種流體的相互穿插。在激波第三次沖擊混合區之后，輕重流體穿插效應所貢獻的混合寬度W滿足W∝(t-t0)θ，這與平面工況中混合寬度的自相似增長相類似。對于混合區內部的平均組分剖面，經過混合區寬度歸一化的曲線近似重合。通過統計混合區前沿位置的分布，本文計算了氣泡和尖釘結構各自的直徑演化。數值結果表明，氣泡和尖釘的直徑與其對應的混合區高度之比在演化后期趨近于常數。

圖12 三維質量組分等值面的時間演化. (a) t=0.006 s為二次激波沖擊后，三次激波沖擊前和(b) t=0.012 s為三次激波作用后；(c) t=0.020 s為激波強度衰減至可以忽略后

從以上的算例可以發現，實際工程問題往往都是耦合多種物理因素和過程的多尺度多物理復雜流動問題。因此為了更好地模擬實際工程問題，我們需要進一步考慮黏性[58-59]、熱傳導[60-61]等耦合物理效應(更多物理效應見Saurel的系列工作及其綜述[62])，以及彈塑性[63-64]、輻射擴散[65]等多物理耦合過程。即需進一步研究上述可壓縮多介質數值方法在這類多物理效應/過程情況下的適用性問題。但是到目前為止，這一方面的研究還遠不夠成熟，我們的研究工作也正在進一步開展。

3 結 論

建立一套高效的可壓縮多介質流動問題的高精度數值模擬方案是一個系統工程，其中涉及的因素較多，主要有：數值框架選擇、非守恒方程相容離散、高精度有界格式構造、界面壓縮、界面-單介質分區計算方法等。本文從上述維度出發，綜述了本研究方向以及課題組的一些工作。主要結論有：

1) 不同數值方案可能導致過界面出現非物理振蕩現象。因此，不管選用何種類型的數值方案，必需對此方案在可壓縮多介質流動問題中的適用性先進行分析。為此，我們針對有限差分方法在此類問題中的適用性進行了分析。

2) 非守恒類型方程必需現相容離散。當守恒方程采用確定的雙曲守恒離散方法(如有限差分方法、有限體積方法)后，非守恒方程如何進行相容離散。這類方法如果離散不恰當，同樣會導致過界面的非物理振蕩。針對此問題，我們提出了一套相對普適的離散方案。

3) 保界/保正是構造適用于可壓縮多介質流動問題的高精度格式的重要要求。因為在對可壓縮多介質流動問題進行模擬時，不僅需要考慮熱力學量的保正問題，還得同時考慮如質量分數這類物理量的保界問題。為此，我們改進了經典的保單調限制器，進一步提升了此算法的保界性能。

4) 在擴散界面模型框架下，界面會一直數值彌散變寬，并且隨著計算時間的增加，這種彌散變寬現象會越來越嚴重。為了避免或抑制界面的持續彌散變寬現象，我們提出了一種界面壓縮方法。

5) 界面-單介質分區計算方案。為了避免模型的退化，必須在純流體中添加一個小的體積分數。我們發現，文獻上常用的這種方法會帶來非物理的波。因此，我們建議采用區域分解的這種多尺度計算方法，開展界面-單介質的分區計算。但如何保證邊界處不產生非物理振蕩需要進一步研究。

通過上述多個維度的工作，我們建立了一套適用于可壓縮多介質流動的低耗散、基本無數值振蕩的高精度歐拉方法，相關數值方法成果已集成至工程問題數值模擬軟件，為相關工程任務提供了重要技術支撐。

由于實際工程問題往往都是耦合其他效應和過程的多尺度多物理問題。因此，為了進一步提高模擬的置信度，需要進一步研究的工作有：

1)多物理效應和多物理過程的建模與計算方案。即需進一步研究“在上述多尺度多物理問題情況下的多介質數值方法的適用性”問題。比如當存在物理黏性時，界面壓縮能否與物理黏性自洽識別、可溶與不可溶界面的能否自動區分等問題。 典型的進展如最近Xiao等最近提出了Boundary Variation Diminishing (BVD)算法[66]，在此方面做了較為成功的嘗試。

2)極端條件下的數值算法研究。相比單介質問題，多介質問題往往更容易遇到各種極端條件(如復雜狀態方程、非牛頓流體、固體的彈塑性效應與超低馬赫數效應、物性差異巨大等)，目前的計算流體數值方法能否適用于這些極端情況、以及如何優化改進現有算法值得進一步研究。

3)高效的大規模計算軟件研發。可壓縮多介質問題中流場每一空間點可能為單介質或多介質，如何構造一種高效表征多介質流動問題的數據結構，以及大規模并行計算情況的下的如何荷載平衡等問題，均對軟件研發提出了挑戰。此外，多介質問題的各類模型并沒有經過大量大規模并行計算的測試與檢驗，進一步對研發高效的大規模計算軟件提出了需求。