第十九章有限元法

第十九章 有限元法

有限元法是一套求解微分方程的系统化数值计算方法，它比传统解法具有理论完整可靠，物理意义直观明确，解题效能强等优点.特别是由于这种方法适应性强，形式单纯、规范，所以近年来在电子计算机的配合下，已推广应用到很多工程技术部门和某些科学领域.

本章是从应用的角度来介绍有限元法的基本知识，首先通过典型的位移法阐述有限元法的一般原理与解算过程，然后叙述了剖分单元的技巧，最后介绍与有限元法有关的弹性力学问题.

本章常用符号规定如下(括号内为力学术语或释例):

Ω,表示区域及其边界.

表示区域Ω的单元及其边界.

表示单元的第i个顶点，简记作节点i.

表示系数(刚度)矩阵.

()表示单元的系数(刚度)矩阵.

(x，y，z)表示总体的直角坐标.

(）表示单元的局部坐标.

(，，），（，，，）等表示单元的自然坐标.

(x，y，)表示节点i的直角坐标.

(u，v，w）表示一组待定函数（分别为沿x，y，z方向的位移分量），其列矢量表示为u .

(u，v，w)表示（u，v，w）在单元上的插值函数，其列矢量表示为u.

(u，v，w）表示节点i的函数（位移）值.

{u，v，w}表示节点i的一组参数值，即函数直到某阶导数在节点i 上的值按一定次序排成的列矢量{u}.例如

{u}= {u，v，w}=（u ，u， u ，u ，v ，v， v ，v ，w ，w， w ，w）式中τ表示转置.

{u，v，w}表示{u，v，w}按单元的节点序号排成的列矢量，表示为{u}.

等表示单元的型函数.

{R}表示n次多项式中含变量x，y，z各项按一定次序排成的列矢量，并以表示其中第k个分量.例如二元二次多项式

{}表示n次多项式中各项相应的系数排成的列矢量，并以表示其中第k个分量.例如对于{}，

{}=

{f，g，h}表示与节点参数相对应的一组已知函数及其导数在节点i 上的值排成的列矢量.

{f，g，h}表示{f，g，h}按单元的节点序号排成的列矢量.

,或放在定义或公式之后表示其中函数u，v，w，变量x，y，z或下标i，j，k作循环替换后，该定义或公式仍然成立.

x=x-x_j，,

表示积分（一维，二维或三维）.

diagΦ表示对角线分块矩阵

（Φ为矩阵，0表示其余元素为0）

§1 一般原理与解算步骤

一、变分原理与有限元法

许多物理、力学问题既可以化为微分方程的定解问题，也可以归结为变分问题，即某物理量的极值问题.相应的变分原理（如能量原理）指出两种问题是等价的.对于前者可采用差分法近似求解；对于后者可采用一种直接解法（能量法）.两种方法各有利弊.有限元法正是吸取它们的特点而发展起来的一种新解法.

在变分问题中，物理量一般可表示为待定函数（如位移函数，温度分布）及其导数的积分式.对其积分区域Ω，可仿照差分法的离散化办法（依一定方式划分网格并取节点）把它划分为有限个子区域.（在这里正由于积分表达式代替了定解问题的微分表达式，其区域的划分比差分法中网格的划分要灵活得多）如可取单元为三角形、四边形、四面体、甚至是曲边的）.待定函数及其导数在子区域的某些节点上的数值，称为节点参数值，显然它是待定的.这些节点的选取也比差分法自由.另一方面，再参考直接解法的逼近方式（一般取一类解析函数的线性组合作为近似解，其系数是待定的），而采用在各个内解析的插值函数（一般取多项式，其次数视计算的精度而定，其系数基本上由节点参数值来确定）来逼近.只要适当扩大节点参数值的范围，一般就能满足对插值函数的光滑性与精度的要求.这一点又比直接解法中选取满足一定边界条件的解析函数类来得自由而易于实现.

总之，有限元法是以变分原理为基础吸取差分格式的思想而发展起来的一种有效的数值解法，它把求解无限自由度的待定函数归结为求解有限个自由度（Ω中待定节点参数值的总个数）的待定值问题.具有按一定分布形式的节点及其一定类型的节点参数值的子区域称为单元.

二、在弹性力学问题上的应用（位移法）

弹性体Ω在荷载f（作用于Ω内的体力），q（作用于边界上的面力）等作用下引起小变形，其变形能可表示为

式中{ε}，{σ}分别表示应变、应力各分量排成的列矢量（§5），它们与位移u的线性关系形式上可写成

{ε}=Bu {σ}=DBu

B是微分算子矩阵，D是与弹性系数E,有关的对称矩阵，视问题的性质而定.于是，弹性体的总势能即变形能与外力势能之和，为u的二次泛函：

最小势能原理指出，对满足边界上一定约束条件的所有位移中，以保持力的平衡状态的位移所造成的总势能达到最小值.利用矩阵D的对称性，可知位移函数u就是变分方程

的解,把Ω划分为有限个单元，其节点i的参数值为{u}，并假定在上u（各分量）的插值函数可表示为同类型的多项式：

（1）

它称为位移模式，式中表示对单元的各节点求和，为坐标变量的多项式，随单元的形状与插值的方式而定（即§2~§4中所列的各种型函数），而对角线分块矩阵的行、列数与u的分量个数有关，例如，对空间问题

为33r（r表示每个分量的节点参数值的个数）矩阵.

把（1）代入，得

(2)

设

(3)

(4)

(5)

(6)

式中表示对所有含节点i的单元求和.表示对所有含i，j两个节点的单元求和，当i，j两个节点不在同一个单元时，；j 也可取i，这时=；其余情况，可能对一个或几个单元求和，视区域划分方式而定.

这些系数的力学意义是很明显的：表示在位移模式（1）的假定下，由于节点j(包括i节点本身)具有位移、转动等等变形值通过弹性体单元的作用而在节点i产生的反力*.就是通常所谓单元的刚度矩阵；虽然包括i、j两个节点的单元个数有各种可能，但依定义(4)可知，就是由于节点j的变形值而引起节点i上的反力.对于整个区域的节点j（实际上只有与节点i邻近的几个节点）求和，就得出由于变形而产生节点i上的总反力.同样，(5)(6)右端的积分分别表示体力、面力等外荷载按位移模式(1)通过单元而分配给第i节点的相应的外力^*，简称为在单元上荷载的等价节点力，而其和则是整个区域荷载分配给节点i的总外力.从(2)提出得到

(7)

除了在给定位移的部分边界(在这部分就不给定面力q)上有关的节点参数值(其变分已等于零!)外，由于各个变分的相互独立性，圆括号内各分量(变分为零的部分除外)都应等于零.这正反映出处于平衡状态时，任一非约束节点的反力与外力之和等于零,亦即力在各节点应取得平衡这一客观事实.

最后，按整个区域的节点编号,依序排列待定的节点参数值,除去已加约束的那一部分，从而构成一个总的节点参数值列矢量,对,也划掉相应的行、列,从而构成总的系数矩阵,即所谓总刚度矩阵(撇号表示已做划行、划列的处理)与右端列矢量,于是(7)可写成

由于定义(3)与D的对称性，单元的，而依(4),，即,把其中第m行、第m列划掉后也是对称的，因此，经过划行、划列处理的K是对称的;当区域细分后,大量的第j节点与第i节点不在同一个单元上，则大量的,这表明K是稀疏的.这种对称性与稀疏性正是改进计算方法与提高解题速度的主要根据.

三、有限元解法的主要步骤

(1) 确定变分方程或虚功泛函方程 §5介绍了弹性力学问题位移法中常用的虚功泛函方程的表达式,应用时要明确哪些是待定的位移函数.例如空间问题有三个位移分量是待定函数,平面问题有二个,薄板弯曲或柱体扭转问题只有一个等等.如果不取位移而取内力分布为待定函数,则应根据最小势能原理建立相应的变分方程.对于其他的微分方程定解问题(例如特征值问题,热导问题),则应参照§5末的论述,乘上待定解的变分并利用格林公式把方程改造成为正定的泛函形式,从而把问题归结为等价的变分问题,作为应用有限元解法的出发点.

在虚功泛函方程中,主要是确定微分算子矩阵B、弹性系数矩阵D以及各种荷载分布(包括体力、面力、线力与集中点力),以便列出有关的积分表达式.

(2) 选择单元并划分区域首先要根据区域形状与计算精度选定单元的形态,并用网格把区域划分为大小不一的单元组合.对于应力集中或变化大的部分,网格要划得密些,变化小的部分可疏些,此外还可以利用问题的对称性把定解区域缩小或退化为低维的情况,以减少节点数与计算量.另一方面,对于不均匀的介质,在系数间断或断裂的地方要放置网格,使在每单元内介质系数保持光滑性;而对集中线力、点力以及检测的地方,则要放置节点,以突出其作用.

其次,要按问题对解的光滑性要求,安排单元的节点与节点参数值.一般地说,对空间或平面问题用基本单元较简便,也可用等参数单元,使其边界能适当接近区域的弯曲边界,以减少边界的扰动误差;对薄板弯曲问题一般应取§4所介绍的拟协调单元.总之,要妥善处理好减少总自由度(即节点参数值总个数)与提高计算精度这一对矛盾,做到既减少计算量,又保证解的收敛性.

对单元形状除了§2,§3提出的要求外,还要注意其顶点不能作为相邻单元边上的内点,而且网格要尽量规则些,以减少计算量.

最后,参考§2~§4,确定同节点参数值对应的型函数与位移模式,对复杂的单元还可以利用待定系数法直接把多项式的系数作为广义节点参数进行单元分析.

(3) 单元分析以插值多项式代替原待定函数,每个单元的总势能二次泛函就离散化为节点参数值的二次函数,也就是把原变分方程离散化为的代数方程组,其系数矩阵就是单元的刚度矩阵.这工作是问题离散化的关键,因为整个区域的总势能可以近似地看作各单元总势能之和,显然它也随着离散化为总节点参数值的二次函数.因此,单元分析的主要内容就是根据单元的几何与力学信息计算其系数(刚度)矩阵与等价函数值(节点力)等,为了程序标准化,其节点序号应统一取局部的.

刚度矩阵的元素依(2),(3)式的定义实际上是×的系数块,其中就是第i节点参数值的个数.例如对于空间问题,可取等参数单元,这时B为6×3的一阶偏微分算子矩阵(见()式),,diag为3×3的对角线分块矩阵,于是

为3×6矩阵，而D为6×6的弹性系数矩阵(见(31)式).代入(3)式不难得出3×3的系数块

(k,m=1,2,3)

其中

只要把(1,2,3)与(x,y,z)同时作循环代换就可得出其他的.显然单元刚度矩阵是对称的,因为.注意,对薄板弯曲问题,由于取拟协调单元，各节点参数值个数不一定等同(见§4).对此情况待定函数只有一个挠度,diag=为1×矩阵,而B改为3×1二阶偏微分算子矩阵,于是是×3的矩阵,以节点i,j的型函数代入(3)式也可得系数块的积分公式.例如,节点i,j的型函数为,,这时=3,=1而按相应的B,D公式则有

式中

这是3×1的系数块,型函数一经确定就可以代入而求出积分值.如果以插值多项式系数为广义节点参数,则单元的各个型函数要全部改为多项式的各项求出的积分值不是系数块而是整个单元的广义系数矩阵与等价节点力等,对此,还得按§4所述,经过变换转到上述的形式才能依(4),(5),(6)式求和.

对体力f，面力q同样可按(5),(6)式求出它们的等价节点力.如果在某线段S或某点P上作用有线力分布p(s)或集中力Q,则要补上相应的等价节点力:

以上的这些数值积分都可以通过标准化程序来计算(参看附录).但要注意,其中型函数多取局部坐标,而微分运算又是对整体坐标系(x,y,z),因此在求积前必须通过坐标变换把它们化为局部坐标系的统一表达式(参看§2~§4).此外,在求单元的各项积分值时,节点序号一般可取局部的;但在参加总体合成时,各单元节点的局部序号都必须改成总体的统一标号,因为在(4),(5),(6)的和式中,i,j只能是各节点的总体标号.

（4）边界条件的处理与总体合成前一段只考虑纯荷载支承的情况,如果在部分边界给定弹性支承与几何约束,在参加总体合成之前还要进行必要的处理.实际上,正是这些条件才保证系数矩阵的正定性与解的唯一性.

假定在空间问题中Ω的部分边界取弹性支承,其弹性耦合系数为3×3的正定矩阵,而对另一部分边界给定位移值.这样,在变分方程(2)中多出一式

而少了面力q沿的积分项，因为在上q对总势能已无贡献,由于空间问题中,,因此对弹性支承部分应把系数块改为，其中

至于约束条件则应对系数矩阵作划行划列的处理,其原理如下：

假定区域剖分后落在上的节点的总体标号为,,…,,依约束条件,(n=1,2,…,).从(2)式可以看出,关于这个节点的变分方程已自然满足,不必列出.这样就把(7)式的系数矩阵划去第,,…,行的系数块(n=1,2,…,）或相应的3行;再把已知的位移值代入(7）式中的部分而移到方程的右边,这样又划去第，，…，列的系数块或相应的3列.经过补块与划行划列处理后的系数矩阵仍然对称,而方程组(7)即可具体列式为