高温气冷堆用碳毡材料导热系数测量及反问题计算

李聪新，任 成，杨星团，姜胜耀，孙艳飞

（清华大学核能与新能源技术研究院先进反应堆工程与安全教育部重点实验室，北京 100084）

20世纪60年代日本Shindo发明以聚丙烯腈（PAN）为原料制造碳纤维，经过50多年的发展，碳纤维材料行业已发展成独立完整的工业体系［1］，碳纤维毡作为绝热保温材料获得了广泛应用［2］。列入《国家中长期科学和技术发展规划纲要》16个重大专项之一的高温气冷堆核电站项目，由于采用了石墨堆芯结构和石墨基体燃料元件，相关设计和实验采用了大量性质兼容的碳纤维材料。高温气冷堆运行温度为750℃［3］，极限事故工况下燃料元件的安全温度限值为1 600℃［4］，相应的热工、余热排出等高温试验需要高性能的耐高温保温材料。碳毡中90%为空隙，保温隔热性能优越、稳定［5］，成为高温气冷堆相关实验优先选择的保温材料。

国内对碳毡材料力学［6］、热学、电学［7］等方面的性能、参数的研究和测定不足，尤其是在其更具优势的高温领域（1 000℃以上）。碳毡的导热性能受温度、密度、石墨化度、碳纤维取向等多种因素影响［8］，并非常量［9］，导热系数仍以实验测量为主。高温条件下通常采用激光闪射法［9－13］和热线法［14］。这两种方法的测量精度受两个因素影响较大，一是小尺寸试样选取不具有代表性，二是测试环境与使用环境不一致。

现代工业越来越多的要求实时精确的测量热物性，因此基于导热反问题的材料物性参数识别方法取得了广泛应用［15］。自60年代Stolz发表简单形状的导热反问题数值计算方法［16］以来，导热反问题的计算方法极大丰富，如正则化方法、LM方法、神经网络算法等［17－19］。本文以实测的200多万组实验数据进行导热系数反问题计算，详细介绍采用非线性导热反问题方法确定材料温度相关导热系数在物理模型简化、非线性导热正问题求解、导热系数反问题识别等过程的具体算法和步骤，为用导热反问题方法测量物性参数提供一个基本思路。

1 实验

1．1 材料性能测试装置

为了给高温气冷堆工程及相关实验提供可靠材料验证和技术支持，填补高温气冷堆在材料相关配套工业方面的不足，清华大学核能与新能源技术研究院建立了一套模拟高温气冷堆温度、环境氛围的材料性能测试装置，以进行相关材料1 600℃以下的性能测试，温度涵盖高温气冷堆从正常运行到极限事故工况下的全部温度范围。该装置设计为石墨电阻炉形式，中心采用高纯度石墨加热元件，上、下和四周为碳毡保温材料。石墨发热体和碳毡保温材料之间形成模拟高温气冷堆高温、碳还原环境的环形材料测试区，该区域可模拟真空（＜30Pa）、氮气、氦气环境。高温区采用特殊设计的可在碳还原环境中耐受高温的钨铼热电偶稳定测量温度场，精度±0．5%。测试装置的内部结构和部分测点位置如图1所示。

通常，保温材料导热系数越小越好，但考虑到承重和强度因素，在材料性能测试装置中，上保温层采用软质碳毡，底部和侧面保温层采用导热系数稍大的硬质碳毡。保温材料的导热系数一般随温度的升高而增加［20］，在设计时，选取的导热系数为厂家给出的设计参考值（表1）。要在很宽温度范围内测量保温材料的导热系数有很大难度，这也是本文采用导热反问题方法的原因之一。

图1 材料测试装置内部结构与测点布置

表1 保温材料热物性参考值

为使材料测试区近似柱坐标一维传热，上下保温层厚度远大于侧保温层厚度。环形区高度为1m，测点T1布置在环形区轴向和径向中心，测点T9、T8距上下保温层分别为125mm，三者温差很小。连同上下保温层中测点，在1 400℃和1 600℃保温阶段，轴向温度分布如图2所示，可认为炉膛区域轴向温度恒定，用一维径向传热模型近似误差可忽略。

侧保温层厚度为200mm，内半径R1为250mm，外半径R3为450mm，在其中心和外边界分别布置测点T2、T3。侧保温层内侧无测点，考虑到环形区温度分布较为均匀，可用测点T1温度作近似代替保温层内壁面T′1处温度来进行导热计算。近似替代引入的误差和修正方法将在后文介绍。

图2 轴向温度分布

1．2 实验过程及结果

一次典型的材料测试实验包括以下过程：1）对炉体小功率预热，同时抽真空，防止碳毡和石墨材料被氧化；2）当炉膛内压力降到100Pa以下，调节功率升温；3）当材料测试区温度上升到1 400℃后，保温24h；4）再次缓慢升温到材料测试区温度为1 600℃，保温22h；5）停止加热、保持冷却水流量不变，缓慢降温；6）降温到165℃时，停止抽真空，充入氮气冷却；7）温度降到120℃时，停止冷却水循环，实验完成。

整个实验过程持续时间约为76h，数据采集系统对各测点的温度采集速率为8次／s，共得到219万组数据，3个用于碳毡导热系数计算的测点T1、T2、T3温度如图3所示。足够多的测量数据对采用灵活的反问题计算方法是有利的。

图3 温度历史曲线

2 正问题

导热正问题求解是反问题计算的基础，正问题求解的精度和速度对反问题计算影响很大，须首先建立合理的正问题的计算模型。

2．1 数学描述与简化

由于轴向导热可忽略，正问题可近似为一维柱坐标下变物性非稳态导热，将求解区域确定为R1～R3之间的碳毡保温层区域，其数学描述如下。

控制方程：

边界条件：

T（R3，τ）＝Y3（τ）

（3）

补充条件：

T（R2，τ）＝Y2（τ）

（4）

目标函数：

其中：r为空间径向坐标；τ为时间；τi为测量时间点；T（r，τ）为r处τ时刻温度计算值；Y1（τ）、Y2（τ）、Y3（τ）为测点T1、T2、T3在τ时刻的温度测量值。待识别的导热系数在满足式（1）～（4）的前提下使得由式（5）定义的误差函数值最小。

为了简化正问题计算，用测点T1的温度近似替代内壁面T′1处温度，将内壁面的辐射边界条件转化为第一类边界条件，式（2）变为：

T（R1，τ）＝Y1（τ）

（6）

假设初始条件不存在，主要考虑一长时间的测量过程，可选取任一时间范围内的数据进行反问题计算，其起始时刻并不一定是温度为室温的绝对实验开始时刻。可用R1、R2、R3处温度测量值的抛物插值作为初始条件，其与此时真实初始条件之间的差异将随计算时间的推移而消除。

至此，由方程（1）、（3）、（6）构成的正问题已满足定解条件。

物性参数ρT为密度，cT为比热容，λT为导热系数，下标T表示其为温度的函数。对于工程中常见的大多数材料，从严格意义上讲，物性参数均非定值，但对于固体，可认为密度不随温度变化，即方程（1）等价于如下方程：

比热容一般随温度变化，但由于两个边界条件均为一类边界条件，不能同时反求导热系数和比热容［21］，为确定导热系数，必须假定比热容为定值，取为厂家给出的参考值。此时控制方程简化为：

将式（8）右端偏导项展开，得到偏导数可直接用差分格式代替的形式：

2．2 数值计算方法

虽然目前主流的商业CFD软件如FLUENT、CFX可精确求解三维固体热传导方程，但是导热反问题计算过程中，需根据反问题算法不断修改待识别参数，进行导热正问题计算，所以高效的做法是根据具体问题统一编写正、反问题计算程序。

如图4所示，将侧保温层沿径向等分为20个网格，Δr＝（R3－R1）／20＝10mm。

图4 离散格式

将式（9）中的偏导项用如下差分形式代替：

则可得到：

其中：Tn为节点n的当前温度；T0n为节点n前一时刻温度；Tn－1、Tn＋1为该节点相邻节点温度；Δτ为时间步长；Δr为空间步长；rn为当前节点半径。将式（11）整理成计算形式为：

需注意，此离散格式有别于常物性一维柱坐标的离散格式：

写成Ax＝d的矩阵形式为：

其中：

系数矩阵为标准的三对角矩阵。由于导热系数为非定值，系数矩阵中元素为温度的函数，可用追赶法等三对角常系数矩阵算法迭代求解。

2．3 辐射边界处理

在上述计算过程中，左边界条件用环形材料测试区中的测点温度Y1代替，该测点温度不同于真实碳毡保温层内壁面温度。测试区为真空氛围，辐射是热传递的主要方式，两者差值与温度和热流密度有关，可采用如下“预估－校正”方法。

1）预估过程

用测点温度Y1作为计算的边界条件，用上述方法计算得到整个实验段内的温度场，将节点2的温度记为T02。假设壁面与温度为T预的灰体处于辐射平衡状态，则有：

其中：ε为灰体表面发射率，实测值为0．8左右；σ为斯特潘－波尔兹曼常量。图5为表面发射率分别取0．6、0．8、1．0时的预估温差比较，可见，相对于1 600℃的温度变化范围，预估温度对表面发射率变化并不十分敏感。由式（16）得到测点与壁面温度的预估温度差值为：

2）校正过程

将预估差值作为测点和壁面温度差值，得到以后计算所用内壁面温度：

用新壁面温度作为以后计算的边界条件。该过程并不需要迭代，只需在反问题计算过程中，对导热系数值变化较大的试算过程，重新执行一次校正过程。

用“预估－校正”方法是为了克服直接处理辐射边界条件带来的计算量增大和发散。预估差值与测点温度Y1的关系示于图5a，由于该差值不是温度的单值函数，不能用测点温度做简单补偿。与直接用测点温度Y1作为内壁面温度相比，使用该方法更接近实际传热过程，可使由式（5）定义的误差值减小，尤其是实验刚开始的升温过程中，其原因可从图5b看出，在刚开始的升温过程中，壁温和测点温差较以后过程大很多。

图5 预估温度差值

3 导热系数反问题计算

3．1 一般反问题计算方法

对导热系数进行识别，预先对其满足的函数关系式进行假设可简化反问题计算。由于大多数保温材料的导热系数可认为是温度的线性函数［20］，可假设其满足如下关系式：

λ＝a＋bT

（19）

对导热系数的函数识别变为对参数a、b的识别。反问题计算即通过不同的参数a、b组合求解导热正问题，使由式（5）定义的误差函数最小，这是一个二维最优化问题。

J（a，b）＝min（Y－T）T（Y－T）

（20）

由多元函数极值条件知，上式取最小值的条件为：

其中，系数矩阵称为敏感度矩阵，敏感度系数为未知参数的函数，可采用如下方法计算：

其中，ε为一个小量。上述问题可采用任何一种基于梯度的迭代求解方法，如最速下降方法、共轭梯度法［22］、Levenberg－Marquardt方法［23］等。这些方法都经过严格的数学推导，具有一定的普适性，但物理概念并不明确，初值选取不当会使迭代过程发散。

3．2 动静导热过程法

根据实验过程的特点和固体热传导规律，给出一种简便的优化方法。本次试验涉及稳态和非稳态导热过程，本文将利用这一特点分别识别a和b，简称为动静导热过程法。

假设导热系数为常数，无论其大小，则稳态过程中测点T1、T2、T3的温度差值之比（Y3－Y2）／（Y2－Y1）将保持不变，即通常认为在热流未知时，用稳态方法不能测量导热系数。要改变其比值，只有改变其导热系数随温度的变化率，据此可确定导热系数随温度的变化率。变化率确定后，改变导热系数零点参考值，即绝对大小，使非稳态过程中温度计算值与测量值偏差减小。这样便将一个二维搜索变为两个一维搜索。

为应用上述方法，将导热系数的表达式（式（19））改写为以下形式：

记a／b为c，分别计算温度对a、c的敏感度系数，结果如所图6所示，其中对应a每变化0．001，c变化0．000 01，6次等差变化的值。可见在稳态阶段，a的变化对温度场几乎无影响，所以当确定c后用非稳态过程优化a不会改变稳态过程的优化效果。

在假设一个a值前提下，可使用任何一种一维最小值搜索算法，如二分法，来确定c。在a取不同值时，使稳态过程误差函数最小的c值几乎不变，为0．000 403左右。测点T2处不同a值假设下搜索到最优c值时的温度历史曲线如图7所示。

图6 温度对a、c的敏感度系数

图7 测点T2处不同a值假设下最优c对应的温度历史

c值确定后，可使用同样的一维最小值搜索算法，得到在确定c值下使误差函数最小的a值，为0．099。

a和c确定后，为保证其组合为最优解，可在极小范围内调整a和c的值，也可以其为共轭梯度法等反问题算法提供迭代初始值，这种良好初值可使反问题算法迅速收敛。

最终导热系数的表达式为：为方便工程应用，实验得出的不同温度下硬质碳毡的导热系数列于表2。

λ＝0．098 2（1＋0．000 404 96T）

（24）

采用这种反问题方法进行物性参数识别是直接在碳毡材料的使用环境中对其整体导热性能的衡量，在工程实践中具有重要的实用价值。

表2 1 600℃以下硬质碳毡的导热系数

3．3 结果分析

在反问题计算过程中，导热系数的优化目标为测点T2位置处温度的测量值与计算值在不同时刻差值的二范数最小，最终二者绘于图8。由图8可见，识别的导热系数使测点T2处温度计算值与测量值在整个实验时间范围内基本一致。

图8 测点T2处温度计算值与测量值比较

T2位置处温度计算值的相对误差如图9所示。温度场恒定和平稳变化的过程中，相对误差基本控制在±2．5%以内。在功率瞬间调整和工况改变的实验开始和结束时，计算温度值的相对误差较大，最大为15%。在这段时间内，传热过程更为复杂，与计算模型之间存在更大的误差。例如时间坐标为80h处，装置内温度已降到200℃以下，辐射传热过程减弱，再充入氮气，气体导热作用不可忽略，再将内壁面作为辐射边界条件已不合适，必然引起更大的误差。

图9 计算温度相对误差

4 结论

随着国内需求量日益增长，碳纤维材料已被列为国家化纤行业重点扶持的新产品，成为国内新材料行业研发的热点。鉴于目前碳纤维材料在高温条件下性能参数数据不足，本文提供了清华大学核能与新能源技术研究院高温气冷堆材料性能测试装置对碳毡材料76h的验证温度数据，并据此采用非线性导热反问题方法确定了1 600℃以下硬质碳毡材料的导热系数，为高温气冷堆相关实验和其他工程应用提供了高温条件下导热性能设计依据。文中详细介绍了采用非线性导热反问题方法确定材料温度相关导热系数的完整过程和具体算法，在正问题计算中提出了处理辐射边界条件的“预估－校正”算法，减少了用测量数据作为第一类边界条件代替实际辐射边界条件引入的误差；在反问题处理中提出了依据物理过程的动静导热过程法进行导热反问题计算，物理概念清晰且简便可行。该方法既可单独使用，又可为其他反问题算法提供良好的迭代初值。

［1］ 赵稼祥．碳纤维的发展与应用［J］．纤维复合材料，1996，46（4）：46－50．

ZHAO Jiaxiang．Development and application of carbon fiber［J］．Fiber Composites，1996，46（4）：46－50（in Chinese）．

［2］ 孙志坚，孙玮，傅加林，等．国内绝热保温材料现状及发展趋势［J］．能源工程，2001，26（4）：26－28．SUN Zhijian，SUN Wei，FU Jialin，et al．Current status and development of thermal insulating materials in China［J］．Energy Engineering，2001，26（4）：26－28（in Chinese）．

［3］ 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．

［4］ NIESSEN H，BALL S．Heat transport and afterheat removal for gas cooled reactors under accident conditions，IAEA－TECDOC－1163［R］．Vienna：IAEA，2001．

［5］ 黎小平，张小平，王红伟．碳纤维的发展及其应用现状［J］．高科技纤维与应用，2005，30（5）：28－34．LI Xiaoping，ZHANG Xiaoping，WANG Hongwei．Modification and application of ultra－high molecular weight polyethylene fiber［J］．Hi－Tech Fiber ＆Application，2005，30（5）：28－34（in Chinese）．

［6］ 易法军，孟松鹤，韩杰才，等．碳毡／碳复合材料超高温力学性能实验研究［J］．无机材料学报，2001，16（6）：1 229－1 234．YI Fajun，MENG Songhe，HAN Jiecai，et al．Mechanical behavior of carbon felt－carbon composites under ultra－high temperature［J］．Journal of Inorganic Materials，2001，16（6）：1 229－1 234（in Chinese）．

［7］ 邱广玮，刘平，韩金铎，等．热处理温度对碳纤维毡导电性能的影响［J］．机械工程材料，2012，36（11）：66－68．QIU Guangwei，LIU Ping，HAN Jinduo，et al．Effect of heat－treatment temperature on conductivity of carbon felt［J］．Materials for Mechanical Engineering，2012，36（11）：66－68（in Chinese）．

［8］ 赵建国，李克智，李贺军，等．碳／碳复合材料导热性能的研究［J］．航空学报，2005，26（4）：501－504．ZHAO Jianguo，LI Kezhi，LI Hejun，et al．Research on the thermal conductivity of C／C composites［J］．Acta Aeronautica et Astronautica Sinica，2005，26（4）：501－504（in Chinese）．

［9］ 陈洁，熊翔，肖鹏．单向C／C复合材料导热系数的计算［J］．炭素技术，2008，27（2）：1－4．CHEN Jie，XIONG Xiang，XIAO Peng．Calculation of thermal conductivity coefficient in unidirectional carbon／carbon composite［J］．Carbon Techniques，2008，27（2）：1－4（in Chinese）．

［10］童文辉，沈峰满，王文忠，等．高炉常用耐火材料导热系数的测定［J］．金属学报，2002，38（9）：983－988．TONG Wenhui，SHEN Fengman，WANG Wenzhong，et al．Measurement of heat－transfer coefficient of refractory for BF［J］．Acta Metallurgica Sinica，2002，38（9）：983－988（in Chinese）．

［11］王东，刘宗奎，张斌，等．激光闪射法测试氮化硅结合碳化硅材料导热系数的研究［J］．耐火与石灰，2012，37（6）：3－6．WANG Dong，LIU Zongkui，ZHANG Bin，et al．Study on coefficient of thermal conductivity of silicon carbide material bonded with silicon nitride determined with laser flash method［J］．Refractories and Lime，2012，37（6）：3－6（in Chinese）．

［12］易法军，张巍，孟松鹤，等．C／C复合材料高温热物理性能实验研究［J］．宇航学报，2002，23（5）：85－88．YI Fajun，ZHANG Wei，MENG Songhe，et al．An experimental study on thermo physical properties of C／C composites at elevated temperature［J］．Journal of Astronautics，2002，23（5）：85－88（in Chinese）．

［13］邹林华，黄伯云，黄启忠，等．C／C复合材料的导热系数［J］．中国有色金属学报，1997，7（4）：135－138．ZOU Linhua，HUANG Boyun，HUANG Qizhong，et al．Thermal conductivity for C／C composites［J］．The Chinese Journal of Nonferrous Metals，1997，7（4）：135－138（in Chinese）．

［14］程远贵，董岱峰，周勇，等．耐火纤维毡的高温导热系数研究［J］．化工设计，2001，11（2）：22－23．CHENG Yuangui，DONG Daifeng，ZHOU Yong，et al．A study on measuring thermal conductivities of Al2O3fibrous refractories under high temperature［J］．Chemical Engineering Design，2001，11（2）：22－23（in Chinese）．

［15］胡国俊．传热中的多宗量反演研究［D］．大连：大连理工大学，2002．

［16］STOLZ G，Jr．Numerical solutions to an inverse problem of heat conduction for simple shapes［J］．Journal of Heat Transfer，1960，82：20－25．

［17］ZUECO J，ALHAMA F，FERN A，et al．An inverse problem to estimate temperature dependent heat capacity under convection processes［J］．Heat and Mass Transfer，2003，39（7）：599－604．

［18］HUANG C，CHEN C．A three－dimensional inverse geometry problem in estimating the space and time－dependent shape of an irregular internal cavity［J］．International Journal of Heat and Mass Transfer，2009，52（7）：2 079－2 091．

［19］王一博，杨海天，邬瑞锋．基于时域精细积分算法的瞬态传热多宗量反演［J］．应用数学和力学，2005，26（5）：512－518．WANG Yibo，YANG Haitian，WU Ruifeng．Precise integral algorithm based solution for transient inverse heat conduction problems with multi－variables［J］．Applied Mathematics and Mechanics，2005，26（5）：512－518（in Chinese）．

［20］姜尧忠．工业电炉［M］．北京：清华大学出版社，1993．

［21］杨晨，Ulrich Gross．基于热传导逆问题方法预测材料热物性参数［J］．化工学报，2005，56（12）：2 415－2 420．YANG Chen，Ulrich Gross．Estimation of thermal properties based on inverse heat conduction method［J］．Journal of Chinese Industry and Engineering（China），2005，56（12）：2 415－2 420（in Chinese）．

［22］张有为，李辉，姜培学，等．采用共轭梯度法的管内壁温度导热反问题求解［J］．工程热物理学报，2009，30（7）：1 188－1 190．ZHANG Youwei，LI Hui，JIANG Peixue，et al．Inverse heat conduction problem of deducing inner wall temperature by using conjugate gradient method［J］．Journal of Engineering Thermophysics，2009，30（7）：1 188－1 190（in Chinese）．

［23］LUGON J，Jr，Neto SILVA A J，SANTANA C C．A hybrid approach with artificial neural networks，Levenberg－Marquardt and simulated annealing methods for the solution of gas－liquid adsorption inverse problems［J］．Inverse Problems in Science and Engineering，2009，17（1）：85－96．