采用序列函数法在身管导热反问题中的应用＊

周渊键 陈 刚 徐春雨

（1.解放军陆军军官学院五系41队 合肥 230031）（2.解放军陆军军官学院机械工程教研室 合肥 230031）（3.解放军陆军军官学院五系42队 合肥 230031）

1 引言

导热反问题从上个世纪60年代起就受到了人们的重视。在军事科学、工业生产研究领域的许多场合都需要知道管内流动的内壁温度。对于管外壁温度，有许多种实验测量方法；而对管内壁温度的实验测量比较困难。利用导热反问题方法，可通过数值求解由管外壁温度得到管内壁温度。导热反问题（IHCP）是相对于导热正问题而言的。对于导热正问题，物体的热物性参数和边界条件为已知，可以由此得出物体的温度分布。而当物体的热物性参数或部分边界条件必须通过已知参数来确定时，就形成了导热反问题。

导热反问题在各个领域中有着广泛的应用背景，但其非适定性、非线性、计算量大等特点，使得求解比较困难。导热反问题作为一个优化问题，已有一些求解方法，如基于梯度的优化方法（高斯牛顿法、共扼梯度法、蚁群算法）等，但都存在一定的局限性。它们不能收敛到全局最优解，并且收敛情况与初始值有关［1］。

文中以对火炮发射时内膛热交换过程的分析，基于导热反问题的研究方法，根据Beck提出的序列函数法建立由结构内部一个或多个温度测点的测量值计算内膛壁温度与热流密度的数学模型，并用FORTRAN语言编写通用计算程序。

2 内壁热交换分析及数学模型的确定

在火炮射击过程中，由于身管内膛壁受到高温高压火药燃气及弹丸摩擦的瞬时热作用。而热量是通过三种基本方式传递，分别是导热方式、对流换热方式、热辐射方式［2］。膛内热交换过程虽然复杂，但是在发射过程结束时，通过内膛面传入身管的热量随时间积累的总热量是一定的。考察边界条件的内涵，是指导热物体在其边界面上与外部环境之间在热交换方面的联系或相互作用。虽然导热问题的三类边界条件在数学表达形式上有差别，但无论导热物体是以哪一种边界条件形式与外部环境发生相互作用，热交换结果在本质上却是一致的，即热量的转移，在形式上均表现为通过某一换热面积的热流量（热流密度）［3］。

导热反问题是由导热微分方程，初始条件，边界条件，即导热微分方程定解条件的三个组成部分再加上作为补充条件或输入条件的内部温度测点信息组成。写成数学表达式为［4］

其中I、B、A分别为初始算子、边界算子和附加算子，A（x）为面积因子，对于直角坐标系A（x）＝1，对于圆柱坐标系A（x）＝r；φ（t），φ（t），Y（x，t）分别为初始条件，边界条件和附加条件。其中只有一个边界条件φ（t），而另一边界条件未知，即为待求未知热流密度。

定义敏感系数X为

根据敏感系数X定义和导热正问题微分方程，可得到关于X的微分方程［5］

由此可知，对线性及准线性问题［6］，X与qM无关，只与位置相关。若tM时刻节点k（位置xk）处热电偶温度测量值为YK，M（℃），则温度可表达为

式中，T＊（k，M）为tM－1≤t≤tM期间热流密度为q＊时节点k处温度值，单位℃，Xk为边界热流对节点k处温度的敏感系数。因此可采用如下计算步骤：任意假设一个迭代初值q＊。

由导热微分方程计算T（x，t，tM－1，qM－1，qM）；对照导热微分方程，由式（3）计算敏感系数 X（x，t，tM－1，qM－1）；利用实验测得的Y（k，M）由式（4）求得qM；再一次求解导热正问题离散方程，即可求出所有节点温度（i＝1，2，3，4，5）。

对于计算中取r个未来时间步的情况，唯一不同的是，假设热流密度q（t）＝q＊在tM－1≤t≤tM内均成立，如果在固体内部不同深度布置J个热电偶，则可采用r个未来时间步的二乘计算方法。最小二乘误差S定义为

对于任一假定热流值q＊（tM－1≤t≤tM），根据温度泰勒展开式

综合式（5）和式（6）得到热流计算表达式

定义测点k在第j未来时间步的加权系数为

以上求解过程为迭代过程，初始时刻给定的假定热流值q＊与真实值之间的差异将影响初始段的热流计算误差，而其后所有时刻的假定热流值均取上一时刻的最终收敛热流。

3 数值实验及结果分析

为验证序列函数法的有效性，针对前面提出的问题，通过计算图1所示铬钢炮管材料组成结构的内膛壁温度、热流密度来验证上述数学计算模型的正确性。

为检验序列函数法的计算效果，本文给定了一组管内膛 壁 在 0℃、45℃、90℃、135℃、180℃节点温度随时间的变化，并作为内膛壁温度的准实验值，节点温度变化如图2所示。

计算检验过程分为求解导热正问题和求解导热反问题两步进行，两个过程都利用FORTRAN语言编制程序进行计算。首先求解导热正问题，即通过求解导热微分方程（1）得到管外壁温度随时间变化。身管外径为20mm，管内径为15mm。由于计算区域的对称性，只计算半个圆周，采用四边形网格，沿圆周向为13个节点，径向为7个节点。身管内膛壁上未知的节点温度由已知的管内膛壁5个节点温度通过插值得到。

计算中取管外侧流体温度为T＝20℃，对流换热系数h＝5w／（m·K），身管的热扩散系数a＝4.5×10－6m2／s，导热系数入＝17W／（m·K）。求解导热微分方程（1），得到管外壁不同位置温度随时间的变化，并做为管外壁温度的准实验值。图3为管内膛壁和管外壁0℃节点温度变化比较。

图1 炮管截面结构示意图

图2 管内膛壁不同位置温度变化

图3 管内膛壁与管外壁0℃位置温度变化比

图4为0℃位置计算得到的管内壁温度与准实验值比较。从图中可以看到，计算值与准实验值的曲线几乎重合在一起，计算值与准实验值相符很好。序列函数法不仅能准确地捕捉到管内壁温度值的较小波动，对于30s～50s之间的较大波动也能精确地捕捉到。比较其它几个节点计算值与准实验值发现，序列函数法也准确捕捉到了其它位置的温度波动，计算值与准实验值很好地吻合。图5为管内壁45℃位置序列函数法计算值与准实验值的比较，计算值与准实验值也符合得非常好。

图4 管内膛壁0℃位置计算值与准实验值比较

图5 管内膛壁45℃位置计算值与准实验值比较

4 结语

通过建立的数学模型是针对二维非线性瞬态导热反问题的，根据Beck's序列函数思想，在计算中采用了多温度测点、多未来时间步长的计算方法。针对算例中内膛壁热流波形曲线变化情况，得到理论值和实际测量值误差不大的理想结果。结果表明，该方法拥有较强的可行性和实用性。

文中建立的方法可以应用在火炮射击特殊状态下内膛壁温度及热流难以准确测量而外部温度测量相对容易的情况下的结构内膛壁热响应特性研究。该方法不需在内壁面布置热电偶和热流密度计，也不需要考虑热传递方式的影响，只需在结构外部或外壁布置温度热电偶，简化了测量难度，并提高了内膛壁热流及温度的精度。

［1］管屏，朱刚，马良.生长竞争蚁群算法求解导热反问题［J］.上海第二工业大学学报，2010（9）：20-21.

［2］陶文铨.传热学［M］.西安：西安工业大学出版社，2006，12：78-103.

［3］吴斌，夏伟，汤勇，等.射击过程中热影响及身管热控制措施综述［J］.兵工学报，2003，24（4）：525-529.

［4］J V Beck，B Blackwell，C R St clair inverse heat conduction-Ill posed problems［M］.New York：Wiley-Interscience，1985（1）：218-243.

［5］J V Beck.Methodology for comparison of inverse heat conduction methods［J］.J.of heat transfer，1988，110：30-37.

［6］Lawton B.Thermal-chemical erosion in gun barrels［J］.Wear，2001，（251）：827-838.

［7］杨冬，陈听宽.导热反问题方法在瞬态传热过程中的应用［J］.核动力工程，1997，18（6）：553-558.

［8］李明海.钢－木组合结构的火烧热响应模拟与导热反问题在火烧试验中的应用［D］.重庆：重庆大学，1998：42-60.

［9］薛齐文，杨海天.应用共轭梯度法求解非线性多宗量稳态热传导反问题［J］.计算力学学报，2005（22）：51-54.

［10］杨海天，薛齐文.两级敏度分析求解非线性稳态多宗量热传导反问题［J］.工程热物理学报，2003，24（3）：463-465.