三维有限元法计算过程

三维有限元法计算过程三维有限元法的计算过程：1）网格单元剖分；2）线性插值；3）单元分析；4）总体刚度矩阵合成；5）求解线性方程组等部分组成。

一、偏微分方程对应泛函的极值问题矿井稳恒电流场分布示意图主要任务是分析在给定边界条件下，求解稳定电流场的Laplace 方程或Poisson方程的数值解，即三维椭圆型微分方程的边值问题：)()((0)(0)()()(000z z y y x x I F u n un u F z u z y u y x u x Lu w D ---=⎪⎪⎪⎩⎪⎪⎪⎨⎧=+∂∂=∂∂=∂∂∂∂+∂∂∂∂+∂∂∂∂≡ΓΓ+Γδδδγσσσ 上述微分方程边值问题等价于下面泛函的极小值问题：dS U dxdydz fU z U y U x U U J w D ⎰⎰⎰⎰⎰Γ+Γ+ΓΩ+-∂∂+∂∂+∂∂=222221}])()()[(2{][γσσ二、网格剖分∞1ρiih ρ......1、网格单元的类型图2-5 网格单元类型2、网格单元剖分原则及其步长选择 因此，网格内的单元剖分应按以下剖分原则1）、各单元节点（顶点）只能与相邻单元节点（顶点）重合，而不能成为其它单元内点；2）、如果求解区域对称，那么单元剖分也应该对称；3）、在场变化剧烈的区域网格剖分单元要密一些，在场变化平缓的区域单元密度应小。

4）、网格单元体的大小变化应逐步过渡。

根据上述剖分原则，以x 、y 、z 坐标轴原点o 为中心，分别向x 、y 、z 方向的两侧作对称变步长剖分，距o 越远，步长应越大。

常用的变步长方法有：c i x x i i )1(1+=∆-∆+ c x x i i =∆∆+/1(i ≠0)c x x i i =∆-∆+111(i ≠0) 以上各式中c 为常数，1+∆i x 、i x ∆为同一坐标轴上相邻步长值。

以x 方向为例，可知，x 正方向与负方向对称，只相差一负号。

若令00=∆x ，只要给出距原点最近节点的坐标1x ∆，由上式即可求出其它相应的步长i x ∆。

同理可求得y 、z 方向上的变步长i y ∆、i z ∆。

3、网格剖分方法图2-6 平面内节点编号示意图图2-7 长方体单元编号示意图在程序设计的过程中，将采用四面体单元对研究区域进行分析。

因此对上面的长方体网格单元还需作进一步的处理。

任取一编号为m 的长方体单元，显然根据m 的值可求出其八个节点的节点编号。

为讨论方便，将长方体单元的各顶点的节点号，将其重新编号为0、1、2、3、4、5、6、7，如图2-8所示。

图2-8 长方体单元节点编号示意图可知长方体单元可划分为六个四面体单元，其方法为用1、2、5、6顶点组成的平面将长方体分成二个三棱柱体，每个三棱柱体又可分成三个四面体单元，其结果如下：（0，1，2，4）、（1，2，4，5）、（2，4，5，6）、（1，2，3，5）、（2，3，5，6）、（3，5，6，7）。

其组合规律为：（j i +，2/)1(1+++j j i ，2/)5(2j j i -++，j i ++4），1,0=i ，2,1,0=j 。

在具体的计算过程中，按长方体网格编号逐一计算，在每个长方体单元内再按以上的组合规律依次对六个四面体计算。

这样处理后四面体的总数为1)(L *L_n *6Y -。

三、线性插值分析任取一四面体单元如图2-9所示，设其顶点相应的节点序号为),,,(m l j i ，并设各节点电位为k u ，m l j i k ,,,=。

对应坐标为),,(i i i z y x ，),,(j j j z y x ，),,(l l l z y x ，),,(m m m z y x ，当单元足够小时，设四面体内部电导率为常数、其电位呈现线性变化。

即：z y x U 4321αααα+++= （2-29）图2-9 四面体单元示意图对m l j i ,,,四个节点有：⎪⎪⎩⎪⎪⎨⎧+++=+++=+++=+++=mm mm l l l l jj j j i i ii z y x u z y x u z y x u z y x u 4321432143214321αααααααααααααααα （2-30）整理式(2-29)、（2-30）得：m m l l j j i i u N u N u N u N U +++= （2-31）式中u 和k N 都是的),,(z y x 坐标函数，其中k N 的值为：)(1z d y c x b a DN k k k k k +++=，m l j i k ,,,= （2-32） 可知：k k k k d c b a ,,,（m l j i k ,,,=）、D 只与节点坐标有关。

为了分析式（2-32）所对应k N 的几何意义，现引入面积坐标和体积坐标的概念。

如图2-10所示为ijl ∆，),(i i y x 、),(j j y x 、),(l l y x 为其顶点坐标，则三角形ijl ∆的面积为：图2-10 三角形平面图ll j j iiijl y x y x y x 11121=∆ 设),(y x p 点为三角形中任一点，则 ljl j l lj ji i i x x zy y x y x y xz c x b a 1111+-=++llj j y x y x yx111= （2-33）上式值为图2-10中三角形pjl ∆面积的两倍，所以：ijlpjlijl i i i i y c x b a y x N ∆∆=∆++=2)(),( （2-34） 同理可得：ijlpil j y x N ∆∆=),(，ijlpij l y x N ∆∆=),(定义),(y x N i ，),(y x N j ，),(y x N l 为),(y x p 在三角形ijl ∆内的面积坐标，具有以下性质：1） 1),(),(),(=++y x N y x N y x N l j i （2-35） 2） 1),(≤y x N i 3） ⎩⎨⎧≠==nk n k y x N n n k 01),(4） 与节点i 对边平行的直线上各点的),(y x N i 值不变。

利用面积坐标的定义，可以方便地确定三角单元中任一点的面积坐标，且可建立一个面积坐标网，如图2-11所示。

因为面积坐标与通常使用的直角坐标系的选择无关，可利用该坐标网构造任意性线、非线性的插值函数。

图2-11 三角形面积坐标网采用与上述相似的分析方法，可知),,(z y x N k (m l j i k ,,,=)为任一点),,(z y x p 在四面体ijlm 内的体积坐标，具有与三角单元面积坐标相似的性质。

ijlm pjlm i D D z y x N =),,( （2-36）ijlm pilm j D D z y x N =),,( （2-37）ijlm pijm l D D z y x N =),,( （2-38）ijlmpijl m D D z y x N =),,( （2-39）上式中ijlm D 为其对应节点所构成四面体的体积。

在式（2-31）假设四面体内电位为线性变化，P 点的电位值是通过该点的体积坐标将各节点的电位值线性化计算得出。

如果四面体内的电位值为非线性变化，则只要改变相应点体积坐标的系数即可。

形如：m m l l j j i i u QN u PN u ON u SN U +++= （2-40）式中S 、O 、P 、Q 为P 点体积坐标对应的系数，由于四面体内电位的非线性变化规律比较复杂，系数值的确定比较困难，因此在以下各章节的讨论中只考虑电位值的线性变化，对于非线性的变化规律将以后讨论。

四、网格单元分析整个求解区域Ω被剖分为一系列小四面体单元后，泛函J （u ）也相应地被离散为各单元上的泛函)(u J e ，故有∑=me u J u J )()(。

m 为四面体单元总个数。

在任一四面体单元上泛函形式为：⎰⎰⎰⎰⎰ΓΩ+-∂∂+∂∂+∂∂=eedS u dxdydz fu z u y u x uu J e 222221}])()()[(2{)(σγσ对任意节点的微分得：⎰⎰⎰⎰⎰ΓΩ∂∂+∂∂-∂∂∂∂⋅∂∂+∂∂∂∂⋅∂∂+∂∂∂∂⋅∂∂=∂∂eedsu uu dxdydz u u f z u u z u y u u y u x u u x u u u J kk k k k k e σγσ})]()()([{)( （2-42）其中m l j i k ,,,= 令 dxdydz zuu z u y u u y u x u u x u k k k e)]()()([∂∂∂∂⋅∂∂+∂∂∂∂⋅∂∂+∂∂∂∂⋅∂∂=Φ⎰⎰⎰Ωσ 则：⎰⎰⎰⎰⎰ΓΩ∂∂+∂∂-Φ=∂∂eeds u uu dxdydz u u f u u J k k k e σγ)( （2-43）现对各项分别进行讨论： 1、对于Φ项由前面的讨论可知：∑=kk k U N U ，)(1z d y c x b a DN k k k k k +++=所以m m l l j j i iU x N U x N U x N U x N x U ∂∂+∂∂+∂∂+∂∂=∂∂ （2-44） xN x UU k k ∂∂=∂∂∂∂)( （2-45） 又因为：Db x N k k =∂∂，k k N U U=∂∂ 所以：)(1)(2m m k l l k j j k i i k k U b b U b b U b b U b b Dx U U x U +++=∂∂∂∂⋅∂∂ （2-46） 同理可得：)(1)(2m m k l l k j j k i i k k U c c U c c U c c U c c D y U U y U +++=∂∂∂∂⋅∂∂ （2-47） )(1)(2m m k l l k j j k i i k k U d d U d d U d d U d d Dz U U z U +++=∂∂∂∂⋅∂∂ （2-48） 则：)()()(zUU z U y U U y U x U U x U k k k ∂∂∂∂⋅∂∂+∂∂∂∂⋅∂∂+∂∂∂∂⋅∂∂ )]()()[(12m k m l k l j k j i k i m k m l k l j k j i k i m k m l k l j k j i k i U d d U d d U d d U d d U c c U c c U c c U c c U b b U b b U b b U b b D+++++++++++=])()()()[(12m k m k m k m l k l k l k l j k j k j k j i k i k i k i U d d c c b b U d d c c b b U d d c c b b U d d c c b b D+++++++++++=（2-49） 记：mk m k m k m l k l k l k l jk j k j k j i k i k i k i k U d d c c b b U d d c c b b U d d c c b b U d d c c b b F )()()()(+++++++++++=所以：kek F V dv F De362σσ==Φ⎰⎰⎰Ω （2-50）2、对于dxdydz u ufek⎰⎰⎰Ω∂∂项由前面分析可知：)()()(000z z y y x x I f ---=δδδ，k kN U U=∂∂所以有： dxdydz z z y y x x z y x N I dxdydz u ufVek Vek)()()(),,(000---=∂∂⎰⎰⎰⎰⎰⎰δδδ （2-51） 由δ函数的性质可知：I z y x N dxdydz u ufk Vek⋅=∂∂⎰⎰⎰),,(000 ⎩⎨⎧=),,(0),,(k k k k k k z y x z y x I 供电点不在点供电点在点 （2-52）3、对于⎰⎰Γ∂∂eds u uukσγ项由式（2-25）分析可知，该项只与研究区域Ω的第三类边界条件有关，所以当四面体完全位于求解区域内部时（即非边界单元），该项值为零。