CN103778294A - 一种热传导线源强度识别反问题的数值通解方法 - Google Patents

Abstract

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

Description

一种热传导线源强度识别反问题的数值通解方法

技术领域

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

背景技术

热传导反问题是基础传热学研究的热点之一，工程上许多问题都可以通过反问题数学模型的描述，将问题的因果关系表达得更清晰，从而找到更有效的解决方法，因此在动力工程、冶金与模具等工业领域有着广泛的应用背景。其中，源项识别反问题成为近年来对扩散问题和热传导反问题的研究热点之一，如模具加热问题等。

目前已有关于反问题的专利，中国发明专利申请号201210350077公开了一种“求解亚音速流动的反问题的数值方法”，针对一类在给定固体壁面压力情况下进行固体壁面几何形状设计的反问题的数值方法。中国发明专利申请号201210366939公开了一种“用拉格朗日形式的欧拉方程求解一类反问题的数值方法”，提供并求解一种新的拉格朗日形式的二维欧拉方程来求解固体壁面几何形状设计的反问题。上述专利均属于几何形状识别反问题的方法，而非源项识别反问题的方法。

中国发明专利申请号200710051566公开了一种“水利水电工程水力学反问题的研究方法”，提出了脉冲谱法、离散优化法、摄动法，控制论方法，构件了水力水电工程水利学反问题的框架。该发明属于分布参数系统反问题的方法，而非源项识别反问题的方法。

目前的模具加热系统设计，通常是采用经验设计，然后再通过温度试验反复调试、修改，直到符合相关标准要求，会造成设计的费时费力，且加工成本高。

发明内容

针对上述不足，本发明提出一种热传导线源强度识别反问题的数值通解方法，所述方法可以反演线源强度，且不受求解域形状的限制，因此具有应用面广、适应性强，工程应用前景好等特点。

为了实现上述目的，本发明的解决方案是：

一种热传导线源强度识别反问题的数值通解方法，包括以下步骤：

步骤1：描述热传导源项识别反问题，建立费用函数；

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

稳态热传导问题可描述为：

$\begin{array}{c}\mathrm{\λ}{\▿}^2\mathrm{\θ}+q_s=0(x,y,z\∈\mathrm{\ϵ\Ω})\\ \mathrm{\θ}{(x,y,z)}_b={\mathrm{\θ}}_b(x,y,z)\\ -\mathrm{\λ}{\left(\frac{\∂\mathrm{\θ}}{\∂n}\right)}_v=h(\mathrm{\θ}{|}_v-{\mathrm{\θ}}_f)\end{array}---\left(1\right)$

式中：θ(x,y,z)为温度；

为拉普拉斯微分算子；Ω为问题的定义域；b为第一类边界条件(指定温度)；v为第三类边界条件(对流换热)；h为表面传热系数或对流换热系数；θ_f为换热介质温度；q_s为热源强度；λ为导热系数，n为边界法向。不失一般性，式（1）不考虑第二类边界条件，反求边界热流可以用本方法求解。

其中，θ(x,y,z)的最小化残差平方费用函数为：

$g=\underset{j}{\overset{m}{\mathrm{\Σ}}}{(\mathrm{\θ}(x_j,y_j,z_j)-{\mathrm{\θ}}_{d,j})}^2---\left(2\right)$

式中：m为测量点数；θ_d,j为测量点j的测量温度。

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

步骤3：解齐次解和特解，构造数值通解；

设有k个不同位置的点源，则热源可表达为：

$q_s=q(x,y,z)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\left(\mathrm{\δ}(x_i,y_i,z_i){\mathrm{\η}}_i\right)---\left(3\right)$

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

首先令式（1）中q_s=0，解得齐次解θ=θ₁；然后在式（1）中分别计算k个不同位置点源强度为1（W/m³）的数值解，解得特解θ=θ_s,i，这样可得到温度场通解θ为热源q_s的函数表达式：

$\mathrm{\θ}={\mathrm{\θ}}_1+\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\η}}_i{\mathrm{\θ}}_{s,i}---\left(4\right)$

式中：η_i为待求未知量；θ_s,i为点源i的特解，是通过有限元数值计算得到的强度为1W/m³的数值解。

将式（4）代入式（2）进行计算，得到：

$g\left(\mathrm{\η}\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{(\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\η}}_i{\mathrm{\θ}}_{s,\mathrm{ij}}+{\mathrm{\θ}}_{1,j}-{\mathrm{\θ}}_{d,j})}^2---\left(5\right)$

式中：θ_s,ij为点源i在测量点j的特解，θ_1,j为测量点j的齐次解。

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

步骤4：求解线性方程组，得到热源强度参数；

令 $\∂g/\∂{\mathrm{\η}}_{\mathrm{\α}}=0,\mathrm{\α}=\mathrm{1,2},...,k$

得到计算点源强度参数的线性方程组：

A·η=B    （6）

式中：η为由k个点源组成的热源强度参数向量。

上述步骤2所述的将线源反问题转化为点源反问题的转换算法，具体包括如下步骤：

1）线源以点源为插值节点作插值，即用点源的插值函数描述线源，得到线源强度的分布函数，从而将其转化成点源反演；

2）设线源分成k-1段，每段用线性插值得其分段表达式为：

$q_s=\left\{\begin{array}{c}{s}_1{\mathrm{\η}}_1(1-s_1){\mathrm{\η}}_2,0\≤s_1\≤1\\ .\\ .\\ .\\ s_{k-1}{\mathrm{\η}}_{k-1}+(1-s_{k-1}){\mathrm{\η}}_k,0\≤s_{k-1}\≤1\\ s_r=(x_{r+1}-x)/(x_{r+1}-x_r),r=\mathrm{1,2},...,k-1\end{array}\right.---\left(8\right)$

特别的，分段常数的线源，可表示为：

$q_s=\left\{\begin{array}{c}{\mathrm{\η}}_{\mathrm{1,2}},x_1\≤x\≤x_2\\ {\mathrm{\η}}_{\mathrm{2,3}},x_2\≤x\≤x_3\\ .\\ .\\ .\\ {\mathrm{\η}}_{k-1,k},x_{k-1}\≤x\≤x_k\end{array}\right.---\left(9\right)$

式中：x₁、x₂、...、x_k为线源分成k-1段后的各小段的端点坐标；η₁为对应x₁端点的强度参数，η_1,2为对应x₁x₂小段的强度参数，余同；

3）通过上述公式（8）或（9），得到点源的位置及个数。

上述热传导线源强度识别反问题的数值通解方法，基于有限元数值解、构造出满足热传导微分方程的、以热源参数为变量的数值通解，将热传导源项识别反问题转化为多元函数极值问题，快速反演出热源参数。该方法可以反演源项的强度，且不受求解域形状的限制，因此应用面广、适应性强，具有很好的工程应用前景。

以下结合附图及具体实施例对本发明做进一步详细描述。

附图说明

图1为本实施例热传导线源强度识别反问题的数值通解方法流程图；

图2为本实施例橡注机模具加热板结构图；

图3为本实施例热管的强度参数分布；

图4为本实施例热板测量点位置示意图；

图5a为本实施例热板表面的点的位置；

图5b为本实施例热板表面的温度场（℃）；

图6为实施例热板表面3个点的实验升温曲线。

具体实施方式

下面结合实例进一步说明本发明所述的热传导线源强度识别反问题的数值通解方法，具体操作步骤如下：

步骤1：描述热传导源项识别反问题，建立费用函数；

橡注机模具加热板热管的线源强度识别反问题可描述为：橡注机模具加热板1，以下简称热板，如图2所示，中间设置有注射孔2，并按一定规律分布有8根发热元件热管，分别安放在8个通孔3中。热板标准要求为加热一定时间后，其表面温度均匀，各点的温差小于±2℃。热管的强度设计是保证表面温度均匀的关键。

上述问题参数代入方程式（1），

$\begin{array}{c}\mathrm{\λ}{\▿}^2\mathrm{\θ}+q_s=0(x,y,\mathrm{z\ϵ\Ω})\\ \mathrm{\θ}{(x,y,z)}_b={\mathrm{\θ}}_b(x,y,z)\\ -\mathrm{\λ}{\left(\frac{\∂\mathrm{\θ}}{\∂n}\right)}_v=h(\mathrm{\θ}{|}_v-{\mathrm{\θ}}_f)\end{array}---\left(1\right)$

得到： $48{\▿}^2\mathrm{\θ}+q_s=0(x,y,z\∈\mathrm{\Ω})---\left(10\right)$

其中边界条件如下：

$\begin{array}{cc}-48\left(\frac{\∂\mathrm{\θ}}{\∂y}\right)=12(\mathrm{\θ}{|}_{\mathrm{xv}}-20)& x_v=\±250\\ -48\left(\frac{\∂\mathrm{\θ}}{\∂x}\right)=12(\mathrm{\θ}{|}_{\mathrm{yv}}-20)& y_v=\±250\end{array}---\left(11\right)$

步骤2：将线源反演转化成点源反演问题

建立热管热源的数学模型，将热管简化成线源，其强度分布为沿热管长度方向的强度变化。受热管制造工艺的限制，将强度设计为分段变化，且相对热管中心对称，如图3所示。按照公式（9），

$q_s=\left\{\begin{array}{c}{\mathrm{\η}}_{\mathrm{1,2}},x_1\≤x\≤x_2\\ {\mathrm{\η}}_{\mathrm{2,3}},x_2\≤x\≤x_3\\ .\\ .\\ .\\ {\mathrm{\η}}_{k-1,k},x_{k-1}\≤x\≤x_k\end{array}\right.---\left(9\right)$

得到热管强度的分段表达式为：

$q_s=\left\{\begin{array}{c}{\mathrm{\η}}_1^{\′},0\≤x\≤l/4;3l/4\≤x\≤l\\ n_2^{\′},l/4\≤x\≤3l/4\end{array}\right.---\left(12\right)$

式中：η′₁为热管两端第一段和第三段的强度参数，η′₂为中间第二段的强度参数；l为热管长度。

上述公式（12）将线源转化为点源表达，即每根热管的强度可以用4个点源描述，其位置分别是l/8、3l/8、5l/8和7l/8。

根据热板和热管的对称性，取1/4进行分析，如图4所示，分析对象有4个线源，用8个点源来描述，编号为a～h。根据温度均匀性要求、热板几何形状和热源分布特点，在热板上选取14个测量点，编号1～14。

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

根据本发明的数值通解方法，对本热板问题的数学模型（10）、（11），令q_s=0，得到齐次解θ₁=20。

对8个点热源分别作用的情况采用有限元方法计算数值特解，并提取数值特解在各测量点的温度值，如表1所示。表中，θ_a列表示在a点热载荷下14个测量点的数值解，其余同。

采用有限元方法计算，即使用目前通用的有限元分析软件，如ansys、Cosmos等，对本热板问题的数学模型式（10）与式（11）进行数值计算。以第1个点热源为例，使用有限元分析软件计算仅在图5中的a点处有1个热源情况下的热板温度场，可以得到表1中第2列，即为对应第1个点热源的特解。依次类推，得到8个点热源的特解，即表1的第2列～第9列。

表1：8个特解的测量点温度数值特解

测量点	θa/℃	θb/℃	θc/℃	θd/℃	θe/℃	θf/℃	θg/℃	θh/℃
									1	0.377	2.902	0.307	1.580	0.218	0.823	0.166	0.498
2	0.979	1.141	0.569	0.865	0.328	0.547	0.227	0.370
									3	1.180	0.486	0.714	0.421	0.408	0.328	0.277	0.257
4	0.318	1.675	0.307	2.314	0.251	1.158	0.204	0.692
									5	0.608	0.865	0.802	0.870	0.447	0.652	0.307	0.467
6	0.763	0.431	0.921	0.406	0.568	0.355	0.392	0.304
									7	0.210	0.798	0.247	1.196	0.289	2.320	0.290	1.392
8	0.315	0.536	0.444	0.683	0.802	0.860	0.546	0.781
									9	0.410	0.341	0.585	0.369	0.936	0.398	0.764	0.404
10	0.159	0.516	0.202	0.763	0.281	1.418	0.350	2.416

11	0.223	0.382	0.314	0.509	0.545	0.762	0.896	0.968
									12	0.317	0.308	0.454	0.365	0.795	0.425	1.262	0.472
13	0.125	0.360	0.167	0.533	0.256	0.965	0.376	1.770
									14	0.166	0.298	0.236	0.412	0.405	0.664	0.708	1.028

步骤4：求解线性方程组，得到热源强度参数；

根据表1的在14个测量点上的温度数值特解，按式（7）进行计算，其中k=8,m=14，下标α=1,2,...,8；j=1,2,...,14，得到（6）式，即计算8个点源强度参数的8阶线性方程组。求解该线性方程组，得到热管参数，如表2所示。

表2：强度参数

ηa	ηb	ηc	ηd	ηe	ηf	ηg	ηh
								51.6	17.3	46.6	16.0	29.6	6.3	55.9	39.3

根据表2的热管功率方案，加工热板，进行加热实验。在热板加热过程中，对热板表面的6、8和9这3个点的温度变化进行监测，得到这3个点的升温曲线，如图6所示。由图6可以看到，当热板达到动平衡的稳态阶段后，热板温度达到了所需的温度值，且3个点的温度值非常接近，比较合理地反映了对热板温度均匀性的要求，说明热板设计符合标准要求，也说明本发明所提出的方法合理有效。

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

Claims (2)

1.一种热传导线源强度识别反问题的数值通解方法，其特征在于包括以下步骤：

步骤1：描述热传导源项识别反问题，建立费用函数；

在热源q_s作用下有温度场θ(x,y,z)，求热源q_s的参数，其中给定补充条件为测量点上给定测量温度θ_d；

稳态热传导问题为：

$\begin{array}{c}\mathrm{\λ}{\▿}^2\mathrm{\θ}+q_s=0(x,y,z\∈\mathrm{\ϵ\Ω})\\ \mathrm{\θ}{(x,y,z)}_b={\mathrm{\θ}}_b(x,y,z)\\ -\mathrm{\λ}{\left(\frac{\∂\mathrm{\θ}}{\∂n}\right)}_v=h(\mathrm{\θ}{|}_v-{\mathrm{\θ}}_f)\end{array}---\left(1\right)$

式中：θ(x,y,z)为温度；

为拉普拉斯微分算子；Ω为问题的定义域；b为第一类边界条件；v为第三类边界条件；h为表面传热系数或对流换热系数；θ_f为换热介质温度；q_s为热源强度；λ为导热系数，n为边界法向；

其中，θ(x,y,z)的最小化残差平方费用函数为：

$g=\underset{j}{\overset{m}{\mathrm{\Σ}}}{(\mathrm{\θ}(x_j,y_j,z_j)-{\mathrm{\θ}}_{d,j})}^2---\left(2\right)$

式中：m为测量点数；θ_d,j为测量点j的测量温度；

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

步骤3：解齐次解和特解，构造数值通解；

热源表达为：

$q_s=q(x,y,z)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\left(\mathrm{\δ}(x_i,y_i,z_i){\mathrm{\η}}_i\right)---\left(3\right)$

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

首先令式（1）中q_s=0，解得齐次解θ=θ₁；然后在式（1）中分别计算k个不同位置点源强度为1（W/m³）的数值解，解得特解θ=θ_s,i，可得到温度场通解θ为热源q_s的函数表达式：

$\mathrm{\θ}={\mathrm{\θ}}_1+\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\η}}_i{\mathrm{\θ}}_{s,i}---\left(4\right)$

式中：η_i为待求未知量；θ_s,i为点源i的特解，是通过有限元数值计算得到的强度为1W/m³的数值解。

将式（4）代入式（2）进行计算，得到：

$g\left(\mathrm{\η}\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{(\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\η}}_i{\mathrm{\θ}}_{s,\mathrm{ij}}+{\mathrm{\θ}}_{1,j}-{\mathrm{\θ}}_{d,j})}^2---\left(5\right)$

式中：θ_s,ij为点源i在测量点j的特解，θ_1,j为测量点j的齐次解；

步骤4：求解线性方程组，得到热源强度参数；

令 $\∂g/\∂{\mathrm{\η}}_{\mathrm{\α}}=0,\mathrm{\α}=\mathrm{1,2},...,k$

得到计算点源强度参数的线性方程组：

A·η=B    （6）

式中：η为由k个点源组成的热源强度参数向量。

2.如权利要求1所述的一种热传导线源强度识别反问题的数值通解方法，其特征在于：上述步骤2所述的将线源反问题转化为点源反问题的转换算法，具体包括如下步骤：

1）线源以点源为插值节点作插值，得到线源强度的分布函数；

2）设线源分成k-1段，每段用线性插值得其分段表达式为：

$q_s=\left\{\begin{array}{c}{s}_1{\mathrm{\η}}_1(1-s_1){\mathrm{\η}}_2,0\≤s_1\≤1\\ .\\ .\\ .\\ s_{k-1}{\mathrm{\η}}_{k-1}+(1-s_{k-1}){\mathrm{\η}}_k,0\≤s_{k-1}\≤1\\ s_r=(x_{r+1}-x)/(x_{r+1}-x_r),r=\mathrm{1,2},...,k-1\end{array}\right.---\left(8\right)$

特别的，分段常数的线源，可表示为：

$q_s=\left\{\begin{array}{c}{\mathrm{\η}}_{\mathrm{1,2}},x_1\≤x\≤x_2\\ {\mathrm{\η}}_{\mathrm{2,3}},x_2\≤x\≤x_3\\ .\\ .\\ .\\ {\mathrm{\η}}_{k-1,k},x_{k-1}\≤x\≤x_k\end{array}\right.---\left(9\right)$

式中：x₁、x₂、...、x_k为线源分成k-1段后的各小段的端点坐标；η₁为对应x₁端点的强度参数，η_1,2为对应x₁x₂小段的强度参数，余同；

3）通过上述公式（8）或（9），得到点源的位置及个数。

Images