飞行器在再入大气层过程中，高速气流流过飞行器表面，由于气体粘性的阻滞作用，会带来气动加热问题，使飞行器表面温度明显升高，以美国阿波罗飞船指令舱和前苏联联盟号飞船返回舱为例，再入时的最高表面温度分别达到了2950K和1900K，这会给飞行安全带来影响。因此，再入飞行器需进行防热系统设计。烧蚀防热是目前主要采用的防热系统设计方法，即通过防热材料的热解、烧蚀后退来吸收并带走气动加热的热量，从而达到热防护的目的[1, 2]。
表面热流的辨识，或称为表面热流的反演，是指在飞行器表面温度较高或表面有烧蚀后退的情况下，不适宜在表面直接安装传感器来进行热流和温度测量，只能在飞行器的防热层内某些点上安装温度传感器进行温度测量，然后根据这些点上的温度测量历程来反算飞行器的表面热流。这类问题是一类典型的热传导逆问题IHCP(Inverse Heat Conduction Problem)，是传热学研究的一个重要领域。目前，针对固定几何域的热传导逆问题研究较多，但由于防热层的烧蚀后退对应着防热层几何外形的变化，因而此时防热层内的热传导逆问题是一个变几何域的热传导逆问题，研究工作相对较少，尤其是在防热层烧蚀后退过程中还伴有材料热解反应、热解气体流动等复杂物理化学过程，使得此时的逆问题更为复杂。文献[3]中，Kanevce利用地面烧蚀试验的温度测量数据对酚醛材料的热物性参数和材料热解反应系数进行了辨识；文献[4, 5]建立了考虑防热材料烧蚀情况下的表面热流反演辨识算法，但文献[4]采用的是一层的简化烧蚀模型，并且通过给定烧蚀表面温度来确定外表面位置，这样处理会导致辨识结果出现非物理振荡；文献[5]采用了两层烧蚀热解面模型，建立了相应的表面热流辨识算法，但未用实测数据进行验证。本文在前期变几何域表面热流辨识方法[6]研究的基础上，借鉴文献[5]的方法，采用两层烧蚀热解面模型，实现了表面热流辨识算法，并用典型试验数据对辨识算法进行了验证。 1 烧蚀模型
防热材料的烧蚀是一个复杂的物理化学过程。本文针对目前工程上通常采用的碳化复合材料开展研究。碳化复合材料的烧蚀模型主要有热解层模型、热解面模型和简化模型，本项工作中采用的是热解面模型[1]，如图 1示，左边界为材料受热边界，右边界为防热材料内部边界，即图中的s3。热解面模型是将防热材料分为两层：碳化层和原始材料层，热解层则简化为一个热解面进行处理。图 1中给定碳化边界和外边界烧蚀后退规律，此时的传热区域包含两个运动边界：一个是表面烧蚀后退的边界，即图中的s1；一个是碳化边界，即图中的s2，对应的数学模型为：
 |
| 图 1 热解面模型Fig. 1 Pyrogeneration-plane ablation model |
。其中，下标“c”表示碳化材料的参数；“v”表示原始材料的参数；“g”表示裂解气体的参数；
为裂解气流动影响；
cΔHc为碳化层烧蚀潜热；pΔHp为裂解化学反应吸收热；
本文采用变几何域的有限控制体积法(FCV: Finite Control Volume method)[7]来对式(1)进行数值求解，式中的s1(t)和s2(t)可通过两种方法来给出，一种是给定为测量结果；另一种是给定烧蚀碳化温度后通过迭代计算来确定其位置。同时，此时由于计算域几何外形随时间变化，计算网格相应随时间变化，因而在离散方程中需考虑由于控制体边界运动而引起的控制体内能变化项。该计算方法处理变几何域传热问题的有效性在文献中已得到了有解析解的一维半无限区域烧蚀表面恒温恒速后退算例的验证。 2 考虑烧蚀情况下的表面热流辨识方法
通常的烧蚀问题计算是在表面热流Q(t)已知的情况下，求解式(1)得出材料表面的后退率和内部温度分布。而表面热流辨识问题则是在表面热流Q(t)未知的情况下，根据材料端面或内部测点的温度历程来反演计算出表面热流Q(t)的值。不失一般性，下面将测量点取在右端面处，建立由测量信息来辨识表面热流Q(t)的方法。首先，利用拉格朗日(Lagrange)乘数法将辨识问题转化为使选取表面热流Q(t)值使如下目标函数达极小的优化问题：
式中的λ为伴随变量。对此式做分部积分，代入式(1)中的相应边界条件后，再做变分，可知为使式(2)中的目标函数达极小，伴随变量需满足如下伴随方程： 而目标函数对Q的梯度为：
在此基础上即可采用梯度类优化算法，如最速下降法和共轭梯度法[6, 8]，来对Q(t)进行辨识计算，以最速下降法为例，优化迭代计算公式为：
式中上标“l”和“l+1”表示迭代步数；下标“i”表示时间方向上的离散；β为步长。在优化计算过程中，为克服不适定性的影响，优化计算设置了如下停止准则：J≤δ； δ=σ2tf； σ为温度测量结果的标准差。 3 验证算例
下面进行算法验证：设烧蚀材料的原始材料层ρ=1300kg/m3，比热Cp=1507J/kgK，kv= 1W/(mK)，初始厚度为L=20mm；碳化后的ρ=406.3kg/m3，比热Cp=2000J/kgK，kc= 3W/(mK)；裂解化学反应吸收热ΔHp= 831500 J/kg；热解面温度为873K；裂解气比热Cg=1300J/kgK。给定表面热流如图 2中的“Qexact”，给定烧蚀端面s1(t)的后退历程如图 3中的“s1 without noise”，采用第1节中的数值方法对温度场进行求解，s2(t)根据热解面温度通过迭代计算来确定，图 4示出了计算出的热解面位置。
 |
| 图 2 表面热流给定值与辨识结果对比Fig. 2 Comparison of exact and estimated value of heat flux |
 |
| 图 3 烧蚀后退量S1随时间变化历程Fig. 3 Time history of ablated surface recession |
 |
| 图 4 热解面时间变化历程Fig. 4 Time history ofcalculated pyrogeneration-plane |
以此时计算出的内壁测点温度历程为实测值，采用第2节介绍的方法对表面热流进行辨识，图 2中的虚线“Qestimated(σ=0K)”示出了此时的辨识结果，与给定值符合较好，但在末时刻有差异，这是由于传热存在一定的滞后性，末段时刻的表面热流信息在内壁温度测量结果中还没有体现出来，所以不能有效辨识。
接下来分析测量噪声的影响，在计算出内壁温度历程上叠加标准差σ=6K(对应测量的相对误差约为1%)的测量白噪声来作为测量结果(图 5中“exp.”)，对表面热流进行辨识，图 2中的虚线“Qestimated(σ=6K)”示出了此时的辨识结果，与给定值仍符合较好；图 5中的“Fitted”给出了利用辨识结果计算出的测点温度历程，和测量值符合较好。进一步将叠加的测量白噪声标准差由σ=6K增大到σ=12K(对应2%的测量误差)，对表面热流进行辨识，图 2中的虚线“Qestimated(σ=12K)”示出了此时的辨识结果，可以看到：在相对误差2%的情况下，表面热流的辨识结果与给定值进一步偏离，但仍较好地反映出了表面热流的变化趋势与强度，辨识算法具有较好的鲁棒性。
 |
| 图 5 测点温度对比Fig. 5 Comparison of temperature measurement |
在实际应用中，烧蚀后退量的测量通常有较大的测量误差，因而在图 3中s1(t)上叠加σ=0.11mm(相当于烧蚀后退量最大值的10%)的白噪声作为测量噪声，得到的测量结果如图 3中的“s1 with noise(σ=0.00011m)”示。利用不考虑s1(t)测量噪声时的测点温度计算值作为实测值(不叠加测量噪声)，对表面热流进行辨识，得到的辨识结果与给定值的比较如图 6示，从图中可以看到，由于烧蚀后退量测量结果存在明显波动，表面热流的辨识结果中也出现了明显波动，二者之间存在较直接的关联关系。
 |
| 图 6 表面热流给定值与辨识结果对比Fig. 6 Comparison of exact and estimated value of heat flux |
下面通过一地面试验的实例来对辨识方法的有效性做进一步分析验证。钝头型碳酚醛材料(Narmco4028)试件在Langley中心11英寸陶瓷加热风洞中的烧蚀试验示意图[9]如图 7所示，材料的几何参数在图中标注，来流中氧的质量分数为0.13，驻点压力5.88atm，驻点焓值2.55MJ/kg。图 8为文献[9]中给出的烧蚀表面、碳化边界的后退历程计算值(图中分别记为"front surface"和"interface")和表面温度测量结果(图中记为"front surface temperature")，其中烧蚀表面和碳化边界在试验结束时刻的计算结果与实测值吻合。基于这组结果，利用前面建立的辨识算法，可辨识出表面加热热流的时间变化历程，如图 9中“Qestimated”示。可以看到，热流辨识结果和陶瓷加热风洞的加热热流3.43MW/m2是基本符合的，偏差约20%。这一结果表明：本文所建立的考虑烧蚀情况下表面热流辨识方法是基本有效的，经过更多算例的检验校核后有望在工程实际中得到更进一步的应用。
 |
| 图 7 钝头型碳酚醛材料地面烧蚀试验Fig. 7 Ablation experiment of blunt Carbon-phenonic material |
 |
| 图 8 烧蚀表面、碳化边界计算结果和表面温度试验结果Fig. 8 Calculated front surface and interface of ablation and experimental surface temperature |
 |
| 图 9 表面热流时间历程辨识结果和加热热流值对比Fig. 9 Comparison of estimated heat flux and heating power |
本文针对烧蚀传热问题，在热解面模型的基础上，推导了伴随方程，建立了基于测点温度辨识表面热流的方法，并对算法的有效性进行了算例考核。结果表明：在不考虑测量噪声时，热流辨识结果与真值符合较好；随着测量误差的增大，辨识结果与真值之间的偏差增大；烧蚀后退量测量结果的误差对辨识结果有较为显著的影响。此后，将该辨识方法用于钝头型碳酚醛材料Narmco4028试件在陶瓷加热风洞中的烧蚀试验结果的分析，结果表明辨识出的表面热流与加热功率基本符合，辨识方法是有效的，在工程实际中有较好的应用前景。
| [1] | 姜贵庆, 刘连元. 高速气流传热与烧蚀热防护. 北京: 国防工业出版社, 2003. |
| [2] | 王希季, 林华宝, 李颐黎, 等. 航天器进入与返回技术(下). 北京: 宇航出版社, 1991. |
| [3] | KANEVCE L P, KANEVCE G H, ANGELEVSKI Z Z. Comparison of two kinds of experiments for estimation of thermal properties of ablative composite. In: Inverse Problems in Engineering: Theory and Practice. 3rd Int. Conference on Inverse Problems in Engineering. Port Ludlow, WA, USA, 1999. |
| [4] | OLIVEIRA A P D, ORLANDE H R B. Estimation of the heat flux at the surface of ablating materials by using temperature and surface position measurements. Inverse Problems in Science and Engineering, 2004, 12(5): 563-577. |
| [5] | HAKKAKI-FARD A, KOWSARY F. Heat flux estimation in a charring ablator. Numerical Heat Transfer, Part A, 2008, 53: 543-560. |
| [6] | 邵元培, 钱炜祺, 周宇, 等. 变几何域传热的表面热流反演方法. 计算力学学报, 2013, 30(2): 296-301. |
| [7] | AMAR A J, BLACKWELL B F, EDWARDS J R. One-dimensional ablation using a full Newton’s method and finite control volume procedure. Journal of Thermophysics and Heat transfer, 2008, 22(1): 71-82. |
| [8] | 钱炜祺, 周宇, 何开锋, 等. 非线性热传导逆问题的表面热流辨识方法. 空气动力学学报, 2012, 30(2): 145-150. |
| [9] | SUTTON K. An experimental study of a carbon-phenolic ablation material. NASA TND-5930, 1970. |