等参数六面体单元在求解胸腔三维恒定电场中的应用
作者：王海滨 王衍文 程敬之
单位：(西安交通大学 生物工程及仪器系，西安 710049)
关键词：等参数单元；胸腔；电场
生物医学工程学杂志980313 内容摘要 等参数任意六面体单元在有限元分析中有着重要的作用。本文目的在于建立该等参数单元的数学模型，并进行了形状函数偏导数的组合；此外，还具体讨论了应用这种等参数8节点单元求解胸腔三维恒定电场分布的方法。
The Use of Isoparametric Element of Arbitrary Hexahedron
for Solving 3-D Thoracic Potential Distribution
Wang Haibin Wang Yanwen Cheng Jingzhi
, http://www.100md.com
(Biomedical Engineering and Instrument Department Xi'an Jiaotong University,Xi'an 710049)
Abstract Isoparametric element of arbitrary hexahedron plays an important role in finite element analysis.This article devotes itself to establishment of mathematical model and combination of derivative of shape function of such element type.In addition,the article also discusses the usage of isoparametric element for solving 3-D thoracic potential distribution in electrical field.
, 百拇医药
Key words Isoparametric element Thorax Electrical field
1 概 述
采用有限元方法研究胸腔三维恒定电场分布是对心阻抗血流图技术中长期存在的阻抗变化(ΔZ)波形成因及心脏每搏量(SV)的无创估算公式精确性等问题所进行的系列基础性理论研究的一个方面。在这方面，国外代表性的研究工作是Kim[1](1988)建立的人体胸腔阻抗心动图的三维有限元模型。研究中，作者明确指出胸腔阻抗改变与主动脉中血液容积改变之间存在的某种定量的关系，也就是说通过对心阻抗血流图的三维建模研究，一定能够合理的解释心阻抗波形图的成因及推导出SV的精确计算公式。但Kim没给出这个定量关系，且Kim是以解剖图谱为依据来选取与他本人几何尺寸大体相当的人体胸腔CT图片，这使得所计算出的模型数值解有一定的误差，即使这样，他的研究成果也给我们提供了一条修正Kubicek[2]均匀园柱并联膨胀物理模型的新途径，但之后的这十几年，有关人体胸腔阻抗心动图三维有限元模型的研究几乎未见报道，偶尔有一些有关人体胸腔阻抗心动图三维有限差分模型的研究报道[3]。而国内对心阻抗血流图的三维有限元的建模研究尚属空白。本文讨论利用等参数六面体单元求解胸腔三维恒定电场分布的问题。在对人体胸腔从颈根部到剑突处所限区域中截取的任意横截面CT图片进行的三维有限元模型研究中，常常需要用较少量的单元来表示边界区域较复杂的各组织器官，从这个角度看，等参数六面体单元较之四面体或正六面体单元显示了它的优越性，它不但有着四面体的灵活性，又保持了正六面体单元较好地反映计算精度的特点。等参数任意六面体单元与其它形式的常用单元如四面体、正六面体单元在建立数学模型的方法上是不同的。例如，为了将任意六面体映射为对应的正六面体，必须运用坐标变换以建立整体坐标与局部坐标之间的对应关系，但在单元变分计算以后所进行的总体合成、边界约束条件的处理以及线性代数方程组的求解步骤与常用单元是相同的。因此本文只侧重由拉普拉斯方程所描述的人体胸腔三维电场分布问题来讨论等参数任意六面体单元的数学模型的建立和形状函数偏导数的组合等问题，而略去总体合成后各步骤的叙述，同时给出了采用该方法对均匀圆柱形模型所进行的计算机仿真结果。
, 百拇医药
2 数学模型的建立
2.1 整体坐标与局部坐标
图1为一用三维直角坐标xyz表示的基本单元，即八节点任意六面体单元e，八个顶点为1，2，3，4，5，6，7，8。可以用形状函数N′来确定单元中任意点在直角坐标中的位置如下[4]：
式中：N′为形状函数；
N′=[N1′，N2′，N3′，N4′，N5′，N6′，N7′，N3′]；
, 百拇医药
xi,yi,zi(i=1，2，3，4，5，6，7，8)分别为八个节点的坐标。
对于任意六面体单元，假设电场中电压在单元内的分布为包含有乘项xyz的三重一次多项式函数，即：
φ=α1+α2x+α3y+α4y+α5xy+α6yz+α7xz+α8xyz(2)
将xyz坐标中的基本单元e看作是ξηζ坐标中的正六面体(标准单元)12345678的映像，标准单元重心为坐标的原点，坐标是无因次的，标准单元六个面的方程为：ξ=±1、η=±1、ζ=±1，如图2所示。这种映射是针对各单元分别进行的。ξηζ坐标仅适用于单个元素，称为局部坐标：而xyz坐标则包含需整体求解电场的区域，称为整体坐标。对于局部坐标中的标准单元我们又可以用形状函数N表示电压φ在单元内的分布，即：
, http://www.100md.com
式中：
—由有限个节点电压表示的电场分布；
N—电场的形状函数，N=[N1，N2，N3，N4，N5，N6，N7，N8]。
如果N′=N，即确定几何特性的形状函数与确定电场的形状函数有相同的数值，这种单元就是本文所要讨论的等参数单元，并且在以后只用N表示。
图1 整体坐标中的八节点任意六面体单元
Fig 1 8-node arbitrary hexahedron element in the whole coordinate
, 百拇医药
对于图2所示的标准单元，其形状函数可表示如下[5]：
Ni=1/8(1+ξξi)(1+ηηi)(1+ζζi)(i=1，2，3，4，5，6，7，8) (4)
式中：ξi,ηi,ζi为i点(i=1，2，3，4，5，6，7，8)所对应的坐标。
由式(4)可看出，形状函数在对应点处的数值为1，而在其它节点处为零。例如N1在节点1上为1，在节点2，3，4，5，6，7，8处的数值都为零。根据坐标变换的原理，图1所示的任意六面体单元在ξηζ坐标上对应地映射出一个标准单元，而且当雅可比行列式的符号在映射区域内的各点都保持不变时，这种对应关系是一一确定的[4]。即：
, http://www.100md.com
将(4)式代入(1)，则可看出：任意六面体单元的分界面或分界线的几何位置在坐标变换中仍具有连续性，从而保证单元之间不互相重叠，也不会彼此分离，这样相邻单元边界面上的电压由四个端点的电压唯一确定，从而保证了边界电压的协调连续。
图2 局部坐标中的标准单元
Fig 2 Standard element in the local coordinate
2.2 从坐标变换的角度导出形状函数对整体坐标的偏导数[5]
令整体坐标与局部坐标的关系为：
, 百拇医药 根据复合函数求偏微分的规则，可以写出Ni(i=1，2，3，4，5，6，7，8)对ξ，η，ζ的偏导数如下：
其中：
式中：J称为雅可比矩阵，而该矩阵所构成的行列式det J,称为雅可比行列式。让式(7)的两端左乘以J的逆矩阵，则有：
从(8)式根据逆矩阵的运算规则可求得J-1。
在(9)式中，局部坐标内形状函数的偏导数
(i=1，2，3，4，5，6，7，8)可由(4)式求导得出，而(12)式中x,y,z分别对ξ，η，ζ所求偏导数可由(1)式得出，从而从坐标换的角度求得了形状函数对整体坐标偏导数的表达式。
, http://www.100md.com
2.3 泛函在坐标变换中的变化
从以上所述可以看到，经过坐标变换，由拉普拉斯方程所描述的人体胸腔三维电场的泛函的形式也发生了变化，由
式(13)表示对应于整体坐标xyz的泛函；而(14)则表示经过坐标变换后对应于局部坐标ξηζ的泛函。从式(14)可看出，经过坐标变换后，积分的上下限变得很简单，分别为-1及+1，这对于简化变分计算是十分有利的。另外在式(13)、(14)的泛函式子中，微元体积的表示方法也发生了变化，即dxdydz=det Jdξdηdζ。
2.4 单元的变分计算[5]
让Πe代表定义在单元e上的泛函，单元变分计算的目的，就是求泛函数值，并建立起单元的矩阵方程式。
, 百拇医药
让上式分别对φ1，φ2，φ3，φ4，φ5，φ6，φ7，φ8求极值，并将(12)式相应的各种偏导数组合代入积分式子进行积分演算，最后可将泛函对φ1，φ2，φ3，φ4，φ5，φ6，φ7，φ8的极值综合用矩阵方程式表示为：
Keφe=0(17)
式中：Ke为单元系数矩阵：
φe为单元电压列向量， φe=[φ1，φ2，φ3，φ4，φ5，φ6，φ7，φ8]T
, http://www.100md.com
逐个对每个单元进行变分计算，并将每个单元的贡献加到主矩阵的相应位置上去。
表1 圆柱分析仿真模型中电压数值解与理论解的比较
Table 1 Comparison between potential simulated values and theoretical values in the cylindrical model
层数
Z方向坐标
(cm)
电压数值解
(V)
电压理论解
, 百拇医药
(V)
相对误差
(%)
1
0.0
0.0(边界)
0.0
0.0
2
1.44
0.002146
0.002149
0.139601
, 百拇医药
3
2.16
0.00322
0.003223
0.093082
4
2.88
0.004309
0.004297
0.58018
5
4.32
0.006497
, 百拇医药
0.006446
0.791183
6
5.04
0.007569
0.007520
0.651598
7
6.48
0.09714
0.009669
0.465401
, 百拇医药
8
7.52
0.011856
0.011818
0.346925
9
9.68
0.014004
0.013967
0.264914
10
12.36
0.016179
, 百拇医药
0.016135
0.272695
11
13.28
0.017549
0.017508
0.234182
12
14.74
0.021541
0.021507
0.158083
, 百拇医药
13
16.18
0.024758
0.024730
0.11322
14
17.62
0.026307
0.026282
0.095124
15
19.06
, 百拇医药 0.028452
0.028431
0.073861
16
19.78
0.029524
0.029505
0.064399
17
21.22
0.031669
0.031654
, 百拇医药
0.047381
18
21.94
0.032742
0.032728
0.042776
19
22.66
0.033814
0.03381
0.011834
20
, 百拇医药 24.10
0.035959(边界)
0.035959
0.0
3 数学模型的应用及讨论
任意六面体作为基本单元在有限元分析中有着重要的地位，它即具备四面体单元灵活性又具有正立方体单元较精确的特点。任意六面体单元与其它常用单元的主要区别在于建立数学模型的方法和解题途径上，其它步骤大体相同。借助于等参数和坐标变换的概念，六面体单元的变分计算变得较为简单可行了。作者应用六面体单元所编的计算机程序[6]对人体胸腔均匀介质圆柱体模型的三维恒定电场分布进行了仿真数值计算，具体为：进行仿真的均匀介质圆柱体模型的高度为24.10 cm，半径为8.0 cm，电阻率设为ρ=150 Ω*cm；圆柱体上下两个面之间加入2 mA的边界驱动恒定电流；圆柱体共分19层，每层的间距可不相等，每一层由30个单元组成，因此仿真模型共有570个单元，740个结点。把每一层上所有结点电压的平均值结果与理论计算值进行了比较，从表1中可知，数值计算精度的相对误差均少于0.8%，完全满足实际应用的要求，从而表明应用等参数六面体单元来解决电场问题是稳定可靠的，为进一步深入研究心阻抗血流技术中存在的问题打下一个良好的基础。
, http://www.100md.com
参 考 文 献
1 Kim DW,Baker LE,Pearce JA et al.Origins of the impedance change in impedance cardiography by a three-dimensional finite element model.IEEE Trans BME,1988;35(12)∶993
2 Kubicek WG,Karnegis JN,Paterson RP et al.Development and evaluation of an impedance cardiac output system.Aerospace Med,1966;37∶1208
3 Wang L,Paterson RP.Multiple sources of the impedance cardiogram based on 3-D finite difference human thorax models.IEEE Trans BME.1995;42(2)∶141
4 盛剑霓.工程电磁场数值分析.西安：西安交通大学出版社，1986∶87
5 李嘉珩.有限单元法及程序实例.长沙：人民铁道出版社，1979∶63
6 王海滨，王衍文，程敬之.人体胸腔三维恒定电场分布的有限元程序设计方法.北京生物医学工程杂志，1997；17(1)∶13
(收稿：1997-04-08 修回1997-10-22), 百拇医药
单位：(西安交通大学 生物工程及仪器系，西安 710049)
关键词：等参数单元；胸腔；电场
生物医学工程学杂志980313 内容摘要 等参数任意六面体单元在有限元分析中有着重要的作用。本文目的在于建立该等参数单元的数学模型，并进行了形状函数偏导数的组合；此外，还具体讨论了应用这种等参数8节点单元求解胸腔三维恒定电场分布的方法。
The Use of Isoparametric Element of Arbitrary Hexahedron
for Solving 3-D Thoracic Potential Distribution
Wang Haibin Wang Yanwen Cheng Jingzhi
, http://www.100md.com
(Biomedical Engineering and Instrument Department Xi'an Jiaotong University,Xi'an 710049)
Abstract Isoparametric element of arbitrary hexahedron plays an important role in finite element analysis.This article devotes itself to establishment of mathematical model and combination of derivative of shape function of such element type.In addition,the article also discusses the usage of isoparametric element for solving 3-D thoracic potential distribution in electrical field.
, 百拇医药
Key words Isoparametric element Thorax Electrical field
1 概 述
采用有限元方法研究胸腔三维恒定电场分布是对心阻抗血流图技术中长期存在的阻抗变化(ΔZ)波形成因及心脏每搏量(SV)的无创估算公式精确性等问题所进行的系列基础性理论研究的一个方面。在这方面，国外代表性的研究工作是Kim[1](1988)建立的人体胸腔阻抗心动图的三维有限元模型。研究中，作者明确指出胸腔阻抗改变与主动脉中血液容积改变之间存在的某种定量的关系，也就是说通过对心阻抗血流图的三维建模研究，一定能够合理的解释心阻抗波形图的成因及推导出SV的精确计算公式。但Kim没给出这个定量关系，且Kim是以解剖图谱为依据来选取与他本人几何尺寸大体相当的人体胸腔CT图片，这使得所计算出的模型数值解有一定的误差，即使这样，他的研究成果也给我们提供了一条修正Kubicek[2]均匀园柱并联膨胀物理模型的新途径，但之后的这十几年，有关人体胸腔阻抗心动图三维有限元模型的研究几乎未见报道，偶尔有一些有关人体胸腔阻抗心动图三维有限差分模型的研究报道[3]。而国内对心阻抗血流图的三维有限元的建模研究尚属空白。本文讨论利用等参数六面体单元求解胸腔三维恒定电场分布的问题。在对人体胸腔从颈根部到剑突处所限区域中截取的任意横截面CT图片进行的三维有限元模型研究中，常常需要用较少量的单元来表示边界区域较复杂的各组织器官，从这个角度看，等参数六面体单元较之四面体或正六面体单元显示了它的优越性，它不但有着四面体的灵活性，又保持了正六面体单元较好地反映计算精度的特点。等参数任意六面体单元与其它形式的常用单元如四面体、正六面体单元在建立数学模型的方法上是不同的。例如，为了将任意六面体映射为对应的正六面体，必须运用坐标变换以建立整体坐标与局部坐标之间的对应关系，但在单元变分计算以后所进行的总体合成、边界约束条件的处理以及线性代数方程组的求解步骤与常用单元是相同的。因此本文只侧重由拉普拉斯方程所描述的人体胸腔三维电场分布问题来讨论等参数任意六面体单元的数学模型的建立和形状函数偏导数的组合等问题，而略去总体合成后各步骤的叙述，同时给出了采用该方法对均匀圆柱形模型所进行的计算机仿真结果。
, 百拇医药
2 数学模型的建立
2.1 整体坐标与局部坐标
图1为一用三维直角坐标xyz表示的基本单元，即八节点任意六面体单元e，八个顶点为1，2，3，4，5，6，7，8。可以用形状函数N′来确定单元中任意点在直角坐标中的位置如下[4]：
式中：N′为形状函数；
N′=[N1′，N2′，N3′，N4′，N5′，N6′，N7′，N3′]；
, 百拇医药
xi,yi,zi(i=1，2，3，4，5，6，7，8)分别为八个节点的坐标。
对于任意六面体单元，假设电场中电压在单元内的分布为包含有乘项xyz的三重一次多项式函数，即：
φ=α1+α2x+α3y+α4y+α5xy+α6yz+α7xz+α8xyz(2)
将xyz坐标中的基本单元e看作是ξηζ坐标中的正六面体(标准单元)12345678的映像，标准单元重心为坐标的原点，坐标是无因次的，标准单元六个面的方程为：ξ=±1、η=±1、ζ=±1，如图2所示。这种映射是针对各单元分别进行的。ξηζ坐标仅适用于单个元素，称为局部坐标：而xyz坐标则包含需整体求解电场的区域，称为整体坐标。对于局部坐标中的标准单元我们又可以用形状函数N表示电压φ在单元内的分布，即：
, http://www.100md.com
式中：
—由有限个节点电压表示的电场分布；N—电场的形状函数，N=[N1，N2，N3，N4，N5，N6，N7，N8]。
如果N′=N，即确定几何特性的形状函数与确定电场的形状函数有相同的数值，这种单元就是本文所要讨论的等参数单元，并且在以后只用N表示。
图1 整体坐标中的八节点任意六面体单元
Fig 1 8-node arbitrary hexahedron element in the whole coordinate
, 百拇医药
对于图2所示的标准单元，其形状函数可表示如下[5]：
Ni=1/8(1+ξξi)(1+ηηi)(1+ζζi)(i=1，2，3，4，5，6，7，8) (4)
式中：ξi,ηi,ζi为i点(i=1，2，3，4，5，6，7，8)所对应的坐标。
由式(4)可看出，形状函数在对应点处的数值为1，而在其它节点处为零。例如N1在节点1上为1，在节点2，3，4，5，6，7，8处的数值都为零。根据坐标变换的原理，图1所示的任意六面体单元在ξηζ坐标上对应地映射出一个标准单元，而且当雅可比行列式的符号在映射区域内的各点都保持不变时，这种对应关系是一一确定的[4]。即：
, http://www.100md.com
将(4)式代入(1)，则可看出：任意六面体单元的分界面或分界线的几何位置在坐标变换中仍具有连续性，从而保证单元之间不互相重叠，也不会彼此分离，这样相邻单元边界面上的电压由四个端点的电压唯一确定，从而保证了边界电压的协调连续。
图2 局部坐标中的标准单元
Fig 2 Standard element in the local coordinate
2.2 从坐标变换的角度导出形状函数对整体坐标的偏导数[5]
令整体坐标与局部坐标的关系为：
, 百拇医药 根据复合函数求偏微分的规则，可以写出Ni(i=1，2，3，4，5，6，7，8)对ξ，η，ζ的偏导数如下：
其中：
式中：J称为雅可比矩阵，而该矩阵所构成的行列式det J,称为雅可比行列式。让式(7)的两端左乘以J的逆矩阵，则有：
从(8)式根据逆矩阵的运算规则可求得J-1。
在(9)式中，局部坐标内形状函数的偏导数
(i=1，2，3，4，5，6，7，8)可由(4)式求导得出，而(12)式中x,y,z分别对ξ，η，ζ所求偏导数可由(1)式得出，从而从坐标换的角度求得了形状函数对整体坐标偏导数的表达式。, http://www.100md.com
2.3 泛函在坐标变换中的变化
从以上所述可以看到，经过坐标变换，由拉普拉斯方程所描述的人体胸腔三维电场的泛函的形式也发生了变化，由
式(13)表示对应于整体坐标xyz的泛函；而(14)则表示经过坐标变换后对应于局部坐标ξηζ的泛函。从式(14)可看出，经过坐标变换后，积分的上下限变得很简单，分别为-1及+1，这对于简化变分计算是十分有利的。另外在式(13)、(14)的泛函式子中，微元体积的表示方法也发生了变化，即dxdydz=det Jdξdηdζ。
2.4 单元的变分计算[5]
让Πe代表定义在单元e上的泛函，单元变分计算的目的，就是求泛函数值，并建立起单元的矩阵方程式。
, 百拇医药
让上式分别对φ1，φ2，φ3，φ4，φ5，φ6，φ7，φ8求极值，并将(12)式相应的各种偏导数组合代入积分式子进行积分演算，最后可将泛函对φ1，φ2，φ3，φ4，φ5，φ6，φ7，φ8的极值综合用矩阵方程式表示为：
Keφe=0(17)
式中：Ke为单元系数矩阵：
φe为单元电压列向量， φe=[φ1，φ2，φ3，φ4，φ5，φ6，φ7，φ8]T
, http://www.100md.com
逐个对每个单元进行变分计算，并将每个单元的贡献加到主矩阵的相应位置上去。
表1 圆柱分析仿真模型中电压数值解与理论解的比较
Table 1 Comparison between potential simulated values and theoretical values in the cylindrical model
层数
Z方向坐标
(cm)
电压数值解
(V)
电压理论解
, 百拇医药
(V)
相对误差
(%)
1
0.0
0.0(边界)
0.0
0.0
2
1.44
0.002146
0.002149
0.139601
, 百拇医药
3
2.16
0.00322
0.003223
0.093082
4
2.88
0.004309
0.004297
0.58018
5
4.32
0.006497
, 百拇医药
0.006446
0.791183
6
5.04
0.007569
0.007520
0.651598
7
6.48
0.09714
0.009669
0.465401
, 百拇医药
8
7.52
0.011856
0.011818
0.346925
9
9.68
0.014004
0.013967
0.264914
10
12.36
0.016179
, 百拇医药
0.016135
0.272695
11
13.28
0.017549
0.017508
0.234182
12
14.74
0.021541
0.021507
0.158083
, 百拇医药
13
16.18
0.024758
0.024730
0.11322
14
17.62
0.026307
0.026282
0.095124
15
19.06
, 百拇医药 0.028452
0.028431
0.073861
16
19.78
0.029524
0.029505
0.064399
17
21.22
0.031669
0.031654
, 百拇医药
0.047381
18
21.94
0.032742
0.032728
0.042776
19
22.66
0.033814
0.03381
0.011834
20
, 百拇医药 24.10
0.035959(边界)
0.035959
0.0
3 数学模型的应用及讨论
任意六面体作为基本单元在有限元分析中有着重要的地位，它即具备四面体单元灵活性又具有正立方体单元较精确的特点。任意六面体单元与其它常用单元的主要区别在于建立数学模型的方法和解题途径上，其它步骤大体相同。借助于等参数和坐标变换的概念，六面体单元的变分计算变得较为简单可行了。作者应用六面体单元所编的计算机程序[6]对人体胸腔均匀介质圆柱体模型的三维恒定电场分布进行了仿真数值计算，具体为：进行仿真的均匀介质圆柱体模型的高度为24.10 cm，半径为8.0 cm，电阻率设为ρ=150 Ω*cm；圆柱体上下两个面之间加入2 mA的边界驱动恒定电流；圆柱体共分19层，每层的间距可不相等，每一层由30个单元组成，因此仿真模型共有570个单元，740个结点。把每一层上所有结点电压的平均值结果与理论计算值进行了比较，从表1中可知，数值计算精度的相对误差均少于0.8%，完全满足实际应用的要求，从而表明应用等参数六面体单元来解决电场问题是稳定可靠的，为进一步深入研究心阻抗血流技术中存在的问题打下一个良好的基础。
, http://www.100md.com
参 考 文 献
1 Kim DW,Baker LE,Pearce JA et al.Origins of the impedance change in impedance cardiography by a three-dimensional finite element model.IEEE Trans BME,1988;35(12)∶993
2 Kubicek WG,Karnegis JN,Paterson RP et al.Development and evaluation of an impedance cardiac output system.Aerospace Med,1966;37∶1208
3 Wang L,Paterson RP.Multiple sources of the impedance cardiogram based on 3-D finite difference human thorax models.IEEE Trans BME.1995;42(2)∶141
4 盛剑霓.工程电磁场数值分析.西安：西安交通大学出版社，1986∶87
5 李嘉珩.有限单元法及程序实例.长沙：人民铁道出版社，1979∶63
6 王海滨，王衍文，程敬之.人体胸腔三维恒定电场分布的有限元程序设计方法.北京生物医学工程杂志，1997；17(1)∶13
(收稿：1997-04-08 修回1997-10-22), 百拇医药