程亞娟，曹 蕊，李祥貴

(北京信息科技大學(xué) 理學(xué)院，北京 100192)

0 引言

在武器物理、慣性約束聚變、天體物理等國(guó)防和前沿科學(xué)工程領(lǐng)域的物理問(wèn)題研究中，要面對(duì)復(fù)雜流體的數(shù)值模擬。由于復(fù)雜流體的控制方程具有高度非線性、強(qiáng)耦合等特點(diǎn)，其物理量在計(jì)算區(qū)域內(nèi)出現(xiàn)突變，高溫、高壓、高速導(dǎo)致流體界面形變較大，產(chǎn)生界面滑移、強(qiáng)激波、流體不穩(wěn)定性等問(wèn)題。要得到這些復(fù)雜流體控制方程的解析解是不可能的，傳統(tǒng)的數(shù)值方法處理這些問(wèn)題的數(shù)值模擬往往力不從心，因此發(fā)展高效數(shù)值方法相當(dāng)重要。

求解復(fù)雜流體控制方程首要的是將方程離散化，即用合適的網(wǎng)格模型劃分求解區(qū)域。根據(jù)方程的離散方式，數(shù)值方法可分為有限差分方法，有限體積方法和有限元方法。其中間斷Galerkin (DG)方法最早是1973年由Reed和Hill提出的。20世紀(jì)80年代后期和90年代初，Cockburn和Shu等結(jié)合龍格—庫(kù)塔(Runge-Kutta)方法[4]，將間斷有限元方法推廣到非線性一維守恒率方程和方程組、高維守恒率方程和方程組。無(wú)論是在歐拉網(wǎng)格[5]還是拉格朗日網(wǎng)格[6]計(jì)算中，DG方法都可以精確地捕捉激波等間斷現(xiàn)象，且易于并行計(jì)算，在光滑解區(qū)域達(dá)到高階精度。

根據(jù)流體與網(wǎng)格的相對(duì)運(yùn)動(dòng)方式，數(shù)值模擬方法可以分為歐拉方法、拉格朗日方法和網(wǎng)格移動(dòng)技術(shù)。歐拉方法的網(wǎng)格固定不動(dòng)，比較容易使用高精度格式[7]，并且在求解流體大變形問(wèn)題的時(shí)候不受網(wǎng)格變形影響[8]，但是對(duì)物質(zhì)界面的計(jì)算誤差較大。拉格朗日方法的網(wǎng)格隨流體移動(dòng)，可以精確分辨物質(zhì)界面，但網(wǎng)格嚴(yán)重變形扭曲時(shí)會(huì)導(dǎo)致計(jì)算精度不高甚至計(jì)算失敗。移動(dòng)網(wǎng)格方法根據(jù)物理解的變化，隨時(shí)調(diào)整網(wǎng)格的疏密和形狀，可以精確分辨物質(zhì)界面性，不易出現(xiàn)計(jì)算中斷的問(wèn)題。網(wǎng)格移動(dòng)技術(shù)結(jié)合拉格朗日方法和歐拉方法的優(yōu)點(diǎn)，可以較好地解決兩者在求解復(fù)雜流體問(wèn)題中的困難[9]，可以讓物質(zhì)界面附近的網(wǎng)格根據(jù)計(jì)算需要選擇移動(dòng)方式。這樣既可以精確分辨物質(zhì)界面，避免不必要的混合網(wǎng)格單元產(chǎn)生，又可以利用DG方法易于保持守恒性的特點(diǎn)，更有利于求解復(fù)雜流體問(wèn)題。

本文主要基于移動(dòng)網(wǎng)格在空間和時(shí)間離散上分別采用DG和龍格—庫(kù)塔方法求解復(fù)雜流體控制方程。第1節(jié)介紹復(fù)雜流體的控制方程。第2節(jié)給出了求解方程的移動(dòng)網(wǎng)格龍格—庫(kù)塔DG方法，包括移動(dòng)網(wǎng)格的物理量及節(jié)點(diǎn)、單元坐標(biāo)的更新、DG方法的空間離散和龍格—庫(kù)塔時(shí)間離散方法3部分內(nèi)容，并討論了詳細(xì)的求解算法。第3節(jié)通過(guò)數(shù)值實(shí)驗(yàn)去驗(yàn)證提出的移動(dòng)網(wǎng)格DG方法的有效性，并對(duì)比速度平均、密度、能量的加權(quán)Roe平均處理單元節(jié)點(diǎn)速度的效果。

1 控制方程

流體流動(dòng)過(guò)程都要受質(zhì)量、動(dòng)量、能量守恒等基本物理定律的支配，即要滿足一定的守恒控制方程組[10]。控制方程組可以寫成如下形式：

(1)

本文研究一維多相反應(yīng)流方程。多相流是復(fù)雜流體中的一種，即有不同的流體成分，它們以其各自的比熱為特征，互不相容。此外，忽略了擴(kuò)散效應(yīng)、表面張力和氣穴現(xiàn)象。一般地，通過(guò)增加附加方程來(lái)完成多相反應(yīng)流的描述，使用參數(shù)μ來(lái)描述反應(yīng)的強(qiáng)度。根據(jù)模型假設(shè)，控制方程可以用守恒的形式寫成

(2)

式中：對(duì)于不同的反應(yīng)，g有不同的表達(dá)形式；λ為質(zhì)量分?jǐn)?shù)；μ為調(diào)控反應(yīng)的常參數(shù)，μ的取值大小由反應(yīng)確定，反應(yīng)強(qiáng)度越大，非線性反應(yīng)強(qiáng)度越大，μ的取值越大。方程式(2)僅限于考慮理想氣體和低溫模型。理想流體力學(xué)的狀態(tài)方程為

p=(γ-1)ρe

(3)

式中γ為流體混合物的比熱。

我們將使用移動(dòng)網(wǎng)格DG方法求解式(1)和式(2)，式(1)和式(2)可以簡(jiǎn)化為

Wt+F(W)x=K

(4)

式中

2 移動(dòng)網(wǎng)格上的高階DG方法

移動(dòng)網(wǎng)格的高階DG方法主要包括3個(gè)方面:移動(dòng)網(wǎng)格單元位置及節(jié)點(diǎn)坐標(biāo)的更新、高階DG的空間離散和同階的TVB龍格—庫(kù)塔時(shí)間離散。

2.1 網(wǎng)格的更新

設(shè)初始時(shí)刻t=0的計(jì)算區(qū)域ΩC是固定的，ΩC被剖分為J個(gè)單元，即ΩC=∑Ii，Ii=(xi,xi+1)(i=1,2,…,J)。利用網(wǎng)格的運(yùn)動(dòng)速度計(jì)算在某一時(shí)刻t=a(a>0)網(wǎng)格移動(dòng)后新的計(jì)算域及相應(yīng)的剖分單元Ii=(xi,xi+1)，此時(shí)網(wǎng)格不再是均勻的。

下面給出式(4)的弱解形式。在區(qū)間Ii=(xi,xi+1)上用測(cè)試函數(shù)ζ(x)∈L2(ΩC)乘以式(4)得到

(5)

式(5)稱為式(4)的弱形式。再對(duì)式(5)使用分部積分有

(6)

(7)

將式(7)代入式(6)得到移動(dòng)網(wǎng)格下一維多相反應(yīng)流方程的弱形式：

(8)

式中

每個(gè)單元的節(jié)點(diǎn)速度ui分別按照速度平均以及節(jié)點(diǎn)左右兩側(cè)的密度、能量加權(quán)來(lái)取Roe平均[12]：

本文的移動(dòng)網(wǎng)格是指在計(jì)算過(guò)程中，讓物質(zhì)界面xs附近的網(wǎng)格和其他區(qū)域的網(wǎng)格按照各自定義的速度移動(dòng)，即

式中θ是一個(gè)正參數(shù)，可以根據(jù)計(jì)算進(jìn)行調(diào)控。這樣既可以精確分辨物質(zhì)界面，又可以避免不必要的混合網(wǎng)格單元產(chǎn)生。

由于網(wǎng)格是不斷運(yùn)動(dòng)的，那么每個(gè)時(shí)刻的單元位置、尺度、節(jié)點(diǎn)坐標(biāo)和物理量，以及一系列與節(jié)點(diǎn)有關(guān)的量都需要進(jìn)行更新。

2.2 DG方法對(duì)空間離散

Wh=(w1,w2,w3,w4)Τ=(ρh,(ρu)h,Eh,(ρλ)h)Τ

使得下述方程成立：

(9)

利用高斯積分公式計(jì)算線性積分

可以采用：

1)Lax-Friedrichs(L-F)數(shù)值通量

這里βi取第i個(gè)節(jié)點(diǎn)左右兩邊的最大聲速。

2)HLLC數(shù)值通量，詳見文獻(xiàn)[13]。

A1i(x)=1

其中xi+1/2=(xi+1+xi)/2。

本文以三階精度DG方法為例給出，對(duì)于高階精度格式可以相似處理。設(shè)在單元i上的自由度為Di(x)，T表示轉(zhuǎn)置。由自由度和基函數(shù)可以確定Wh。令

那么

(10)

x∈Iij=1,2,3；l=1,2,3,4

(Di)tMi=Z(Di)

式中：Z(Di)為右邊的空間離散算子；Mi被稱作質(zhì)量矩陣，表示為

本文的Mi是一個(gè)特征值為1、1/3、4/45的(3×3)對(duì)角矩陣。

2.3 TVB龍格—庫(kù)塔時(shí)間離散

三階TVB龍格—庫(kù)塔方法對(duì)時(shí)間離散的詳細(xì)步驟如下。

1)更新基函數(shù)。由于網(wǎng)格是不斷移動(dòng)的，所以在每個(gè)時(shí)間步上都需要先更新單元尺度、基函數(shù)和節(jié)點(diǎn)坐標(biāo)等。先更新單元尺度和節(jié)點(diǎn)位置，Δt是固定的時(shí)間步長(zhǎng)，它的取值由下式?jīng)Q定:

(11)

最后得到更新的基函數(shù)：

(12)

式中

(j,r=1,2,3；m=1,2,n+1)

為了保證數(shù)值解的穩(wěn)定性和抑制解的虛假震蕩，在每個(gè)時(shí)間步長(zhǎng)都引入了一個(gè)限制器，此處引入TVB限制器[15]。對(duì)于標(biāo)量方程，TVB限制器是利用單元平均的差值作用在Dilj(j=1,2,3)上，保證總網(wǎng)格函數(shù)值變差減小的性質(zhì)。令

h=max(xi+1-xi)

使得

不論針對(duì)方程組還是標(biāo)量方程，j=1時(shí)限制器均不起作用。以方程組為例，需要在每個(gè)時(shí)間步上對(duì)特征變量作出限制，例如對(duì)Dil2做出限制：

式中

Δ+Di2=Di+1,l1-Di,l1

Δ-Dil2=Di,l1-Di-1,l1

m是minmod函數(shù)，定義為

m(a1,a2，…,an)=

上述格式的精度在極值點(diǎn)附近退化為一階。為了克服這一缺陷，采用修正的minmod函數(shù)

Wh=Dil1A1i+Dil2A2i

TVB龍格—庫(kù)塔時(shí)間離散，不但保證在每一時(shí)間離散中間步是TVB的，而且有效地抑制了震蕩。

2.4 算法設(shè)計(jì)

求解過(guò)程基于兩個(gè)獨(dú)立的部分:網(wǎng)格更新和算法求解。第一部分是基于式(7)和式(11)的更新過(guò)程。第二部分在第一部分的基礎(chǔ)上進(jìn)行。求解過(guò)程如下。

3 數(shù)值算例

激波管問(wèn)題算例如下。

給出左右常數(shù)狀態(tài)的初值條件：

此問(wèn)題又稱為L(zhǎng)ax的Riemann問(wèn)題。在[0，10]進(jìn)行計(jì)算，取均勻的200個(gè)網(wǎng)格為初始計(jì)算網(wǎng)格，γ=1.4。給出t=1.0 s時(shí)的數(shù)值結(jié)果。圖1和圖2是分別采用HLLC數(shù)值通量和L-F數(shù)值通量進(jìn)行計(jì)算的結(jié)果，兩種數(shù)值通量在高階格式下均能夠近似精確解。由于L-F數(shù)值通量的使用需要人工粘性的引入，所以相比HLLC數(shù)值通量計(jì)算誤差更大。兩者在高階精度格式下都接近精確解，區(qū)別不大，所以使用高階精度的格式計(jì)算時(shí)，兩種通量幾乎能達(dá)到一樣的效果，而相較于HLLC數(shù)值通量，L-F數(shù)值通量程序簡(jiǎn)單，計(jì)算量小。

圖1 采用HLLC數(shù)值通量的一階格式和二階格式

圖2 L-F數(shù)值通量的一階格式和二階格式

另外，使用HLLC數(shù)值通量，在處理節(jié)點(diǎn)速度時(shí)發(fā)現(xiàn)采用3種加權(quán)Roe平均的圖像幾乎達(dá)到了重合，由此可知，對(duì)于一階和二階格式的激波管問(wèn)題，不論取何種Roe平均差別都不大。圖2為3種加權(quán)Roe對(duì)比。

圖3 速度平均、密度、能量加權(quán)的Roe平均

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

本文給出了一種高階的求解一維多相反應(yīng)流方程的間斷有限元方法，并通過(guò)數(shù)值實(shí)驗(yàn)比較L-F、HLLC兩種數(shù)值通量在精度方面的誤差，發(fā)現(xiàn)使用一階格式時(shí)HLLC的精度更好，但對(duì)于高階精度格式兩者之間的差別不大，均能很好地近似精確解。另外通過(guò)比較不同權(quán)重的Roe平均處理速度的效果可知，權(quán)重對(duì)誤差影響甚小?；谝苿?dòng)網(wǎng)格技術(shù)的DG方法求解一維多相流反應(yīng)方程能夠達(dá)到不錯(cuò)的效果，未來(lái)我們將進(jìn)行二維多相反應(yīng)流方程的求解。