一种热传导热源位置识别反问题的数值通解方法

摘要

本发明涉及一种热传导热源位置识别反问题的数值通解方法，该数值通解方法包括如下步骤：热传导源项位置识别反问题的描述，若热传导源为点源，则直接进入下一步骤，若热传导源为非点源，则采用转换算法，将非点源反问题转化为点源反问题，再进入下一步骤；解齐次解和特解，构造数值通解；求解线性方程组，得到热源位置参数。本发明基于有限元数值解、构造出满足热传导微分方程的、以热源参数为变量的数值通解，将热传导位置识别反问题转化为多元函数极值问题，快速反演出热源参数。该方法不仅可以反演点热源位置，且可以反演任意形状热源位置，因此应用面广、适应性强，具有很好的工程应用前景。

说明书

技术领域

本发明涉及一种热传导源项位置识别反问题的方法，特别涉及一种热传导热源位置识别反问题的数值通解方法。

背景技术

热传导反问题应用是近年来工程上的研究热点，它是传热学基础研究的重点。近年来热传导反问题在工程上正展现出良好的应用前景，在动力工程、冶金与模具等工业领域有着广泛的应用背景。

热传导反问题是指对于热传导问题，已知某一时间后的结构温度分布，反演结构的热物性参数、热源函数、或几何形状、或边界条件等，这就是热传导反问题。很显然，热传导反问题相比热传导正问题，其求解的难度更大、也更费时。

目前关于反问题求解算法的专利已有不少，应用面也比较广。如中国发明专利申请号200710051566公开了一种“水利水电工程水力学反问题的研究方法”，提出了脉冲谱法、离散优化法、摄动法，控制论方法，构件了水力水电工程水利学反问题的框架。该发明属于分布参数系统反问题的方法，而非源项识别反问题的方法。

中国发明专利申请号201210350077公开了一种“求解亚音速流动的反问题的数值方法”，针对一类在给定固体壁面压力情况下进行固体壁面几何形状设计的反问题的数值方法。该发明属于几何形状识别反问题的方法，而非源项识别反问题的方法。

中国发明专利申请号201210366939公开了一种“用拉格朗日形式的欧拉方程求解一类反问题的数值方法”，提出一种新的拉格朗日形式的二维欧拉方程来求解固体壁面几何形状设计的反问题。该发明也属于几何形状识别反问题的方法，而非源项识别反问题的方法。

中国发明专利申请号201410095196.7给出了“基于全局优化算法的农田组分温度反演方法”，采用全局优化算法模拟退火法对目标函数进行极小化实现了农田组分温度的反演。

中国发明专利申请号201410032593.X给出了“一种热传导线源强度识别反问题的数值通解方法”，这是关于热传导热源强度参数识别的方法。

发明内容

针对现有技术中存在的上述问题，本发明的目的在于提出一种热传导热源位置识别反问题的数值通解方法。该方法不仅可以反演点热源位置，且可以反演任意形状热源位置，因此应用面广、适应性强，具有很好的工程应用前景。

所述的一种热传导热源位置识别反问题的数值通解方法，其特征在于包括如下步骤：

1)热传导源项位置识别反问题的描述

该热传导源项位置识别反问题描述如下：在热源q_s作用下有温度场θ(x,y,z)，求热源q_s的参数，其中给定补充条件为测量点上给定测量温度θ_d；

稳态热传导问题描述为如式(1)所示：

$>\left(\begin{array}{cc}\begin{array}{c}\begin{array}{cc}\lambda {▿}^2\theta +q_s=0& (x,y,z\in \Omega )\end{array}\\ \theta {(x,y,z)}_b={\theta }_b(x,y,z)\\ -\lambda {\left(\frac{\partial \theta }{\partial n}\right)}_v=h(\theta {|}_v-{\theta }_f)\end{array}& \left(1\right)\\ -\lambda {\left(\frac{\partial \theta }{\partial n}\right)}_w=f_2& \left(2\right)\end{array}\right)$>

式中：θ(x,y,z)为温度；▽为拉普拉斯微分算子；Ω为问题的定义域；b为第一类边界条件；v为第三类边界条件；w为第二类边界条件，h为表面传热系数或对流换热系数；θ_f为换热介质温度；q_s为热源函数；λ为导热系数，n为边界法向，f₂为热流密度，其中第一类边界条件为指定温度；第三类边界条件为对流换热；第二类边界条件为热流密度，热源函数包含了强度与位置；

2)若步骤1)所述的热传导源为点源，则直接进入步骤3)，若步骤1)所述的热传导源为非点源，则采用转换算法，将非点源反问题转化为点源反问题，再进入步骤3)；

3)解齐次解和特解，构造数值通解

根据微分方程解的基本理论，式(1)通解由齐次解与特解组成，其中齐次解是在式(1)中令q_s＝0求解得，特解则令q_s＝1求解得；

设问题(1)有k个不同位置的点源，热源表达为：

$>q_s=q(x,y,z)=\underset{i=1}{\overset{k}{\Sigma }}\left(\delta \right(x_i,y_i,z_i\left){\eta }_i\right)---\left(3\right)$>

式中：δ(x_i,y_i,z_i)为位置函数，(x_i,y_i,z_i)是位置参数；η_i为第i个点源的强度参数，稳态问题，η_i为常数；k为点源的个数；

a根据问题的性质或工程意义，对源项的位置给出其可能的位置范围,

式(3)中第i个源项q_si的位置变量是(x_i,y_i,z_i)，该位置变量设其变化范围为：x_i1≤x_i≤x_i2、y_i1≤y_i≤y_i2、z_i1≤z_i≤z_i2，其中x_i1,x_i2,y_i1,y_i2,z_i1,z_i2为已知值；引进无量纲位置参数变量有助于后续公式的简化：$>\begin{array}{ccc}{\alpha }_i^1=(x_{i1}-x_i)/(x_{i1}-x_{i2}),& {\alpha }_i^2=(y_{i1}-y_i)/(y_{i1}-y_{i2}),& {\alpha }_i^3=(z_{i1}-z_i)/(z_{i1}-z_{i2})\end{array},$>且有$>0\le {\alpha }_i^k\le 1,k=1,2,3;$>

b只考虑一个变量，如x_i，计算特解；

分别计算k个点源强度为η_i的有限元数值解，这里有限元数值解的含义是：对式(1)给出的问题定义域Ω，利用计算机通过有限元方法对第i个点源计算在2个端点(x_i1,x_i2)有强度为η_i的点热源作用下定义域Ω内的温度场，记为数值特解

令式(1)中q_s＝0，解得齐次解θ＝θ₁；采用无量纲位置参数变量，从而构造温度场数值通解θ表示为：

$>\theta ={\theta }_1+\underset{i=1}{\overset{k}{\Sigma }}\left(\right(1-{\alpha }_i^1){\theta }_{qi}^1+{\alpha }_i^1{\theta }_{qi}^2)---\left(4\right)$>

式中：为对应坐标x_i的无量纲位置参数，是待求未知量；为点源i的特解；

通过给定补充条件，即在若干测量点上考虑数值通解(4)与给定测量温度θ_d的误差，得到问题的残差平方费用函数，如下：

$>g\left({\alpha }^1\right)=\underset{j=1}{\overset{m}{\Sigma }}{({\theta }_{1j}+\underset{i=1}{\overset{k}{\Sigma }}(\left(1-{\alpha }_i^1\right){\theta }_{qij}^1+{\alpha }_i^1{\theta }_{qij}^2)-{\theta }_{dj})}^2---\left(5\right)$>

其中：m为测量点数；θ_dj为测量点j的测量温度、θ_1j是测量点j的齐次解、是第i个点源在j点的特解；

源项识别反问题式(1)转化为式(5)所表示的一个以热源位置参数为变量的多元函数的极值问题，解得极值，求得反问题的解；

c对位置变量y_i,z_i，分别重复b步，得到类似(5)式的分别以为变量的表达式，如(6)式给出，然后转到步骤4)求解；

$>\begin{array}{c}g\left({\alpha }^2\right)=\underset{j=1}{\overset{m}{\Sigma }}{({\theta }_{1j}+\underset{i=1}{\overset{k}{\Sigma }}(\left(1-{\alpha }_i^2\right){\theta }_{qij}^1+{\alpha }_i^2{\theta }_{qij}^2)-{\theta }_{dj})}^2\\ g\left({\alpha }^3\right)=\underset{j=1}{\overset{m}{\Sigma }}{({\theta }_{1j}+\underset{i=1}{\overset{k}{\Sigma }}(\left(1-{\alpha }_i^3\right){\theta }_{qij}^1+{\alpha }_i^3{\theta }_{qij}^2)-{\theta }_{dj})}^2\end{array}---\left(6\right)$>

4)求解线性方程组，得到热源位置参数

求式(5)、式(6)的极值问题，对位置参数变量求导数，即令：

$>\partial g/\partial {\alpha }_i^p=0,i=1,2,...,k;p=1,2,3$>

式中：i对应第i个未知位置变量参数，p对应坐标的三个方向，其中1对应式(5)，2对应式(6)的第一个表达式，3对应式(6)的第二个表达式，由此得到计算点源位置参数的线性方程组：

A·α＝B(7)

$>\begin{array}{c}B={\left(\underset{j}{\overset{m}{\Sigma }}\right({\theta }_{dj}-{\theta }_{1j}-\underset{i}{\overset{k}{\Sigma }}{\theta }_{qij}^1\left)\right({\theta }_{qpj}^2-{\theta }_{qpj}^1\left)\right)}_{k\times 1}\\ =\underset{j}{\overset{m}{\Sigma }}\left(\begin{array}{c}({\theta }_{dj}-{\theta }_{1j}-\underset{i}{\overset{k}{\Sigma }}{\theta }_{qij}^1)({\theta }_{q1j}^2-{\theta }_{q1j}^1)\\ .\\ .\\ .\\ ({\theta }_{dj}-{\theta }_{1j}-\underset{i}{\overset{k}{\Sigma }}{\theta }_{qij}^1)({\theta }_{qkj}^2-{\theta }_{qkj}^1)\end{array}\right)\end{array}---\left(9\right)$>

其中：α为由k个点源组成的无量纲热源位置参数向量，每一次求解只对应着某一个坐标方向；

求解式(7)后，根据求得的无量纲α得到热源的位置参数，如下：

$>\begin{array}{c}{x}_i=x_{i1}-(x_{i1}-x_{i2}){\alpha }_i^1\\ y_i=y_{i1}-(y_{i1}-y_{i2}){\alpha }_i^2\\ z_i=z_{i1}-(z_{i1}-z_{i2}){\alpha }_i^3\end{array},(i=1,2,\mathrm{...},k)---\left(10\right).$>

所述的一种热传导热源位置识别反问题的数值通解方法，其特征在于第二类边界条件中的热流密度f₂视为边界热源，其参数的反演与热源q_s相同。

所述的一种热传导热源位置识别反问题的数值通解方法，其特征在于步骤2)所述非点源转化为点源反问题的转换算法，具体包括如下步骤：将非点源反问题借助插值方法将其转化成点源反问题，然后采用步骤3)和步骤4)，进行求解。

所述的一种热传导热源位置识别反问题的数值通解方法，其特征在于对空间域的热源，转化成点源反问题的方法如下：

将空间热源离散成若干六面体或四面体热源单元的组合，即：

$>F=\underset{i}{\Sigma }{F}_i---\left(11\right)$>

以点热源作为插值节点，对六面体或四面体单元进行插值，用点源为插值节点的插值函数来描述空间热源，从而将其转化成点热源的反演问题，得点源的位置参数，从而得到空间热源的位置，即将热源单元用单元插值函数表达：

$>F_i=\underset{j}{\overset{n}{\Sigma }}{N}_j(x,y,z)P_j---\left(12\right)$>

其中：n是热源单元的插值节点数，对六面体单元而言一般有8个节点，四面体单元一般有4个节点。

使用本方法时须注意各个物理量的量纲，即要量纲统一。如采用国际单位制SI，基本量纲为：米m，千克kg，秒s，则热源强度为W(瓦)、材料导热系数的量纲为W/(m·℃)、对流换热系数的量纲为W/(m²·℃)；而温度的量纲采用℃(摄氏度)。

通过采用上述技术，与现有技术相比，本发明的有益效果如下：

1)本发明通过采用限定的构造方法，根据问题的性质或工程意义，对源项的位置给出其可能的位置范围，施加热源，求解方程(1)的齐次解和特解；

2)本发明提出的一种热传导热源位置识别反问题的数值通解方法，基于有限元数值解、构造出满足热传导微分方程的、以热源参数为变量的数值通解，将热传导位置识别反问题转化为多元函数极值问题，快速反演出热源参数。该方法不仅可以反演点热源位置，且可以反演任意形状热源位置，因此应用面广、适应性强，具有很好的工程应用前景。

附图说明

图1为求解热传导点源位置识别反问题的流程图；

图2为本发明实施案例的x_i1作用点热源的有限元解；

图3为本发明实施案例的x_i2作用点热源的有限元解。

具体实施方式

下面结合说明书附图及实施例对本发明作进一步说明，但本发明的保护范围并不仅限于此：

如图1所示，本发明所述的热传导热源位置识别反问题的数值通解方法，包括如下步骤：

1)热传导源项位置识别反问题的描述

该热传导源项位置识别反问题描述如下：在热源q_s作用下有温度场θ(x,y,z)，求热源q_s的参数，其中给定补充条件为测量点上给定测量温度θ_d；

稳态热传导问题描述为如方程(1)所示：

$>\left(\begin{array}{cc}\begin{array}{c}\begin{array}{cc}\lambda {▿}^2\theta +q_s=0& (x,y,z\in \Omega )\end{array}\\ \theta {(x,y,z)}_b={\theta }_b(x,y,z)\\ -\lambda {\left(\frac{\partial \theta }{\partial n}\right)}_v=h(\theta {|}_v-{\theta }_f)\end{array}& \left(1\right)\\ -\lambda {\left(\frac{\partial \theta }{\partial n}\right)}_w=f_2& \left(2\right)\end{array}\right)$>

式中：θ(x,y,z)为温度；▽为拉普拉斯微分算子；Ω为问题的定义域；b为第一类边界条件，即指定温度；v为第三类边界条件，即对流换热；w为第二类边界条件，即热流密度，h为表面传热系数或对流换热系数；θ_f为换热介质温度；q_s为热源函数(包含了强度与位置)；λ为导热系数，n为边界法向，f₂为热流密度；本发明中，对第二类边界条件中的热流密度f₂，可视为边界热源，其参数的反演与热源q_s相同，可以用本方法求解，因此本方法对稳态热传导源项位置识别反问题的描述可用方程(1)表达，不考虑第二类边界条件；

2)若步骤1)所述的热传导源为点源，则直接进入步骤3)，若步骤1)所述的热传导源为非点源，则采用转换算法，将非点源反问题转化为点源反问题，再进入步骤3)；

3)解齐次解和特解，构造数值通解

根据微分方程解的基本理论，方程(1)通解由齐次解与特解组成，其中齐次解是在(1)式令q_s＝0求解得、特解则令q_s＝1求解；

热源q_s的形式一般可以用点、线、面及体的热源来描述，其中点是基本的，其它形式的热源可以通过点源一定规律的组合来描述。因此以下先讨论点热源。设问题(1)有k个不同位置的点源，热源可表达为：

$>q_s=q(x,y,z)=\underset{i=1}{\overset{k}{\Sigma }}\left(\delta \right(x_i,y_i,z_i\left){\eta }_i\right)---\left(3\right)$>

式中：δ(x_i,y_i,z_i)为位置函数，(x_i,y_i,z_i)是位置参数；η_i为第i个点源的强度参数，稳态问题，η_i为常数；k为点源的个数。

数值通解方法的核心是采用数值求解方法如有限元方法数值求解方程(1)的齐次解和特解，从而构造以热源强度参数或位置参数为变量的数值通解。在强度反演问题中各点源的位置参数是已知的，因此求解方程(1)的齐次解和特解时很容易施加热源。但在位置反演问题中各点源的强度参数是已知的，而位置是待定，这样无法施加热源，也就无法求解方程(1)的齐次解和特解。

因此本发明的特征在于，给出如下构造方法：

a根据问题的性质(或工程意义)，对源项的位置给出其可能的位置范围。

式(3)中第i个源项q_si的位置变量是(x_i,y_i,z_i)，该位置变量可设其变化范围为：x_i1≤x_i≤x_i2、y_i1≤y_i≤y_i2、z_i1≤z_i≤z_i2，其中x_i1,x_i2,y_i1,y_i2,z_i1,z_i2为已知值。引进无量纲位置参数变量有助于后续公式的简化：$>\left(\begin{array}{cc}{\alpha }_i^1=(x_{i1}-x_i)/(x_{i1}-x_{i2}),& {\alpha }_i^2=(y_{i1}-y_i)/(y_{i1}-y_{i2}),\end{array}\right)$>$>{\alpha }_i^3=(z_{i1}-z_i)/(z_{i1}-z_{i2}),$>且有$>0\le {\alpha }_i^k\le 1,k=1,2,3.$>

b只考虑一个变量，如x_i，计算特解。

分别计算k个点源强度为η_i的有限元数值解。这里有限元数值解的含义是：对(1)式给出的问题定义域Ω，利用计算机通过有限元方法对第i个点源计算在2个端点(x_i1,x_i2)有强度为η_i的点热源作用下定义域Ω内的温度场，记为数值特解

令式(1)中q_s＝0，解得齐次解θ＝θ₁；采用无量纲位置参数变量，从而构造温度场数值通解θ表示为：

$>\theta ={\theta }_1+\underset{i=1}{\overset{k}{\Sigma }}\left(\right(1-{\alpha }_i^1){\theta }_{qi}^1+{\alpha }_i^1{\theta }_{qi}^2)---\left(4\right)$>

式中：为对应坐标x_i的无量纲位置参数，是待求未知量；为点源i的特解。

通过给定补充条件，即在若干测量点上考虑数值通解(4)与给定测量温度θ_d的误差，得到问题的残差平方费用函数，如下：

$>g\left({\alpha }^1\right)=\underset{j=1}{\overset{m}{\Sigma }}{({\theta }_{1j}+\underset{i=1}{\overset{k}{\Sigma }}(\left(1-{\alpha }_i^1\right){\theta }_{qij}^1+{\alpha }_i^1{\theta }_{qij}^2)-{\theta }_{dj})}^2---\left(5\right)$>

其中：m为测量点数；θ_dj为测量点j的测量温度、θ_1j是测量点j的齐次解、是第i个点源在j点的特解。

这样，源项识别反问题式(1)就转化为式(5)所表示的一个以热源位置参数为变量的多元函数的极值问题，很容易解得极值，从而求得反问题的解。

c对位置变量y_i,z_i，分别重复2.2步，得到类似(5)式的分别以为变量的表达式，如(6)式给出，然后转到步骤4)求解。

$>\begin{array}{c}g\left({\alpha }^2\right)=\underset{j=1}{\overset{m}{\Sigma }}{({\theta }_{1j}+\underset{i=1}{\overset{k}{\Sigma }}(\left(1-{\alpha }_i^2\right){\theta }_{qij}^1+{\alpha }_i^2{\theta }_{qij}^2)-{\theta }_{dj})}^2\\ g\left({\alpha }^3\right)=\underset{j=1}{\overset{m}{\Sigma }}{({\theta }_{1j}+\underset{i=1}{\overset{k}{\Sigma }}(\left(1-{\alpha }_i^3\right){\theta }_{qij}^1+{\alpha }_i^3{\theta }_{qij}^2)-{\theta }_{dj})}^2\end{array}---\left(6\right)$>

4)求解线性方程组，得到热源位置参数

求(5)、(6)式的极值问题，可对位置参数变量求导数，即令：

$>\partial g/\partial {\alpha }_i^p=0,i=1,2,...,k;p=1,2,3$>

式中：i对应第i个未知位置变量参数，p对应坐标的三个方向，1对应(5)式，2对应(6)式的第一个表达式，3对应(6)式的第二个表达式。由此得到计算点源位置参数的线性方程组：

A·α＝B(7)

$>\begin{array}{c}B={\left(\underset{j}{\overset{m}{\Sigma }}\right({\theta }_{dj}-{\theta }_{1j}-\underset{i}{\overset{k}{\Sigma }}{\theta }_{qij}^1\left)\right({\theta }_{qpj}^2-{\theta }_{qpj}^1\left)\right)}_{k\times 1}\\ =\underset{j}{\overset{m}{\Sigma }}\left(\begin{array}{c}({\theta }_{dj}-{\theta }_{1j}-\underset{i}{\overset{k}{\Sigma }}{\theta }_{qij}^1)({\theta }_{q1j}^2-{\theta }_{q1j}^1)\\ .\\ .\\ .\\ ({\theta }_{dj}-{\theta }_{1j}-\underset{i}{\overset{k}{\Sigma }}{\theta }_{qij}^1)({\theta }_{qkj}^2-{\theta }_{qkj}^1)\end{array}\right)\end{array}---\left(9\right)$>

其中：α为由k个点源组成的无量纲热源位置参数向量(每一次求解只对应着某一个坐标方向)。

求解方程(7)后，根据求得的无量纲α可以得到热源的位置参数，如下：

$>\begin{array}{c}{x}_i=x_{i1}-(x_{i1}-x_{i2}){\alpha }_i^1\\ y_i=y_{i1}-(y_{i1}-y_{i2}){\alpha }_i^2\\ z_i=z_{i1}-(z_{i1}-z_{i2}){\alpha }_i^3\end{array},(i=1,2,\mathrm{...},k)---\left(10\right)$>

对于非点源的热传导位置识别反问题，可以借助插值方法将其转化成点源反问题，然后采用上述步骤2和3，进行求解。如对空间域的热源，转化成点源反问题的方法思路如下：

将空间热源离散成若干六面体或四面体热源单元的组合，即：

$>F=\underset{i}{\Sigma }{F}_i---\left(11\right)$>

以点热源作为插值节点，对六面体或四面体单元进行插值，即将热源单元用单元插值函数表达：

$>F_i=\underset{j}{\overset{n}{\Sigma }}{N}_j(x,y,z)P_j---\left(12\right)$>

其中：n是热源单元的插值节点数，对六面体单元而言一般有8个节点，四面体单元一般有4个节点。

这样用点源为插值节点的插值函数来描述空间热源，从而将其转化成点热源的反演问题。通过上述步骤，可得点源的位置参数，从而得到空间热源的位置。

实施例：以点源为例

第1步：热传导源项位置识别反问题描述

该一维热传导源项位置识别反问题的例子可描述为：一维问题，长为1m、中点作用5W点热源，一端指定温度10℃，另一端放置在环境温度为25℃的空气中；材料导热系数20，对流换热系数25。要求反演点热源的位置。

本实例中涉及的有限元分析，可以使用目前通用的有限元分析软件，如Ansys、Cosmos等完成；

第2步：确定点源反演的相关参数

根据本发明提出的反演算法，首先确定点源反演的相关参数如下：

1)设定点源位置可能的范围，即对x_i1≤x_i≤x_i2，取l/4≤x_i≤3l/4；

2)在x＝1位置取一个测量点，根据精确解考虑偏差后给出补充条件θ₁₁＝138(单位是摄氏度℃，下同，略)；

第3步：解齐次解和特解，构造通解

根据本发明的数值通解方法，采用有限元法在测量点分别计算齐次解、x_i1、x_i1点作用5W点热源的特解，见图1与图2；提取测量点的温度值，有：$>{\theta }_{11}=10,{\theta }_{q11}^1=72.5,{\theta }_{q11}^2=\mathrm{197.5.}$>

第4步：求解方程，得到热源位置参数

根据公式(7)～(9)，可得点源位置参数为：

$>\begin{array}{c}\alpha =({\theta }_{d1}-{\theta }_{11}-{\theta }_{q11}^1)/({\theta }_{q11}^2-{\theta }_{q11}^1)\\ =(138-10-72.5)/(197.5-72.5)=0.444\end{array}---\left(13\right)$>

根据公式(10)，已知点源位置参数的无量纲量是0.5，因此本算法给出的反演结果是比较理想的。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。