基于格子玻尔兹曼方法的多孔材料导热模拟

Simulation of Thermal Conductivity of Porous Materials Based on Lattice Boltzmann Method

摘要: 多孔材料作为一种新型材料，因其轻质、高导热等特性应用于多个领域。对于多孔材料进行数值重构，并研究其传热特性对工程应用具有重要意义。本文根据多孔材料内部的复杂结构进行建模，应用格子玻尔兹曼方法计算多孔材料的等效导热系数，进行算法验证，并分析网格数目和两相导热系数比对等效导热系数计算结果的影响。其中，等效导热系数模拟结果会在一定网格区间内趋于准确，并且两相导热系数比越大，等效导热系数变化越平稳，趋于定值。

Abstract: As a new kind of material, porous material is applied in many fields because of its light weight and high thermal conductivity. The numerical reconstruction of porous materials and the study of their heat transfer characteristics are of great significance for engineering applications. In this paper, the complex structure of porous material is modeled, and the lattice Ludwig Boltzmann method is used to calculate the effective thermal conductivity of porous material, validate the algorithm, and analyze the influence of mesh number and two-phase thermal conductivity ratio on the calculation results of equivalent thermal conductivity. The simulation results of equivalent thermal conductivity tend to be accurate in a certain grid range, and the higher the two-phase thermal conductivity ratio, the more stable the change of equivalent thermal conductivity, tending to a fixed value.

1. 引言

核能是未来太空开发和探索的重要能源之一。以核能为核心的空间核动力系统，主要由核反应堆、热电转换系统、辐射屏蔽系统、热排放系统、控制系统、任务模块和推进装置几大系统组成[1]。其中，热排放系统由散热器和辅助设备组成，目前使用的散热器主要采用金属材料制成，如果采用具有高导热性的多孔材料，更有利于减小部件尺寸、减轻设备重量，实现结构紧凑化、运行高效化[2]。多孔材料作为一种新型材料，其热传导过程十分复杂，包括实体骨架内部的热传导、实体骨架与填充相之间的对流换热以及实体骨架之间的辐射换热。对于不考虑对流、辐射的多孔材料导热研究，Calmidi和Mahajan [3]分别用空气和水作为填充相，测量了能源研究发电公司所生产泡沫铝的有效导热系数，并根据泡沫铝自身的六角形结构推导了相关的二维解析模型，与实验结果吻合良好。之后，在Calmidi和Mahajan研究工作的基础上，Boomsma和Poulikakos [4]将金属多孔材料结构用圆柱形状的韧带表示，基于金属多孔材料理想的三维基本单元几何结构建立相关的等效导热系数模型，并与实验数据进行对比，吻合程度较为理想。美国橡树岭国家实验室所属的Klett等[5]采用光学图像分析对石墨泡沫的结构进行了分析，建立了一个简单的双参数模型，根据中间相沥青基石墨泡沫有效导热系数实测结果，拟合得到石墨泡沫有效导热系数的半经验公式，可以应用于具有相似加工参数的石墨泡沫。

国内张赛等[6]研究多孔材料时根据其物理组成，将其分为具有分形结构的团聚体集合和具有非分形特性的固相孔相，建立了简化的单元模型。李仁民等[7]基于LBM (Lattice Boltzmann Method，格子玻尔兹曼方法)采用四数法对黏土微观结构进行了理论分析。马强等[8]应用LBM来计算多孔介质内的传热过程，讨论了不同的分形特性对多孔介质导热过程的影响。

采用传统实验研究的方法了解多孔材料的传热特性，所需成本高[9]、实验周期长且无法准确获得多孔材料内部详细的温度场分布。除此之外，多孔材料的传热特性不仅与固体骨架及填充相自身的热物理特性有关，也与材料内部复杂的微观结构有关，而实验所得的结果难以适用于不同种类的多孔材料。由此来看，基于数值分析的多孔材料温度场建模仿真分析具有显著优势。而对此的研究大多采用近似模型或理想模型的建模方法，本文基于随机四参数法，以及根据多孔材料的特点、性质以及结构参数，编写程序进行多孔材料的二维、三维数值重构。在此基础上，本文进一步完成模型优化，并应用格子玻尔兹曼方法计算多孔材料的等效导热系数[10]。

2. 重构模型

本文基于介观统计的数值方法，研究多孔材料的传热特性，采用的多孔材料重构方法便是随机四参数法。对于随机四参数法来说，多孔材料建模时所需要关注的孔隙率、连通性、孔隙尺寸等结构参数都可以通过调整骨架或孔隙的生长参数来控制，随机生成具有不同孔隙形状的多孔材料，用于后续的研究设计。

在整个计算域中随机分布生长核，以各个生长核为起点，按给定的与各向同性/各向异性相关的概率向各个方向生长，等到符合设定的孔隙率ε，停止生长，步骤如下[11]：1) 计算域内随机分布初始生长核；2) 根据材料的各向同性/各向异性确定八个方向(见图1)的生长概率；3) 遍历每个生长核，通过各方向生长概率进行生长；4) 重复步骤3)，直至模型孔隙率达到设定值。

Figure 1. Schematic diagram of the growth direction of the growth phase

图1. 生长相生长方向示意图

在重构过程中，可以自由选择固相或者孔相作为生长相，另一相即可自动填充模型所占计算域。考虑到骨架的连续性以及孔隙的连通性，孔隙率低于0.5时用孔相作为生长相，在整个固相中随机生长；孔隙率高于0.5时，将固相作为生长相，不过需要注意初始生长核的分布，尽量在整个模型区域内均匀分散，且数量足够，不然很容易出现大量固相孤岛，影响后续程序计算稳定性。因此，本文将孔隙率小于0.5的多孔材料作为研究对象，对于孔隙率大于0.5的多孔材料不作分析。当然，如果忽视多孔材料结构的合理性，仅考虑对其传热性质影响较大的孔隙率、各向同性、连通性等性质，孔隙率大于0.5的多孔材料可看作是孔隙率小于0.5的材料的反相。

在初始生长核概率和各个方向生长概率给定条件下，可以得到孔隙率为0.4的各向同性多孔材料的二维建模如图2所示，三维建模如图3所示，其中黑色代表生长的固体材料，白色代表孔洞，横纵坐标代表网格数量。可以看出，图2中出现了固相“孤岛”的现象，“孤岛”具体指生成的固相在整体中属于独立存在，在实际材料中是不存在的。而三维模型实际上是二维模型堆叠而成，增加了一个生长坐标轴，基本上不会出现“孤岛”现象，更接近坐标轴，基本上不会出现“孤岛”现象，更接近实际模型。但考虑到二维模型生成时间远低于三维模型生成时间，从研究的效率上来考虑，使用二维模型进行后续的研究分析。

Figure 2. Two-dimensional model of porous materials

图2. 多孔材料二维模型

Figure 3. Three-dimensional model of porous materials

图3. 多孔材料三维模型

二维建模快速高效，但会产生很多微小甚至单网格的固体相“孤岛”，影响后续计算程序收敛性，考虑优化重构方法，过滤微小特别是单网格的固相颗粒。

本文中进行模型优化时选择过滤单个网格的固相颗粒，对于团聚在一起的固相不作处理。对于孔隙率为0.4的模型(图2)优化效果如图4所示。优化模型时需要检测孔隙率是否达标，可以通过“生成–优化”的循环来解决问题。较为节省计算空间的方式是在生成模型时将设定的孔隙率减去一个小数(大约是孔隙率的十分之一)，为优化模型提供裕量。

Figure 4. Optimization of porous material model

图4. 多孔材料模型优化

3. 数学模型

如果只考虑导热，不考虑对流及辐射，可将能量方程离散为[12]：

$f_i\left(r+c_i\Delta t,t+\Delta t\right)=f_i\left(r,t\right)+\frac{\Delta t}{\tau }\left[f_i^{eq}\left(r,t\right)-f_i\left(r,t\right)\right]$ (1)

式中，f——分布函数；τ——无量纲松弛因子；i——格子速度方向；$f_i^{eq}$ ——平衡分布函数；c——离散速度；r——离散位置；t——时间，Δt——时间步长，由温度求得：

$f_i^{eq}=T{w}_i$ (2)

式中，T——宏观温度；$w_i$ ——分布函数在格子速度方向所占权重，由格子玻尔兹曼离散模型所确定，满足${\displaystyle {\sum }_i{w}_i}=1$ 。

应用D2Q9模型(见图5)进行导热计算，考虑了多孔材料两相的共轭热[13]效应，f₀、f₁、f₂、f₃、f₄、f₅、f₆、f₇、f₈的权重因子$w_i$ 分别为0、1/6、1/6、1/6、1/6、1/12、1/12、1/12、1/12，对应的速度为c₀(0, 0)、c₁(1, 0)、c₂(−1, 0)、c₃(0, 1)、c₄(0, −1)、c₅(1, 1)、c₆(−1, 1)、c₇(−1, −1)、c₈(1, −1)，选取格子数据∆x = 1、∆t = 1，可将能量方程分为碰撞步骤：

$f_i\left(x,t+\Delta t\right)=f_i\left(x,t\right)+\frac{\Delta t}{{\tau }_n}\left[f_i^{eq}\left(x,t\right)-f_i\left(x,t\right)\right]$ (3)

和迁移步骤：

$f_i\left(x+\Delta x,t\right)=f_i\left(x,t\right)$ (4)

Figure 5. D2Q9 model

图5. D2Q9模型

碰撞步骤中的松弛因子[13] [14]求解为：

${\tau }_n=\frac{3}{2}\frac{k_n}{{\left(\rho c_p\right)}_n{c}^2}+\frac{1}{2}$ (5)

式中，n——多孔材料固相/孔相；k——各相导热系数；$\rho c_p$ ——各相体积热容，为保证界面温度和热流连续性一般将其设为1，不影响等效导热系数计算结果；c——伪声速，理论上可为任意正值，主要用于控制程序收敛性，需要根据实际应用调整，尽量使τ值在0.5~2范围内。

由此，节点热流及等效导热系数求解为：

$T={\displaystyle \sum _i{w}_i{f}_i}$ (6)

$q=\left({\displaystyle \sum _i{c}_i{f}_i}\right)\frac{{\tau }_n-0.5}{{\tau }_n}$ (7)

等温边界采用非平衡外推方程[15]，即：

$f_i\left(r,t\right)-f_i^{eq}\left(r,t\right)=-\left[f_{i^{\prime }}\left(r,t\right)-f_{i^{\prime }}^{eq}\left(r,t\right)\right]$ (8)

式中，i、i'——格子链上的相反方向，如图5所示的1和3、2和4、5和7、6和8。

绝热边界的温度梯度为0，参考文献[15]中认为粒子分布函数控制热流方向，为了保证绝热边界热流不外溢，提供绝热边界条件：

$f_i\left(m,t\right)=f_{i^{\prime }}\left(m,t\right)$ (9)

式中，m——边界节点。

4. 结果分析

计算模型如图6所示。

4.1. 算法验证

对于导热程序验证，将程序运行结果与参考文献[16]的实验数据进行对比验证。实验采用Cu/solder复合材料，其中Cu颗粒均匀地分散在solder中，由于两相皆为固体，不用考虑对流情况。利用随机四参数法再现该两相颗粒复合材料的二维微观结构，solder作为非生长相，Cu颗粒为生长相，该结构被认为是各向同性的[17] [18]。各组分导热系数分别为k_Cu = 398.0 W/(m∙K)，k_solder = 78.1 W/(m∙K)。图7显示了数值模拟的等效导热系数，作为Cu体积分数ε (可视作孔隙率)的函数，并进行了比较。拟合结果与实际数据符合良好，证明程序的可靠性。

Figure 6. Calculation model

图6. 计算模型

Figure 7. Comparison of Cu/solder material simulation and test data

图7. Cu/solder材料的模拟与实验数据之间的比较

事实上，由于随机四参数法生成多孔材料结构具有随机性，实验数据有一定偏差。为了尽量降低随机性的影响，图7中模拟结果采用十个随机生成多孔材料结果数值的平均值。

4.2. 网格数的影响

为了尽量降低随机性的影响，图8中数值模拟结果采用十个随机生成多孔材料结果数值的平均值。图8中表示出不同导热比条件下，选取不同网格数对模拟结果的影响。可以明显看出，网格数目在3600~14,400范围，等效导热系数趋于平稳；低于3600网格数目，等效导热系数有较为明显的变化；在14,400~62,500范围，等效导热系数曲线有一定的起伏波动，但两相导热比为$1:\frac{1}{2}$ 的等效导热系数曲线变化细微，这是因为两相导热性质接近，整个计算域数据稳定，但是两相导热比的增加会影响程序稳定性，在数值模拟中需要注意这一点。

4.3. 孔隙率的影响

考虑到适用性更全面的模拟研究，后续研究采用固相孔相导热系数比的形式进行研究分析。为了尽量降低随机性的影响，图9中数值模拟结果采用十个随机生成多孔材料结果数值的平均值。图9中表示出不同孔隙率对模拟结果的影响。对各项数据进行拟合，得到拟合关系式：

Figure 8. Effect of grid number on simulation results

图8. 网格数对模拟结果的影响

Figure 9. Effect of porosity on simulation results

图9. 孔隙率对模拟结果的影响

$\{\begin{array}{l}{k}_{\text{eff}}=-\text{0}\text{.51094}\epsilon +\text{0}\text{.95605},两相导热比1:\frac{1}{2}\\ k_{\text{eff}}=-\text{1}\text{.16566}\epsilon +\text{0}\text{.8509},两相导热比\text{1}:\frac{\text{1}}{10}\\ k_{\text{eff}}=-\text{1}\text{.50674}\epsilon +\text{0}\text{.7784},两相导热比\text{1}:\frac{\text{1}}{100}\\ k_{\text{eff}}=-\text{1}\text{.60183}\epsilon +\text{0}\text{.76957},两相导热比\text{1}:\frac{\text{1}}{500}\\ k_{\text{eff}}=-\text{1}\text{.62186}\epsilon +\text{0}\text{.76876},两相导热比\text{1}:\frac{\text{1}}{1000}\end{array}$ (10)

可以明显看出，随着孔隙率的增加，等效导热系数呈线性下降趋势，这是由于孔相所选取的导热系数小于固相所选取的导热系数，孔相占比的增加会使整个多孔材料的等效导热系数减小。孔隙率的增加会出现更多的孔相，增加了多孔材料两相交界面，降低了程序稳定性。

4.4. 两相导热比的影响

为了尽量降低随机性的影响，图10中数值模拟结果采用十个随机生成多孔材料结果数值的平均值。图10中表示在不同孔隙率条件下，不同两相导热比对模拟结果的影响。可以明显看出，随着两相导热比的增加，等效导热系数也在逐渐降低，并趋于平缓，这是由于固相所选取的导热系数系数大于孔相所选取的导热系数，不同孔隙率条件下，孔相导热系数的减少很容易导致整个多孔材料等效导热系数的减少，

注：(1) 两相导热比为$1:\frac{1}{2}$ ；(2) 两相导热比为$1:\frac{1}{10}$ ；(3) 两相导热比为$1:\frac{1}{100}$ ；(4) 两相导热比为$1:\frac{1}{500}$ ；(5) 两相导热比为$1:\frac{1}{1000}$ ；(6) 两相导热比为$1:\frac{1}{1500}$ ；(7) 两相导热比为$1:\frac{1}{2000}$ 。

Figure 10. Effect of two-phase thermal conductivity ratio on simulation results

图10. 两相导热比对模拟结果的影响

并且在固相导热系数与孔相导热系数比值达到$1:\frac{1}{500}$ 时，孔相导热系数降低对整个多孔材料等效导热系数的影响已经微乎其微了，此时可认为孔相中无填充相，处于真空环境。

5. 结论

本文介绍了格子玻尔兹曼方法在多孔材料导热计算中的应用，并且分析了计算域网格数目、多孔材料两相导热比等对计算多孔材料等效导热系数的影响，并得出结论：

1) 计算域网格数目对模拟结果的影响不是网格数目越大，结果越精确，而是在3600~14,400范围内，结果较为稳定，当网格数量较少时，虽然会加快计算速度，但由于数量较少，尺度过大，计算结果偏小。当网格数量较多时，不仅影响计算速度，还会使计算结果出现偏差，尺度过小会引入更多不稳定性，考虑网格质量和网格数目的影响，选取100 × 100网格进行计算。

2) 通过生成不同孔隙率两相多孔介质材料，计算其等效导热系数来观察两相导热比对等效导热系数的影响。多孔材料两相导热比对等效导热系数的影响主要在于，两相导热比越大，整体的等效导热系数变化越平缓，趋于定值。其中，等效导热系数随着导热比越来越大，其受孔隙率的影响越来越大。

NOTES

^*通讯作者。