酚醛树脂气凝胶复合材料热物性参数预测方法

张春云，陈雄斌，刘健，崔苗，*

1.大连理工大学 工业装备结构分析优化与CAE 软件全国重点实验室，大连 116024

2.大连理工大学 辽宁省空天飞行器前沿技术重点实验室，大连 116024

3.北京空天技术研究所，北京 100074

酚醛树脂基热防护材料［1-3］应用比较成熟，但是密度较大，导热系数仍相对较高，对于长时间热防护应用，厚度和重量难以承受。酚醛树脂气凝胶复合材料（Inorganic-Organic Hybrid Phenolic Aerogel Composites，IPC）是近年来发展的一种新型轻质高效隔热材料［4］，其抗氧化能力强，残碳率高，非常适合长时间热防护应用，发展前景广阔。但是与酚醛树脂基热防护材料相似，IPC 在400 ℃以上出现热解碳化，伴随这一复杂的物理化学过程，材料的导热系数等热物性参数会发生明显的变化。由于裂解产物（油气）对设备的影响，很难通过直接测量［5-6］获取材料的高温热物性，只能进行粗略估计，这降低了温度场预测的精度，进而影响热防护结构的设计与优化。因此，获取较为准确的IPC 随温度变化的热物性参数对于热防护结构的精细化设计具有重要意义。

酚醛树脂气凝胶复合材料在碳化过程中发生热解反应，密度、导热系数和比热容等物性参数随材料成分的变化而变化；并且其传热等物理化学模型的建立与求解非常具有挑战性。因此，在工程分析中通常采用等效导热系数［7］等概念，将IPC 碳化裂解以及热物性参数随温度变化等影响都考虑到等效热物性参数中。在不建立化学反应方程的条件下，通过在热传导方程中增加源项等方式模拟化学反应的影响，将微观量转化为宏观量，以满足精度需要。

反问题［8］为酚醛树脂气凝胶复合材料热物性参数的准确确定提供了新思路。反问题一般是基于可测/易测的物理量测量信息来确定系统的结构特征或输入信息，在航空航天、损伤识别、冶金工程、核反应堆等领域广泛应用［9-10］。周焕林等［11］结合边界单元法和梯度正则化算法确定了随温度变化的导热系数，并且证明了该算法具有较好的精度和稳定性。Tahmasbi 和Noori［12］采用Levenberg-Marquardt（LM）算法确定了正交碳化复合材料纤维方向和厚度方向的导热系数。Liu 等［13］采用遗传算法反演纳米多孔介质二氧化硅气凝胶复合材料随温度变化的等效导热系数，并与实验测得的结果对比，证明了该方法的准确性和鲁棒性。姜贵庆等［14］采用Newton-Raphson 算法通过预测3 个待定参数确定了防热涂层材料碳化层的导热系数。Xie 等［15］采用Levenberg-Marquardt 算法预测了硅气凝胶材料的等效导热系数。张红军等［16］基于共轭梯度法，通过求解烧蚀传热方程，确定了纳米酚醛气凝胶材料原始层和碳化层的导热系数。林旭文等［17］基于动态贝叶斯网格对烧蚀热防护材料的导热系数进行了预测。Wang 等［18］基于蒙特卡罗法确定了烧蚀材料原始层和碳化层随温度变化的导热系数。

通过分析前人的研究工作，可以看出，对于近些年新开发的酚醛树脂气凝胶复合材料，关于其热物性参数确定的报道较少，缺乏数据库，没有通用性，不利于工程应用；对于烧蚀型复合材料，现有的研究存在部分温度区间热物性参数预测精度低和计算效率低的问题，不便于高超声速飞行器热防护结构的整体性能评估或者优化设计。

酚醛树脂气凝胶复合材料热物性参数的预测可以采用随机法［19］和梯度法［20-21］。随机类算法的优点是全局收敛性好，但是计算效率比较低，而梯度算法不需要进行全局搜索，计算效率和精度高，更适用于工程应用。然而，对于梯度法而言，灵敏度矩阵的准确确定至关重要。复变量求导 法（Complex Variable-Differentiation Method，CVDM）［22］在计算一阶偏导数时具有高的精度和好的通用性，因此，本文采用复变量求导法计算灵敏度系数，保证预测结果的准确性。同时，在求解不适定问题时，灵敏度矩阵的准确确定避免了迭代过程中目标函数的震荡，提高了计算效率。

基于上述分析，本文首先建立求解瞬态非线性热传导反问题的改进Levenberg-Marquardt 梯度算法；然后，采用该算法预测随温度变化的热物性参数，验证该算法的有效性，并研究测量误差和初始猜测值对反演结果的影响；最后，开展酚醛树脂气凝胶复合材料测温试验，并基于试验测量数据预测IPC 碳化过程随温度变化的等效导热系数。

1 IPC 热物性参数预测的理论基础

1.1 瞬态传热问题

预测酚醛树脂气凝胶复合材料随温度变化的热物性参数需要求解瞬态传热问题。瞬态非线性热传导控制方程的差分表达式［23］为

初始条件为

边界条件表示为

以对流换热边界条件为例，边界节点的离散方程为

式中：上标i为节点编号；下标t为时刻；Δt为时间步长；λ为导热系数，W·m-1·℃-1；T为温度，℃；q为热流密度，W·m-2；Δx为空间步长；T0为初始温度；Q为体热源，W·m-3；ρ为密度，kg·m-3；c为比热容，J·kg-1·℃-1；t为时间，s；Tw为壁面 温度，℃；Tf为流体温度，℃；h为对流 换热系数，W·m-2·℃-1；a为热扩散率，m2·s-1；n为法线方向。

1.2 基于梯度算法的热物性参数预测

基于试验测点温度信息，采用改进的Levenberg-Marquardt 梯度算法，预测酚醛树脂气凝胶复合材料随温度变化的热物性参数。在LM 算法中，瞬态热传导反问题的无量纲优化目标函数［24］表示为

式中：Ti和分别为测点物理量的计算值和测量值；M为测量物理量的个数；Y为预测的热物性参数，可为随温度变化的导热系数、比热容或热扩散率，Y=[ ]Y1，Y2，…，YN，其中N为预测物性参数的数量。

预测的热物性参数根据式（8）进行迭代计算：

式中：w为收敛因子，设置为1；P为迭代次数。此外，热物性参数的更新值δ由式（9）确定：

式中：μ为阻尼因子，阻尼因子是Levenberg-Marquardt 算法中的关键参数，采用无量纲目标函数的幂指数形式，随着迭代次数增加不断更新［24］：

其中：指数n=1.5。J为灵敏度矩阵，具体形式为

为了准确计算灵敏度矩阵，保证算法的精度，采用复变量求导法［25］计算灵敏度系数。其偏导数的计算与空间步长大小无关，因此求解精度高。在复变量求导法中，以x+ih代替实函数f(x)中的变量x，h的取值很小，一般取10-30。将f(x+ih)展开成泰勒级数的形式为

因为h的取值很小，所以f（x）的一阶导数为

式中：Im 代表虚部。

当无量纲优化目标函数值很小或者目标函数的变化很小时，迭代停止，即

式中：收敛容差ζ为一个小的正整数，一般为1×10-6。

该算法求解非线性热传导反问题的流程如下：

1）读取待预测参数的初始猜测值，对于提出的算法，由于其是一种梯度类算法，算法的收敛性依赖于初始猜测值，因此，初始猜测值应尽量接近真实值，然后设置收敛容差、收敛因子以及测量温度信息等。

2）基于几何条件、边界条件以及初始猜测值，采用有限差分法求解瞬态非线性热传导问题，得到测点的计算温度值，同时得到由复变量求导法计算的灵敏度矩阵系数。

3）根据式（7）计算无量纲优化目标函数S，并通过式（14）判断目标函数是否收敛。如果满足，则迭代停止，输出参数预测值和优化目标函数值；如果不满足，则继续执行下一步计算。

4）基于目标函数值计算阻尼因子μ，并根据式（9）计算热物性参数的更新值δ。

5）基于式（8）更新随温度变化的热物性参数，并更新迭代步数P，然后返回步骤2）进行循环。

热物性参数预测的流程图如图1 所示。

图1 预测算法流程图Fig.1 Flowchart of prediction algorithm

2 导热系数预测验证算例

为了确保该算法能够有效地预测IPC 随温度变化的等效导热系数，验证该梯度算法求解瞬态非线性热传导正反问题的有效性和精度。

2.1 瞬态非线性热传导问题

以一个一维模型为例，其几何模型和边界条件如图2 所示。该模型长10.0 mm，两端为随时间变化的温度边界条件，如表1 所示。材料热物性参数随温度变化，如表2 所示。初始温度为20 ℃，计算时间为100 s。该瞬态热传导问题采用有限差分法进行求解。

表1 验证算例的边界条件Table 1 Boundary conditions of validation model

表2 验证算例的材料热物性参数Table 2 Thermo-physical properties of validation model

图2 验证算例的几何模型和边界条件Fig.2 Geometry and boundary conditions of validation model

分别采用有限单元法（Finite Element Method，FEM）和有限差分法（Finite Difference Method，FDM）对瞬态非线性热传导问题进行求解，计算结果如图3 所示。

图3 节点数量对计算结果的影响Fig.3 Effects of number of nodes on calculation results

由图3 可以看到，采用有限差分法计算的3.0、6.4、6.6 mm 3 个位置处的温度值与有限单元法的计算结果吻合较好，51、101、201 和3014 组节点的计算温度值与有限单元法计算结果的平均相对误差分别为 0.026%、0.010%、0.017%和0.018%。可以看出，对于该瞬态非线性热传导问题，只要满足收敛稳定性条件，节点数量对结果精度影响较小。

在保证有限差分法满足收敛稳定性条件的前提下，选取101 个节点，分析差分步长对计算精度的影响。计算结果如图4 所示。可以看出，采用不同差分步长计算的温度值与有限单元法计算的结果吻合较好，当差分步长为0.000 2、0.001 0 和0.005 0 时，其计算温度与FEM 计算结果的平均相对误差分别为0.009%、0.010%和0.017%。可以得出，计算误差随着差分步长的增加而增加。综合考虑计算精度和效率，选取0.001 0 的差分步长进行计算。

图4 差分步长对计算结果的影响Fig.4 Effects of difference step on calculation results

2.2 随温度变化导热系数的预测

为了验证提出的预测算法的优势，采用该方法预测随温度变化的导热系数，并与传统LM 算法进行比较。对于此类问题，由于无法推导出灵敏度系数的解析表达式，在传统LM 算法中，通常采用差分格式来确定灵敏度矩阵：

以该模型3.0、6.4、6.6 mm 处的正问题温度计算值作为测量信息，每2 s 读取一次温度，则一共有150 个测量信息。求解瞬态热传导反问题时，需要预测7 个温度处的导热系数，其他计算条件同正问题一致。数值算例均在CPU Intel®CoreTMi5-8400、2.80 GHz、内存16 GB 的计算机上完成。计算得到的导热系数收敛曲线如图5 所示。本文算法与传统Levenberg-Marquardt 算法的对比如表3 所示。

表3 不同预测算法计算结果比较Table 3 Comparison of computational results with different prediction methods

图5 验证算例导热系数的收敛过程Fig.5 Convergence process for thermal conductivities of validation model

由图5 可以看到，5 次迭代后，7 个参数收敛至各自的真实值，这说明该方法在预测热物性参数方面具有较高的效率和精度。

由表3 可以看出，当步长较小时，如1×10-16，传统的梯度算法不收敛，并且步长影响参数预测的迭代步数和目标函数。而本文算法在步长为1×10-30时仍能收敛，迭代步数不随步长的变化而变化，目标函数随步长的变化较小，说明该算法受步长影响小，具有很好的稳定性和精度。在实际应用中，该算法不需要调整步长，通用性好，具有可推广性。此外，由表3 还可以看到，本文算法的计算时间均比传统算法少，说明该算法效率高。这是因为采用差分格式计算灵敏度系数时，反复计算温度场，增加了计算量。

2.3 测量误差的影响

在试验过程和实际工程应用中，温度的测量必然存在着误差。反问题具有不适定性，很小的测量误差都可能使预测参数产生很大的波动。为了研究测量误差对算法的影响，定义平均相对误差Ere为

式中：N为预测参数的个数；为预测的热物性参数的真实值。

存在随机测量误差的测量温度［24］可以表示为

式中：Tex为真实的测点温度；Tes为存在测量误差的温度；γ为-1～1 之间的随机数；ξ为测量误差。

表4 给出了施加不同测量误差时，随温度变化的导热系数的预测结果。可以看到，测量误差影响预测算法的精度，测量误差越大，该算法预测精度越低。但是，当测量误差为3%时，导热系数预测的平均相对误差为1.448%，仍在合理范围之内。这说明该算法受测量误差影响较小，采用该方法确定随温度变化的导热系数具有高的精度。此外，由表4 还可以看到，平均相对误差均小于各自的测量误差，说明了该预测算法鲁棒性好。

表4 测量误差对预测结果的影响Table 4 Effects of measurement errors on prediction results

2.4 初始猜测值的影响

相比于随机算法，梯度算法的计算结果受初始猜测值影响较大。为了研究初始猜测值对随温度变化的导热系数预测的影响，选取3 组不同的初始猜测值进行反演计算，其他计算条件均相同。随温度变化的导热系数的预测结果如表5 和图6 所示。可以看出：当导热系数的初始猜测值取0.1、0.5、1.0 时，预测的导热系数均可以收敛到真实值，这说明该预测算法受初始猜测值影响较小，具有很好的精度和稳定性。此外，由表5 进一步分析可以看出，初始猜测值与导热系数的真实值越接近，迭代次数越少，预测算法的收敛速度越快。

表5 初始猜测值对预测结果的影响Table 5 Effects of initial guessed values on prediction results

图6 初始猜测值不同时导热系数的预测值Fig.6 Predicted thermal conductivities with different initial guessed values

3 试验结果与分析

开展了IPC 测温实验，并基于试验测量温度值预测酚醛树脂气凝胶复合材料随温度变化的等效导热系数。

3.1 IPC 测温试验

试验件为15 mm 厚的IPC 防热层和8.6 mm厚的铝合金背板的双层平板，尺寸为300 mm×300 mm。采用石英灯对酚醛树脂气凝胶复合材料表面进行加热。虽然未能考虑到实际工程应用中边界气体对流及热阻塞效应的影响，但是，IPC 是微烧蚀材料，热失重非常小，质量引射对降热的效果较弱。此外，热阻塞效应主要发生在边界层，使作用到表面的冷壁热流减小，而不改变材料内部的导热系数，因此，其对IPC 等效导热系数的预测结果几乎无影响。

在该试验件四周包裹隔热材料，以减少试验件对环境的辐射散热和对流换热，如图7 所示。石英灯加热热流如图8 所示，铝合金板背面为自然对流边界条件。图9 为该测温试验的物理模型。

图7 试验示意图Fig.7 Schematic diagram of experiment

图8 热流边界条件Fig.8 Heat flux boundary condition

图9 试验模型图Fig.9 Experiment model diagram

试验持续时间为750 s，初始温度为25 ℃。该试验酚醛树脂气凝胶复合材料的密度为640 kg·m-3，室温比热容为950 J·kg-1·℃-1。在试验中，该试验件表面温度采用多个热电偶和红外测温仪测量。同时，为了获取IPC 内部测点随时间变化的温度数据，在该试验件内部布置4 处热电偶，热电偶距加热面的距离分别为3.0、6.4、6.6、10.0 mm。处理后的测点随时间变化的温度如图10 所示。可以看到受热面附近的温度曲线中间部分下降后上升，这个现象的主要原因是石英灯加热热流的变化，如图8 所示，受热面附近的温度变化趋势与热流变化趋势一致。当节点向外输出的热量大于吸收的热量时，温度下降；当节点向外输出的热量小于吸收的热量时，温度上升。

图10 测点温度分布Fig.10 Temperature distributions of measured points

3.2 IPC 碳化过程等效导热系数预测

在该IPC 随温度变化的等效导热系数预测时，以受热面已知测点温度的平均值作为受热面的边界条件，10.0 mm 处的测量温度值作为下表面边界条件，3.0、6.4、6.6 mm 处的温度值作为测量数据。正问题计算过程中，差分步长设置为0.001 0，总节点数为101，满足稳定性条件。

预测酚醛树脂气凝胶复合材料随温度变化的等效导热系数，在整个温度区间范围内采取分段线性插值，选取100～1 100 ℃的温度点，间隔为100 ℃。11 个参数的初始猜测值均为0.03 W·m-1·℃-1。在3 个测点位置分别选取150 个测量数据，则一共有450 个测量数据。IPC 随温度变化的等效导热系数的预测结果如表6 和图11 所示。可以看到，该试验IPC 的等效导热系数随着温度的升高而增加，并且在800 ℃以下，等效导热系数均低于0.1 W·m-1·℃-1，预测结果符合物理规律。

表6 预测的典型温度点的等效导热系数Table 6 Inverted effective thermal conductivities for typical temperature points

图11 等效导热系数随温度的变化规律Fig.11 Variation of effective thermal conductivity with temperature

预测过程中，优化目标函数的变化过程如图12 所示，可以看到，经过11 次迭代，该预测算法收敛，并且优化目标函数下降很快。这说明该预测算法在确定IPC 碳化过程等效导热系数时效率高。

图12 目标函数的变化过程Fig.12 Change process for objective function

图13 比较了试验测量温度和采用预测的等效导热系数计算的温度值。可以看到，3 个测点处温度的测量值与计算值吻合较好，测量数据与计算温度值的平均相对误差为3.279%，说明了本文算法精度高，预测结果的可信度高。

图13 IPC 测量温度与计算温度的对比Fig.13 Comparison of measurement values with calculated temperatures for IPC

4 结论

建立了酚醛树脂气凝胶复合材料（IPC）碳化过程热物性参数的预测算法。通过与传统梯度算法比较，表明了本文算法的良好性能。开展了IPC 测温实验，基于试验测量数据，对IPC 随温度变化的等效导热系数进行了预测。主要研究结论如下：

1）该算法受步长影响较小，具有很好的稳定性和精度。对于所研究的算例，该算法在步长为1×10-30时仍能收敛，且迭代步数和目标函数受步长影响小，而传统的梯度算法在步长为1×10-16时就无法收敛。

2）相比于传统的梯度算法，该算法的计算时间明显缩短，针对所研究的算例，计算时间由75 s下降至35 s，说明该算法效率高。

3）采用本文算法预测的酚醛树脂气凝胶复合材料的等效导热系数随着温度的增加而增加，试验测量值和计算温度吻合较好，平均相对误差为3.279%。预测结果符合物理规律，且精度高。

该工作为获取热防护材料随温度变化的高温热物性参数提供了一种有效方法，为热防护结构传热分析的整体优化设计计算提供了关键参数。