基于Lattice－Boltzmann方法的SiO2多孔绝热材料传热分析

孙 泽，杨自春，苏高辉

(海军工程大学 船舶动力学院，湖北 武汉430044)

近几年，随着新型多孔介质绝热材料的不断涌现，与此相关的微观传热机理分析越来越受到研究人员的重视，希望能从传热机理上为材料的设计和制备提供指导［1－2］。微米、纳米尺度的传热问题使得传统的宏观模型不再有效，目前比较适合于分析微尺度传热问题的主要方法有如下几类［3－4］:Boltzmann方程方法、分子动力学方法、直接Monte－Carlo模拟方法和量子分子动力学方法等。Boltzmann方程方法被普遍认为是现有方法中用来分析微尺度能量输运现象的最具有普遍适用性和强有力的工具。几乎所有的宏观输运方程如Fourier定理、Ohm定理、Fick定理及双曲线型热传导方程，均可由该方程导出［5］。Lattice－Boltzmann(LBM)方法是基于Boltzmann输运方程得到的，对于微纳米尺度材料的输运性质如热导率、电导率等能进行很好的模拟。此外LBM方法还有如下优点:

(1)使用范围更广，易应用到各种复杂情形中;

(2)适用于大规模并行计算;

(3)边界条件容易设定;

(4)与传统的数值方法相比，达到相同的精度所需要的计算网格数及计算时间明显减少。

由于LBM方法所具有的独特的优点，使得该方法越来越受到国内外学者的重视。特别是，在最近的几年里，LBM方法已经被国内外研究者用来研究微尺度流动和传热问题。

微观研究的另一个主要困难是如何获取多孔介质的孔隙结构数据。对于一般的多孔介质可以通过数字扫描的方法得到灰度图像，并用图像处理方法获得孔隙形貌。但这种方法不但成本较高，而且孔隙结构的精度依赖仪器的分辨率。因此理论研究中人们经常使用人工构造的多孔介质模型。目前，很多学者提出了多种能够反映多孔介质统计信息的方法，如随机配置法、球体沉降法、硬球Monte－Carlo方法、分数布朗运动方法［6－8］。Wang［9］等采用与LBM方法紧密结合的四参数随机生成方法，该方法结合多孔介质构造的孔隙生长模型和丛生理论，构造过程更接近于多孔介质的生成过程，能够通过参数调整控制生成介质的形貌特征。

本文基于多孔介质四参数随机生成法的基本原理，构造SiO2多孔绝热材料微观结构形态，编制相应的程序，采用LBM方法分析构造的SiO2多孔绝热材料的热学性能。

1 Lattice－Boltzmann导热模型

1.1 D2Q9不可压格子多相LBM热模型

本文采用二维九速度(D2Q9)不可压格子LBM热模型，在多孔材料中无相变无热源无对流的每一相的温度演化方程表示为［10］

式中r——坐标矢量;

t——时间;

n——第n相;

δt——时间步长;

ei(i=0，1，…8)——离散速度矢量(参见图1);

图1 D2Q9模型速度矢量示意图

gi(r，t)——与i方向相对应的粒子温度分布函数;

离散速度矢量ei为

无量纲松弛时间为

其中k为导热系数，c=δx/δt，设定时需将 τn的值控制在0.5～2.0。

宏观参数温度及热流密度由以下公式给出

导热系数为

1.2 边界条件的处理

LBM方法中，边界条件起着十分重要的作用，其处理方法会对数值计算的精度、计算的稳定性以及计算效率产生很大的影响。这里采用通用性较好的由Guo［11］等人提出的非平衡外推方法，其基本思想是，将边界节点上的分布函数分解为平衡态和非平衡态两部分，其中，平衡态部分由边界条件的定义近似获得，而非平衡态部分则用非平衡外推确定。

平衡态部分

其中下标b，f分别表示边界上的格点和边界相邻的格点。

非平衡态部分

(1)恒温边界(温度已知)

边界上的温度分布函数为

(2)定热流边界(温度梯度已知)

边界上的温度分布函数为

其中

特殊情况，对应绝热边界即热流密度q=0，

2 四参数随机生成法生成多孔介质的微观结构

Wang［9］等提出的四参数随机生成法(quartet structure generation set，QSGS)可以通过参数调整控制生成介质的形貌特征。QSGS首先假设系统中包含一个非生长相和若干(设为N)个生长相，其生长步骤如下:

(1)在构造网络上以概率pc随机布置第一个生长相的生成核，pc要小于该相的体积分数;

(2)对于每一个生长核，以概率pi向第i个方向生长(i=1～8)，即第i个方向的相邻单元的状态以概率pi变为第一相态;

(3)重复(2)直到该相达到预设的体积分数;

(4)对第n＞1个生长相，其生长过程与相之间的相互作用有关。如果该相与其他已生成的相无关，则按照上述步骤进行生长;否则，该相以概率

(5)重复(4)直至该相达到其体积分数;

(6)重复(4)～(5)直至n=N;

(7)最后，网络中没有被生长相占据的为非生长相单元。

显然SiO2多孔介质微观结构由孔隙和固体颗粒两相构成，即N=2，令固体颗粒生长相，孔隙为非生长相，初始相全为孔隙。也就是只需要进行(1)～(3)步骤直至非生长相(孔隙)达到给定的孔隙率。

图2 二维节点8个生长方向

3 程序实现与算例验证

本文的主要编程思路为利用Matlab编程生成多孔介质图像已及相函数矩阵，并将相函数矩阵导入C++编的LBM算法中作为固液相判定矩阵，进而对多孔介质进行LBM数值传热分析。

3.1 编程实现

3.1.1 生成结构的编程实现

根据QSGS的基本思想和实现步骤，采用Matlab编制结构生成程序，程序的流程图如图3所示，整个计算域分为网络，网络中固相用1表示，液相用0表示，通过Matlab的散点图可以生成多孔介质图像，并导出相函数矩阵。

图3 QSCS流程图

从而为LBM传热分析提供相分布的数值依据。

3.1.2 多相LBM算法的程序流程

LBM在运行前需要先读取QSGS生成的相函数矩阵，赋值给程序中的相判定变量。在演化过程中先判定为固相或孔隙，再进行相应的演化。LBM算法程序的流程图如图4所示。

图4 LBM流程图

3.2 算列验证

为了验证模型的正确性，以及所编程序的可行性，利用D2Q9不可压格子多相LBM热模型计算二维两种不同材料的串联结构的温度分布(见图5)。

图5 串联传热结构模型

两种材料的导热系数分别为k0、k1，对于如图所示的串联结构，其稳态温度场的解析解为

算例中取L=6 mm，k0=6 W/m·K，k1=236 W/m·K，格子数取100×100。用D2Q9模型和解析方法计算得到y轴中点的温度分别2.474 5 K和2.479 3 K，误差仅为0.2%。图6给出了D2Q9模型计算的串联结构沿y轴方向上的分布曲线与解析解比较的曲线，可以看出计算的结果与解析解吻合的比较好，说明用该模型计算的温度场有着相当高的精度。

图6 LBM计算值与解析解的比较

4 二维多孔绝热材料的多相导热过程的模拟

SiO2多孔绝热材料具有优异的绝热性能，是一种新型轻质保温隔热的理想材料［2］。

对于多孔介质内部传热的研究，首先需要构造出相应的多孔结构。按照前文所提到的QSGS方法，构造SiO2多孔绝热材料的结构特征，取生长核构造概率分别为pc=0.1、pc=0.1，1～4方向生长概率p1=p2=p3=p4=0.1，5～8方向生长概率p5=p6=p7=p8=0.025，孔隙率 ε=0.85。构造的二维结构分别如图7所示。

图7 QSGS法构造的多孔介质

本文算列中，样本材料为SiO2多孔绝热材料。骨架为SiO2，孔隙为空气。物性参数如表1所示。综合考虑计算机的运算速度和运算时间，计算区域采用边长为5×10－4m，划分为200×200的网络，格子长度d x=2.5×10－6m，时间步长为d t=2.5×10－7s。

表1 SiO2多孔绝热材料物性参数

边界条件取上下定温，上边界为285 K，下边界为275 K，左右两边界为绝热。

图8和图9分别表示两种不同结构(pc分别等于0.1和0.01)的多孔材料达到稳态传热后(演化50 000步)的温度分布图，图10为pc=0.1的结构演化过程中(演化5 000步)的温度分布图。对照图8，可明显看出模拟是按照真实的温度传递演化的。

图8 p c=0.1生成的多孔材料稳态温度场

图9 p c=0.01生成的多孔材料稳态温度场

图10 p c=0.1，演化5 000步生成的多孔材料温度场

对照结构图，明显得到该多孔材料能量主要沿骨架向下传递，高温的空气再向周围的孔隙传递能量，最后达到平衡。比较图8与图9的等温线，图8的比图9的更平直，因为虽然都是采用随机生成法生成的结构，但是前者pc大，固相生长核更均匀些，相应的材料结构也更均匀，局部导热系数的差异较小，而后者局部导热系数差异较大，但整体两侧温差很小，所以当选择更大的计算区域时，整体表现为各向同性，所以该随机生成的结构的坐标轴方向仍然是其导热系数主轴方向。

图11是两种不同结构在y方向上的有效导热系数随孔隙率的变化关系。

图11 有效导热系数与孔隙率的关系

计算结果显示，两种结构的斜率偏差明显。所以骨架结构对有效导热系数的影响明显。pc越小的结构有效导热系数也明显偏低，由图7可得在相同孔隙率的情况下，pc越小的结构平均孔径越大，也就是说在相同孔隙率下孔径越大，有效导热系数越低。

为了讨论以上结论，我们可以观察图8的结构。当上下边界为等温边界且材料为均匀材料时，等温线必然都为相互平行的直线，而图8中pc偏大的结构的等温线更接近于均匀材料的等温线(此时孔径偏小)。那么根据温度梯度矢量统一原则［12］，便可得到在孔隙率相同的情况下，孔径越小该多孔材料有效系数越高。当然这是在不考虑传热的尺度效应的前提下，而以上算列中Knudsen数小于10－1，属于分子自由流区域，并不需要考虑传热的尺度效应［13］。这更验证了本算列的正确性。

随着孔隙率的减小，有效导热系数明显增加，也就是说增加气相导热比重能够增强该绝热材料的绝热性能。

5 结论

本文采用四参数随机生成法构造了SiO2多孔绝热材料微观结构，该方法较于传统构造多孔介质结构，更真实的反应了材料结构的随机性。引入D2Q9不可压格子多相LBM热模型，该模型能够方便的计算格子微尺度并具有复杂边界或多相材料中的导热问题。给出了从结构的构造到具体的LBM传热分析的编程实现过程，以及程序流程图。对于生成的不同结构的SiO2多孔绝热材料进行了传热的计算和分析。

对于SiO2多孔绝热材料中的导热过程模拟结果表明，在Knudsen数小于10－1时，相同孔隙率下孔径越大，有效导热系数越低。随着孔隙率的减小，有效导热系数明显增加，增加气相导热比重是增强该绝热材料绝热性能的有效途径，骨架结构对有效导热系数的影响显著。

［1］张贺新，赫晓东，何飞.气凝胶隔热性能及复合气凝胶隔热材料研究进展［J］.材料工程，2007(s1):94－97.

［2］杨自春，陈德平.SiO2纳米多孔绝热材料的制备与绝热性能研究［J］.硅酸盐学报，2009，37(10):1740－1743.

［3］刘静.微米纳米尺度传热学［M］.北京:科学出版社，2001.

［4］Z.M.Zhang，Nano/Microscale heat transfer［M］.New York:The McGraw－Hill，2007.

［5］C.L.Tien，A.Majamdar，F.M.Gemer，Microscale energy transport［M］.Washington:Taylor＆Francis，1998.

［6］M.Madadi and M.Sahimi.Lattice Boltzmann simulation of flow in fracture networks with rough，self－affline surfaes［J］.Phys.Rev.E，6302－6309.

［7］H.F.Zhang，X.S.Ge，and H.Ye.Randomly mixed model for predicting the effective thermal conductivity of moist porous media［J］.J.Phys.D，2006(39):220－226.

［8］R.S.Maiera，D.M.Kroll，R.S.Bnard et al.Porescale simulation of dispersion［J］.Phys.Fluids，2000(12):2062－2079.

［9］M.Wang，J.Wang，N.Pan，S.Y.Chen，Mesoscopic predictions of the effective thermal conductivity for microscale random porous media［J］.Physical Review E，2007(311):562－570.

［10］J.Wang，M.Wang，Z.Li，A lattice Boltzmann algorithm for fluid－solid conjugate heat transfer［J］.International Journal of Thermal Sciences，doi:10.1016/j.ijthermalsci，2007(46):228－234.

［11］Guo Z L，Shi B C，Zheng C G.A coupled lattice BGK model for Boussinesq equations［J］.Int J Numer Methods Fluids，2002，39(4):325－342.

［12］X.G Cheng，ZX Li，ZY Guo，Constructs of the highly effective heat transport paths by bionic optimization［J］.Sciences in China(Series E)，2003，43(3):296－302.

［13］郑林.微尺度流动与传热传质的格子Boltzmann方法［D］.武汉:华中科技大学，2010.