摘要本文介绍了一种用于求解人体胸腔三维恒定电场有限元方程组的算法。该算法采用直接法中克劳特分解法与分块技术相结合来求解已合成的总体有限元方程组，从而解决了计算机因内存有限而使程序经常死机的问题。该算法已成功地应用于心阻抗血流图三维有限元建模的理论研究中。

An Algorithm to Solve the Finite Element Eguations for

Thoracic 3-Dimentional Steady Electric Field Analysis

Wang Jianqi

(Dept.of BME, Fourth Military Medical University, Xi′an 710032)

Zhu xinya

(Dept.of BME, Fourth Military Medical University, Xi′an 710032)

Wang Haibin

(Dept.of BME, Fourth Military Medical University, Xi′an 710032)

Dong Xuozhen

(Dept.of BME, Fourth Military Medical University, Xi′an 710032)

Qi Jiaxue

(Dept.of BME, Fourth Military Medical University, Xi′an 710032)

Abstract

In this paper is introduced an algorithm to solve the finite element equations for thoracic 3-dimentional steady electric field analysis. The algorithm can give the solution of resultant finite element equations by means of the combination of direct Crout decomposition method and blocking technology. Thereby, we have solved the problem that the computer stops running the program because of not enough memory. In has been applied to theoretical study of modeling of 3-dimentional finite element of impedance cardiograph successfully.

Key words:Finite element; Algorithm; Computer; Memory

0引言

我们采用三维有限元建模技术进行心阻抗血流图基础理论研究时，需对求解的三维电场有限元方程组的若干问题进行讨论，这主要包括总体系数阵元素的存储技术，强加的第一类边界条件处理，有限元方程求解算法的选择等问题。本文在此重点就有限元方程求解中算法的选择问题进行研究。提出一种适合于有限内存的有限元方程组求解的新算法，以解决程序经常死机的问题。

在选择有限元方程组求解算法时，必须考虑的两个主要因素是计算机的储存量和计算时间，对586/133/8而言，前者更为重要。有限元方程组的解法可以分为迭代法与直接法两大类。迭代法的优点是所需的存贮量较少，程序编制较容易，总体系数矩阵中非零元素的分布不受限制。但是，在许多实际工程问题中，它的收敛速度较慢，计算时间难以估计。直接法虽然需要较多的存贮量与较复杂的程序编制技巧，但它的计算时间可以估计，所需的时间一般也比较短，而且文中把分块技术引入直接法中^［1］，使其存贮量大的问题得到改善，并且随着计算机硬件的迅速发展，有的缺点已不成为决定性的问题，因而大多数有限元方程都采用直接解法。

我们知道人体胸腔是一个容积导体，而采用Kubicek恒流式环形四电极法所限人体颈根部到胸部剑突处该空间区域内的电场问题可表示成边值问题，即用偏微分方程和定解条件描述，最后形成的待求解的有限元方程为［K］［φ］=0，其中，φ=φ(x,y,z)表示结点电压，［K］为合成的三维总体系数矩阵，详细内容请参考文献［2］。

1直接法中的克劳特(Crout)三角分解算法

由于［K］是N阶对称、正定矩阵，它可以唯一分解为：

［K］=［L］［U］(1)

其中［L］是下三角矩阵，［U］是主元素为1的上三角矩阵。并且［L］与［U］中的元素有下列关系：

u_pi=l_ip/l_pp (i=1,2,…,n; p=1,2,…,i-1)(2)

即［U］中第i列的第一行至第i-1行各元素分别等于［L］中第i行的第一列至第i-1列各元素除以各列的主元素。由于(2)式，可知［L］与［U］两个矩阵中独立元素的个数与对称矩阵［K］中独立元素的个数一样都是n(n+1)/2个。根据(1)式的等号两边矩阵中的对应元素应相等可得：

K_ij=l_i1u_1j+l_i2u_2j+…+l_i,j-1u_j-1,j+l_ij×1

K_ii=l_i1u_1i+l_i2u_2i+…+l_i,i-1u_i-1,i+l_ii×1

上述两式可改写为：

(3)

(4)

将由(2)式得到的l_ip=u_pil_pp代入(4)式得

(5)

由(2)、(3)、(5)三式实现了(1)式所示的［K］的三角分解，即由(3)计算［L］的副元素，由(5)计算［L］的主元素，由(2)计算［U］的副元素。

现将(1)式代入加入强加边界条件之后的有限元方程［K］［φ］=［R］中，得

［L］［U］［φ］=［R］(6)

令［U］［φ］=［G］(7)

将(7)代入(6)式得

［L］［G］=［R］(8)

于是，解一个对称系数矩阵线性代数方程组的问题化成为解两个三角系数矩阵线性代数方程组(7)和(8)的问题。

又由于［L］是个下三角矩阵，因此利用(8)式便可以自上而下逐个解出中间未知量［G］的各个元素。一般地，由(8)式第i个方程

l_i1g₁+l_i2g₂+…+l_i,j-1g_i-1+l_ijg_i=R_i(9)

便可以求得l_iig_i为

(10)

在求l_iig_i时，g_p(p=1,2,…,i-1)已在此之前解得。

将由(2)式得到的l_ip=u_pil_pp代入(10)式，并令

l_iig_i=f_i, (i=1,2,…,n)(11)

即得

(12)

由此式便可将f₁至f_n逐个解出，再利用由(11)式得来的

g_i=f_i/l_ii, (i=1,2,…,n)(13)

即可解得g₁至g_n。

由于［U］是个上三角矩阵，而［G］又已由(12)和(13)两式解得，便可利用(7)式自下而上逐个解出未知量［φ］的各个元素。一般地，由(7)式中第i个方程

φ_i+u_i,i+1φ_i+1+…+u_i,nφ_n=g_i(14)

便可以求得φ_i为

(15)

在求φ_i时，φ_p(p=n,n-1,…,i+1)已在此之前解得。利用(15)式，便可将φ_n至φ₁逐个解出。

为了便于编制出质量较高的程序，使得所需存贮量较少，运算时间较短，还需对(1)、(3)、(5)、(12)、(13)、(15)各式作进一步的分析。由这些公式可见，［K］中任何一个元素K_ij(或K_ii)仅由(3)、(5)两式计算l_ij(或l_ii)时要用到，在算得l_ij(或l_ii)后就用不着K_ij(或K_ii)了。因此，可以将由这两个公式算得的l_ij(或l_ii)就存放在原来存放K_ij(或K_ii)的位置上，不必另开新的存贮单元。类似地，由(12)式算得f_i后可以存入在R_i的位置上；由(13)式算得g