基于孔隙尺度的方腔內(nèi)多孔介質(zhì)融化模擬研究

王猛，朱衛(wèi)兵

（哈爾濱工程大學(xué)航天與建筑工程學(xué)院，黑龍江哈爾濱150001）

基于孔隙尺度的方腔內(nèi)多孔介質(zhì)融化模擬研究

王猛，朱衛(wèi)兵

（哈爾濱工程大學(xué)航天與建筑工程學(xué)院，黑龍江哈爾濱150001）

由于單松弛（LBGK）格子Boltzmann模型在用反彈格式處理無滑移邊界時存在缺陷。基于孔隙尺度，采用多松弛（MRT）格子Boltzmann模型研究封閉方腔內(nèi)多孔介質(zhì)的自然對流融化過程，其中，通過焓方法考慮相變潛熱。分析了Rayleigh數(shù)和Prandtl數(shù)對融化的影響。結(jié)果表明：采用的多松弛模型能很好的預(yù)測導(dǎo)熱和對流融化過程；多孔介質(zhì)的導(dǎo)熱融化界面不再與垂直壁面平行，自然對流融化界面呈現(xiàn)不規(guī)則形狀；Rayleigh數(shù)和Prandtl數(shù)對多孔介質(zhì)的融化有較大影響。

多孔介質(zhì)；孔隙尺度；融化；自然對流；格子Boltzmann方法

多孔介質(zhì)的融化存在于很多領(lǐng)域：如地質(zhì)學(xué)中，固相和液相共存于在巖漿庫中，因此而具有復(fù)雜的動力學(xué)行為，巖漿系統(tǒng)經(jīng)歷再加熱和冷卻的演化過程，其融化或凝固導(dǎo)致復(fù)雜的局部融化分?jǐn)?shù)的變化；在航空航天領(lǐng)域中，碳／陶復(fù)合材料在燒蝕過程中表現(xiàn)出明顯的多孔介質(zhì)特性，SiC／ZrC氧化生成的SiO2／ZrO2，當(dāng)溫度達(dá)到其熔點(diǎn)時，會發(fā)生融化。因此，有必要對多孔介質(zhì)的融化做深入的探究。

格子Boltzmann方法（lattice Boltzmann method，LBM）具有模型簡單、處理復(fù)雜邊界方便及并行性好等獨(dú)特優(yōu)勢，在多孔介質(zhì)、多組分、多相流、化學(xué)反應(yīng)以及微尺度流動等領(lǐng)域有廣泛的應(yīng)用［1］。利用LBM研究固液相變問題的方法主要是相場法（phase?field method）和焓方法（enthalpy?based meth?od）。相場法為了得到精細(xì)的相界面要求足夠小的網(wǎng)格尺度，這就需要巨大的計(jì)算量。焓方法則對此沒有要求，Jiaung等［2］首次將焓方法引入LBM中，研究了導(dǎo)熱融化和凝固問題。Chakraborty和Chat?terjee［3］將焓方法與流動耦合研究了晶體生長等固液相變問題。Johann等［4］利用多松弛（multiple re?laxation time，MRT）模型［5］求解流場，有限差分法求解能量方程研究了自然對流融化，結(jié)果表明所用模型具有二階精度。Gao等［6］基于表征體元尺度，采用焓方法模擬了填充了多孔介質(zhì)的方腔內(nèi)相變材料的自然對流融化過程。Chen等［7］基于孔隙尺度研究了金屬泡沫內(nèi)相變材料的融化過程。Huber等［8］基于孔隙尺度，對多孔介質(zhì)融化的Rayleigh?Be'nard對流做了模擬，但是該文僅給出了融化界面和溫度的演化過程，未做深入的分析。上述研究采用的模型大多為單松弛模型（single relaxation time，SRT）（也稱為LBGK（lattice Bhatnagar?Gross?Krook））［9］，LBGK是最簡單也是目前應(yīng)用最廣泛的LBM模型，但是該模型的數(shù)值穩(wěn)定性較差，且當(dāng)用反彈格式處理無滑移邊界時，邊界處的滑移速度依賴于松弛時間［10］。而多松弛模型可以克服LBGK的這一固有缺陷，MRT有多個的松弛時間，可通過調(diào)節(jié)自由參數(shù)提高數(shù)值穩(wěn)定性和精度。

本文采用MRT模型模擬流場，LBGK模擬溫度場，利用焓方法考慮相變潛熱，首先對純物質(zhì)的導(dǎo)熱融化和自然對流融化進(jìn)行模擬，以驗(yàn)證模型的正確性；然后基于孔隙尺度研究了封閉方腔內(nèi)多孔介質(zhì)的自然對流融化過程，考慮了Rayleigh數(shù)和Prandtl數(shù)對融化的影響。

1 格子Boltzmann模型

1.1 流動的多松弛格子Boltzmann模型

考慮到MRT具有較好的數(shù)值穩(wěn)定性及可消除LBGK的邊界處速度依賴于粘性的缺陷，本文流場采用考慮外力的MRT模型［11］，其在矩空間執(zhí)行碰撞過程，而在速度空間執(zhí)行遷移過程，其演化方程如下：

式中：f為粒子的密度分布函數(shù)的矢量形式，對于本文采用的二維九速（D2Q9）模型中，f＝（f0，f1，…，f8）T；δt為時間步長；ei為各方向的粒子速度，由下式確定：

式中：c＝δx／δt為格子速度，δx為格子間距。

M為變換矩陣，具體形式見文獻(xiàn)［5］；m和meq分別為矩空間的分布函數(shù)和平衡態(tài)分布函數(shù)：

式中：ρ為宏觀密度，ux和uy分別為x和y方向的宏觀速度。

F＝（F0，F(xiàn)1，···，F(xiàn)8）T為外力項(xiàng)，定義如下：

式中：F'＝（F'0，F(xiàn)'1，···，F(xiàn)'8）T，各分量具體形式如下：

式中：ωi為權(quán)系數(shù)，ω0＝4／9，ωi＝1／9（i＝1，2，3，4），ωi＝1／36（i＝5，6，7，8）；為格子聲速；G＝-β（T-T0）ρg，β為熱膨脹系數(shù)，T為溫度，T0為參考溫度，g為重力加速度。

S為松弛矩陣，確定如下：

當(dāng)si＝1／τ（i＝0，1，2，…，8）時，MRT模型退化為LB? GK模型。本文中，取se＝sε＝sv，sq＝8（2-sν）／（8-sν）［12］。sv用于確定剪切粘性系數(shù)：

宏觀密度ρ和速度u由下式計(jì)算：

1.2 熱格子Boltzmann模型

溫度場采用LBGK二維五速（D2Q5）模型，基于焓方法，通過在溫度方程中引入一個源項(xiàng)模擬相變，利用Huber等［8］的處理方法，溫度場的演化方程為

式中：gi和分別為粒子的溫度分布函數(shù)和平衡態(tài)溫度分布函數(shù)；τh為溫度場松弛時間；ω'i為權(quán)系數(shù)，ω'0＝1／3，ω'i＝1／6（i＝1，2，3，4）；L為相變潛熱；C為熱容；γfl為液相分?jǐn)?shù)。gieq確定為

熱擴(kuò)散系數(shù)和宏觀溫度分別為

將總焓H分為顯焓和潛熱2部分

通過上式確定H后，計(jì)算時間步t的液相分?jǐn)?shù)為

式中：Hs和Hl分別為固相線和液相線的焓值，Tm為融化溫度。

1.3 簡化與假設(shè)

本文的模擬主要基于以下幾個條件：1）融化后的液體為牛頓流體且不可壓；2）流體流動基于Boussinesq假設(shè)；3）材料融化前后的熱物性相同且保持常數(shù)。

2 模型驗(yàn)證

2.1 導(dǎo)熱融化

首先選取純物質(zhì)的導(dǎo)熱融化問題驗(yàn)證模型及程序。對于此問題，Neumann得到了半無限大空間的解析解，溫度分布和固液界面的位置分別為

式中：T1為過熱溫度，k為導(dǎo)熱系數(shù)。λ的值由以下超越方程確定：

式中：St＝ρc（T1-Tm）／L為Stefan數(shù)。

初始固相溫度為相變溫度Tm，左壁面溫度突然升高到T1，開始融化。在模擬中，上下壁面為周期性邊界，網(wǎng)格數(shù)取為100×10。此處，Stefan數(shù)取為0.013，對應(yīng)的λ為0.25。圖1和圖2分別為融化前沿位置（固液界面）隨時間的變化和不同時刻溫度的分布曲線，可以看到，模擬結(jié)果與解析解吻合較好。

圖1 融化前沿位置隨時間的變化Fig.1 The evolution of melting front position with time

圖2 不同時刻溫度的分布Fig.2 The temperature distribution at different mon?ents

2.2 對流融化

為了進(jìn)一步驗(yàn)證模型，本節(jié)對封閉方腔內(nèi)純物質(zhì)的自然對流融化進(jìn)行模擬。初始固相溫度為相變溫度Tm，左壁面溫度突然升高到T1，并保持不變；右壁面仍然保持在相變溫度Tm，上下壁面為絕熱壁面。對于速度場，4個壁面均采用二階精度的Half-Way反彈格式處理無滑移邊界。長寬均為l，網(wǎng)格數(shù)為200×200。對于對流融化問題，由幾個基本的無量綱數(shù)確定融化特性：Rayleigh數(shù)Ra、Prandtl數(shù)Pr和Ste?fan數(shù)。Rayleigh數(shù)和Prandtl數(shù)分別定義為

其中，ΔT為左右壁面溫差ΔT＝T1-Tm。此處，取Ra＝2.5×104，Pr＝0.02，St＝0.01。

Nusselt數(shù)定義為

式中：θ＝T-Tm／T1-Tm為無量綱溫度。

Mencinger［13］采用自適應(yīng)網(wǎng)格研究了自然對流融化過程，下面將本文結(jié)果與Mencinger的結(jié)論進(jìn)行對比。圖3為總液相分?jǐn)?shù)隨Fourier數(shù)的變化，可以看到，本文結(jié)果與Mencinger的結(jié)果吻合很好。圖4給出了左壁面（熱壁面）平均Nusselt數(shù)隨Fou?rier數(shù)的變化，本文結(jié)果與Mencinger結(jié)果差別不大，其差異可能是由數(shù)值方法不同引起的。通過以上算例的驗(yàn)證，證明了模型的正確性。

圖3 總液相分?jǐn)?shù)隨Fourier數(shù)的變化Fig.3 The evolution of the total liquid fraction with the Fourier number

圖4 左壁面平均Nusselt數(shù)隨Fourier數(shù)的變化Fig.4 The evolution of the average Nusselt number a?long the left wall with the Fourier number

3 方腔內(nèi)多孔介質(zhì)的融化

上節(jié)計(jì)算了純物質(zhì)的導(dǎo)熱融化和自然對流融化，本節(jié)對封閉方腔內(nèi)多孔介質(zhì)的融化進(jìn)行模擬。模擬所用多孔介質(zhì)結(jié)構(gòu)如圖5所示，圖中黑色區(qū)域?yàn)楣腆w骨架，白色區(qū)域?yàn)榭紫丁３跏技斑吔鐥l件與2.2節(jié)純物質(zhì)的自然對流融化相同，不再贅述。網(wǎng)格數(shù)為300× 300，Prandtl數(shù)和Stefan數(shù)均取為1.0，Rayleigh數(shù)在1.0×105～1.0×108取4個不同的值進(jìn)行討論。圖6為多孔介質(zhì)融化界面隨Fourier數(shù)的演化，與純物質(zhì)自然對流融化類似，初始導(dǎo)熱作用占主導(dǎo)，融化界面基本與垂直壁面平行，隨著融化的進(jìn)行，由于對流作用的增強(qiáng)，多孔介質(zhì)上部融化快于下部，且界面呈現(xiàn)不規(guī)則形狀，這是由多孔介質(zhì)的不均勻性引起的。對比圖6（a）和（b）可以看到，Rayleigh數(shù)越大（對流作用越強(qiáng)），界面形狀越不規(guī)則。在Ra＝1.0×108（見圖6（b）），F(xiàn)o＝0.003時，由于介質(zhì)存在一個孔隙連通區(qū)域，使流動超過了融化前沿位置，但是在Rayleigh數(shù)較小時（見圖6（a）），未出現(xiàn)此現(xiàn)象。從以上分析可知，Rayleigh數(shù)對多孔介質(zhì)的融化過程影響較大。這里也給出了方腔內(nèi)多孔介質(zhì)導(dǎo)熱融化過程，見圖6（c），可以看到，多孔介質(zhì)與純物質(zhì)的導(dǎo)熱融化也存在較大差別，在純物質(zhì)導(dǎo)熱融化中，融化界面始終與垂直壁面平行，但由于多孔介質(zhì)的不均勻性，導(dǎo)致多孔介質(zhì)導(dǎo)熱融化過程中融化前沿的上、下不再同步，當(dāng)融化界面接近右壁面時，下部先融化完（這取決于具體的多孔介質(zhì)結(jié)構(gòu)）。數(shù)的變化，為了說明對流傳熱對融化的影響，圖7也給出了相應(yīng)的導(dǎo)熱融化過程，以便比較。初始融化過程靠導(dǎo)熱控制，不同Rayleigh數(shù)下液相分?jǐn)?shù)線重合，隨液相分?jǐn)?shù)的增大，對流作用增強(qiáng)，Rayleigh數(shù)越大，對流區(qū)出現(xiàn)的越早，且融化越快。在對流主導(dǎo)區(qū)，融化分?jǐn)?shù)隨Fourier數(shù)基本呈線性增加，當(dāng)融化界面接近右壁面時，融化速率減慢。

圖5 多孔介質(zhì)結(jié)構(gòu)Fig.5 The structure of porous media

圖6 多孔介質(zhì)融化界面隨Fourier數(shù)的變化Fig.6 The evolution of melting interface of porous media with time

圖7 不同Rayleigh數(shù)下總液相分?jǐn)?shù)隨Fourier數(shù)的變化Fig.7 The evolution of total liquid fraction with the Fourier number at different Rayleigh numbers

圖8 不同Rayleigh數(shù)下Nusselt數(shù)隨Fourier數(shù)的變化Fig.8 The evolution of average Nusselt number with the Fourier number at different Rayleigh numbers

不同Rayleigh數(shù)下熱壁面平均Nusselt數(shù)隨Fourier數(shù)的變化趨勢見圖8，可以看到，由于開始融化時，導(dǎo)熱為主要傳熱機(jī)制，使Nusselt數(shù)急劇減小，隨著對流作用的增強(qiáng)，限制了Nusselt數(shù)的減小，使其達(dá)到一個最小值，然后趨于穩(wěn)定。Rayleigh數(shù)越大，穩(wěn)定后的Nusselt數(shù)越大。圖中還可以看到，當(dāng)Ra＝1.0×108時，Nusselt數(shù)存在一個峰值，對比圖6（b）Fo＝0.003，這是因?yàn)槎嗫捉橘|(zhì)為非均勻的，此時流動超過了融化前沿位置，在流動后面未融化的骨架區(qū)會對流動產(chǎn)生擾動作用，從而增強(qiáng)對流，使Nusselt數(shù)上升，當(dāng)此部分固體融化后，Nusselt數(shù)又有所下降。Gobin等［14］研究了Prandtl數(shù)對純物質(zhì)自然對流融化的影響（0.01＜Pr＜0.1），認(rèn)為Prandtl數(shù)對融化有較大影響。本文考慮大范圍Prandtl數(shù)對多孔介質(zhì)自然對流融化的影響，計(jì)算中，Rayleigh數(shù)和Stefan數(shù)分別取為1.0×105為和1.0，而Prandtl數(shù)取0.01、0.1、1.0和10.0 4組做對比。圖9和圖10分別為Prandtl數(shù)對總液相分?jǐn)?shù)和熱壁面Nusselt數(shù)的影響，可以看到，在低Prandtl數(shù)范圍內(nèi)，Prandtl數(shù)對多孔介質(zhì)的融化影響較大，但是Pr＝1.0和Pr＝10.0得到的總液相分?jǐn)?shù)和Nusselt數(shù)曲線差別很小，也就是說，當(dāng)Pr＞1.0時，繼續(xù)增大Prandtl數(shù)對融化（傳熱）基本沒有影響。在開始融化階段，主要為純導(dǎo)熱融化，不同Prandtl數(shù)下液相分?jǐn)?shù)曲線重合（見圖9），一旦液相分?jǐn)?shù)達(dá)到一定比例，對流開始起作用，在低Prandtl數(shù)范圍內(nèi)，Prandtl數(shù)越大（相對動量擴(kuò)散率越大，或者相對熱擴(kuò)散率越小），對流越強(qiáng)，融化更快，且達(dá)到穩(wěn)定時的Nusselt數(shù)也更大（見圖10）。

圖9 不同Prandtl數(shù)下總液相分?jǐn)?shù)隨Fourier數(shù)的變化Fig.9 The evolution of total liquid fraction with the Fourier number at different Prandtl numbers

圖10 不同Prandtl數(shù)下Nusselt數(shù)隨Fourier數(shù)的變化Fig.10 The evolution of average Nusselt number with the Fourier number at different Prandtl numbers

4 結(jié)論

1）通過驗(yàn)證，本文采用的多松弛模型能夠很好的預(yù)測導(dǎo)熱和自然對流融化；

2）方腔內(nèi)多孔介質(zhì)的導(dǎo)熱融化界面不再與垂直壁面平行，自然對流融化出現(xiàn)不規(guī)則形狀的固液界面；

3）增大Rayleigh數(shù)，多孔介質(zhì)融化加快，融化界面更不規(guī)則，熱壁面Nusselt數(shù)增大；在低Prandtl數(shù)范圍內(nèi)，Prandtl數(shù)越大，融化越快，Nusselt數(shù)越大，但到一定值后，繼續(xù)增大Prandtl數(shù)基本不再影響融化過程。

［1］NOURGALIEV R.The lattice Boltzmann equation method：the?oretical interpretation，numerics and implications［J］.Interna?tional Journal of Multiphase Flow，2003，29（1）：117?169.

［2］JIAUNG W S，CHUN J R H.Lattice Boltzmann method for the heat conduction problem with phase change［J］.Numeri?cal Heat Transfer，Part B：An International Journal of Com?putation and Methodology，2001，39（2）：167?187.

［3］CHAKRABORTY S，CHATTERJEE D.An enthalpy?based hybrid lattice?Boltzmann method for modelling solid?liquid phase transition in the presence of convective transport［J］.Journal of Fluid Mechanics，2007，592：155?175.

［4］JOHANN M，KUZNIK F，JOHANNES K，et al.Develop?ment and validation of a new LBM?MRT hybrid model with enthalpy formulation for melting with natural convection［J］.Physics Letters A，2014，378（4）：374?381.

［5］LALLEMAND P，LUO L S.Theory of the lattice Boltzmann method：Dispersion，dissipation，isotropy，Galilean invari?ance，and stability［J］.Physical Review E，2000，61（6）：6546?6562.

［6］GAO D Y，CHEN Z Q.Lattice Boltzmann simulation of nat?ural convection dominated melting in a rectangular cavity filled with porous media［J］.International Journal of Thermal Sciences，2011，50（4）：493?501.

［7］CHEN Z，GAO D，SHI J.Experimental and numerical study on melting of phase change materials in metal foams at pore scale［J］.International Journal of Heat and Mass Transfer，2014，72：646?655.

［8］HUBER C，PARMIGIANI A，CHOPARD B，et al.Lattice Boltzmann model for melting with natural convection［J］.In?ternational Journal of Heat and Fluid Flow，2008，29（5）：1469?1480.

［9］QIAN Y H，D'HUMIèRES D，LALLEMAND P.Lattice BGK models for Navier?Stokes equation［J］.EPL（Euro?physics Letters），1992，17（6）：479.

［10］HE X，ZOU Q，LUO L S，et al.Analytic solutions of sim?ple flows and analysis of nonslip boundary conditions for the lattice Boltzmann BGK model［J］.Journal of Statistical Physics，1997，87（1?2）：115?136.

［11］GUO Z，ZHENG C.Analysis of lattice Boltzmann equation for microscale gas flows：Relaxation times，boundary condi?tions and the Knudsen layer［J］.International Journal of Computational Fluid Dynamics，2008，22（7）：465?473.

［12］朱衛(wèi)兵，王猛，陳宏，等.多孔介質(zhì)內(nèi)流體流動的格子Boltzmann模擬［J］.化工學(xué)報，2013，64（S1）：33?40.ZHU Weibing，WANG Meng，CHEN Hong，et al.Lattice Boltzmann simulations for fluid flow through porous media［J］.CIESC Journal，2013，64（S1）：33?40.

［13］MENCINGER J.Numerical simulation of melting in two?di?mensional cavity using adaptive grid［J］.Journal of Compu?tational Physics，2004，198（1）：243?264.

［14］GOBIN D，BENARD C.Melting of metals driven by natural convection in the melt：influence of Prandtl and Rayleigh numbers［J］.Journal of Heat Transfer，1992，114（2）：521?524.

Simulation study for melting of porous media in a square cavity based on the pore scale

WANG Meng，ZHU Weibing

（College of Aerospace and Civil Engineering，Harbin Engineering University，Harbin 150001，China）

：When a bounce?back scheme is used to handle no?slip boundary conditions，the lattice Bhatnagar?Gross?Krook（LBGK）model in lattice Boltzmann method（LBM）will result in viscosity?dependence of boundary locations.So in this paper，the multiple relaxation time（MRT）model was chosen to research the naturally convec?ted melting process of porous media in a closed square cavity based on the pore scale.In particular，the enthalpy?based method was used to simulate solid?liquid phase change.The effects of Rayleigh and Prandtl numbers on melt?ing were also analyzed.The results indicate that the MRT model can well predict conduction?dominated and convec?tion?dominated melting.The solid?liquid interface is no longer parallel to vertical walls for conduction melting of porous media.The interface appears irregular shapes in the natural convection?dominated melting of porous media.The Rayleigh and Prandtl numbers have larger influence on the melting of porous media.

porous media；pore scale；melting；natural convection；lattice Boltzmann method

10.3969／j.issn.1006?7043.201404044

V257

：A

：1006?7043（2015）06?0774?05

http：／／www.cnki.net／kcms／detail／23.1390.u.20150428.1116.019.html

2014?04?13.網(wǎng)絡(luò)出版時間：2015?04?26.

高超聲速沖壓發(fā)動機(jī)技術(shù)重點(diǎn)實(shí)驗(yàn)室開放基金資助項(xiàng)目（20140103008）；哈爾濱市科技創(chuàng)新人才研究專項(xiàng)資助項(xiàng)目（2014RFXXJ062）.

王猛（1988?），男，博士研究生；朱衛(wèi)兵（1961?），男，教授，博士生導(dǎo)師.

朱衛(wèi)兵，E?mail：zhuweibing＠hrbeu.edu.cn.