1、九九九九 有限元法有限元法是一套求解微分方程的系统化数值计算方法，它比传统解法具有理论完整可靠，物理意义直观明确，解题效能 强等优点.特别是由于这种方法适应性强，形式单纯、 规范，所以近年来在电子计算机的配合下，已推广应用到很多工程技术部门和某些科学领域.本章是从应用的角度来介绍有限元法的基本知识，首先通过典型的位移法阐述有限元法的一般原理与解算过程，然后叙述了剖分单元的技巧，最后介绍与有限元法有关的弹性力学问题.本章常用符号规定如下(括号内为力学术语或释例):, 表示区域及其边 界.表示区域 的 单元及其边界.eiP表示单元 的第 i 个顶点，简记作节点 i.e)(ijK 表示系数(刚度)矩
2、阵.( )表示 单元 的系数(刚度) 矩阵.eij e(x，y，z)表示总体的直角坐标.( ,）表示单元的局部坐标.( 1， 2， 3），（ 1， 2， 3， 4）等表示单元的自然坐标.(x i，y ， )表示节点 i 的直角坐标.iz(u，v，w）表示一组待定函数（分别为沿 x，y，z 方向的位移分量），其列矢量表示 为 u .(u e，v ，w e)表示（u ，v，w）在单元 上的插值函数，其列矢量表示为 u e.e(u i，v ，w i）表示节点 i 的函数（位移）值.u，v，w 表示节点 i 的一组参数值，即函数直到某阶导数在节点 i 上的值按一定次序排成的列矢量u i.例如u i=
3、u，v，w i= （u i ，u xi， u yi ，u xyi ，v i ，v xi， v yi ，v xyi ，w i ，w xi， w yi ，w xyi） 式中 表示转置.u，v，w e表示u ，v，w i按单元 的节点序号排成的列矢量，表示为u e.e)(,Pi等表示单元的型函数.R )(n表示 n 次多项式中含变量 x，y，z 各项按一定次序排成的列矢量，并以 表示其)(nkR中第 k 个分量 .例如二元二次多 项式 Rxy)2(32)2( ),1 )(n表示 n 次多项 式中各项相应的系数排成的列矢量，并以 )(nk表示其中第 k 个分量.例如对于 2R， )2(= , )2(6
4、5)2(43)2(1f，g，h i表示与节点参数相 对应的一组已知函数及其导数在节点 i 上的值排成的列矢量.f，g，h表示f ，g，h i按单元 的节点序号排成的列矢量 .e, 或 放在定义或公式之后表示其中函数 u，v，w，变量 x，y，z 或下标i，j，k 作循环替换后，该定义或公式仍然成立.x ij=x -xj， ,表示积分（一维，二维或三维）.ddiag表示对角线分块矩阵 0（为矩阵，0 表示其余元素 为 0）1 一般原理与解算步骤一、变分原理与有限元法许多物理、力学问题既可以化为微分方程的定解问题，也可以归结为变分问题，即某物理量的极值问题.相应的变分原理（如能量原理）指出两种问题
5、是等价的.对于前者可采用差分法近似求解；对于后者可采用一种直接解法（能量法）.两种方法各有利弊.有限元法正是吸取它们的特点而发展起来的一种新解法.在变分问题中，物理量一般可表示为待定函数（如位移函数，温度分布）及其导数的积分式.对其积分区域 ，可仿照差分法的离散化办法（依一定方式划分网格并取节点）把它划分为有限个子区域 .（在这里正由于积分表达式代替了定解问题的微分表达式，其区域的划分比e差分法中网格的划分要灵活得多）如可取单元为三角形、四边形、四面体、甚至是曲 边的）.待定函数及其导数在子区域的某些节点上的数值，称为节点参数值，显然它是待定的.这些节点的选取也比差分法自由.另一方面，再参考直
6、接解法的逼近方式（一般取一类解析函数的线性组合作为近似解，其系数是待定的），而采用在各个 内解析的插值函数（一般取多项式，其e次数视计算的精度而定，其系数基本上由节点参数值来确定）来逼近.只要适当扩大节点参数值的范围，一般就能满足对插值函数的光滑性与精度的要求.这一点又比直接解法中选取满足一定边界条件的解析函数类来得自由而易于实现.总之，有限元法是以变分原理为基础吸取差分格式的思想而发展起来的一种有效的数值解法，它把求解无限自由度的待定函数归结为求解有限个自由度（ 中待定节点参数值的总个数）的待定值问题.具有按一定分布形式的节点及其一定类型的节点参数值的子区域 称为e单元.二、在弹性力学问题上
7、的应用（位移法）弹性体 在荷载 f（作用于 内的体力）， q（作用于 边界 上的面力）等作用下引起小变形，其变形能可表示为 d21d21V式中， 分别表示应变、应力各分量排成的列矢量（5），它们与位移 u 的线性关系形式上可写成=Bu = DBu)(DB 是微分算子矩阵，D 是与弹性系数 E, 有关的对 称矩阵，视问题的性质而定.于是，弹性体的v总势能即变形能与外力势能之和，为 u 的二次泛函： dd)(21)( qufuuDBFVP最小势能原理指出，对满足边界上一定约束条件的所有位移中，以保持力的平衡状态的位移所造成的总势能达到最小值.利用矩阵 D 的对称性，可知位移函数 u 就是变分方程
8、0)()()() fuBP的解,把 划分为有限个单 元 ，其 节点 i 的参数 值为u i，并假定在 上 u（各分量）的插e e值函数可表示为同类型的多项式：（1）ei iidag它称为位移模式，式中 ei表示对单元 的各节点求和， 为坐标变量的多项式，随 单元的形状与插值的方式而定（即24 中所列的各种型函数），而对角线分块矩阵 的行、列数idag与 u 的分量个数有关，例如，对空间问题 iiii 0dag为 33r（r 表示每个分量的节点参数值的个数）矩阵.把（1）代入 0)(uP，得 e j jii ii ee BDB ddag)dag( uei iie f(2)0diagei iequ
9、设(3)e BDKjiij d)iag()()(4)ijij)(5)iei ieii da)( fff(6) iei iii g)( qq式中 ie表示对所有含节点 i 的单元求和. ije表示 对所有含 i，j 两个节点的单元求和，当 i，j两个节点不在同一个单元时， 0ijK；j 也可取 i，这时 ije= ；其余情况，可能 对一个或几l个单元求和，视区域划分方式而定.这些系数的力学意义是很明显的： jeij)(u表示在位移模式（1）的假定下，由于节点 j(包括 i 节点本身)具有位移、 转动等等变形值 j通过弹性体单元 的作用而在节点 i 产生的反e力 *. )(eijK就是通常所谓单元
10、的刚度矩阵；虽然包括 i、j 两个节点的单元个数有各种可能，但依定义(4 )可知 jeij)(u，就是由于节点 j 的变形值而引起 节点 i 上的反力.对于整个区域的节点 j（实际上只有与节点 i 邻近的几个节点）求和，就得出由于变形而产生节点 i 上的总反力.同样，(5)(6)右端的积分分别表示体力、面力等外荷载按位移模式(1)通过单元 而分配给第 i 节e点的相应的外力 *，简称为在单元 上荷载的等价节点力，而其和 ii,qf则是整个区域荷载e分配给节点 i 的总外力.从(2)提出 iu得到 0)( j iijiiKf(7)除了在给定位移的部分边界(在这部分就不给定面力 q)上有关的节点参
11、数值(其变分已等于零!)外，由于各个变分 iu的相互独立性，圆括号内各分量(变分为零的部分除外)都应等于零.这正反映出处于平衡状态时，任一非约束节点的反力与外力之和等于零,亦即力在各节点应取得平衡这一客观事实.最后，按整个区域的节点编号,依序排列待定的节点参数值,除去已加约束的那一部分，从而构成一个总的节点参数值列矢量 u,对 ijK, iiqf也划掉相应的行、列,从而构成总的系数矩阵,即所谓总刚度矩阵 )(ijK(撇号表示已做划行、划列的处理)与右端列矢量 ,qf,于是(7)可写成 f由于定义(3) 与 D 的对称性，单元的 )()(ejiij，而依(4), jiejieijij KKjii
12、j )()( ，即jiijK,把其中第 m 行、第 m 列划掉后也是对称的，因此， 经过划行、划列处理的 K 是对称的; 当区域细分后,大量的第 j 节点与第 i 节点不在同一个 单元上，则大量的 0jiij,这表明 K 是稀疏的.这种对称性与稀疏性正是改进计算方法与提高解题速度的主要根据.三、有限元解法的主要步骤(1) 确定 变分方程或虚功泛函方程 5 介绍了弹性力学问题位移法中常用的虚功泛函方程的表达式,应用时要明确哪些是待定的位移函数.例如空间问题有三个位移分量是待定函数,平面问题有二个,薄板弯曲或柱体扭转问题只有一个等等.如果不取位移而取内力分布为待定函数,则应根据最小势能原理建立相应
13、的变分方程.对于其他的微分方程定解问题(例如特征值问题,热导问题),则应参照5 末的论述,乘上待定解的变分并利用格林公式把方程改造成为正定的泛函形式,从而把问题归结为等价的变分问题,作为应用有限元解法的出发点.在虚功泛函方程中,主要是确定微分算子矩阵 B、弹性系数矩阵 D 以及各种荷载分布(包括体力、面力、线力与集中点力),以便列出有关的积分表达式 .(2) 选择单元并划分区域 首先要根据区域形状与计算精度选定单元的形态,并用网格把区域划分为大小不一的单元组合.对于应力集中或变化大的部分,网格要划得密些,变化小的部分可疏些,此外还可以利用问题的对称性把定解区域缩小或退化为低维的情况,以减少节点
14、数与计算量.另一方面,对于不均匀的介质,在系数间断或断裂的地方要放置网格,使在每单元内介质系数保持光滑性;而对集中线力、点力以及 检测的地方 ,则要放置节点,以突出其作用.其次,要按问题对解的光滑性要求,安排单元的节点与节点参数值.一般地说,对空间或平面问题用基本单元较简便,也可用等参数单元,使其边界能适当接近区域的弯曲边界,以减少边界的扰动误差;对薄板弯曲问题一般应取4 所介绍的 拟协调单元.总之,要妥善处理好减少总自由度(即节 点参数值总个数) 与提高计算精度这一对矛盾,做到既减少计算量,又保证解的收敛性.对单元形状除了2,3 提出的要求外,还要注意其顶点不能作为相邻单元边上的内点,而且*
15、 反力与外力在这里都是指广义的,即包括弯矩、扭矩等 ,与 节点参数值中的转角、扭矩等广 义位移相对应.网格要尽量规则些,以减少计算量.最后,参考24,确定同节点参数值对应的型函数与位移模式,对复杂的单元还可以利用待定系数法直接把多项式的系数作为广义节点参数进行单元分析.(3) 单 元分析 以插值多项式 eu代替原待定函数,每个单元的总势能二次泛函就离散化为节点参数值 eu的二次函数,也就是把原变分方程离散化为 eu的代数方程组,其系数矩阵就是单元的刚度矩阵 )(eijK.这工作是问题离散化的关 键,因为整个区域的总势能可以近似地看作各单元总势能之和,显然它也随着离散化为总节点参数值 e的二次函
16、数.因此,单元分析的主要内容就是根据单元的几何与力学信息计算其系数(刚度)矩阵 )(ijK与等价函数值(节点力)iei,)(qf等,为了程序标准化,其节点序号应统一取局部的.刚度矩阵的元素 依(2),(3)式的定义实际上是 ir j的系数块,其中 ir就是第 i 节点参数)eij值的个数.例如对于空间问题,可取等参数单元,这时 B 为 63 的一阶偏微分算子矩阵(见( 92)式), ii,diag i为 33 的对角线分块矩阵,于是 00)diag( xyzyyzxBiiiii iiii 为 36 矩阵，而 D 为 66 的弹性系数矩阵(见( 31)式).代入(3)式不难得出 33 的系数块(
17、k,m=1,2,3)yxrKejikmij d),()其中yixjyjxi zjijijiGr Gjiji ),(12),()只要把(1,2,3)与(x,y,z)同时作循环代换就可得出其他的 ),(jikmr.显然单元刚度矩阵 )(eijK是对称的,因为 .注意,对薄板弯曲问题,由于取拟协调单元，各节点参数值个数 ir不一定等同(见),=ijkmjikr4).对此情况待定函数只有一个挠度,diag i= i为 1 i矩阵,而 B 改为 31 二阶偏微分算子矩阵,于是 是 ir3 的矩阵,以节点 i,j 的型函数代入(3)式也可得系数块 )(eij的)diag(iB积分公式.例如,节点 i,j
18、的型函数为 (iii, ,这时 ir=3, j=1 而按相应的0jB,D 公式 则有 eyxvEhK yxxByxiiiiyiixxeij ji d)1(2 2diag,2)dag(3)( 0202222式中 )(12)(2)(2 )(1()( 0221120 000 02222220 yxvxvyxvxy yyxvxvyvx iii iii2 iii 这是 31 的系数 块,型函数一经确定就可以代入而求出积分值.如果以插值多项式系数为广义节点参数,则单元的各个型函数要全部改为多项式的各项 )(nkR求出的积分值不是系数块而是整个单元的广义系数矩阵 )(eK与等价节点力 )(,eqf等,对此,
19、还得按4 所述,经过变换转到上述的形式才能依(4),(5),( 6)式求和.对体力 f，面力 q 同样可按(5),(6)式求出它们的等价节点力 iei,)()(qf.如果在某线段 S或某点 P 上作用有线力分布 p(s)或集中力 Q,则要补 上相应的等价节点力: psiieSii)(dagd)()(以上的这些数值积分都可以通过标准化程序来计算(参看附录).但要注意,其中型函数多取局部坐标,而微分运算又是对整体坐标系(x,y,z ),因此在求积前必须通过坐标变换把它们化为局部坐标系的统一表达式(参看24). 此外,在求单 元的各项积分值时,节点序号一般可取局部的;但在参加总体合成时,各单元节点的
20、局部序号都必 须改成总体的统一标号,因为在( 4),(5),(6)的和式中,i,j 只能是各节点的总体标号.（4） 边 界条件的处理与总体合成 前一段只考虑纯荷载支承的情况,如果在部分边界给定弹性支承与几何约束,在参加总体合成之前还要进行必要的处理.实际上,正是这些条件才保证系数矩阵的正定性与解的唯一性.假定在空间问题中 的部分 边界 取弹性支承,其 弹性耦合系数为 33 的正定矩阵3)(ijcC,而对另一部分边界 给定位移值 ),(wvu.这样,在变分方程(2) 中多出一式1 ei jiiie C3 ddag) u而少了面力 q 沿 的积分项，因为在 上 q 对总势能已无贡献,由于空间问题中
21、,1 1ii,因此对弹性支承部分应把系数块 eijK改为 )()(eijij，其中 CCjijieij ee )ddagd33)( 至于约束条件则应对系数矩阵 )ij作划行划列的处理,其原理如下：假定区域剖分后落在 上的节点的总体标号为 1i, 2, 0Ni,依约束条件,10,iniinu(n=1,2, 0N).从(2)式可以看出 ,关于这 个节点的变分方程已自然满足,不必列出.这样就把(7) 式的系数矩阵 ijK划去第 1i, 2, 0i行的系数块 (n=1,2, 0N）ijK或相应的 3 行;再把已知的位移值 代入(7）式中 的部分而移到方程的右边,这样又0Ninuju划去 第 1i， 2
22、， 0Ni列的系数块或相应的 3 列.经过补块与划行划列处理后的系数矩ijK0阵 仍然对称,而方程组(7)即可具体列式为),(021Nnii01,Nniniiiijij KKn qf其中 仍按(5) 计算,而(6 )式应改为ifi ieeiii 1)( dagqq对薄板弯曲问题，边界条件处理较为复杂，但原理一样.如果还有热效应，则可根据问题的性质把热当量荷载作为体力的一部分加到单元的等价节点力 中去.ie)(f经过这些处理,再按(4),(5),(6)等式求和就可得出关于全部节点参数值的线性方程组,总体合成就算结束.它同单元分析一样,也可以用标准化程序来实现.（5） 方程组的数值求解 方程组的系数矩阵一般是正定的,这就保证了解的唯一性.此外,它还有对称性与稀疏性,只要适当调整节点标号可以使非零的系数集中在矩阵主对角线的附近形成“带状矩 阵”, 利用这 些性质可以大大节省求解的 计算量与存储量.对这类型线性方程组的数值解法主要有两类,即迭代法与直接法.最后,还要对计算成果进行综合与分析.例如根据解出的参数位移值求出单元或节点的应力分布,画图列表进行分析等等.
