高温气冷堆球床接触导热计算方法

，，，∗，，

（1．环境保护部核与辐射安全中心，北京100082；2．清华大学核能与新能源技术研究院，北京100084）

高温气冷堆球床接触导热计算方法

李聪新1，任成2，许超1，∗，温丽晶1，刘宇生1

（1．环境保护部核与辐射安全中心，北京100082；2．清华大学核能与新能源技术研究院，北京100084）

接触导热是低温条件下高温堆球床堆芯传热的主要方式，接触导热的计算精度决定了低温条件下高温气冷堆球床传热计算的精度。本文建立了类似热离散单元法的球床接触导热的计算模型，通过计算流体力学方法进行模型中两个球之间理想接触的导热量计算，给出了所编写的大规模随机堆积球床接触导热计算程序的详细步骤。通过与计算流体力学方法对球床局部堆积结构接触导热计算结果的比较，验证了该计算程序的准确性。该方法与球床局部结构计算流体力学方法精度接近，并可用于颗粒尺度大规模球床计算。

高温堆；球床；接触导热；颗粒尺度；理性接触；热离散单元法

高温气冷堆球床堆芯是由球形燃料元件随机堆积而成的［1］，在低温条件下，球床的主要传热方式为接触导热［2］，球床接触导热的计算精度决定了低温条件下高温气冷堆球床传热计算的精度。目前，从颗粒尺度进行球床传热计算的主要方法为热离散单元法［3］（TDEM，Thermal Discrete Element Method）。TDEM是在颗粒流动计算的离散单元法 （DEM，Discrete Element Method）基础上增加了传热模型［4］，将每个球用一个平均温度表示，从颗粒尺度进行球床传热计算［5］，从而得到每个球的温度和传热量，既可以研究球床的局部传热特性［6］，又可以采用统计平均方法研究球床等效导热系数［7］、热流量等宏观参数，是目前球床流动传热计算研究的热点［8，9］。热离散单元法的主要问题是缺乏严格的数学证明，模型中存在大量的工程假设和经验参数，所以热离散单元法可以认为是一种工程方法，而不是严格的理论计算方法。目前开源的DEM软件有LIGGGHTS［10］和OpenFOAM［11］，商业软件有EDEM和PFC3D等。在这些软件中，颗粒间的接触导热量通常表示为：

ΔT为两个接触球的温差，hc为接触面传热系数，其表达式为：

λs1、λs2为两个接触颗粒各自的材料导热系数；FN为重力；r∗为两个接触颗粒的几何平均半径；E∗为等效杨氏模量；括号中的部分表示颗粒间的接触面积。从公式 （2）可以看出，接触面传热系数表示为了导热系数和接触面积的一个简单函数，而接触面积则是按照弹性力学假设的一个经验关系式。该公式将传热和受力变形的经验关系式结合在一起，导致接触传热系数的计算精度难以控制，用热离散单元法进行的球床传热计算结果也将难以验证。

本文在接触导热模型中将受力变形假设和传热计算分离，只进行理想接触条件下的传热计算，通过与计算流体力学方法对球床局部堆积结构接触导热计算结果的比较，验证了计算程序的准确性。该方法可用于颗粒尺度大规模球床导热计算，并得到与球床局部结构计算流体力学方法接近的精度。在进行实际球床传热计算时，可以通过实际接触面积与理想接触面积的比值耦合力学参数。

l 接触导热计算模型

在低温条件下，球表面的辐射传热可以忽略，则球之间通过接触面的导热将成为主要传热方式。将每一个球定义为一个单元，每个球的温度用球的平均温度Ti代表，定义流入球的热流方向为正，相邻两个球之间的导热量为：Ti为当前球的平均温度；Tj为与当前球接触的球的温度；hci定义为接触面传热系数。

根据能量守恒，传入球内的总热量等于球的内能变化：

其中，n为与当前球接触的球的个数。Cs为球的热容，定义为球的质量与比热容的乘积：

当处于稳态时，球的净热流量为零，则：

对于靠近壁面的球，还需要相应的边界条件才能使方程组闭合。根据目前球床导热计算的需要，只处理给定壁温和绝热两类边界条件。

当球处于给定壁温边界时：hci为接触半径相同的球之间的接触面传热系数，Tw为接触壁面的温度。即在同等接触面积时，球和壁面之间的接触面传热系数为两个球之间的接触面传热系数的两倍。

当处于绝热边界条件时，相当于减少了一个热流通道：

本文的接触导热模型与一般热离散单元法的导热模型类似，区别在于接触面传热系数的处理上。热离散单元法将传热和受力变形的经验关系式结合在一起，导致接触传热系数计算精度难以控制。本文将受力变形假设和传热计算分离，进行单纯的接触导热计算，受力变形关系由离散单元法在确定堆积结构时确定。球之间的接触导热采用理想接触面积计算，将由表面粗糙度引起的

Rc是半个球的接触热阻。研究两个球之间接触热阻的文章很多，早期多采用解析方法和简化假设得到近似解析解［12，13］，但不同解析解的差异较大，且都存在较大误差。文献［14］分析了导致差异和误差的原因，并通过计算流体力学方法精细计算了不同接触半径的球之间的导热量，得到了一个较为精确的接触面传热面系数表达式：

其中r为球的半径，λs为球的导热系数，Ks是一个由相对接触半径确定的形状函数，其通过计算流体力学方法计算半径为1m、导热系数为1W·m-1·K-1、相对接触半径为γ的两个半球在一定温差下的导热量得到：实际接触面积减小用一个0和1之间的系数Kn加以考虑。

在球床接触导热计算中，关键在于获得两个球之间的接触面传热系数hci。接触面传热系数的引入是为了计算方便，其定义为两个半球之间接触热阻的倒数，即：

形状函数Ks确定以后，理想接触的两球之间的导热量可以通过公式 （10）确定。当考虑实际物体的表面粗糙度时，实际导热面积小于理想接触面积，因此引入一个表示实际接触面积和理想接触面积的比例系数Kn，其取值范围为0到1之间。因此，表示两个球之间接触面传热系数的完整公式为：

2 计算程序开发

在球床的接触导热计算中，将每个球的温度用一个平均温度表示，每个球的能量守恒可以确定一个线性方程，所有球的线性方程组成一个闭合的线性方程组。由于球床的随机堆积状态，每个球所接触的球的个数和方向不同，导热模型的计算过程相当于一个非结构化网格的固体导热计算。球床的接触导热计算的完整计算程序包括前处理、方程求解、后处理三个部分，以下对本文所开发的程序的这三个部分进行简要说明。

2.l前处理

前处理用于获得球所接触的球和壁面信息，相当于流体力学计算中的网格信息。由于离散单元法程序只是给出了一定堆积结构内的球的位置坐标，需要从球的坐标中提取球间的接触和边界信息。以立方体随机堆积为例，首先对球的坐标沿某一方向排序 （可取主热流方向），然后对每个球按排序结果编号，球的编号是求解计算的基础，其与方程和节点相对应，对于由一万个球堆积成的球床，球的编号即为1～10 000。完成编号以后，对球床中任意两个球的球心距离做一次比较，当小于球直径时即认为两个球接触。在对球心距离进行比较的过程中，如图1所示，需要存储当前球的接触球总数、每个接触球的编号和相对接触半径。最后进行一次壁面边界信息搜索，比如对于球心Z方向坐标小于一倍球半径的球认为与下壁面接触，确定与边界接触的球的总个数以及每个边界球的编号和相对接触半径。在立方体堆积球床中，有六个壁面，在环形球床中有四个壁面，壁面边界信息中还应存储边界条件信息，如温度边界条件或绝热边界条件。在获得了每一个球和壁面的接触信息之后，前处理过程完成。

图l 每个球和壁面的存储信息Fig.l Stored information of each pebble and wall

2.2 方程求解

方程求解的关键是通过球和壁面信息组成一个线性方程组，形式为：

A为系数矩阵，维度为N×N，N为球床中球的个数，x是N×1矩阵，即未知温度，b方程中与未知温度无关的常数，由边界条件确定。由于上述方程可采用任何线性方程组的数值解法进行求解，因此计算的难点在于组成矩阵A和b。

根据前处理得到的球和壁面信息，可以通过两个大的循环完成A和b的计算。初始时将A和b的全部元素值赋零。首先是对球的循环，对于每个球，根据公式 （12）确定当前球i与邻球j的接触面传热系数，j为邻球的编号。矩阵A的主对角线 （i，i）的值为其与所有邻球的接触面传热系数和的负值，矩阵A的 （i，j）元素的值为当前球与邻球j的接触面传热系数。

其次是对壁面的循环，给定温度和绝热边界条件的处理方法不同。对于给定温度的边界条件，每个壁面球按公式 （7）确定其与壁面的接触面传热系数。对于编号为j的壁面球：

hwj为当前球和壁面的接触面传热系数，Tj为壁面的温度。A（j，j）n-1为在此次循环之前确定的系数矩阵中主对角线上元素的值，b（j）n-1为当前循环之前b（j）的值。

对于给定的绝热边界条件，其相当于对该接触点对球能量守恒方程无贡献，直接忽略该循环即可。

在完成了上述两个循环之后，系数矩阵A和右端项b即计算完成。在一般流体力学计算中，系数矩阵的非零元素在二维情况下为五对角矩阵，在三维情况下为七对角矩阵。与一般流体力学固体导热计算的系数矩阵不同，由于球床中球接触的随机性，系数矩阵也具有一定的随机性。以前文模拟的立方体随机堆积球床为例，系数矩阵A（左上角的部分）的非零元素如图2所示，形成一个宽度约为600的随机对对角线，每一行中非零元素为2～13个，由于壁面效应造成的配位数变化也可以在一定程度上反应在带状边界宽度上。由于系数矩阵A为严格对角占优矩阵，因此方程组是稳定可解的。为提高计算效率，可采用一些大型稀疏矩阵算法，如multifrontal算法［15］。

图2 系数矩阵A中非零元素位置Fig.2 Positions of non-zero elements in the coefficient matrix A

2.3 后处理

后处理主要是根据计算得到的每个球的温度和热流量进行相关物理参数的计算和统计，如球床温度分布、边界热流量、球床等效导热系数等。由于球床内每个球的温度都已经确定，这些物理参数不过是一个求和、平均的统计过程。

以立方体内的球在上下两侧给定温度、四周绝热条件下的热流量和等效导热系数为例。给定温度的壁面导热量的计算与第一类边界条件设置算法类似，每个与壁面接触的球与壁面之间传递的热量为：

其中hwi为球与壁面的接触面传热系数，Tw为设定的壁面温度，Twi为计算得到的与壁面接触的球的温度。通过该壁面的总的热流量为：

其中ΔT为球床两次温度差，ΔL为球床的高度，A为球床的横截面积。

3 模型验证

n为前处理阶段得到的与该壁面接触的球个数。球床的等效导热系数为：

3.l一维验证

将10个直径为60 mm，导热系数为100 W· m-1·K-1的石墨球等间距连成一排，球四周绝热，左右两侧所接触壁面的温度分别为0℃和100℃，通过改变球心距离计算在不同相对接触半径下通过单个球的热流量。球心距离在12～58 mm之间以1 mm为步长递增，在58．1～59．9 mm之间以0．1 mm为步长递增。

采用计算流体力学软件进行计算， （如CFX、Fluent等）计算结果作为比较对象，在不同接触半径下，球床接触导热模型和计算流体力学软件得到的通过每球的热流量如图3所示，在计算范围内。两者计算得到的热流量的相对误差最大值不超过2%，导热模型和计算流体力学方法的计算结果应保持一致，但可以看出，二者之间仍然存在一定误差。误差产生的原因是球接触热阻形状函数是通过数值方法得到的，其计算结果与网格疏密有关，在计算形状函数时，单一球的等效网格数目为千万量级，远远大于10个球计算时每个球的网格数目。若可以进一步增加网格数量，两种计算方法的计算结果应趋于一致。

图3 不同接触半径下通过单个球的热流量比较Fig.3 Comparison of heat flow through one pebble with different contact radius

3.2三维验证

为了验证球场接触导热计算方法可以应用于三维堆积结构中，本文仍然采用与计算流体力学软件直接计算结果作比较的方法。采用60 mm、导热系数为100 W·m-1·K-1的石墨球通过立方体 （SC）堆积而成，X、Y、Z方向球的个数分别为6、5、4。球床六个表面同时施加不同的温度边界条件，其中X方向为0℃和400℃，Y方向为100℃和500℃，Z方向为200℃和600℃，比较球内温度为三维分布时接触导热模型与计算流体力学软件计算得到的球的温度和热流量。由于120个球的计算流体力学软件固体导热计算耗时较长，因此不再进行每一接触半径下的比较。由于随着相对接触半径减小，同样的导热量计算精度所需的网格数据迅速增加，因此只比较球心距离为58 mm，实际接触半径为7．68 mm，相对接触半径为0．256的情况，此时所需网格数目尚在目前16G内存服务器的计算能力范围内。

由于接触导热模型基于一维球接触热阻得到，其在三维条件下是否适用是模型能否用于随机球床导热计算的关键。因此，上述边界条件的设置是为了使球床中的温度处于完全的三维分布状态。随机取出堆积区域中间的两个球，计算流体力学软件计算得到的等温面和热流线如图4所示。由于接触点的温度各不相同，球内部的温度场已经完全处于三维分布状态，从随机取出的两个球看，其等温面已经不存在相似之处。

图4 球等温面和热流线的计算流体利力学软件计算结果Fig.4 The isothermic surface and Heat flux lines of CFD software Calculated result

采用计算流体力学方法和球床接触导热方法计算得到的温度分布如图5所示，从整体看，两种方法得到的计算结果都反映了所设置的边界条件特征，温度最小值出现在X=0的平面，温度最大值出现在最上层。由于堆积结构的六个表面设置的温度各不相同，球床中的球温度均显示出缓慢渐变的特点。导热模型相当于计算了球的平均温度，不能反映球在接触点附近的温度变化，例如在X=0的平面上设置的温度为0℃，可以从计算流体力学软件计算结果中明显看出X=0的平面上球接触点的温度均为0℃，而导热模型计算结果则为求平均后的温度。

图5 三维球床温度分布比较Fig.5 Comparison of 3D temperature distribution of pebble bed

为了定量比较导热模型和计算流体力学软件计算得到的球温度，取Z方向最上层的30个球，编号如图5（b）所示。对计算流体力学软件计算结果采取体积平均的方式计算每个球的平均温度，与导热模型得到的球的温度作比较，结果如图6所示。可见，两种方法得到的球温度几乎相同，最大相对误差不超过1．7%。从图6可以看到两种方法计算得到的温度都存在5个周期变化，每个周期有6个温度点，正好与球的排列一致。这30个球的温度变化幅度较大，这是由于所设置的边界条件剧烈变化所致，可见采用接触导热模型并不限于小温差假设。

图6 Z方向最上层球温度比较Fig.6 Temperature comparison of the top layer pebble in Z direction

在上述堆积方式中，六个边界上的总热流量见表1，其中热流以流入球内的方向为正。在计算流体力学软件计算中，边界总热流量采用对每个接触面的热流密度的积分得到；在球床接触导热计算方法中，采用公式 （16）计算，即对与相应壁面接触的每个球的热流量求和。从表1可以看出，边界热流量最大相对误差为4．7%，考虑到每个球内部的三维温度分布和非常大的温度梯度，这一精度完全可以接受。

通过接触导热模型与计算流体力学方法在热流量和温度计算方面的比较，可见基于一维热阻推导的接触导热模型在三维球床导热计算中同样适用。接触导热模型可以进行大规模球床导热的计算，与计算流体力学方法的相对误差控制在5%左右。这一精度是目前热离散单元法难以达到的。由于现有的计算流体力学软件无法进行大规模的随机堆积球床的导热计算，因此不再进行大规模球床的接触导热验证。

表l 边界总热流量比较Table l Comparison of boundary total heat flux

4 结论

本文建立了类似热离散单元法的球床接触导热的计算方法，将导热计算和力学经验关系式区分离，相当于一种理想条件下球床导热的高精度计算方法。通过计算流体力学方法进行了模型中两个球之间理想接触的导热量计算，提高了接触面传热系数计算的净多。本文给出了所编写的大规模随机堆积球床接触导热计算程序的详细步骤。通过与计算流体力学方法对球床局部堆积结构接触导热计算结果的比较，验证了本文计算程序的准确性。通过引入一个表示实际接触面积和球心距离确定的理想接触面积的比例系数，可以进行实际球床接触导热的工程计算。

［1］Zhang Z，Wu Z，Sun Y，et al．Design aspects of the Chinese modular high-temperature gas-cooled reactor HTR-PM［J］．Nuclear Engineering and Design，2006，236（5）：485-490．

［2］Nishino K，Yamashita S，Torii K．Thermal contact conductance under low applied load in a vacuum environment［J］．Experimental thermal and fluid science，1995，10（2）：258-271．

［3］张品．基于DEM方法的散体颗粒热传导模拟 ［D］．湖南：中南大学，2011．

［4］Nguyen V D，Cogn E C，Guessasma M，et al．Discrete modeling of granular flow with thermal transfer：application to the discharge of silos［J］．Applied thermal engineering，2009，29（8）：1846-1853．

［5］Bu C，Liu D，Chen X，et al．Modeling and Coupling Particle Scale Heat Transfer with DEM through Heat Transfer Mechanisms［J］．Numerical Heat Transfer，Part A：Applications，2013，64（1）：56-71．

［6］Malone K F，Xu B H．Particle-scale simulation of heat transfer in liquid-fluidised beds［J］．Powder Technology，2008，184（2）：189-204．

［7］Zhang H W，Zhou Q，Xing H L，et al．A DEM study on the effective thermal conductivity of granular assemblies［J］．Powder Technology，2011，205（1）：172-183．

［8］Nguyen V，Benhabib K，Marie C，et al．Heat Transfer in a Granular Media Modeled by a Coupled DEM-Finite Difference Method：Application to Fluidized Bed Processes［J］．Procedia Engineering，2012，42：895-904．

［9］Brosh T，Levy A．Modeling of Heat Transfer in Pneumatic Conveyer Using a Combined DEM-CFD Numerical Code［J］．Drying Technology，2010，28（2）：155-164．

［10］Vedachalam V，Virdee D．Discrete Element Modelling Of Granular Snow Particles Using LIGGGHTS［D］．Ph．D Dissertation．Edinburgh Parallel Computing Centre．The University of Edinburgh，UK，2011．

［11］Jacobsen N G，Fuhrman D R，Fredsoe J．A wave generation toolbox for the open-source CFD library：OpenFoamⒸ［J］．International Journal for Numerical Methods in Fluids，2012，70（9）：1073-1088．

［12］Carslaw H S，Jaeger J C．Conduction of heat in solids［J］．Oxford：Clarendon Press，1959，2nd ed．，1959，1．

［13］Yovanovich M M．Thermal contact resistance across elastically deformed spheres［J］．Journal of Spacecraft and Rockets，1967，4（1）：119-122．

［14］李聪新．高温气冷堆球床等效导热系数研究 ［D］．清华大学，2014．

［15］Davis T A，Duff I S．A combined unifrontal/multifrontal meth

od for unsymmetric sparse matrices［N］．ACM Transactions on Mathematical Software（TOMS），1999，25（1）：1-20．

Calculating the Contact Heat Conduction in a High Temperature Gas-cooled Reactor

LI Congxin1，REN Cheng2，XU Chao1∗，LIU Yusheng1，ZHANG Pan1
（1．Nuclear and Radiation Safety Center，MEP，Beijing 100082，China；2．Institute of Nuclear and New Energy Technology，Tsinghua University，Beijing，100084，China）

At low temperatures，heat transfer within the core of the high temperature gas-cooled reactor（HTGR）is dominated by contact heat conduction，the prediction of which thereby determining how accurate the heat transfer calculation for any HTGR could be．On the basis of the thermal discrete element method，this paper establishes a model to calculate the contact heat conduction occurring in the pebble bed．In doing so，we first resolve the conduction flux between two ideally-contacting balls via CFD approach，and then develop the detailed procedure of our code which calculates the contact heat conduction for large-scale，randomly-packing pebble bed．Our code has been validated by evaluations against predictions obtained by CFD approach for the contact heat conduction within local structure of pebble bed．Our method enjoys accuracy comparable to CFD simulation of local pebble bed structure，and it is further applicable to large-scale pebble bed simulation in the particle level．

HTGR；pebble bed；contact heat conduction；particle scale；ideal contact；thermal discrete element method

TL331

：A

：1672-5360（2016）04-0052-07

2016-06-24

2016-08-03

中国科学院先导专项，项目编号 XDA02050500

李聪新 （1984—），男，河北衡水人，工程师，核科学与工程，现主要从事核安全相关热工水力试验研究工作

∗通讯作者：许 超，E-mail：seanwillian＠126．cn