刘曌阳，王婷婷，王发杰，张耀明

(山东理工大学 数学与统计学院， 山东 淄博 255049)

近年来，许多学者对热传导反问题进行了深入研究[1]，并广泛应用于冶炼、航空、化学工程等方面.

边界元法[2-3]将问题的维数降低一维，且只需对物体边界划分单元，既降低了工作难度，又提高了计算精度，是目前很有潜质的一种数值算法，已广泛应用于各种工程问题.然而，对于复杂几何边界的问题，划分单元仍然是非常耗费精力的.此外，边界元法需要计算大量复杂边界积分，这些积分甚至是奇异的和几乎奇异的，这就给解决实际问题带来了很大的困难.张耀明提出了一种新的边界型无网格方法——平均源边界节点法[4-5](Average Source Boundary Node Method, ASBNM)，该方法基于完全规则化边界积分方程和平均源技术[6]，具有无网格、无积分、半解析的特征.

本文是平均源边界节点法在模拟三维热传导反问题的一次尝试.采用截断奇异值分解(truncated singular value decomposition TSVD)和Tikhonov正则化技术[7-8]求解病态线性方程组，利用广义交叉校验准则(generalized cross validation criterion GCV)来确定正则化参数[9].

1 三维热传导反问题的平均源边界节点法

1.1 控制方程

本文假定问题的区域为有界区域Ω⊂R3，其边界为Γ=∂Ω.边界Γ由两部分组成，即Γ=Γ1∪Γ2，这里Γ1∩Γ2=Ø，函数u(x)满足Laplace控制方程：

(1)

边界条件为

(x1,x2,x3)∈Γ1

(2)

(3)

控制方程(1)的基本解为

(4)

式中|x-y|表示两点x和y之间的欧几里得距离.

1.2 平均源边界节点法

为了避免直接计算强、弱奇异积分，本文基于文献[4-5]的平均源边界节点法：

(5)

式中：Gij,Hij为系数矩阵；N为总的边界节点数；uj,qj是节点j处的温度和法向通量.

Gij,Hij如下：

(6)

(7)

远离边界的内部点y温度和温度梯度表示为

(8)

(9)

k=1,2,3

(10)

式中rk=xk-yk.

2 正则化方法

2.1 截断奇异值分解

考虑如下线性方程组

Ax=b

(11)

式中A∈Rm×n，x∈Rn，b∈Rm，且m≥n.

截断奇异值分解的基本思想[7]是用K阶矩阵AK来逼近M×N阶矩阵A.其中，AK可以表示为

(12)

这里U=(u1,…,uM)和V=(v1,…,vN)分别满足UTU=IM;VTV=IN；Σ=diag(σ1,…σN)表示非负对角矩阵；K为正则化参数，与之对应的截断奇异值为

(13)

2.2 Tikhonov正则化方法

Tikhonov法的基本思想[8]如下：

把正则化泛函

Jα(x)=‖Ax-b‖2+α2‖x‖2,α>0

(14)

的极小元xα作为方程Ax=b的解，可表示为如下形式

(15)

2.3 广义交叉校验准则

广义交叉校验准则是由Golub.G.H.[9]提出，该方法以正则化参数K为参变量，当求得GCV的极小值时，对应的K值为最优值.其计算公式为

(16)

其中AI满足

(17)

可用来产生正则化参数.

3 数值算例

数值算例中采用以下公式[8]来施加边界数据扰动：

(18)

这里，bi是精确解，rand是一个-1到1之间的随机数，δ表示扰动幅度.为了对比不同方法数值解的有效性，引入下列平均相对误差

Average relative error =

(19)

算例1球型区域

u(x1,x2,x3)=x1x2+x1x3+x2x3+

x1+x2+x3+2

(20)

图1 球型区域Fig.1 Problem sketch

计算时，球域外壳上配置200个边界节点.分别用TSVD和Tikhonov正则化方法对问题进行求解，并用GCV法选取正则化参数.如图2所示，两种方法分别在k=155和l=0.007处取得最优正则化参数.图3和图4分别描述了结合TR-GCV、TSVD-GCV法求解未知边界的温度和通量的误差曲面图.从图中可以看出，两种方法均能对三维Cauchy反问题进行有效求解.

(a)TSVD正则化参数 (b)Tikhonov正则化参数图2 GCV法选取的正则化参数Fig.2 GCV Curvers for TSVD and TR

(a)温度 (b) 通量图3 未知边界上温度和通量的误差曲面图(结合TR-GCV法)Fig.3 Relative error surfaces of the temperature and the flux on the unavailable boundary obtained using the TR-GCV method

(a) 温度 (b) 通量图4 未知边界上温度和通量的误差曲面图(结合TSVD-GCV法)Fig.4 Relative error surfaces of the temperature and the flux on the unavailable boundary obtained using the TSVD-GCV method

为了验证方法的收敛性，计算时分别在边界配置200、625、840、960、1 200个边界节点.表1给出了不同节点数下，结合TSVD和TR技术，GCV法选择的最优正则化参数.图5给出了未知边界温度和通量的平均相对误差收敛曲线.从图中可以看出，随着边界节点数的增加，边界温度和通量的平均相对误差呈现明显下降的趋势，表明该方法具有很好的收敛性，从而说明此方法对三维Cauchy反问题可以进行有效地求解.

表1 不同单元数下正则化参数的选取

Tab.1 Optimal regularization parameters for different boundary nodes

节点数GCV法参数kTR法参数λ2001550.0076254820.0018406420.0029607340.0011 2009140.001

算例2空心球区域

考虑空心球区域的稳态温度场问题.内边界Γ2的温度和通量都未知，外边界Γ1的温度和通量已知，且温度场精确解为

u(x1,x2,x3)=x1x2x3+10x1+10x2+10x3

(21)

(a) 温度 (b) 通量图5 不同边界节点数下未知边界温度和通量的平均相对误差Fig.5 Average relative error of temperatures and fluxes on underspecified boundary of different boundary nodes using the TSVD technique and the TR technique in conjunction with the GCV criterion

表2给出了GCV法选取的正则化参数，在k=924和l=1.00×10-3处获得最优参数.利用两种正则化方法及其相应的最优化参数，图6给出了在3%扰动程度下温度和通量的精确解与本文计算结果的比较.从图中可以看出，在已知条件存在随机扰动的情况下，边界温度和通量与精确解能够较好地吻合，两种方法均能求得较好的结果.

表2 不同扰动下正则化参数的选取

Tab.2Optimalregularizationparametersforinvestigatedregularizationmethodatvariousnoiselevels

正则化方法0%1%3%5%TSVD法924924924924Tikhonov法1.00×10-31.00×10-31.00×10-31.00×10-3

(a) 温度:TSVD-GCV (b) 温度：TR-GCV

(c)温度：TSVD-GCV-3%扰动 (d) 温度：TR-GCV-3%扰动 图6 未知边界上温度的数值解与解析解的对比图Fig.6 Relative error surfaces of the temperature on the unavailable boundary compared with the analytical solution

4 结束语

本文采用平均源边界节点法(ASBNM)求解三维Cauchy热传导反问题，对求解中产生的不适定线性系统，采用TSVD和Tikhonov技术求解结合GCV的正则化方法来求解.与其他现有的求解病态的Cauchy热传导反问题的方法相比，该方法无积分计算、无网格，具有计算精度高、收敛速度快、程序实现简单，适合于高维问题等优点.为Cauchy热传导反问题提供了新的求解思路，且拓宽了平均源边界节点法的应用范围.