劉曉川，王存海，黃 勇，朱克勇

(1.北京航空航天大學(xué)航空科學(xué)與工程學(xué)院，北京 100191，2.北京科技大學(xué)能源與環(huán)境工程學(xué)院，北京100083)

輻射傳輸方程描述了輻射能量在介質(zhì)中的傳遞，在許多科學(xué)和工程領(lǐng)域具有重要作用，例如大氣輻射傳輸[1]、光學(xué)層析成像[2]、天體物理學(xué)[3]及核工程[4]等.輻射傳輸方程是一個(gè)高維、復(fù)雜的積分微分方程，輻射強(qiáng)度涉及波長、時(shí)間、空間和角度等，求其解析解十分困難.學(xué)者們提出發(fā)展了很多種數(shù)值方法來求解輻射傳輸方程，如蒙特卡洛法[5]，離散坐標(biāo)法[6]，有限體積法[7]，有限元法[8]等.

近年來，利用格子-Boltzmann方法(LBM)來求解輻射傳輸方程吸引了許多學(xué)者的興趣.LBM起源于格子氣自動(dòng)機(jī)，已經(jīng)發(fā)展成為了一種計(jì)算流體力學(xué)的有力數(shù)值工具[9].并且，LBM已經(jīng)被拓展到求解許多線性和非線性系統(tǒng)問題，例如聲子輸運(yùn)[10]，波傳播[11]，反應(yīng)擴(kuò)散[12]，對(duì)流擴(kuò)散[13]等.相比于其他的求解輻射傳輸方程的數(shù)值方法，LBM不需要計(jì)算大量的光線軌跡，也不需要離散復(fù)雜的偏微分方程.LBM具有容易實(shí)現(xiàn)，高并行效率等優(yōu)點(diǎn).目前，對(duì)于利用LBM來求解輻射方程還不完善，發(fā)展完善的LBM用于求解輻射傳輸方程是必要的.

Mishra等[14]假定了可調(diào)節(jié)的虛擬光速和輻射平衡條件，將LBM推廣到分析參與性介質(zhì)中的輻射問題.Ma等[15]基于輻射流體力學(xué)，提出了一維輻射的格子-Boltzmann模型.Zhang等[16]通過采用全隱式后項(xiàng)差分格式處理輻射方程中的瞬態(tài)項(xiàng)，將LBM擴(kuò)展到求解參與性介質(zhì)中的一維瞬態(tài)輻射傳輸.Mink等[17]在將P1近似應(yīng)用輻射傳輸方程的基礎(chǔ)上提出了一種三維的格子-Boltzmann模型，然而此模型僅適用于光學(xué)厚介質(zhì).Yi等[18]通過引入虛擬的擴(kuò)散項(xiàng)，將輻射傳輸方程視為一種特殊的對(duì)流擴(kuò)散方程，從而提出了一種二維穩(wěn)態(tài)輻射傳輸方程的格子-Boltzmann模型.Wang等[19]將瞬態(tài)輻射傳輸方程處理為雙曲守恒方程，然后提出了一種求解瞬態(tài)輻射和中子輸運(yùn)的格子-Boltzmann模型.

目前，求解輻射方程的多松弛的格子-Boltzmann模型還未見報(bào)道.本文提出了一種多松弛格子-Boltzmann模型(multiple-relaxation-time lattice Boltzmann model).基于擴(kuò)散尺度下的Maxwell迭代，輻射傳輸方程可以嚴(yán)格地從格子Boltzmann方程推導(dǎo)得出，并且不引入任何限制和近似.本文發(fā)展的多松弛格子-Boltzmann模型可以精確地求解參與性介質(zhì)內(nèi)的多維瞬態(tài)及穩(wěn)態(tài)輻射傳輸問題.數(shù)值結(jié)果表明該模型具有二階精度和收斂速率.并且，相比于單松弛模型，多松弛模型具有更好的穩(wěn)定性.該模型可以進(jìn)一步推廣到求解參與性介質(zhì)內(nèi)的輻射傳輸問題.

1 輻射傳輸方程的多松弛格子-Boltzmann模型

1.1 輻射傳輸方程

考慮吸收、發(fā)射和散射介質(zhì)內(nèi)的輻射傳輸方程，其離散坐標(biāo)形式可以寫為[20]

(1)

公式中：cL為介質(zhì)內(nèi)的光速；I為輻射強(qiáng)度；r為位置坐標(biāo)；β=ka+ks為衰減系數(shù)；Ωm=μmi+ηmj+ξmk為離散方向，源項(xiàng)S可以表示為

(2)

公式中：N為總的離散方向，m=1，2，…，N，m′=1，2，…，N；wm′為對(duì)應(yīng)方向的權(quán)重.

考慮漫發(fā)射和反射壁面，邊界條件可以寫為

(3)

公式中：εw為發(fā)射率；ρw為反射率；Iext為外部入射輻射強(qiáng)度.

1.2 多松弛格子-Boltzmann模型

瞬態(tài)輻射常常發(fā)生于極短的時(shí)間內(nèi)，在瞬態(tài)輻射的模擬中，通常引入無量綱時(shí)間來避免過小的時(shí)間步長.將無量綱時(shí)間t*=cLt/LR代入方程(1)中，得到時(shí)間無量綱形式的輻射傳遞方程[21]

(4)

公式中

F(r，Ωm，t*)=LRS(r，Ωm，t*)-LRβI(r，Ωm，t*)

(5)

公式中：LR為介質(zhì)的參考長度.

本文提出的時(shí)間無量綱形式的輻射傳輸方程的格子Boltzmann方程如下

(6)

公式中：fi(r，t*)為分布函數(shù)；M為變換矩陣；S=diag(s0，s1，…，sn)為松弛參數(shù)矩陣，平衡函數(shù)的表達(dá)式為

(7)

輻射強(qiáng)度可以由平衡函數(shù)給出，關(guān)系如下

(8)

LBM方法中采用DmQn格子模型，對(duì)于一維和二維問題，本文分別采用D1Q3和D2Q9模型.對(duì)于D1Q3模型，其格子信息為

[c0，c1，c2]=eic=[0 1 -1]c

(9)

(10)

(11)

對(duì)于D2Q9模型，其格子信息為

(12)

(13)

(14)

1.3 從格子Boltzmann方程到輻射傳輸方程

本節(jié)基于擴(kuò)散尺度Δt*=γ(Δx)2下的Maxwell迭代，不引入任何限制和假設(shè)，從多松弛格子-Boltzmann模型嚴(yán)格推導(dǎo)得出輻射傳輸方程.這種擴(kuò)散尺度是針對(duì)模型中的無量綱時(shí)間步長和空間步長的尺度.

首先，令f(r，t*)=(f0(r，t*)，f1(r，t*)，…，f8(r，t*))T，ω=(ω0，ω1，…，ω8)T，時(shí)間無量綱形式的輻射傳遞方程(6)可以寫成矢量形式

f(r+ciΔt*，t*+Δt*)-f(r，t*)=-M-1SM[f(r，t*)-feq(r，t*)]+Δt*ωF(r，Ωm，t*)

(15)

方程(15)左邊應(yīng)用Taylor展開，

(16)

其中微分算子Ds

(17)

矩陣

Ex=diag(e0，x，e1，x，…，e8，x)Ey=diag(e0，y，e1，y，…，e8，y)

(18)

公式中：p和q均為非負(fù)整數(shù).

令m=M·f，meq=M·feq，將Taylor展開形式代入方程(15)并整理得到

(19)

其中

(20)

(21)

m=meq-S-1Lm+γ(Δx)2FS-1Mω

(22)

基于擴(kuò)散尺度下的Maxwell[22]迭代，從m0=meq開始，方程(19)經(jīng)過三次迭代得到：

(23)

根據(jù)矢量方程(23)的第零項(xiàng)及各算子作用結(jié)果，可以得到輻射傳遞方程

(24)

至此，我們從多松弛格子-Boltzmann模型出發(fā)，基于擴(kuò)散尺度下的Maxwell迭代，嚴(yán)格推導(dǎo)得出了輻射傳輸方程，并且可以從方程(24)理論上得出該模型具有二階的精度.一般而言，對(duì)于對(duì)流擴(kuò)散問題，計(jì)算流體力學(xué)等問題的LB模型，其中的松弛系數(shù)與宏觀方程中的擴(kuò)散系數(shù)，流體黏性系數(shù)等有定量關(guān)系.需要指出的是，根據(jù)從多松弛格子-Boltzmann模型嚴(yán)格推導(dǎo)得出輻射傳輸方程可知，本文提出的多松弛格子-Boltzmann模型中的松弛參數(shù)均是自由的，與其他參數(shù)無關(guān).對(duì)于一維和二維LB模型，我們?nèi)∪缦碌乃沙趨?shù)矩陣

S=diag(1，sr，1)

(25)

S=diag(1，1，1，sr，1，sr，1，1，1)

(26)

對(duì)流擴(kuò)散方程的多松弛LB模型也采用了同樣的處理方法，其中一維模型中的松弛參數(shù)s1，二維模型中的松弛參數(shù)s3和s5與擴(kuò)散系數(shù)有關(guān)，而其他的松弛參數(shù)均取1.由于松弛矩陣中的松弛參數(shù)有無限種組合方式，因此出于通用性考慮，我們選擇了這種處理方法.同時(shí)需要指出的是當(dāng)松弛參數(shù)矩陣中的松弛系數(shù)相同時(shí)，多松弛模型退化到單松弛模型，即松弛矩陣中的松弛參數(shù)均為sr.

2 結(jié)果及分析

2.1 具有高斯型發(fā)射場的一維無限大平板

考慮一充滿吸收發(fā)射性介質(zhì)的一維無限大平板內(nèi)的輻射傳遞問題，平板內(nèi)具有一高斯型發(fā)射場，該問題由如下方程控制

(27)

考慮如下邊界條件

I(0，ξ)=β-1e-b2/α2，ξ>0

(28)

該問題存在解析解形式，其表達(dá)式如下

(29)

考慮方向ξ=1.0，a=0.02，b=0.5，采用LBM來模擬衰減系數(shù)為β=1，10和50 m-1時(shí)介質(zhì)內(nèi)輻射強(qiáng)度的分布，取100個(gè)格子，無量綱時(shí)間步長取Δt*=0.000 1，單松弛模型得到的結(jié)果和解析解對(duì)比，如圖1所示，LBM得到的輻射強(qiáng)度分布和解析解得到的輻射強(qiáng)度分布吻合地很好.

圖1 衰減系數(shù)為β=1，10和50m-1時(shí)LBM得到的輻射強(qiáng)度分布和解析解對(duì)比

接下來，我們進(jìn)一步研究一維多松弛模型的穩(wěn)定性和精度.為了研究穩(wěn)定性，我們考慮衰減系數(shù)為10 m-1的情況，取100個(gè)格子，研究不同松弛參數(shù)下所允許的最大時(shí)間步長.數(shù)值解和解析解的相對(duì)誤差定義為

(30)

穩(wěn)定性標(biāo)準(zhǔn)為數(shù)值解和解析解的相對(duì)誤差小于10-2.表1給出了不同松弛參數(shù)下所允許的最大時(shí)間步長，不同參數(shù)的最大時(shí)間步長得到是根據(jù)我們定義的穩(wěn)定性標(biāo)準(zhǔn)，然后通過數(shù)值實(shí)驗(yàn)得到的，可以發(fā)現(xiàn)多松弛模型允許的最大時(shí)間步長可以隨松弛參數(shù)調(diào)整，尤其當(dāng)松弛參數(shù)小于1時(shí)，所允許的時(shí)間步長大于單松弛模型，結(jié)果表明相比單松弛模型，多松弛模型可以在更大的時(shí)間步長內(nèi)保持穩(wěn)定，具有更好的穩(wěn)定性.多松弛模型的碰撞過程發(fā)生在矩空間，與多個(gè)速度分布函數(shù)相關(guān)聯(lián)，相比單松弛模型發(fā)生在速度空間的碰撞，多松弛模型本身在穩(wěn)定性方面展現(xiàn)了很大的優(yōu)勢，數(shù)值結(jié)果證明了多松弛模型在穩(wěn)定性上的優(yōu)勢.此外，表2給出了不同格子數(shù)下單松弛和多松弛模型的相對(duì)誤差，可以看出多松弛模型相比單松弛模型具有更高的精度.

表1 衰減系數(shù)β=10 m-1，100個(gè)格子下，單松弛(BGK)和多松弛(MRT)模型允許的最大時(shí)間步長

表2 衰減系數(shù)β=10 m-1，不同格子數(shù)下，單松弛(BGK)和多松弛(MRT)模型的相對(duì)誤差

2.2 受高斯型脈沖照射的一維純散射介質(zhì)

考慮厚度為L=1 m的一維半透明平板介質(zhì)內(nèi)的瞬態(tài)輻射傳輸問題.介質(zhì)為各向同性散射，壁面和介質(zhì)溫度均為0 K，無發(fā)射.介質(zhì)邊界為透明邊界，環(huán)境為真空.平板介質(zhì)的衰減系數(shù)為1 m-1，右側(cè)邊界無照射，左側(cè)邊界受到如下法向平行光入射輻射的照射：

(31)

公式中：I0為脈沖輻射強(qiáng)度；H(t)為Heaviside階躍函數(shù).

圖2(a)和圖2(b)中的松弛參數(shù)分別取sr=0.8，LBM得到的計(jì)算結(jié)果和文獻(xiàn)[24]中Liu和Hsu采用間斷有限元方法得到的結(jié)果進(jìn)行對(duì)比.由圖2可知本文的計(jì)算結(jié)果和文獻(xiàn)結(jié)果對(duì)比吻合很好，證明了本文提出的多松弛格子-Boltzmann模型可以穩(wěn)定精確地求解參與性介質(zhì)內(nèi)一維瞬態(tài)輻射傳遞問題.本文的LB模型執(zhí)行基于碰撞和遷徙的顯式瞬態(tài)演化過程.其它基于離散化偏微分方程的數(shù)值方法，如FVM、DOM、無網(wǎng)格法和FEM等，一般都需要在每個(gè)時(shí)間步長內(nèi)進(jìn)行全局迭代收斂，在效率上并不優(yōu)于LB模型.因此，本文提出的LB模型非常適合于求解瞬態(tài)輻射傳輸問題.在處理穩(wěn)態(tài)輻射傳輸問題時(shí)，與基于離散化偏微分方程的數(shù)值方法相比，LB模型的缺點(diǎn)是計(jì)算效率較低.這是因?yàn)長BM通過依賴時(shí)間的演化過程來求解穩(wěn)態(tài)問題.對(duì)于穩(wěn)態(tài)問題，這些基于離散化偏微分方程的數(shù)值方法只需對(duì)穩(wěn)態(tài)輻射傳輸方程進(jìn)行全局迭代收斂即可獲得收斂解，因而對(duì)穩(wěn)態(tài)輻射傳輸問題具有較高的計(jì)算效率.

圖2 高斯型脈沖照射下平板界面處時(shí)域反射率和透射率信號(hào)

2.3 二維方腔內(nèi)輻射傳遞問題

本算例考慮二維方腔內(nèi)的輻射傳遞問題，方腔的邊長為L=1 m，衰減系數(shù)為1 m-1，方腔內(nèi)介質(zhì)為熱介質(zhì)，所有壁面保持0 K，四個(gè)壁面均為黑壁面.方腔內(nèi)充滿吸收散射性介質(zhì)，散射為各向異性散射，散射相函數(shù)F2的不對(duì)稱因子為0.669 72(見文獻(xiàn)[25]).將計(jì)算域劃分為60×60的格子，無量綱時(shí)間步長取Δt*=0.000 5，空間離散采用S8方案.松弛參數(shù)取0.8時(shí)的多松弛模型得到的方腔底部無量綱熱流，如圖3所示.同時(shí)文獻(xiàn)[25]中采用離散坐標(biāo)法得到的計(jì)算結(jié)果也顯示在圖3中作為對(duì)比，結(jié)果顯示LBM得到的計(jì)算結(jié)果和離散坐標(biāo)法得到的計(jì)算結(jié)果吻合得很好，進(jìn)一步驗(yàn)證了本文發(fā)展得多松弛格子-Boltzmann模型可以準(zhǔn)確地求解參與性介質(zhì)內(nèi)多維輻射傳輸問題.

圖3 不同散射反照率下方腔底部的無量綱輻射熱流

3 結(jié) 論

對(duì)于輻射方程，本文提出建立了一種多松弛的格子-Boltzmann模型.對(duì)多松弛格子-Boltzmann模型進(jìn)行Maxwell迭代，可以嚴(yán)格地得出輻射方程.相比于已有的格子-Boltzmann模型，本文提出的模型沒有任何近似和限制.

數(shù)值實(shí)驗(yàn)表明本文發(fā)展的多松弛格子-Boltzmann模型可以穩(wěn)定精確地求解參與性介質(zhì)內(nèi)一維、二維瞬態(tài)和穩(wěn)態(tài)輻射傳輸問題.同時(shí)，理論結(jié)果表明該模型具有二階精度.本文構(gòu)建的多松弛的格子-Boltzmann模型可以退化到單松弛格子-Boltzmann模型，且多松弛模型相比于單松弛模型具有更好的穩(wěn)定性.采用三維的格子模型，多松弛模型也可以用來求解三維輻射傳輸問題，與此同時(shí)，本文提出的多松弛的格子-Boltzmann模型可以進(jìn)一步推廣到求解參與性介質(zhì)內(nèi)復(fù)雜的輻射傳輸問題.