基于孔隙尺度的方腔内多孔介质融化模拟研究

王猛，朱卫兵

（哈尔滨工程大学航天与建筑工程学院，黑龙江哈尔滨150001）

基于孔隙尺度的方腔内多孔介质融化模拟研究

王猛，朱卫兵

（哈尔滨工程大学航天与建筑工程学院，黑龙江哈尔滨150001）

由于单松弛（LBGK）格子Boltzmann模型在用反弹格式处理无滑移边界时存在缺陷。基于孔隙尺度，采用多松弛（MRT）格子Boltzmann模型研究封闭方腔内多孔介质的自然对流融化过程，其中，通过焓方法考虑相变潜热。分析了Rayleigh数和Prandtl数对融化的影响。结果表明：采用的多松弛模型能很好的预测导热和对流融化过程；多孔介质的导热融化界面不再与垂直壁面平行，自然对流融化界面呈现不规则形状；Rayleigh数和Prandtl数对多孔介质的融化有较大影响。

多孔介质；孔隙尺度；融化；自然对流；格子Boltzmann方法

多孔介质的融化存在于很多领域：如地质学中，固相和液相共存于在岩浆库中，因此而具有复杂的动力学行为，岩浆系统经历再加热和冷却的演化过程，其融化或凝固导致复杂的局部融化分数的变化；在航空航天领域中，碳／陶复合材料在烧蚀过程中表现出明显的多孔介质特性，SiC／ZrC氧化生成的SiO2／ZrO2，当温度达到其熔点时，会发生融化。因此，有必要对多孔介质的融化做深入的探究。

格子Boltzmann方法（lattice Boltzmann method，LBM）具有模型简单、处理复杂边界方便及并行性好等独特优势，在多孔介质、多组分、多相流、化学反应以及微尺度流动等领域有广泛的应用［1］。利用LBM研究固液相变问题的方法主要是相场法（phase⁃field method）和焓方法（enthalpy⁃based meth⁃od）。相场法为了得到精细的相界面要求足够小的网格尺度，这就需要巨大的计算量。焓方法则对此没有要求，Jiaung等［2］首次将焓方法引入LBM中，研究了导热融化和凝固问题。Chakraborty和Chat⁃terjee［3］将焓方法与流动耦合研究了晶体生长等固液相变问题。Johann等［4］利用多松弛（multiple re⁃laxation time，MRT）模型［5］求解流场，有限差分法求解能量方程研究了自然对流融化，结果表明所用模型具有二阶精度。Gao等［6］基于表征体元尺度，采用焓方法模拟了填充了多孔介质的方腔内相变材料的自然对流融化过程。Chen等［7］基于孔隙尺度研究了金属泡沫内相变材料的融化过程。Huber等［8］基于孔隙尺度，对多孔介质融化的Rayleigh⁃Be'nard对流做了模拟，但是该文仅给出了融化界面和温度的演化过程，未做深入的分析。上述研究采用的模型大多为单松弛模型（single relaxation time，SRT）（也称为LBGK（lattice Bhatnagar⁃Gross⁃Krook））［9］，LBGK是最简单也是目前应用最广泛的LBM模型，但是该模型的数值稳定性较差，且当用反弹格式处理无滑移边界时，边界处的滑移速度依赖于松弛时间［10］。而多松弛模型可以克服LBGK的这一固有缺陷，MRT有多个的松弛时间，可通过调节自由参数提高数值稳定性和精度。

本文采用MRT模型模拟流场，LBGK模拟温度场，利用焓方法考虑相变潜热，首先对纯物质的导热融化和自然对流融化进行模拟，以验证模型的正确性；然后基于孔隙尺度研究了封闭方腔内多孔介质的自然对流融化过程，考虑了Rayleigh数和Prandtl数对融化的影响。

1 格子Boltzmann模型

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

考虑到MRT具有较好的数值稳定性及可消除LBGK的边界处速度依赖于粘性的缺陷，本文流场采用考虑外力的MRT模型［11］，其在矩空间执行碰撞过程，而在速度空间执行迁移过程，其演化方程如下：

式中：f为粒子的密度分布函数的矢量形式，对于本文采用的二维九速（D2Q9）模型中，f＝（f0，f1，…，f8）T；δt为时间步长；ei为各方向的粒子速度，由下式确定：

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

M为变换矩阵，具体形式见文献［5］；m和meq分别为矩空间的分布函数和平衡态分布函数：

式中：ρ为宏观密度，ux和uy分别为x和y方向的宏观速度。

F＝（F0，F1，···，F8）T为外力项，定义如下：

式中：F'＝（F'0，F'1，···，F'8）T，各分量具体形式如下：

式中：ωi为权系数，ω0＝4／9，ωi＝1／9（i＝1，2，3，4），ωi＝1／36（i＝5，6，7，8）；为格子声速；G＝-β（T-T0）ρg，β为热膨胀系数，T为温度，T0为参考温度，g为重力加速度。

S为松弛矩阵，确定如下：

当si＝1／τ（i＝0，1，2，…，8）时，MRT模型退化为LB⁃ GK模型。本文中，取se＝sε＝sv，sq＝8（2-sν）／（8-sν）［12］。sv用于确定剪切粘性系数：

宏观密度ρ和速度u由下式计算：

1.2 热格子Boltzmann模型

温度场采用LBGK二维五速（D2Q5）模型，基于焓方法，通过在温度方程中引入一个源项模拟相变，利用Huber等［8］的处理方法，温度场的演化方程为

式中：gi和分别为粒子的温度分布函数和平衡态温度分布函数；τh为温度场松弛时间；ω'i为权系数，ω'0＝1／3，ω'i＝1／6（i＝1，2，3，4）；L为相变潜热；C为热容；γfl为液相分数。gieq确定为

热扩散系数和宏观温度分别为

将总焓H分为显焓和潜热2部分

通过上式确定H后，计算时间步t的液相分数为

式中：Hs和Hl分别为固相线和液相线的焓值，Tm为融化温度。

1.3 简化与假设

本文的模拟主要基于以下几个条件：1）融化后的液体为牛顿流体且不可压；2）流体流动基于Boussinesq假设；3）材料融化前后的热物性相同且保持常数。

2 模型验证

2.1 导热融化

首先选取纯物质的导热融化问题验证模型及程序。对于此问题，Neumann得到了半无限大空间的解析解，温度分布和固液界面的位置分别为

式中：T1为过热温度，k为导热系数。λ的值由以下超越方程确定：

式中：St＝ρc（T1-Tm）／L为Stefan数。

初始固相温度为相变温度Tm，左壁面温度突然升高到T1，开始融化。在模拟中，上下壁面为周期性边界，网格数取为100×10。此处，Stefan数取为0.013，对应的λ为0.25。图1和图2分别为融化前沿位置（固液界面）随时间的变化和不同时刻温度的分布曲线，可以看到，模拟结果与解析解吻合较好。

图1 融化前沿位置随时间的变化Fig.1 The evolution of melting front position with time

图2 不同时刻温度的分布Fig.2 The temperature distribution at different mon⁃ents

2.2 对流融化

为了进一步验证模型，本节对封闭方腔内纯物质的自然对流融化进行模拟。初始固相温度为相变温度Tm，左壁面温度突然升高到T1，并保持不变；右壁面仍然保持在相变温度Tm，上下壁面为绝热壁面。对于速度场，4个壁面均采用二阶精度的Half-Way反弹格式处理无滑移边界。长宽均为l，网格数为200×200。对于对流融化问题，由几个基本的无量纲数确定融化特性：Rayleigh数Ra、Prandtl数Pr和Ste⁃fan数。Rayleigh数和Prandtl数分别定义为

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

Nusselt数定义为

式中：θ＝T-Tm／T1-Tm为无量纲温度。

Mencinger［13］采用自适应网格研究了自然对流融化过程，下面将本文结果与Mencinger的结论进行对比。图3为总液相分数随Fourier数的变化，可以看到，本文结果与Mencinger的结果吻合很好。图4给出了左壁面（热壁面）平均Nusselt数随Fou⁃rier数的变化，本文结果与Mencinger结果差别不大，其差异可能是由数值方法不同引起的。通过以上算例的验证，证明了模型的正确性。

图3 总液相分数随Fourier数的变化Fig.3 The evolution of the total liquid fraction with the Fourier number

图4 左壁面平均Nusselt数随Fourier数的变化Fig.4 The evolution of the average Nusselt number a⁃long the left wall with the Fourier number

3 方腔内多孔介质的融化

上节计算了纯物质的导热融化和自然对流融化，本节对封闭方腔内多孔介质的融化进行模拟。模拟所用多孔介质结构如图5所示，图中黑色区域为固体骨架，白色区域为孔隙。初始及边界条件与2.2节纯物质的自然对流融化相同，不再赘述。网格数为300× 300，Prandtl数和Stefan数均取为1.0，Rayleigh数在1.0×105～1.0×108取4个不同的值进行讨论。图6为多孔介质融化界面随Fourier数的演化，与纯物质自然对流融化类似，初始导热作用占主导，融化界面基本与垂直壁面平行，随着融化的进行，由于对流作用的增强，多孔介质上部融化快于下部，且界面呈现不规则形状，这是由多孔介质的不均匀性引起的。对比图6（a）和（b）可以看到，Rayleigh数越大（对流作用越强），界面形状越不规则。在Ra＝1.0×108（见图6（b）），Fo＝0.003时，由于介质存在一个孔隙连通区域，使流动超过了融化前沿位置，但是在Rayleigh数较小时（见图6（a）），未出现此现象。从以上分析可知，Rayleigh数对多孔介质的融化过程影响较大。这里也给出了方腔内多孔介质导热融化过程，见图6（c），可以看到，多孔介质与纯物质的导热融化也存在较大差别，在纯物质导热融化中，融化界面始终与垂直壁面平行，但由于多孔介质的不均匀性，导致多孔介质导热融化过程中融化前沿的上、下不再同步，当融化界面接近右壁面时，下部先融化完（这取决于具体的多孔介质结构）。数的变化，为了说明对流传热对融化的影响，图7也给出了相应的导热融化过程，以便比较。初始融化过程靠导热控制，不同Rayleigh数下液相分数线重合，随液相分数的增大，对流作用增强，Rayleigh数越大，对流区出现的越早，且融化越快。在对流主导区，融化分数随Fourier数基本呈线性增加，当融化界面接近右壁面时，融化速率减慢。

图5 多孔介质结构Fig.5 The structure of porous media

图6 多孔介质融化界面随Fourier数的变化Fig.6 The evolution of melting interface of porous media with time

图7 不同Rayleigh数下总液相分数随Fourier数的变化Fig.7 The evolution of total liquid fraction with the Fourier number at different Rayleigh numbers

图8 不同Rayleigh数下Nusselt数随Fourier数的变化Fig.8 The evolution of average Nusselt number with the Fourier number at different Rayleigh numbers

不同Rayleigh数下热壁面平均Nusselt数随Fourier数的变化趋势见图8，可以看到，由于开始融化时，导热为主要传热机制，使Nusselt数急剧减小，随着对流作用的增强，限制了Nusselt数的减小，使其达到一个最小值，然后趋于稳定。Rayleigh数越大，稳定后的Nusselt数越大。图中还可以看到，当Ra＝1.0×108时，Nusselt数存在一个峰值，对比图6（b）Fo＝0.003，这是因为多孔介质为非均匀的，此时流动超过了融化前沿位置，在流动后面未融化的骨架区会对流动产生扰动作用，从而增强对流，使Nusselt数上升，当此部分固体融化后，Nusselt数又有所下降。Gobin等［14］研究了Prandtl数对纯物质自然对流融化的影响（0.01＜Pr＜0.1），认为Prandtl数对融化有较大影响。本文考虑大范围Prandtl数对多孔介质自然对流融化的影响，计算中，Rayleigh数和Stefan数分别取为1.0×105为和1.0，而Prandtl数取0.01、0.1、1.0和10.0 4组做对比。图9和图10分别为Prandtl数对总液相分数和热壁面Nusselt数的影响，可以看到，在低Prandtl数范围内，Prandtl数对多孔介质的融化影响较大，但是Pr＝1.0和Pr＝10.0得到的总液相分数和Nusselt数曲线差别很小，也就是说，当Pr＞1.0时，继续增大Prandtl数对融化（传热）基本没有影响。在开始融化阶段，主要为纯导热融化，不同Prandtl数下液相分数曲线重合（见图9），一旦液相分数达到一定比例，对流开始起作用，在低Prandtl数范围内，Prandtl数越大（相对动量扩散率越大，或者相对热扩散率越小），对流越强，融化更快，且达到稳定时的Nusselt数也更大（见图10）。

图9 不同Prandtl数下总液相分数随Fourier数的变化Fig.9 The evolution of total liquid fraction with the Fourier number at different Prandtl numbers

图10 不同Prandtl数下Nusselt数随Fourier数的变化Fig.10 The evolution of average Nusselt number with the Fourier number at different Prandtl numbers

4 结论

1）通过验证，本文采用的多松弛模型能够很好的预测导热和自然对流融化；

2）方腔内多孔介质的导热融化界面不再与垂直壁面平行，自然对流融化出现不规则形状的固液界面；

3）增大Rayleigh数，多孔介质融化加快，融化界面更不规则，热壁面Nusselt数增大；在低Prandtl数范围内，Prandtl数越大，融化越快，Nusselt数越大，但到一定值后，继续增大Prandtl数基本不再影响融化过程。

［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］朱卫兵，王猛，陈宏，等.多孔介质内流体流动的格子Boltzmann模拟［J］.化工学报，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.网络出版时间：2015⁃04⁃26.

高超声速冲压发动机技术重点实验室开放基金资助项目（20140103008）；哈尔滨市科技创新人才研究专项资助项目（2014RFXXJ062）.

王猛（1988⁃），男，博士研究生；朱卫兵（1961⁃），男，教授，博士生导师.

朱卫兵，E⁃mail：zhuweibing＠hrbeu.edu.cn.