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

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

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

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

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

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

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

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

1.1 輻射傳輸方程

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

(1)

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

(2)

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

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

(3)

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

1.2 多松弛格子-Boltzmann模型

瞬態輻射常常發生于極短的時間內，在瞬態輻射的模擬中，通常引入無量綱時間來避免過小的時間步長.將無量綱時間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方程到輻射傳輸方程

本節基于擴散尺度Δ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模型出發，基于擴散尺度下的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 具有高斯型發射場的一維無限大平板

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

(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得到的輻射強度分布和解析解對比

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

(30)

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

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

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

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

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

(31)

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

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

圖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得到的計算結果和離散坐標法得到的計算結果吻合得很好，進一步驗證了本文發展得多松弛格子-Boltzmann模型可以準確地求解參與性介質內多維輻射傳輸問題.

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

3 結 論

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

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