当前位置: 首页 > 期刊 > 《北京生物医学工程》 > 2000年第1期
编号:10257679
一种用于求解人体胸腔三维恒定电场有限元方程组的算法
http://www.100md.com 《北京生物医学工程》 2000年第1期
     作者:王建琪 朱新亚 王海滨 董秀珍 漆家学

    单位:第四军医大学生物医学工程系 西安 710032

    关键词:有限元;算法;计算机;内存

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

    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)
, http://www.100md.com
    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.
, http://www.100md.com
    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]中的元素有下列关系:

    upi=lip/lpp (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)式的等号两边矩阵中的对应元素应相等可得:

    Kij=li1u1j+li2u2j+…+li,j-1uj-1,j+lij×1

    Kii=li1u1i+li2u2i+…+li,i-1ui-1,i+lii×1

    上述两式可改写为: (3) (4)
, http://www.100md.com
    将由(2)式得到的lip=upilpp代入(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个方程

    li1g1+li2g2+…+li,j-1gi-1+lijgi=Ri (9)

    便可以求得liigi (10)

    在求liigi时,gp(p=1,2,…,i-1)已在此之前解得。
, 百拇医药
    将由(2)式得到的lip=upilpp代入(10)式,并令

    liigi=fi, (i=1,2,…,n) (11)

    即得 (12)

    由此式便可将f1至fn逐个解出,再利用由(11)式得来的

    gi=fi/lii, (i=1,2,…,n) (13)

    即可解得g1至gn
, 百拇医药
    由于[U]是个上三角矩阵,而[G]又已由(12)和(13)两式解得,便可利用(7)式自下而上逐个解出未知量[φ]的各个元素。一般地,由(7)式中第i个方程

    φi+ui,i+1φi+1+…+ui,nφn=gi (14)

    便可以求得φi (15)

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

    为了便于编制出质量较高的程序,使得所需存贮量较少,运算时间较短,还需对(1)、(3)、(5)、(12)、(13)、(15)各式作进一步的分析。由这些公式可见,[K]中任何一个元素Kij(或Kii)仅由(3)、(5)两式计算lij(或lii)时要用到,在算得lij(或lii)后就用不着Kij(或Kii)了。因此,可以将由这两个公式算得的lij(或lii)就存放在原来存放Kij(或Kii)的位置上,不必另开新的存贮单元。类似地,由(12)式算得fi后可以存入在Ri的位置上;由(13)式算得gi后可以存放在fi的位置上,也就是Ri的位置上;由(15)式算得φi后可以存放在gi的位置上,也就是Ri的位置上,同样不必另开新的存贮单元。此外,在(5)、(12)、(13)、(15)各式中,只要用到[L]的主元素以及[U]的副元素,在(3)式中只要用到[L]中i行的第j列之前的副元素(即lip)以及[U]的副元素,因此,在用(3)式计算lij时,并不需要将[L]中第i行之前各行的副元素存贮起来。根据这个特点,不必另开存贮[U]的存贮单元,而只须将由(1)式算得的[U]中第i列的第1行至第i-1行的副元素存放在[L]中第i行的第1列至第i-1列的副元素位置上,亦即将upi存放在Kip的位置上。采用这种存贮办法后,需要将(1)、(3)、(5)、(12)、(13)、(15)六式中[U]元素的行、列下标对调一下,即
, 百拇医药
    (3)式成为:(i=1,2,…,n; j=1,2,…,i-1) (16)

    (1)式成为:uip=lip/lpp (i=1,2,…,n; p=1,2,…,i-1) (17)

    (5)式成为: (18)

    (12)式成为: (19)

    (13)式成为:gi=fi/lii, (i=1,2,…,n) (20)
, 百拇医药
    (15)式成为: (21)

    上面前三个公式是用来对[K]进行三角分解的,称为三角分解公式,第4、5两式是前代公式,最后一式是回代公式。

    2 分块求解有限元方程组

    采用有限元法求解时,为了达到一定精度要求,往往离散模型划分的单元较多,结点相应增多,得到的求解方程阶数一般都很高,系数矩阵往往不能全部进入计算机内存进行运算,否则造成计算机经常死机的现象。如何针对现有586/133/8机型来解决这一难点,是确保心阻抗血流图三维有限元仿真研究顺利进行的关键。文中在上述算法的基础上,采用分块技术解决了有限的计算机内存对大矩阵容量运算的限制问题。

    现假设已建立的有限元方程的个数为N=12,其总体系数矩阵的非零元素在图1中用.表示,而按一维变带宽存贮的[K]所需的存贮量NH=8,存贮这NH个元素的一维数组记为SK(NH),采用[K]中主对角线元素序号指示向量MA(N)来表征一维数组SK(NH)的元素与三维总系数矩阵[K]的元素之间的关系,具体为MA(12)=[1 3 6 9 13 18 23 28 31 36 39 41]。
, 百拇医药
    图1 有限元方程分块求解示意图

    若计算机的内存资源只给该有限元方程提供的存贮容量为20,即NSPACE=20,显而易见,图1所示的矩阵必须分块,且第一块可存贮第1行的1号元素到第6行的18号元素(文中以行为存贮单位)。将这6行的K元素在内存中形成后,即可自上而下分解再前代,接着将分解后的元素(即副元素为uij,主元素为lii)存入文件中,以便在自下而上回代时再调用。那么第二块是否可存贮从第7行的19号元素到第10行的36号元素呢?回答是否定的。因为在做这些元素的分解运算时,还需要用到这些元素所在行前面有关各行的已分解好的元素。例如,为了将K74改变为l74(此时K73已改变为l73),由(16)式得:

    l74=K74-l73u43
, http://www.100md.com
    其中要用到第4行第3列的元素u43;又例如,为了将l73改变为u73,由(17)式得:

    u73=l73/l33

    其中要用到第3行第3列的元素l33。那么,在第一块中哪些元素仍需保留在第二块中呢?由分解公式可知,应该是从第7行开始到最后一行为止的所有各行中最左边那个非零元素所在列号(第3列)的对应行(第3行)开始到第6行为止的各行中第3列以右(包括第3列)的那些元素保留在第二块中,亦即图1中用○框起来的10个元素需保留在第二块中。于是,第二块只能新增加第7、8两行的10个K元素。在第7、8两行的10个K元素形成后,即可将它们分解再前代,接着将第二块写入另一文件,以便在回代时调用。同样道理,第三块中包含由第二块中保留下来的用△框起来的6个元素以及新形成的第9至第12行的13个元素。在第9至第12行的13个K元素形成后,即可将它们分解再前代,接着自下而上地回代解出第12至第9个未知量,再调出第二块回代解出第8、7两个未知量,再调出第一块回代解出第6至第1个未知量。
, http://www.100md.com
    以据上述的原理即可用程序实现采用克劳特算法分块求解三维有限元方程的过程。

    3 讨 论

    本程序是在586/133/8兆内存组装机上调试通的。由于机器内存所限,在程序调试过程中,又因所形成的整体系数矩阵的总容量非常大,当按32位浮点数开设一维数组的动态内存超过18000个字节时,计算机经常死机。针对这种情况,对8结点等参数单元程序采用了分块直接解法,在分块直接解法中,整体系数矩阵被分成若干块,随后逐块在内存中形成、分解与前代,再逐块在内存中回代,解得全部结点电压的未知值。由于整体系数矩阵中每一块的容量只是总容量的一部分,因此就能以较少的内存容量计算较大规模的题目。当然,分块解比整块解所需的计算时间略有增加,但它能解决上述死机的问题。由于整体系数矩阵是在内存中分块形成的,这可避免重复形成单元系数矩阵以节省更多的计算时间,在程序设计中采用了先将每个单元的系数矩阵形成后顺序存贮在一个数据文件中,当逐块在内存中形成整体系数矩阵时,再从该文件中顺序读出单元系数矩阵进行集合。另外,对于强加边界条件采用了“取大数”处理方法,即在有强加边界条件约束的结点上直接设置一个大数1030的值,而使[K]的形式不发生任何改变,简化了计算过程。同时,运用该算法开发的软件包已成功地应用于胸腔阻抗血流图三维有限元建模的理论研究中[3]
, http://www.100md.com
    作者简介:王建琪(1962—)男,第四军医大学生物医学工程系电子教研室主任。

    4 参考文献

    [1] 王勖成,邵敏.有限单元法基本原理与数值方法.北京:清华大学出版社,1988,第1版:262-292

    [2] 王海滨,董秀珍,漆家学.用于人体胸腔阻抗血流图研究的三维有限元程序设计.第四军医大学学报,1997,18(增刊):23-25

    [3] 王海滨.心阻抗血流图检测技术的理论研究.西安交通大学博士论文,1998.9:42-61

    (1999-01-25收稿), 百拇医药