劉學(xué)哲,林 忠,王瑞利,余云龍

(北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所,北京 100094)

輻射流體力學(xué)在武器物理、慣性約束聚變、磁約束聚變、地質(zhì)學(xué)和天體物理等諸多領(lǐng)域有廣泛的應(yīng)用，描述多物理過程耦合的輻射流體力學(xué)方程組具有強(qiáng)耦合、強(qiáng)間斷、強(qiáng)非線性等特征，很難得到解析解，利用計(jì)算機(jī)進(jìn)行數(shù)值求解是重要的研究手段。

數(shù)值模擬是否正確求解了方程，模擬結(jié)果能否再現(xiàn)真實(shí)的物理過程，都需要對數(shù)值模擬程序進(jìn)行驗(yàn)證與確認(rèn)[1-3]。人為構(gòu)造解方法[4]是程序驗(yàn)證的重要方法之一，對某些特殊需求(正確性驗(yàn)證、收斂性考察，等)具有不可替代性。在求解雙曲型流體力學(xué)方程組程序人為解驗(yàn)證方面已有很好的工作[5-14]，但大多基于歐氏方程，對拋物型輻射擴(kuò)散方程解析解及人為解驗(yàn)證方面也有很多工作[15-16]，而適用于求解輻射與流體耦合方程組及拉氏程序驗(yàn)證的人為解模型尚不多見。

針對拉氏輻射流體力學(xué)程序正確性驗(yàn)證的需要，文獻(xiàn)[17-18]中構(gòu)造了一類一維拉氏流體力學(xué)和輻射流體力學(xué)人為解模型，但一維模型的網(wǎng)格運(yùn)動僅改變網(wǎng)格點(diǎn)的疏密，不涉及到網(wǎng)格的變形，而網(wǎng)格隨流體運(yùn)動而變形是拉氏計(jì)算的特點(diǎn)，因此，構(gòu)造適應(yīng)流體大變形的二維拉氏輻射流體力學(xué)人為解模型對實(shí)際應(yīng)用程序的正確性驗(yàn)證很有意義。本文中基于二維坐標(biāo)變換關(guān)系式，研究了拉氏輻射流體力學(xué)人為解方法，構(gòu)造了適用于輻射流體力學(xué)程序驗(yàn)證的二維人為解模型，并應(yīng)用于非結(jié)構(gòu)拉氏應(yīng)用程序LAD2D[19]輻射流體力學(xué)計(jì)算的正確性考核，取得了很好的效果。

1 二維拉氏輻射流體力學(xué)方程組

考慮如下拉氏形式的二維輻射流體力學(xué)方程組：

(1)

(2)

(3)

(4)

p=pρ,T,e=eρ,T,κ=κρ,T

(5)

和適當(dāng)?shù)某踹呏禇l件，則計(jì)算模型封閉。

2 二維拉氏輻射流體力學(xué)人為解構(gòu)造

2.1 二維坐標(biāo)變換

引入拉氏坐標(biāo)x0=x0,y0,給定拉氏空間x0,y0,τ到歐氏空間x,y,t的坐標(biāo)變換關(guān)系式：

x=xx0,y0,τ,y=yx0,y0,τ,t=τ

(6)

其微分關(guān)系為：

(7)

記Jx0,y0,τ為坐標(biāo)變換的Jacobi矩陣：

(8)

其行列式記為：

(9)

(10)

為流體速度。

函數(shù)ux0,y0,τ、vx0,y0,τ與坐標(biāo)變換的Jacobi矩陣Jx0,y0,τ應(yīng)滿足Cauchy相容關(guān)系：

(11)

下面給出歐氏空間中的空間導(dǎo)數(shù)和拉氏隨體導(dǎo)數(shù)在拉氏空間的表達(dá)式。

設(shè)物理變量z在歐氏空間關(guān)于坐標(biāo)x,y,t的函數(shù)關(guān)系為f，在拉氏空間關(guān)于坐標(biāo)x0,y0,τ的函數(shù)關(guān)系為g，即：

z=fx,y,t=fxx0,y0,τ,yx0,y0,τ,τ≡

gx0,y0,τ=gx0x,y,t,y0x,y,t,t≡fx,y,t

(12)

則有：

(13)

(14)

(15)

而:

即歐氏空間的拉氏隨體導(dǎo)數(shù)就等于拉氏空間的時(shí)間導(dǎo)數(shù)：

(16)

2.2 二維拉氏輻射流體力學(xué)人為解構(gòu)造方法

從拉氏計(jì)算中使用的二維拉氏輻射流體力學(xué)方程(1)～(5)出發(fā)，給定位移解函數(shù)和溫度解函數(shù)，利用二維坐標(biāo)變換關(guān)系式，由無源項(xiàng)的質(zhì)量方程解出密度解函數(shù)，再將已知解函數(shù)代入動量方程和能量方程，殘余項(xiàng)視為方程源項(xiàng)。具體構(gòu)造過程如下:

(1)給定位移解函數(shù)x=xx0,y0,τ，y=yx0,y0,τ和溫度解函數(shù)T=Tx0,y0,τ。

(2)根據(jù)式(9)和(10)，由位移解函數(shù)x=xx0,y0,τ,y=yx0,y0,τ,可計(jì)算出坐標(biāo)變換Jacobi行列式J=Jx0,y0,τ和流體速度u=ux0,y0,τ,v=vx0,y0,τ。

(3)解連續(xù)性方程求出密度解函數(shù)。

利用空間導(dǎo)數(shù)關(guān)系式(15)和Cauchy相容關(guān)系式(11)求出速度散度：

(17)

代入連續(xù)性方程(1)，并利用歐氏空間的拉氏隨體導(dǎo)數(shù)與拉氏空間的時(shí)間導(dǎo)數(shù)之間的關(guān)系式(16)，得：

(18)

式中：ρ0x0,y0=ρx0,y0,0，J0x0,y0=Jx0,y0,0。

(4)根據(jù)溫度解函數(shù)確定動量方程源項(xiàng)和內(nèi)能方程源項(xiàng)。

由動量方程(2)～(3)導(dǎo)出源項(xiàng)：

(19)

根據(jù)時(shí)間導(dǎo)數(shù)關(guān)系式(16)和空間導(dǎo)數(shù)關(guān)系式(15)，可得：

(20)

(21)

由(18)式可得：

(22)

(23)

其中用到了狀態(tài)方程關(guān)系式(5)：

px0,y0,τ=pρx0,y0,τ,Tx0,y0,τ

由內(nèi)能方程(4)，導(dǎo)出源項(xiàng)：

(24)

利用時(shí)間導(dǎo)數(shù)關(guān)系式(16)、速度散度表達(dá)式(17)和空間導(dǎo)數(shù)關(guān)系式(15)，得：

(25)

其中用到了狀態(tài)方程關(guān)系式(5)：

ex0,y0,τ=eρx0,y0,τ,Tx0,y0,τ,κx0,y0,τ=κρx0,y0,τ,Tx0,y0,τ

3 二維拉氏輻射流體力學(xué)程序驗(yàn)證

3.1 拉氏輻射流體二維人為解模型

為方便計(jì)，在不致引起混淆的情況下，時(shí)間變量統(tǒng)一用t表示。

由上節(jié)的人為解方法可知，選擇不同的位移解函數(shù)和溫度解函數(shù)，可以構(gòu)造出不同的人為解模型。

給定計(jì)算區(qū)域:

0≤x0,y0≤1； 0≤x,y≤1； 0≤t≤1

和狀態(tài)方程:

p=0,e=cVT,κ=k0cV=1，κ0=1

選擇位移解函數(shù)：

xx0,y0,t=x0+x01-x0bsin2πtcosπy0,

yx0,y0,t=y0+y01-y0bsin2πtcosπx0b=0.8

和溫度解函數(shù):

T=Tx,y,t=T1+sin2πxcos(2πy)T1=2

初始密度:

ρ0x0,y0=ρ0=1

則可推導(dǎo)出以下。

(1)坐標(biāo)函數(shù)初、邊值：

0≤x0≤1， 0≤y0≤1

xlefty0,t=0,xrighty0,t=1，ytopx0,t=1,ybottomx0,t=0 0≤t≤1

且:

xlefty0,t0=x0left=0，xrighty0,t0=x0right=1，ytopx0,t0=y0top=1,ybottomx0,t0=y0bottom=0

式中：xleft、xright、ytop、ybottom分別為計(jì)算區(qū)域的左、右、上、下邊界。

(2)坐標(biāo)變換Jacobi行列式：

1+(1-2x0)bsin2πtcosπy01+(1-2y0)bsin2πtcosπx0-

πx01-x0bsin2πtsinπy0πy01-y0bsin2πtsinπx0

不難驗(yàn)證，當(dāng)0≤x0≤1，0≤y0≤1，00;當(dāng)t=0時(shí)，J(x0,y0,t)=1。

(3)速度解函數(shù)：

ux0,y0,t=2πx01-x0bcos2πtcosπy0

vx0,y0,t=2πy01-y0bcos2πtcosπx0

(4)密度解函數(shù)：

(5)壓力解函數(shù)：

px0,y0,t=0

(6)比內(nèi)能解函數(shù)：

ex0,y0,t=T

(7)內(nèi)能方程源項(xiàng)：

cos2πxcos2πyx01-x0cosπy0-sin2πxsin2πyy01-y0cosπx0+

8π2k0sin2πxcos2πyJ

(8)動量方程源項(xiàng)：

(9)邊界條件。

流體為固壁邊界條件：

ulefty0,t=0，urighty0,t=0，vtopx0,t=0,vbottomx0,t=0

輻射左右邊界為第一類邊界條件，上下邊界為第二類邊界條件：

不難看出，構(gòu)造的人為解滿足保正的物理?xiàng)l件，即對任意的:

x0,y0,t∈xleft,xright×ytop,ybottom×R+

都有:

Jx0,y0,t>0，ρx0,y0,t>0，Tx0,y0,t>0，ex0,y0,t>0

3.2 二維拉氏輻射流體力學(xué)LAD2D程序正確性驗(yàn)證

將構(gòu)造的人為解模型應(yīng)用于二維非結(jié)構(gòu)網(wǎng)格拉氏輻射流體力學(xué)應(yīng)用程序LAD2D[19]的正確性驗(yàn)證。

初始時(shí)，計(jì)算區(qū)域均勻劃分為32×32的方形網(wǎng)格，人為解中參數(shù)b=0.8，取固定時(shí)間步長Δt=7.812 5×10-4，計(jì)算到t=1.0時(shí)刻。

由構(gòu)造的人為解位移解函數(shù)可知，它為關(guān)于時(shí)間t的周期函數(shù)，這里選取了t=0，0.125，0.25等3個(gè)典型時(shí)刻，圖1給出了不同時(shí)刻流體運(yùn)動網(wǎng)格圖。

圖1 不同時(shí)刻流體運(yùn)動網(wǎng)格Fig.1 Fluid meshes at different times

圖2給出了t=0.125時(shí)刻流體網(wǎng)格上的密度、溫度、x方向速度、y方向速度的數(shù)值解和人為精確解的等值線，紅線為數(shù)值解，綠線為人為精確解，兩者基本一致，驗(yàn)證了程序輻射流體求解的正確性。

圖2 t=0.125時(shí)刻流體運(yùn)動網(wǎng)格上密度、溫度、x方向速度、y方向速度等值線(紅：數(shù)值解 綠：精確解)Fig.2 Contours of density, temperature, x and y components of velocity(red: numerical solution; green: exact solution)

通常的擴(kuò)散方程求解，無論是正交網(wǎng)格還是大變形扭曲網(wǎng)格，求解過程中網(wǎng)格始終保持不變，在二維運(yùn)動網(wǎng)格上考察擴(kuò)散問題目前尚未見文獻(xiàn)報(bào)道。由本文人為解構(gòu)造可知，通過設(shè)定特殊形式的狀態(tài)方程(p=0)，可以在流體運(yùn)動網(wǎng)格上求解擴(kuò)散模型。 圖3給出了不同時(shí)刻流體運(yùn)動網(wǎng)格上溫度數(shù)值解和人為精確解的比較，紅線為數(shù)值解，綠線為人為精確解，從圖中可以看出，從初始時(shí)刻開始，網(wǎng)格在不斷運(yùn)動，網(wǎng)格變形程度逐漸增大，而我們構(gòu)造的溫度解函數(shù)在歐氏空間來看不隨時(shí)間變化，數(shù)值解較好地保持了精確解的這一性質(zhì)。

圖3 不同時(shí)刻流體運(yùn)動網(wǎng)格上溫度等值線(紅：數(shù)值解 綠：精確解)Fig.3 Contours of temperature on fluid meshes at different time (red: numerical solution; green: exact solution)

進(jìn)一步考察人為解的數(shù)值收斂性，在不同的網(wǎng)格規(guī)模下，網(wǎng)格比Δt/Δx2=0.8保持不變，計(jì)算到t=1.0時(shí)刻，表1給出了溫度的誤差收斂性分析。數(shù)值結(jié)果顯示溫度L2模和最大模誤差二階收斂，驗(yàn)證了計(jì)算結(jié)果與理論分析的一致性。

表1 溫度收斂誤差和收斂階Table 1 Convergence errors and orders for temperature

4 結(jié) 論

從拉氏計(jì)算中使用的輻射流體力學(xué)方程組出發(fā)，基于二維坐標(biāo)變換技術(shù)，研究了一類質(zhì)量方程無源項(xiàng)，動量方程和能量方程包含源項(xiàng)的人為解構(gòu)造方法，構(gòu)造了一類適用于輻射流體力學(xué)程序驗(yàn)證的二維人為解模型，并應(yīng)用于非結(jié)構(gòu)拉氏程序LAD2D輻射流體力學(xué)計(jì)算的正確性考核，為流體運(yùn)動網(wǎng)格上的輻射擴(kuò)散計(jì)算提供了新的途徑。