三维有限元法计算过程.doc

PAGE ２６ PAGE ７ 三维有限元法计算过程 三维有限元法的计算过程： 1）网格单元剖分； 2）线性插值； 3）单元分析； 4）总体刚度矩阵合成； 5）求解线性方程组等部分组成。 一、偏微分方程对应泛函的极值问题 矿井稳恒电流场分布示意图 主要任务是分析在给定边界条件下，求解稳定电流场的Laplace方程或Poisson方程的数值解，即三维椭圆型微分方程的边值问题： 上述微分方程边值问题等价于下面泛函的极小值问题： 二、网格剖分 图2-4 单巷道地质模型 1、网格单元的类型 图2-5 网格单元类型 2、网格单元剖分原则及其步长选择 因此，网格内的单元剖分应按以下剖分原则 1）、各单元节点（顶点）只能与相邻单元节点（顶点）重合，而不能成为其它单元内点； 2）、如果求解区域对称，那么单元剖分也应该对称； 3）、在场变化剧烈的区域网格剖分单元要密一些，在场变化平缓的区域单元密度应小。 4）、网格单元体的大小变化应逐步过渡。 根据上述剖分原则，以、、坐标轴原点为中心，分别向、、方向的两侧作对称变步长剖分，距越远，步长应越大。常用的变步长方法有： (i0) (i0) 以上各式中为常数，、为同一坐标轴上相邻步长值。以方向为例，可知，正方向与负方向对称，只相差一负号。若令，只要给出距原点最近节点的坐标，由上式即可求出其它相应的步长。同理可求得、方向上的变步长、。 3、网格剖分方法 图2-6 平面内节点编号示意图 图2-7 长方体单元编号示意图 在程序设计的过程中，将采用四面体单元对研究区域进行分析。因此对上面的长方体网格单元还需作进一步的处理。任取一编号为的长方体单元，显然根据的值可求出其八个节点的节点编号。为讨论方便，将长方体单元的各顶点的节点号，将其重新编号为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）。其组合规律为：（，，，），，。 在具体的计算过程中，按长方体网格编号逐一计算，在每个长方体单元内再按以上的组合规律依次对六个四面体计算。这样处理后四面体的总数为。 三、线性插值分析 任取一四面体单元如图2-9所示，设其顶点相应的节点序号为，并设各节点电位为，。对应坐标为，，，，当单元足够小时，设四面体内部电导率为常数、其电位呈现线性变化。即： （2-29） 图2-9 四面体单元示意图 对四个节点有： （2-30） 整理式(2-29)、（2-30）得： （2-31） 式中和都是的坐标函数，其中的值为： ， （2-32） 可知：（）、D只与节点坐标有关。 为了分析式（2-32）所对应的几何意义，现引入面积坐标和体积坐标的概念。如图2-10所示为，、、为其顶点坐标，则三角形的面积为： 图2-10 三角形平面图 设点为三角形中任一点，则 （2-33） 上式值为图2-10中三角形面积的两倍，所以： （2-34） 同理可得：， 定义，，为在三角形内的面积坐标，具有以下性质： 1） （2-35） 2） 3） 4） 与节点对边平行的直线上各点的值不变。 利用面积坐标的定义，可以方便地确定三角单元中任一点的面积坐标，且可建立一个面积坐标网，如图2-11所示。因为面积坐标与通常使用的直角坐标系的选择无关，可利用该坐标网构造任意性线、非线性的插值函数。 图2-11 三角形面积坐标网 采用与上述相似的分析方法，可知()为任一点在四面体内的体积坐标，具有与三角单元面积坐标相似的性质。 （2-36） （2-37） （2-38） （2-39） 上式中为其对应节点所构成四面体的体积。 在式（2-31）假设四面体内电位为线性变化，P点的电位值是通过该点的体积坐标将各节点的电位值线性化计算得出。如果四面体内的电位值为非线性变化，则只要改变相应点体积坐标的系数即可。形如： （2-40） 式中、、