基于格子波尔兹曼方法的回热器数值模拟

夏宇栋 陈 曦 马诗旻 张 华

(上海理工大学能源与动力学院 上海 200093)

基于格子波尔兹曼方法的回热器数值模拟

夏宇栋 陈 曦 马诗旻 张 华

(上海理工大学能源与动力学院 上海 200093)

利用格子玻尔兹曼方法，直接对蚀刻薄片和层叠丝网回热器的微观结构流场进行了模拟。得到了两种回热器填料的微观流场和两端的压差。模拟结果显示，当回热器的直径、水力直径和填充率相近情况下，不同流速下蚀刻薄片卷裹式回热器的稳态阻力系数均比层叠丝网回热器小。稳态阻力系数的模拟变化趋势与实验一致。

格子波尔兹曼方法 回热器 流阻系数 数值模拟

1 引言

回热器是低温制冷机的关键部件，也是影响制冷机性能的最重要因素。目前，对回热器的研究主要有实验研究、传统计算方法研究和模拟计算方法研究。随着计算机技术的迅速发展，对回热器的模拟计算已经成为了回热器研究领域的一个新的热点。张晓青，郭方中等［1］利用基于网络模型开发了仿真软件对热声驱动的制冷机用回热器的长度、水力半径以及回热器在驻波声场中的位置等参数对其性能的影响进行了模拟分析，并对其进行了优化，为更深地理解回热器的机理和实质提供了数值基础。章忠敏［2］根据豪森的数学模型，利用有限差分法求解对称情况下的理想回热器传热微分方程，并编制了计算稳流回热器效率的通用程序，分析了理想稳流回热器的传热特点和近似的数学描述。M.D.Atrey等［3］运用有限差分法重点分析了回热器对整机性能的影响，分析给出了工质质量流率、平均压力、流速、传热周期等对回热器效率的影响，得出回热器效率是丝网孔隙率的函数。

近年来，作为一种崭新的数值模拟手段格子波尔兹曼法也逐渐被应用到对多孔介质的数值模拟当中。格子波尔兹曼方法(Lattice Boltzmann Method(LBM))是一种基于粒子分布函数演化的数值方法。早期格子波尔兹曼方法是从格子气元胞自动机方法发展而来［4］。后来研究表明，通过离散时间、物理空间和速度空间，该方法可以直接从连续波尔兹曼方程严格导出。与传统CFD方法比较，其控制方程的对流项是线性的，流体的压力也无需求解Possion方程得到，同时该方法具有二阶精度和良好的数值稳定性。

本文就将运用格子波尔兹曼方法对层叠丝网和蚀刻薄片卷裹式回热器填料的微观机构进行数值模拟，结合宏观尺度上的回热器数值计算模型，对回热器性能做出预测。

2 回热器物理模型建立

格子玻尔兹曼方法需将计算空间划分为离散的流体微团。因此首先得建立两种回热器的物理模型。

2.1 层叠丝网回热器物理模型建立

层叠丝网是回热器填料主要是采用高目数的不锈钢或磷青铜层叠而成，其结构形式如图1所示。

图1 层叠丝网回热器填料的结构Fig.1 Structure of mesh screen regenerator

其结构特点是径向为规整的编织网片，轴向为网片的随机旋转叠加。本次模拟选用的层叠丝网的详细几何参数见表1所示。

表1 层叠丝网回热器几何参数Table 1 Structure parameters of mesh screen regenerator

根据层叠丝网的几何参数，利用Wolfram Mathematic 8.0软件，通过编写层叠丝网的几何方程得到层叠丝网的物理计算模型，如图2所示，所生成的层叠丝网为7层，每层丝网厚度Lr1=8.5×10-5m。计算空间为0.6 mm×0.28 mm×0.28 mm。并取格子尺寸δx=5×10-6m。则总共的网络计算节点约为588万个。

图2 三维层叠丝网回热器填料结构Fig.2 Three-dimensional model of meshscreen regenerator

2.2 蚀刻薄片卷裹式回热器物理模型的建立

蚀刻薄片卷裹式回热器是由本课题组自主研发的新型回热器，它是在金属薄片表面采用化学腐蚀或光刻的方式形成一定的图案，通过卷裹形成流道。其几何尺寸如图3所示。

图3 蚀刻薄片的几何参数(mm)Fig.3 Geometric parameter of etched foil

同样利用Wolfram Mathematic 8.0软件，生成最终的物理计算模型如图4所示。

该结构的计算空间为0.275 mm×1.4 mm×2.6 mm，取网格尺寸δx=5×10-6m，则计算的网络节点约为800.8万个。

图4 层叠丝网回热器填料结构Fig.4 Three-dimensional model of etched foil regenerator

3 D3Q19不可压等温格子模型用于填料计算

运用开源的CFD框架Palabos可以对于这两种回热器结构进行模拟。Palabos是由日内瓦大学Jonas Latt博士、Orestis博士等成员联合开发的开源的通用CFD计算库，它基于格子玻尔兹曼法。2011年9月，由FlowKit公司开始运营，为用户提供解决方案。它基于C++的标准模板库(STL)，提供了完整的程序库(Library)框架，用户可以在其框架的基础上开发基于格子玻尔兹曼的计算程序。

3.1 回热器计算模型设定

本次计算的回热器的轴向长度仅为0.6 mm和2.6 mm，其两端温差较小，因此忽略温度变化，即模拟的结构为一个等温的模型。所以计算选用等温D3Q19-BGK模型，其恢复的方程与标准的不可压Navier-Stokes方程组一致，适用于低马赫数的流动。

对计算物理模型，定义其两端的速度边界。此处采用了“湿节点”法中的 Regularized边界条件［5］，这种边界条件处理方法的特征是，边界节点中的所有局部分布函数都会被替换，不管它们是否已知。

3.2 运行结果处理及分析

对层叠丝网结构回热器，为了能有效地与实验进行对比，根据实际的回热器工作参数，模拟选择了入口速度为2、5、8、11 m/s 4 个工况点。

3.2.1 流场分布

图5所示为入口速度为5 m/s时，两种填料结构的速度场向量图。从图中可以清晰的看出工质在进入层叠丝网结构前速度较小，当工质进入填料间的流道，由于流道变窄，工质加速，最大速度可达到入口速度的7—10倍大小。

图6为两种结构回热器填料流道内的速度分布向量图，两种回热器填料结构的差异使其流场有很大的不同，进而影响回热器的回热性能。

图5 入口速度为5 m/s的三维速度场分布Fig.5 Velocity field at inlet velocity of 5 m/s

图6 入口速度为5 m/s的流道内速度场分布Fig.6 Velocity field of flow pass at inlet velocity of 5 m/s

从图中可以分析看出，对于层叠丝网结构的回热器，丝网叠加的随机性带来了流动的随机性，使得其性能很难预测。即使是相同材料、相同孔隙率、相同大小、相同水力直径的回热器，由于丝网叠加形式的不可控，也会造成性能的不同。而对于蚀刻薄片卷裹式回热器来说，其速度分布更为规则和均匀，没有出现图6(a)中部分填料区域没有发生换热的情况。工质在平行于流动方向的流道只受到边界粘应力的影响，其阻力较小。而斜流道的出现，适时地打破热边界层对换热的影响，同时对流阻贡献不大，可以看到，工质流过斜流道后到达直流道，速度迅速达到恢复原来的大小。这种流场换热均匀，流动阻力小，对提高回热器的换热效率，降低流阻教层叠丝网回热器有较大的优势。

3.2.2 稳态压降

根据不同流动速度下，模拟的层叠丝网以及蚀刻薄片式回热器结构密度云图，可以得到回热器结构进出口的格子密度差，参考样机的制冷机的工作压力pa=1.75 MPa，取平均温度为Ta=190 K。查氦气物性后。由式(1)计算得到稳态的流阻系数。两种回热器结构的稳态流阻如图7所示。

图7 层叠丝网回热器和蚀刻薄片卷裹式的流阻系数的对比Fig.7 Resistance coefficient vs Reynolds number of mesh screen and etched foil regenerator

式中:ρ为工质密度，(kg/m3)，qv为回热器内工质体积流量，m3/s;dh为回热器水力直径，m;μ为工质动力粘度，Pa·S;d为回热器直径，m;φ为孔隙率;L为回热器长，m。

由图7分析可得:随着工质雷诺数Re的增加，两种回热器的工质稳态流动阻力系数f均减少，层叠丝网回热器的阻力系数比蚀刻薄片卷裹式回热器大，说明蚀刻薄片卷裹式回热器有流动阻力小的特点。

4 回热器稳态流阻测试及分析

4.1 稳态流动的测试试验系统

如图8所示为稳态测试系统中的实验与测试系统图。测试时压力为0.3 MPa，温度为290 K。测试对象为280目的层叠丝网回热器和蚀刻薄片卷裹式回热器，工质为氦气和氮气。

图8 稳态流动系统测试图1.进口混合腔;2.出口混合腔;3.直流电源转换器;4.差压变送器;5.压力传感器;6.Agilent 34901A数据采集卡;7.T型热电偶;8.Agilent 34970A数据采集仪。Fig.8 Test bench of stable-state flow system

4.2 LBM模拟数据与实验结果比较

将运用LBM的模拟数据与实验所测得的数据进行比较，其比较结果如图9所示。

图9 回热器模拟值和实验值的对比Fig.9 Value of simulation vs value of experiment

从图中可以看出，模拟值在趋势上是符合实验值的，但出现了整体的偏移。模拟值较实验值偏小50%左右。

上述模拟结果与实验值的对比可以看出，模拟结果在趋势上符合实际，但在压降的精确计算上产生了偏差。本文认为，可能是由以下因素造成的:

(1)流体进入填料结构时，流场会有一个发展过程。实际的回热器由于流体发展所占的空间较回热器总体积很小，可以忽略不计。而模拟中的模型体积过小，由于流场的入口效应，流场可能还没进入充分发展的阶段，因此对于微观结构的模拟结果造成了较大影响。

(2)模拟中使用的回热器模型极小，忽略了温度对回热器性能的影响。而回热器实际使用中整体温度、压力跨度大，工质的物性也变化也较大。该模拟无法完整的模拟回热器大温差的特点。所以对模拟的结果有一定的影响。

(3)本文模拟的模型尺寸过小。多孔介质具有复杂微观结构，离散时需要大量的网格来描述，很多时候远远超出可承受的范围，这是以往文献中很少见到关于对回热器微观尺度模拟的原因。以叠层丝网计算模型为例，模拟模型长度为0.6 mm，而实验采用的模型为29.88 mm，实际尺寸约为模拟尺寸的50倍。

5 结论

(1)建立了层叠丝网式和蚀刻薄片卷裹式回热器的物理模型，并使用格子波尔兹曼的等温D3Q19格子模型对其进行模拟计算。计算得到了两种填料结构的速度场分布和两端的密度差。并在此基础上分析了两种回热器结构的对于流场的影响。

(2)通过模拟得到填料两端的密度差，计算了不同雷诺数下的稳态流阻系数。结果显示:雷诺数越大，稳态流阻系数越小且变化幅度区域平缓。在相同雷诺数下，水力直径和填充率相近的卷裹式回热器的稳态流阻系数较层叠丝网式回热器的低。

(3)对模拟结果与实验值进行了对比，模拟结果在趋势上符合实际，但在压降的精确计算上产生了偏差。根据实验结果分析了模拟误差产生的原因。对进一步通过格子波尔兹曼方法模拟回热器结构有一定的指导意义。

(4)下一步工作展望:1)本文计算的模型尺寸在0.1—1 mm之间，相对回热器来说，这个尺度显得太小了，速度边界层和充分发展边界条件会干扰流场内的流动，四周的壁面也带来了额外的沿程阻力损失，使得模拟结果与实际有很大的偏差。所以下一步应该扩大模型的尺寸。此外，格子玻尔兹曼计算过程中每个离散微团的独立性非常适合显卡的GPU运算，因此并行GPU运算是今后格子玻尔兹曼方法求解回热器微观流场的一个可以考虑的方向。2)可压热格子玻尔兹曼方法的引入:本文所计算的模型忽略了回热器沿着轴线的巨大温差，假设其为等温，采用不可压等温格子玻尔兹曼方法计算稳态压降。若增加计算规模，建立大尺寸的回热器填料模型，则可以引入热格子玻尔兹曼模型，将温度的影响列入考虑范围，进一步使模拟结果接近实际工况。

1 张晓青，郭方中，董凯军，等.热声热机回热器特性参数影响的模拟分析及优化［J］.低温工程，2000，(3):7-11.

2 章忠敏.回热器效率的差分计算［J］.制冷学报，1982(3):18-28.

3 Atrey M D，Bapat S L，Narayankhedkar K G.Theoretical analysis and performance investigation of Stirling cycle regenerators［J］.Cryogenics，1991，31(12):1044-1052.

4 He Xiaoyi，Luo Lishi.A priori derivation of the lattice Boltzmann equation［J］.Physical Review E，1997，55(6):R6333-6336.

5 Latt J.Hydrodynamic limit of lattice Boltzmann equations［D］.Genève;Univ.Genève，2007.

6 何雅玲，王 勇，李 庆.格子Boltzmann方法的理论及应用［M］.北京:科学出版社，2009.

7 郭永飞.低温制冷机用径轴向混合填充式回热器的流动与回热特性研究［D］.上海:上海理工大学，2011.

8 王 强.新型蚀刻金属薄片式回热器的流动特性研究［D］.上海:上海理工大学，2012.

9 马诗旻.低温制冷机回热器的数值模拟［D］.上海:上海理工大学，2012.

10 Silva G，Semiao V.First-and second-order forcing expansions in a lattice Boltzmann method reproducing isothermal hydrodynamics in artificial compressibility form［J］.Journal of Fluid Mechanics，2012，698:282-303.

11 Wolff K，Marenduzzo D，Cates M E.Cytoplasmic streaming in plant cells:the role of wall slip［J］.Journal of the Royal Society Interface，2012，9(71):1398-1408.

12 Wu F，Shi W，Liu F.A lattice Boltzmann model for the Fokker-Planck equation ［J］.Communications in Nonlinear Science and Numerical Simulation，2012，17(7):2776-2790.

Numerical simulation of cryocooler regenerator using Lattice Boltzmann Method(LBM)

Xia Yudong Chen XiMa Shimin Zhang Hua

(School of Energy and Power Engineering，University of Shanghai for Science and Technology，Shanghai 200093，China)

Based on Lattice Boltzmann Method(LBM)，the flow field in the microstructure of mesh screen and etched foil regenerators were simulated.The velocity and pressure drop were obtained by LBM.The simulation results show that velocity field in etched foil regenerator is better distributed than that of mesh screen.And the etched foil regenerator has less resistance coefficient than that of mesh screen regenerator.The simulation results basically agree with the experiments.

Lattice Boltzmann Method;regenerator;resistance coefficient;numerical simulation

TB651

A

1000-6516(2012)05-0041-05

2012-07-17;

2012-09-25

国家自然科学基金(50906054)，教育部博士点新教师基金项目(20093120120006)资助。

夏宇栋，男，24岁，硕士研究生。