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

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

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

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

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

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

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

1.1 輻射傳輸方程

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

(1)

公式中：cL為介質內的光速；I為輻射強度；r為位置坐標；β=ka+ks為衰減系數；Ωm=μmi+ηmj+ξmk為離散方向，源項S可以表示為

(2)

公式中：N為總的離散方向，m=1，2，…，N，m′=1，2，…，N；wm′為對應方向的權重.

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

(3)

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

1.2 多松弛格子-Boltzmann模型

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

(4)

公式中

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

(5)

公式中：LR為介質的參考長度.

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

(6)

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

(7)

輻射強度可以由平衡函數給出，關系如下

(8)

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

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

(9)

(10)

(11)

對于D2Q9模型，其格子信息為

(12)

(13)

(14)

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

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

首先，令f(r，t*)=(f0(r，t*)，f1(r，t*)，…，f8(r，t*))T，ω=(ω0，ω1，…，ω8)T，時間無量綱形式的輻射傳遞方程(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)左邊應用Taylor展開，

(16)

其中微分算子Ds

(17)

矩陣

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

(18)

公式中：p和q均為非負整數.

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

(19)

其中

(20)

(21)

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

(22)

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

(23)

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

(24)

至此，我們從多松弛格子-Boltzmann模型出發(fā)，基于擴散尺度下的Maxwell迭代，嚴格推導得出了輻射傳輸方程，并且可以從方程(24)理論上得出該模型具有二階的精度.一般而言，對于對流擴散問題，計算流體力學等問題的LB模型，其中的松弛系數與宏觀方程中的擴散系數，流體黏性系數等有定量關系.需要指出的是，根據從多松弛格子-Boltzmann模型嚴格推導得出輻射傳輸方程可知，本文提出的多松弛格子-Boltzmann模型中的松弛參數均是自由的，與其他參數無關.對于一維和二維LB模型，我們取如下的松弛參數矩陣

S=diag(1，sr，1)

(25)

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

(26)

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

2 結果及分析

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

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

(27)

考慮如下邊界條件

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

(28)

該問題存在解析解形式，其表達式如下

(29)

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

圖1 衰減系數為β=1，10和50m-1時LBM得到的輻射強度分布和解析解對比

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

(30)

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

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

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

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

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

(31)

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

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

圖2 高斯型脈沖照射下平板界面處時域反射率和透射率信號

2.3 二維方腔內輻射傳遞問題

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

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

3 結 論

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

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