在许多化学物理问题中, 常常要考虑有热源的物体的温度分布规律, 而热源的发热量又与该物体温度有关, 例如内部能量的来源, 或者是计算机芯片能源的发电量, 或在微波炉加热过程中, 或在化学反应过程和许多设计与制造领域等都需要探测未知的发热源.由于热源的强度未知, 要精确地描述物体的温度分布规律, 必须获得未知源项的信息来决定热源的强度, 而要直接测定热源的发热量是相当困难的, 进而人们转向利用附加的测量数据, 通过求解相应的数学模型, 来获得热源的信息. Huang ^[1]使用共轭梯度法结合伴随矩阵识别了板内部的热源, Yang ^[2-5]基于线性最小二乘残差使用了逆矩阵法、符号法和数值序贯法求解了源强问题.本文主要对二维热传导方程点源强度问题, 采用插值方法, 将源强识别问题转化为参数反演问题, 通过微分进化算法结合最佳摄动量法给出其求解结果.

考虑如下含有初边值的二维热传导方程:

$ \begin{eqnarray} \rho c{\partial _t}T(x, y, t) = k{\nabla ^2}T(x, y, t) + g(t)\delta (x - {x^*}, y - {y^*}), (x, y) \in D, 0 < t < {t_f}, \end{eqnarray} $

(2.1)

$ \begin{eqnarray} T(x, y, 0) = {T_0}(x, y), (x, y) \in D \cup \partial D, \end{eqnarray} $

(2.2)

$ \begin{eqnarray} {\partial _x}T(0, y, t) = 0, 0 \le y \le L, 0 \le t \le {t_f}, \end{eqnarray} $

(2.3)

$ \begin{eqnarray}{\partial _x}T(L, y, t) = P(y, t), 0 \le y \le L, 0 \le t \le {t_f}, \end{eqnarray} $

(2.4)

$ \begin{eqnarray}{\partial _y}T(x, 0, t) = 0, \le x \le L, 0 \le t \le {t_f}, \end{eqnarray} $

(2.5)

$ \begin{eqnarray} {\partial _y}T(x, L, t) = Q(x, t), 0 \le x \le L, 0 \le t \le {t_f}, \end{eqnarray} $

(2.6)

其中 $D = \left\{ {(x, y)|0 < x < L, 0 < y < L} \right\}$, $\delta ( \cdot )$是狄拉克(Dirac)函数, 且 $({x^*}, {y^*}) \in D, g(t)$为热源强度, $k, ~ \rho, ~ c$分别称为材料的热导率, 密度和比热, 均为常数.

对于问题(2.1)-(2.6), 若已知 $g(t)$, ${T_0}(x, y)$, $P(y, t)$, $Q(x, t)$, 求解 $T(x, y, t)$构成了热传导方程的正问题.现在的问题是 $g(t)$未知, 如何根据附加的测量数据, 例如定点 $({x_0}, {y_0}) \in D$上的温度值 $T({x_0}, {y_0}, t) = E(t)$来反演 $g(t)$, 便构成了热传导方程的源强识别反问题.

由于反演离不开高精度的正问题求解, 本文采用有限元法对其求解.假设 ${{\rm T}_h}$是 $D$的外表面三角剖分的集合, 满足 $ \mathop {\max }\limits_{\tau \in {{\rm T}_h}}{\hbox{ diam}}(\tau ) \le h$, 其中 ${\hbox{diam}}(\tau )$是三角形 $\tau $的顶点之间的最大距离.令 ${S_h}$是 $D$上的连续函数的有限维空间.引入空间 ${S_h}$上的内积

$ (f, g) = \iint\limits_D {fgdxdy}. $

任取 $v(x, y, t) \in {S_h}$且满足 $v{|_{\partial D}} = 0$, 乘以方程(2.1) 两端, 然后在 $D$上积分得

$ \begin{equation}\rho c(T, v) = k({\nabla ^2}T, v) + \left( {g(t)\delta (x - {x^*}, y - {y^*}), v} \right), \end{equation} $

(2.7)

利用格林公式消去其中的二阶偏导数项得

$ \begin{eqnarray} \rho c\iint\limits_D {{\partial _t}T(x, y, t)vdxdy} &=& - k\iint\limits_D \nabla T(x, y, t)\nabla vdxdy+\\ &&\iint\limits_D g(t)\delta (x - {x^*}, y - {y^*})vdxdy \nonumber\\ &&+k\int_{\partial D} {n \cdot \nabla Tvdxdy}, \end{eqnarray} $

(2.8)

把边界条件(2.3), (2.6) 代入(2.8) 式可得

$ \begin{equation}\rho c\iint\limits_D {{\partial _t}T(x, y, t)vdxdy} = - k\iint\limits_D {\nabla T(x, y, t)\nabla vdxdy} + \iint\limits_D {g(t)\delta (x - {x^*}, y - {y^*})vdxdy}, \end{equation} $

(2.9)

由狄拉克(Dirac)函数的性质可得

$ \begin{equation}\rho c\iint\limits_D {{\partial _t}T(x, y, t)vdxdy} + k\iint\limits_D {\nabla T(x, y, t)\nabla vdxdy} - g(t)v({x^ * }, {y^ * }) = 0. \end{equation} $

(2.10)

设方程的近似解为

$ T(x, y, t) \approx \sum\limits_{j = 1}^N {{T_j}(t)} {\Phi _j}(x, y), $

其中 $\{ {\Phi _j}(x, y)\} _{j = 1}^N$是标准的三角剖分的基函数, 代入方程(2.10) 可得

$ \begin{equation}\rho c\sum\limits_{i = 1}^N {{{T'}_i}(t)} \left( {{\Phi _i}, {\Phi _j}} \right) + k\sum\limits_{i = 1}^N {{T_i}(t)\left( {\nabla {\Phi _i}, \nabla {\Phi _j}} \right)} - g(t){\Phi _j}({x^*}, {y^*}) = 0, ~~ j = 1, 2, \cdots, N.\end{equation} $

(2.11)

若令

$ \begin{eqnarray*}&& T = {\left( {{T_1}(t), {T_2}(t), \cdots , {T_N}(t)} \right)^T}, \\ && T' = {\left( {\frac{{d{T_1}(t)}}{{dt}}, \frac{{d{T_2}(t)}}{{dt}}, \cdots , \frac{{d{T_N}(t)}}{{dt}}} \right)^T}, \end{eqnarray*} $

代入初始条件(2.2), 我们得到如下矩阵形式的有限元离散

$ \begin{equation}\left\{ \begin{array}{l} \rho cAT' + kBT = {F_1}, \\ T{|_{t = 0}} = {T_0}({x_i}, {y_i}), \end{array} \right. \end{equation} $

(2.12)

其中

$ {A_{ij}} = \left( {{\Phi _i}, {\Phi _j}} \right) = \iint\limits_D {{\Phi _i}{\Phi _j}dxdy}, {B_{ij}} = \left( {\nabla {\Phi _i}, \nabla {\Phi _j}} \right) =\\ \iint\limits_D {\nabla {\Phi _i}\nabla {\Phi _j}dxdy}, {F_{{1_j}}} = g(t){\Phi _j}({x^*}, {y^*}) . $

对方程(2.12) 的时间导数采用差分近似, 假设在 $t$时刻的值为 ${T_t}$, 则

$ \frac{{dT}}{{dt}} = \frac{{{T_{t + \Delta t}} - {T_t}}}{{\Delta t}}, $

代入(2.12) 式可得

$ \left( {\rho c\frac{A}{{\Delta t}} + kB} \right){T_{t + \Delta t}} = \rho c\frac{A}{{\Delta t}}{T_t} + {F_1}. $

这样便可求出在区域内任一点及任一时刻的数值解.

有了正问题的求解, 下面来确定 $g(t)$.为了确定源强 $g(t)$, 需要附加条件, 本文选取已知点 $({x_0}, {y_0}) \in D$的函数值

$ \begin{equation}T({x_0}, {y_0}, t) = E(t), \end{equation} $

(3.1)

其中 $E(t)$为已知函数.

设未知函数 $g(t)$在时刻 ${t_i}$的值为 $g({t_i})$, $\left( {i = 0, 1, \cdots, n} \right)$, 并记

$ \begin{equation}g({t_i}) = {a_i} t \in [0, {t_f}], i = 0, 1, \cdots, n.\end{equation} $

(3.2)

选取冒状基函数 ${w_i}(t)$, 对 $g(t)$进行插值, 即确定一个 $n+1$维实向量 $a = {({a_0}, {a_1}, \cdots, {a_n})^T} \in {R^{n + 1}}$, 使得

$ \begin{equation} g(t) \approx \tilde g(t) = \sum\limits_{i = 0}^n {{a_i}{w_i}(t)} = {a^T}W(t), \end{equation} $

(3.3)

其中 $W(t) = {\left( {{w_0}(t), {w_1}(t), \cdots, {w_n}(t)} \right)^T}$, ${w_i}(t)$, $\left( {i = 0, 1, \cdots, n} \right)$具有如下形式的线性分段函数

$ {w_i}(t) = \left\{ {\begin{array}{*{20}{c}} {\frac{{t - {t_{i - 1}}}}{h}, ~~~~ t \in \left[{{t_{i-1}}, {t_i}} \right]}, \\ {\frac{{{t_{i + 1}} - t}}{h}, ~~~~ t \in \left[{{t_i}, {t_{i + 1}}}\right]}, \\ 0, ~~~~~~~~{\mbox{其它}}, ~~~~~~~~ \end{array}} \right. ~~~i = 1, 2, \cdots, n - 1, $

当 $i = 0, n$时, 定义

$ w{}_0(t) = \left\{ {\begin{array}{*{20}{c}} {\frac{{{t_1} - t}}{h}, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} t \in \left[{{t_0}, {t_1}} \right], {\kern 1pt} }\\ {0, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt}\mbox{其它} }, \end{array}} \right.~~~ {w_n}(t) = \left\{ {\begin{array}{*{20}{c}} {\frac{{t - {t_{n - 1}}}}{h}, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} t \in \left[{{t_{n-1}}, {t_f}} \right]}, \\ {0, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt}\mbox{其它} }, \end{array}} \right. $

其中 $h = \frac{{{t_f}}}{n}$.则将源强 $g(t)$的反演转换为参数 ${a_i}$的反演, 从而将未知函数的反演转化为具体参数的反演.

设 $T({x_0}, {y_0}, {t_i})$是固定点 $({x_0}, {y_0})$在 ${t_i}$时刻的第 $i$个温度观测值, 即 $T({x_0}, {y_0}, {t_i}) = E({t_i})$, $T({a_0}, {a_1}, \cdots, {a_n}, {x_0}, {y_0}, {t_i})$是以 ${a_0}, {a_1}, \cdots, {a_n}$为参数解正问题(2.1)-(2.6) 所获得的相应计算值, 则反问题转化为如下非线性优化问题:

$ \begin{equation}\mathop {\min }\limits_{({a_1}, {a_2}, \cdots, {a_n})} \sum\limits_{i = 0}^n {{{\left( {E({t_i}) - T({a_0}, {a_1}, \cdots , {a_n}, {x_0}, {y_0}, {t_i})} \right)}^2}} \begin{array}{*{20}{c}} {}&{} \end{array}{\kern 1pt} {\kern 1pt} ({x_0}, {y_0} \in D), \end{equation} $

(3.4)

其中 $T({a_0}, {a_1}, \cdots, {a_n}, {x_0}, {y_0}, {t_i})$满足方程(2.1)-(2.6).

对于非线性优化问题, 传统的优化方法大多数是基于迭代原理的各种数值方法, 函数的局部性质将严重影响算法的收敛性和收敛速度, 若要求非线性问题的全局最优解更是困难.近年来, 基于自然法则的随机搜索算法——演化算法, 在优化领域发挥了一定的优势.假设逆热源问题的源强精确解为 $p(t)$, 对于问题(3.4), 首先要给出未知量 ${a_i}{\kern 1pt} (i = 0, 1, \cdots, n)$的初始猜测值 $a_i^0{\kern 1pt} (i = 0, 1, \cdots, n)$.由于最佳摄动量法^[7-8]具有局部寻优性, 初值的选取决定了能否找到全局最优解, 而微分进化算法是一种基于种群进化的多点搜索算法^[9-10], 是进化算法的一个分支, 与其它进化算法相比, 包含着相似的机理, 如初始化种群, 进行变异、交叉和选择操作, 不断进化更新, 判断停止条件等等.最大的不同之处是, 它的变异算子是从当前种群中选取的两个或多个任意个体做差值运算, 并乘以系数得到的, 因此, 将这两种方法结合将是一种比较理想的反演方法.

其算法步骤如下:

(1) 首先采用微分进化算法反演出全局近似的最优解 $a_0^0, a_1^0, \cdots, a_n^0$;

(2) 应用有限元法求解初边值问题(2.1-(2.6) 的数值解 $T\left( {{g_0}(t);{x_0}, {y_0}, {t_i}{\kern 1pt} } \right)$, 其中

$ {g_0}(t) = \sum\limits_{i = 0}^n {a_i^0{w_i}(t)}; $

(3) 用数值微分计算

$ {a_{i, j}} = \frac{{\partial T({g_0}(x);{x_0}, {y_0}, {t_i})}}{{\partial {k_j}}} = \frac{1}{\tau }\left[{T({g_0}(t) + \tau {w_j}(t);{x_0}, {y_0}, {t_i})- T({g_0}(t);{x_0}, {y_0}, {t_i})} \right], $

即分别对 ${g_0}(t) + \tau {w_{j.}}(t)$ $(i = 0, 1, \cdots, n)$应用数值方法求边值问题(2.1)-(2.6) 的数值解 $T\left( {{g_0}(t) + \tau {w_j}(t);{x_0}, {y_0}, {t_i}} \right) ;$

(4) 求解方程组: $\left( {{A^T}A + \alpha I} \right)\delta {a_0} = {A^T}\left( {{U^*} - U}\right)$, 其中

$ A = {\left( {{a_{i, j}}} \right)_{n \times n}}, {a_{i{\rm{, }}j}}{\rm{ = }}\frac{{\partial T\left( {{g_0}(x);{x_0}, {y_0}, {t_i}} \right)}}{{\partial {k_j}}}, U = \left[{\begin{array}{*{20}{c}} {T({g_0}(t);{x_0}, {y_0}, {t_1})}\\ \cdots \\ {T({g_0}(t);{x_0}, {y_0}, {t_n})} \end{array}} \right], {U^*} = \left[{\begin{array}{*{20}{c}} {E({t_1})}\\ \cdots \\ {E({t_n})} \end{array}} \right], $

并用式 $\delta {g_0}(t) = \sum\limits_{i = 0}^n {\delta a_i^0{w_i}(t)} = \delta a_0^TW(t)$求得扰动量 $\delta {g_0}(t)$, 取新的初始猜测值为 ${g_1}(t) = {g_0}(t) + \delta {g_0}(t)$, 重复上述过程, 到指定迭代次数为止.

在问题(2.1)-(2.6) 中, 设 $D = \left\{ {(x, y)|0 < x < 1, 0 < y < 1} \right\}$, $0 < t \le 1$, $k = \rho = c = 1.$考虑如下方程

$ \begin{eqnarray*}&& {\partial _t}T(x, y, t) = {\nabla ^2}T(x, y, t) + g(t)\delta (x - 0.5, y - 0.5), 0 < t < 1, \\ && T(x, y, 0) = 1, 0 < t < 1, \\&& \frac{{\partial T(x, y, 0)}}{{\partial n}} = 0, 0 \le t \le 1. \end{eqnarray*} $

数值模拟时, 先假设 $g(t)$已知, 取 $g(t) = \left\{ {\begin{array}{*{20}{c}} {2t, ~~ 0 \le t \le 0.5}, \\ {1, ~~ 0.5 \le t \le 1}. \end{array}} \right.$采用本文所提供的有限元算法, 对区域 $D$进行三角剖分, 取时间步长 $\Delta t = 0.1$, 我们得到 $T(x, y, t)$在点 $(0.25, 0.25)$及时刻 ${t_i}$处的函数值如下表所示:

现在假设 $g(t)$未知, 把表 4-1的结果作为一组附加条件, 采用本文的插值方法, 取 ${t_i}$为区间 $\left( {0, 1} \right)$的第 $i$个结点, $n = 10$, 则所求源强 $g(t)$可以近似表示为 $\tilde g(t) = \sum\limits_{i = 0}^{10} {{a_i}{w_i}(t)} {\kern 1pt} $, 其中 ${a_i}({\kern 1pt} {\kern 1pt} i = 0, 1, \cdots, 10)$为所求参数.

微分进化算法的参数设置:群体大小: 50;终止代数: 100;交叉概率: 0.1;交叉因子: 0.5.将得出的最优估计值作为最佳摄动量法的初值.由于初值的选取不再是任意的给出的, 最佳摄动量法只需要很少的迭代也很接近真值, 故指定迭代次数为50代, 反演结果如表 4-2所示, 真解与反演值的比较见图 1.

由表 4-2和图 1可以看出, 本文的反演精度较高, 计算结果和原函数图像吻合度较好, 文献[6]对此类问题也进行了研究, 并给出了数值计算结果.本文结果明显优于文献[6].且和单纯的最佳摄动量法相比, 初值的选择更具合理性, 避免了盲目性.

本文通过有限元理论和插值方法对二维逆热源问题的源强进行了识别, 模拟结果表明算法的可行性和有效性, 也为求解此类反问题提供了一个可行的方法.在利用最佳摄动量法求解时, 采用了微分进化算法选取初值, 将两种方法结合, 发挥了各自优势, 同时有效克服了各自的弊端, 反演效果令人满意.