非线性热传导逆问题的表面热流辨识方法

钱炜祺， 周 宇， 何开锋， 袁军娅， 黄建栋

（1.中国空气动力研究与发展中心 空气动力学国家重点实验室，四川 绵阳 621000；2.中国运载火箭技术研究院，北京 100076）

0 引 言

热传导逆问题，或称为热传导反问题（IHCP：Inverse Heat Conduction Problem），是利用实验手段测得材料内部或边界上某点或某些点上的温度随时间变化历程，通过求解传热微分方程来反演边界热流、材料热传导系数或材料内部热源分布等参数。热传导逆问题在航天、核物理、冶金等工业研究领域中有着广泛的应用背景［1］，是传热学研究的主要领域之一。现有大多数热传导逆问题研究针对的都是线性热传导方程，即材料的热物性参数，包括热传导系数和比热是常数的情况。而在实际工程应用中，很多材料的热物性参数是随温度变化的，此时的热传导方程是非线性的偏微分方程，对应的热传导逆问题也称为非线性热传导逆问题［2－3］。非线性热传导逆问题主要包括两方面研究内容：一是边界热流已知，根据材料内部或边界温度场的测量结果来辨识材料热物性参数随温度的变化规律；二是材料热物性参数随温度的变化规律已知，需要根据温度的测量结果来辨识表面热流。对于第一类问题，其求解方法在文献［4－5］中有详细介绍，在此不做讨论；对于第二类问题，文献［3］采用顺序函数法来进行了表面热流辨识，文献［7］给出了伴随方程的形式，但未给出其具体推导。相比之下，国内对这一问题的研究开展得相对较少，限制了表面热流辨识方法在工程实际中的使用范围。因此，本文工作将以一维非线性热传导逆问题为研究对象，系统地建立起两类表面热流辨识方法：顺序函数法和共轭梯度法，并结合算例来对算法的有效性和鲁棒性进行分析。

1 非线性热传导逆问题的数学模型

典型的一维传热如图1示，厚度为L的材料上部受表面热流Q（t），下部绝热，温度传感器置于测点P处，此时的热传导方程为：

初始条件：t＝0：T＝T0。

观测方程：

其中，ρ为材料密度，k（T）为材料热传导系数，Cp（T）为材料比热，都是随温度T变化的函数，xm为测点位置，v（t）为测量噪声。此时的热传导逆问题就是在k（T）和Cp（T）函数已知、表面热流Q（t）未知的情况下，由观测方程式（2）中的信息来反演辨识表面热流Q（t）。

图1 一维传热示意图Fig.1 Sketch of 1Dheat conduction problem

该逆问题等价于求合适的Q（t）使如下目标函数达极小的优化问题：

其中，t＝［0，tf］表示温度测量的时间段，T（xm，t，Q）表示表面热流为Q时测点温度历程的计算值；为相应的实测值。由于Q（t）是时间的函数，因此式（3）实际上是一个泛函表达式，考虑到数值求解式（1）时对时间的离散，此时的待优化变量为各时刻的热流值Qi，即

其中Nt为数值计算中时间方向的离散层数。则式（3）的优化变成了一参数优化问题，下面介绍两种优化策略：顺序函数法（SFM：Sequential Function specification Method）和共轭梯度法（CGM：Conjugate Gradient Method）。

2 顺序函数法

针对式（1）描述的表面热流辨识问题，当用顺序函数法来辨识未知热流Q（t）时，问题可描述为：已知tM时刻的前M－1个时刻处的热流1，…，－1和后r个时刻的温度值TM，TM＋1，…，TM＋r－1，要求辨识tM时刻的热流QM。例如，当M＝1时，即需已知初始热流Q0及后r个温度测值T1，T2，…，Tr，再由Q0，T1，T2，…，Tr的信息来确定Q1。由于热传导的延迟性，QM决定tM时刻之后时间里的温度值，即TM＋r－1不 仅 取 决 于QM，还 取 决 于QM＋1，…，QM＋r－1，因此，在辨识QM时，还必须建立QM＋i（i＝1，r）与QM的关系。若设QM－1至QM＋r－1是线性变化的，当时间间隔是常数时，则有：

因此，辨识QM的状态方程为：

观测方程为（设采样为等时间间隔）：

辨识的目标函数取为：

视QM为参数，用灵敏度法来辨识，则其迭代的牛顿－拉夫逊（Newton－Raphson）算式为：

式中Q的上标k、k＋1表示迭代层次，∂T／∂QM为灵敏度，其满足的方程可近似写为：

综上所述，用顺序函数法辨识热流时，是在时间方向上按顺序向后辨识。当M＝1时，Q0已知，方程（6）中的T′（x）即为初始时刻已知的全场温度分布，由初估的Q1及式（6）可求出后r个时刻的温度测值T（xm，ti）（i＝1，…，r），再由式（10）求出灵敏度∂T（xm，ti）／∂QM（i＝1，…，r）后，用式（9）来对QM优化，优化得QM的最优估值则取为Q1的辨识值。然后，再沿时间方向推进，用及T2，…，Tr＋1按照同样的方法来辨识Q2，依此类推。

在上述辨识过程中，涉及到状态方程式（6）和灵敏度方程式（10）的求解，本文采用有限控制体积法［1，6］来对其进行数值求解。

3 共轭梯度法

共轭梯度法的基本思想是将传热模型的偏微分方程看作为对未知参数的约束，利用拉格朗日（Lagrange）乘数法将辨识问题转化为使如下目标函数达极小的无约束优化问题［7］：

式中的λ称为伴随变量。为了获得伴随变量满足的方程，首先需推导灵敏度方程。

对式（1）中的Q（t）做扰动ΔQ，设相应的温度场变化为ΔT，热传导系数和比热变为：

显然，扰动后的物理量也满足方程式（1）的形式，与未扰动时的方程相减，即可得ΔT满足的灵敏度方程：

接下来对式（11）中的目标函数做变分，并引入上面灵敏度方程的边界条件后，有：

由ΔT的任意性，可得伴随变量满足的伴随方程为：

边界条件：x＝0：－k（T）

时域条件：t＝tf，λ＝0。

同时可求出目标函数对Q的梯度为：

将此梯度值代入如下优化算法，即为共轭梯度法：

其中，γn

U是由式（12）求得的ΔQ＝Pn引起的温度场变化值；该优化算法的计算停止准则为：

这一辨识方法中涉及到的状态方程式（6）、灵敏度方程式（12）和伴随方程（14），同样都采用有限控制体积法［1，6］来进行数值求解。

4 算例分析

下面给出表面热流辨识的算例，对于如图1所示的传热模型，设材料密度为ρ＝1650kg／m3，厚度为L＝5mm，材料热物性参数随温度变化规律如图2示，表面x＝0处的给定热流如图3中的“Exact”示，温度测点位于x＝L处，利用此给定的热流可以仿真计算出测点温度历程如图4。先不考虑测量噪声，将此温度历程作为实测值，采用顺序函数法和共轭梯度法对表面热流进行辨识，顺序函数法中的r取为30，对应的将来时间长度tM＋r－1－tM－1＝6s。图3中示出了采用顺序函数法（图中“SFM”）和共轭梯度法（图中“CGM”）辨识出的表面热流与给定值的比较，从中可以看到：两种方法的辨识结果与给定值符合很好，两种辨识方法都是可行的。

图2 材料热物性参数随温度变化函数Fig.2 Temperature dependent thermal property of the material

图3 表面热流辨识结果比较Fig.3 Comparison of estimated heat flux

就此问题而言，如果在辨识过程中不考虑材料热物性参数随温度的变化，以往的做法是将材料热物性参数用近似的常数值代替后再用线性热传导逆问题的表面热流辨识算法［8］来处理，这里也给出这一处理方法的结果。从图4可以看到，测点温度区间的中值约为T＝550K，可将材料的热物性参数取为T＝550K时对应的常数值：k＝2.4W／（m·K），Cp＝1800J／（kg·K）；同样以图4中的测点温度计算结果做为测量值，采用顺序函数法来对表面热流进行辨识，得到的表面热流如图5中的“k＝2.4，Cp＝1800”所示，可以看到，由于温度测量值是在考虑了材料热物性参数随温度变化的情况下计算出来的，但此时辨识过程中材料热物性参数采用了常数而未考虑随温度的变化，因而此时的辨识结果与图3中的“SFM”辨识结果有相当大的差异。由此可见，对于当前工况，材料热物性参数用近似常数值代替的辨识方法得到的表面热流会有较大误差，在辨识过程中考虑材料热物性参数随温度的变化是必要的。

图4 测点温度计算结果Fig.4 Calculated temperature history at measurement point

图5 物性参数选取对热流辨识结果影响Fig.5 Influence of thermal property on estimated heat flux

接下来在非线性热传导逆问题的表面热流辨识中考虑测量噪声的影响，在图4的测点温度历程计算结果上叠加标准差σ＝2.5K的白噪声（对应测量的相对误差约为0.5%）来作为模拟的测量结果，如图6中的“Exp.”所示。同样采用顺序函数法和共轭梯度法来对表面热流进行辨识，顺序函数法中的r取为40，对应的将来时间长度tM＋r－1－tM－1＝8s。图7中示出了这两种方法的辨识结果与给定值的比较，图6中的“Fitted”则示出了采用顺序函数法辨识结果计算出的测点温度拟合曲线，可以看到：采用热流辨识结果计算出的测点温度与测量值拟合较好，表面热流的辨识结果也与给定值符合。

图6 测点温度比较（σ＝2.5K）Fig.6 Comparison of temperature histories at measurement point（σ＝2.5K）

图7 表面热流辨识结果比较（σ＝2.5K）Fig.7 Comparison of estimated heat flux（σ＝2.5K）

进一步增大测量噪声的标准差到σ＝10K（对应测量的相对误差约为2%），并将该测量噪声叠加到测点温度历程计算结果上来模拟测量结果，如图8中的“Exp.”示。图9中示出了此时采用顺序函数法和共轭梯度法的辨识结果与给定值的比较，顺序函数法中的r取为45，对应的将来时间长度tM＋r－1－tM－1＝9s。图8中的“Fitted”则示出了采用共轭梯度法辨识结果计算出的测点温度拟合曲线，可以看到：在相对误差约为2%的情况下，采用热流辨识结果计算出的测点温度与测量值拟合较好，表面热流的辨识结果仍与给定值符合较好，但是由于测量噪声的增大，辨识结果与给定值之间的偏差有所增大。

图8 测点温度比较（σ＝10K）Fig.8 Comparison of temperature histories at measurement point（σ＝10K）

图9 表面热流辨识结果比较（σ＝10K）Fig.9 Comparison of estimated heat flux（σ＝10K）

从上述辨识结果可以看到：对于非线性热传导逆问题的表面热流辨识，顺序函数法和共轭梯度法都能给出较好的辨识结果，并且算法都具有较好的抗噪性能。相比较而言：

（1）顺序函数法是对热流按时间顺序估计，算法的推导和实现较为简便；而共轭梯度法是在全时间域内对热流一次进行估计，算法的推导和实现相对较为复杂。

（2）顺序函数法中存在人为确定的参数r，该参数值的选取需遵循一定的原则，文献［1，8］对此进行了专门分析。共轭梯度法中虽然没有自由参数，但需事先知道式（17）中的测量结果标准差σ，以对优化计算的步数进行控制。实际上，顺序函数法中确定r的原则也需要用到测量结果标准差σ的信息［1］，从这一点来看，两种辨识方法是类似的，但具体的要求形式不同。

（3）从计算结果的精度来看，两种辨识方法的计算结果精度大致相当，但在有测量噪声的情况下，t＜5s段共轭梯度法的辨识结果精度明显低于顺序函数法，其原因尚有待进一步深入分析。

（4）从计算效率看，两种辨识方法的计算效率与温度测量的时间段长度tf大小有关，当测量时间较短时，顺序函数法的计算效率较高，共轭梯度法的计算效率较低；而当测量时间较长时，顺序函数法的计算效率则相对较低，共轭梯度法的计算效率较高。

5 小 结

本文建立了非线性热传导逆问题的两种表面热流辨识方法：顺序函数法和共轭梯度法，给出了这两种辨识方法的基本思想和具体算法推导，并针对典型算例进行了仿真辨识，结果表明：（1）对于算例中的工况而言，采用材料热物性参数用近似常数值代替的辨识方法得到的表面热流会有较大误差，在辨识过程中考虑材料热物性参数随温度的变化是必要的。（2）本文针对非线性热传导逆问题建立的两种表面热流辨识方法虽然在算法构造、计算效率方面存在一定的差异，但都能给出较好的辨识结果，并且算法受测量噪声的影响较小，具有较好的鲁棒性。此外，本文工作虽然只讨论了一维非线性热传导逆问题的表面热流辨识算法，但相关的公式推导和求解方法可以较容易地推广应用于高维情况。

［1］BECK J V，BLACKWELL B，Jr CLAIR C R S.Inverse heat conduction ill－posed problems［M］.John Wiley ＆Sons.New York，1985.

［2］OSMAN A M，DOWDING K J，BECK J V.Numerical solution of the general two－dimensional inverse heat conduction problem （IHCP）［J］.JournalofHeat Transfer，1997，119：38－45.

［3］YANG C Y.Estimation of boundary conditions in nonlinear inverse heat conduction problems［J］.Journalof ThermophysicsandHeatTransfer，2003，17（3）：389－395.

［4］HUANG C H，YAN J Y.An inverse problem in simultaneously measuring temperature－dependent thermal conductivity and heat capacity［J］.Int.J.HeatMass Transfer，1995，38（18）：3433－3441.

［5］唐中华，钱国红，钱炜祺.材料热传导系数随温度变化函数的反演方法［J］.计算力学学报，2011，28（3）：377－382.

［6］钱炜祺，蔡金狮.用灵敏度法辨识热传导系数和热流参数［J］.空气动力学学报，1998，16（2）：226－231.

［7］ALIFANOV O M.Inverse heat transfer problems［M］.Springer－Verlag，Berlin，1994.

［8］钱炜祺，何开锋，桂业伟，等.非稳态表面热流反演算法研究［J］.空气动力学学报，2010，28（2）：155－161.