第二章 弹力学问题有限元方法的一般原理和表达格式.doc

第二章 弹性力学问题有限元方法的一般原理和表达格式2.1 引言 本章将讨论通过弹性力学变分原理建立弹性力学问题有限元法列式的基本步骤。最小位能原理的未知场变量是位移，以结点位移为基本未知量，并以最小位能为基础建立的有限单元位移元。它是有限元方法中应用最普遍的单元。对于一个力学或物理问题，在建立其数学模型以后，用有限元方法对它进行分析的首要步骤是选择单元形式。平面问题三结点三角形单元是有限元方法最早采用，而且至今仍经常采用的单元形式。我们将以此作为典型，讨论如何应用广义坐标建立单元位移模式与位移插值函数，以及如何根据最小位能原理建立有限元求解方程的原理、方法与步骤，并进而导出弹性力学问题有限元方法的一般列式。2.2 弹性力学平面问题的有限元列式2.2.1 单元位移模式及插值函数典型的三结点三角形单元结点编码为i,j,m。每个结点有两个位移分量，如图2.2所示。每个结点的位移可用位移矢量表示，即 每个单元有6个结点位移分量（称为6个自由度），于是单元结点的位移向量可表示为为单元结点位移列阵。1单元的位移模式和广义坐标在有限元方法中单元的位移模式，是指在单元内位移的插值函数，其一般形式采用多项式作为近似函数，因为多项式运算简单，并且随着项数的增多，可以逼近任何一段光滑的函数曲线。假设3结点三角形单元位移模式选取一次多项式 （2.2.1）它的矩阵形式是 （2.2.2）其中 ，由于三个结点也在单元内，满足位移模式，于是得 （2.2.3）上式是关于的线性方程组。是待定常数，也称为广义坐标。它可由（2.2.3）式求出。上式的系数行列式是 （2.2.4）上式中当i,j,m的编号顺序与坐标顺序一致时（此处为逆时针），A值为正，其大小为三角形面积，因此为了方便单元的编号一般与坐标的转向一致。 （2.2.5）同理，y也有三个线性方程组 由上面方程组可求得，即 （2.2.6）在（2.2.5）和（2.2.6）式中 （2.2.7）上式下标轮换，可得及。2位移插值函数将求得的广义坐标代入（2.2.1）式中，则位移函数可表示为结点位移的函数，即 （2.2.8）其中 （2.2.9）称为单元的插值函数或形函数，这里它是的一次函数，其中，及是常数，由表达式可知，它完全由单元的大小和方位确定，一旦单元确定了，这些常数也完全确定。 （2.2.9）中的单元面积可表示为 （2.2.10）（2.2.8）式的矩阵形式是 （2.2.11）称为插值函数矩阵或形函数矩阵。插值函数具有如下性质：（1）在结点上插值函数的值有 （2.2.12）也就是说在i结点上，在j,m点上。这一点与位移插值函数的表达式一致。（2）在单元中任一点各插值函数之和应等于1，即 （2.2.3）以上结论可以由下面得到解释，若单元发生刚体位移，如x方向有刚体位移，显然单元内任何一点的位移都应该为，当然结点位移也为，即，由（2.2.8）式得于是有。若插值函数不满足此要求，在不能反映单元的刚体位移，用以求解必然得不到正确的结果。单元内的性质称为规一性。（3）对于现在的单元，插值函数是线性的，在单元内部及单元的边界上位移也是线性的，可由结点上的位移值唯一确定。由于相邻单元公共结点的位移是相等的，因此保证了相邻单元公共边界上位移的连续性。3应变矩阵和应力矩阵 应变确定了单元位移后，可以很方便地利用几何方程和物理方程求得单元的应变和应力。在（1.4.21）式的几何方程中，位移用（2.2.11）式代入，得到单元应变为 （2.2.14）B称为应变矩阵，L是平面的微分算子。应变矩阵的分块矩阵是 （2.2.15）对（2.2.9）式求导得 ， （2.2.16） 代入（2.2.15）式得 （2.2.17）3结点三角形单元的应变矩阵为 （2.2.18）式中，及由单元结点坐标确定，当单元确定了，这些常数也完全确定，因此B市常量矩阵。而应变，当单元结点位移确定后，单元内的应变都是常量，因此3结点三角形单元为常应变单元。在应变变化较大的区域，单元划分应适当加密，否则不能反映应变的真实变化而导致较大误差。 应力单元应力可以根据物理方程求得， 而对于平面问题式中，因此，对于平面应力问题，弹性矩阵D为由（2.2.14）可知，从而可得 （2.2.19）上式中S写成分块矩阵的形式，则 （2.2.20）S称为应力矩阵。写成分块矩阵的形式 （2.2.21）对于平面应变问题：，对于具体计算时只要对不同的赋值，就得到不同的平面问题。 从而可以到，三大物理参量，都可以用单元结点位移向量表示：由于N，B，S都是已知的矩阵，只要求得，则单元内的位移、应变和应力就可以就得，问题是：如何求结点位移向量2.2.2 利用最小位能原理建立有限元方程 最小位能原理的泛函总位能的表达式如（1.4.51） （1.4.51）在平面问题中最小位能原理的总位能与（1.4.51）式的形式一致。现改为矩阵形式： （2.2.24）其中，t是二维体的厚度；f是体积力；T是表面力。积分域也作了相应的改变。对于离散模型，系统位能是各单元位能的和，将（2.2.11）及（2.2.14）式代入（2.2.24）式，即得到离散模型的总位能为 （2.2.5）将结构总位能的各项矩阵表达成各个单元总位能的各对应项矩阵之和，隐含着要求单元各项矩阵的阶数和结构各项矩阵的结束（自由度数）相同。令 （2.2.28）和称为单元的刚度矩阵和单元等效结点载荷列阵。于是 （2.2.25b）由于每个单元的结点不同，相应的结点位移不同，为了统一起见，引入一个单元结点自由度和整体结构结点自由度的转换矩阵G，从而将单元结点位移列阵用整体结构结点位移列阵表示，即 （2.2.26）其中（2.2.27）于是（2.2.26）（2.2.28）式代入（2.2.25b）式，则离散形式的总位能可表示为 （2.2.29）令 （2.2.30）K和P分别称为整体结构刚度矩阵和结构结点载荷列阵，于是总位能可表示为 （2.2.31）由于离散形式的总位能的未知量是结构的结点位移，根据变分原理，泛函取驻值的条件是它的一次变分为零，即由于的任意性，则欲使上式成立，只有 （2.2.32）上式表明：总位能对任一自由度分量求导都为零，若结点数为n，总的自由度数为2n（平面问题），即得2n个方程。其矩阵形式为 （2.2.33）由（2.230）式可以看出结构整体刚度矩阵K和结构结点载荷列阵P都是由单元刚度矩阵和单元等效结点载荷列阵集合而成。2.2.3 单元刚度矩阵1单元刚度矩阵由（2.2.28）式定义的单元刚度矩阵，由于应变矩阵B对于3结点三角形单元是常量矩阵，因此有 （2.2.34）将弹性矩阵D和应变矩阵B代入上式，任意分块矩阵可表示成 （2.2.35）显然 （2.2.37）2单元刚度矩阵的力学意义和性质 为了进一步理解单元刚度矩阵的物理意义，同样可以利用最小位能原理建立一个单元的有限元方程 （2.2.38）是单元等效结点载荷（体力和面力引起的等效结点力）， 是其他单元对该单元的作用力，和统称为单元结点力。 （2.2.40）这是单元结点平衡方程，每个结点在x和y方向上各有一个平衡方程，3个结点共有6个平衡方程。令，而其余结点位移皆等于零，由上式可得 （2.2.41）上式表明：单元刚度矩阵第一列的元素的物理意义是当i结点在x方向有一单位位移，其他结点位移都为零时，在单元各结点位移方向上施加结点力的大小。当然，单元在这些结点力作用下应处于平衡，因此在x和y方向上结点力之和应为零，即在x方向 在y方向 （2.2.42）单元刚度矩阵的性质：（1）对称性这一性质在迦辽金形式的加权余量法和里兹法已加以说明，在3结点三角形单元的刚度矩阵的具体表达式（2.2.37）已加以证明。实际上对于各种形式的单元都普遍具有这种对称性。（2）奇异性 由（2.2.42）可以见到，由于单元刚度矩阵的每一列都代表结构的平衡，各列元素之和等于零。因此单元刚度矩阵是奇异的。因此有限元方程组的解不唯一，而是包含了一系列的刚体位移，为消除结构的刚体位移，结构应有足够的约束。（3）主元恒大于零 单元刚度矩阵对角线上的元素称为主元，由（2.235）及（2.2.36）式可见他们是恒正的。 的意义是要使i结点在x方向的位移1，在x方向所施加的结点力。显然恒为正。例2.1 给定3结点三角形如图2.3 所示。边长厚度，材料常数，。计算单元的刚度矩阵，并验证矩阵的对称性，奇异性和主元恒正。解：按（2.2.7）式计算，得到2.2.4 单元等效结点载荷列阵 由（2.2.28）式得到单元等效结掉载荷是体积力的等效结点载荷：其中 面积力的等效结点载荷： 常见的载荷如下：1均质等候单元的自重单元的单位体积重量为，坐标方向如图2.4所示。按照（2.2.28）式，应有 （2.2.45）其中，每个结点的等效结点载荷为 （2.2.46）自重的等效结点载荷是 （2.2.47）2均布压力 设在ij边上有一压力q，如图2.5所示。ij边长为l，与x轴的夹角为。压力q在x和y方向的分量分别为作用在单元边界上的面力为 （2.2.48）在边界上取局部坐标s，由面积坐标可知，沿ij边插值函数可以写作， （2.2.49）将（2.2.28）代入（2.2.48）及（2.2.49）式代入（2.2.28），得 （2.2.50）3x方向均布力（如图2.6） （2.2.51）4x方向三角形分布载荷 载荷作用在ij边，如图2.7所示。这时边界面积力写成局部坐标s的函数，即则单元等效结点载荷为 （2.2.52）2.2.5 结构刚度矩阵和结构结点载荷列阵的集成 （2.2.30）式给出了结构刚度矩阵和结构结点载荷列阵由单元刚度矩阵和单元等效结点载荷列阵集成的表达式。1单元刚度矩阵的转换刚度矩阵的转换表示为 （2.2.53）其中n为结构结点总数；i,j,m为单元编号。在上式的矩阵中除标明的9个子块外，其它位置上的元素均为零。i,j,m均为双行和双列。单元矩阵的这种变换起到两个作用：（1）将单元刚度矩阵扩大到与整个结构刚度矩阵同阶，便于矩阵相加；（2）将单元刚度矩阵中各子矩阵按照单元结点的实际编码安放在扩大的矩阵中。即对号入座。2单元等效结点载荷列阵的转换单元等效结点载荷列阵的转换表示为 （2.2.54）i,j,m均为双行，单元等效结点载荷包括体积力和面积力等的等效结点载荷。由（2.2.54）可见单元等效结点载荷列阵的转换是将单元结点载荷扩大到与结构载荷列阵同价，并将单元结点载荷按结点自由度顺序入位。即对号入座。3结构刚度矩阵和结构载荷列阵的集成举例：4结构刚度矩阵的性质和特点（1）对称性（2）奇异性；（3）稀疏性（4）非零元素在矩阵对角线附近呈带状分布。半带宽：，D为相邻结点的最大差值2.2.6 引入位移边界条件 对于求解位移场问题，由于最后所得的有限元方程的系数矩阵为奇异阵，所表示的物理含义：满足方程的位移包含刚体位移。为了求解一个具体结构的受力行为，应有足够的约束，以消除结构的刚体位移，从而消除刚度矩阵的奇异性。1直接代入法在方程组（2.2.33）中将已知结点位移的自由度消去，得到一组修正方程，用以解其他待定的结点位移。其原理是按结点位移已知和代定重新组合方程 （2.2.55）其中为待定结点位移， 为已知结点位移；由刚度矩阵的对称性可知。由（2.2.55）式展开，得 （2.2.56）由于已知，最后的求解方程可写成 （2.2.57）其中， （2.2.58）若总体结点位移为n个，其中已知结点位移m个，则得到一组求解n-m个待定结点位移的修正方程组，为n-m阶方阵。修正方程组的意义是在原来方程中，只保留与未知结点位移相应的n-m个方程，并将方程中左端的已知位移和相应刚度系数的乘积移至方程右端作为载荷的修正项。这种方法：组成的新的方程阶数降低了，但结点位移的顺序已被破坏，这给编程带来一定的麻烦。2对角元素改1法（主元置1法） 当给定位移是零位移时，可以在系数矩阵K中与零结点位移相对应的行和列中，将主对角线元素改为1，其它元素改为0，并在载荷列阵中将与零结点位移相对应的元素改为0。例如： （2.2.59）3对角元素乘大数法（主元赋大值） 当有结点位移为给定位移，对第j个方程的方程的对角线元素乘以大值，并将用代替，则 （2.2.60）则第j个方程变为由于，项比其它项要大得多，于是总是能保证对于多个给定位移时则按其所在位置进行修正。2.2.7 线性代数方程组的求解及应力计算有限元求解方程在引入位移边界条件，消除了K的奇异性后，就可以求解了。具体解法将在第六章专门讲解。2.3 广义坐标有限元的一般格式2.3.1 广义坐标有限元位移模式的选择和插值函数的构造 常用的二维单元有三角形或矩形，常用的三维单元有四面体、五面体或平行六面体。同样形状的单元还可有不同的单元结点数，如二维三角形3结点单元、三角形6结点单元、三角形10结点单元，因此单元种类较多。图2.11中列举了一些常用的单元形式。如何选择合适的单元进行计算，涉及求解问题的类型、对精度的要求以及经济性多方面的因素。当然，不同的单元位移模式也不相同。 本节讨论选择广义坐标有限元位移模式的一般原则，以及建立位移插值函数的一般步骤。1选择广义坐标有限元位移模式的一般原则（1）广义坐标是由结点场变量确定的，因此它的个数应与结点自由度数相等。如3结点三角形单元有六个自由度，广义坐标（）应取6个，因此两个方向的位移各取三项多项式。对4结点的矩形单元广义坐标数为8，位移函数可取四项多项式作为近似函数。（2）选择多项式时，位移模式必须完备，即包含可能有的变形。对于类问题，就是位移模式应有常数项和坐标的一次项。位移模式中的常数项和一次项反映了单元刚体位移和常应变特性。当单元数趋于无穷时，单元缩小趋于一点，此时单元应变应趋于常应变。（3）多项式的选取应由低阶到高阶，尽量选取完全多项式以提高单元精度。（4）多项式的选取不应偏惠于某一坐标，即x或y应等同出现。一般可根据巴斯卡三角形（或杨辉三角形）选取。具体见表2.1 2建立广义坐标有限元位移插值函数的一般步骤）假设位移模式；）将各结点坐标代入，得到关于广义坐标的线性方程组，从而求出广义坐标；）将广义坐标回代入一般位移模式中得到由单元结点位移列阵所表示的位移模式；）由位移模式，由矩阵形式的几何方程，求导数可得到应变矩阵B，即）由矩阵形式的物理方程，则弹性矩阵乘以应变向量，得。2.3.2 弹性力学问题有限元分析的执行步骤在根据问题的类型和性质选定了单元的形式，并构造了它的插值函数以后，可按以下步骤对问题进行有限元分析。（1）对结构进行离散。按问题的几何特点和精度要求等因素划分单元并形成网格，既将原来的连续体离散为在结点处互相联结的有限单元组合体。（2）形成单元的刚度矩阵和等效结点载荷列阵。单元刚度矩阵的一般形式为 （2.3.7）单元等效结点载荷的一般形式为 （3）集成结构的刚度矩阵和等效结点载荷列阵 （2.3.10）其中是直接作用于结点上的集中力。（4）引入强制边界条件（给定位移）。（5）求解有限元方程，得到结点位移。 （2.3.12）（6）计算单元应变和应力。（7）进行必要的后处理。如果单元有初应力和初应变，则则单元的位能为对于具有初应力和初应变的情况，由于初应力或初应变不会随真实的位移而变，考虑其应变能时，不必除以2，于是总位能为令由得式中包括边界上具有面力的表面的等效结点力。此时只要对结点力进行修正。2.4 有限元解的性质和收敛准则 有限元法作为求解微分方程的一种数值方法可以认为是里兹法的一种特殊形式，不同之处在于有限元法的试探函数是定义于单元（子域）而不是全域。因此有限元解的收敛性可以与里兹法的收敛性对比进行讨论。里兹法的收敛条件是要求试探函数具有完备性和连续性，即如果试探函数满足完备性和连续性要求，当试探函数的项数时，则里兹法的近似解将趋于微分方程的精确解。现在要研究有限元解的收敛性。在有限元法中，场函数的总体泛函是单元泛函集成的，如果采用完全多项式（无穷多项）作为单元的插值函数，则有限元解在一个有限尺寸的单元内可以精确地和精确解一致。但是实际上有限元的试探函数只能取有限项多项式，因此有限元解只能是精确解的一个近似解答。有限元解的收敛准则需要回答的是，在什么条件下当单元尺寸趋于零时，有限元趋于精确解。下面仍以含有一个待求标量场函数为例，微分方程是 （2.4.1）相应的泛函是 （2.4.2）假定泛函中包含和它的直至m阶的各阶导数是非零的，则近似函数至少必须是m次多项式。若取p次完全多项式为试探函数，则必须满足。假设仅是x的函数，则及其各阶导数在一个单元内的表达式为： （2.4.3）由上式可见，因为是p次完全多项式，所以它的直至m阶导数的表达式中都包含有常数项。当单元尺寸趋于零时，在每一单元内及其m阶导数将趋于精确解，即趋于常数。因此，每一个单元的泛函有可能趋于它的精确解。如果试探函数还满足连续性要求，那么整个系统的泛函将趋于它的精确解。即解是收敛的。收敛准则：准则1完备性要求。如果出现在泛函中场函数的最高阶导数是m阶，则有限元解收敛的条件之一是单元内场函数的试探函数至少是m次多项式。或者说试探函数中必须包括本身和直至m阶导数为常数的项。当单元的插值函数满足上述要求时，称这样的单元是完备的。准则2协调性要求。如果出现在泛函中的最高阶导数是m阶，则试探函数在单元交界面上必须具有连续性，即在相邻单元的交界面上函数应有直至m-1阶的连续导数。当单元的插值函数满足上述要求时，称这样的单元是协调的。当单元为完备的协调单元，则有限元解收敛，即细分单元其解趋于精确解。2.4.2 收敛准则的物理意义 在平面问题中，泛函中出现的是位移u和v的一次导数，即应变，因此。收敛准则1要求插值函数或位移函数至少是x,y的一次完全多项式。我们知道位移及其一阶导数为常数的项是代表与单元的刚体位移和常应变状态相应的位移模式。实际分析中，各单元的变形往往包含着刚体位移，同时单元尺寸趋于无穷小时各单元的应变也趋于常应变。所以完备性要求由插值函数所构成的有限元解必须能反映单元的刚体位移和常应变状态。若不能满足上述要求，那么赋予结点以单元刚体位移（零应变）或常应变的位移模式时，在单元内部将产生非零或非常值的应变，这样有限元解将不可能收敛于精确解。应该指出，在Bazeley等人开始提出上述收敛准则时，是要求在单元尺寸趋于零的极限情况下满足完备性收敛准则，如果将此收敛准则用于有限尺寸时，将使解的精度得到改进。对平面问题，协调性要求是连续性，即要求位移函数u,v的零阶导数，也就是位移函数自身在单元交界面上是连续的。如果在交界面上位移不连续表现为当结构变形时将在相邻单元间产生缝隙或重叠，这意味着将产生无限大的应变，这时应该将发生在交界面上的附加应变能补充到系统的应变能中去。但在建立泛函时，没有考虑到这种情况，只考虑了产生于各个单元内部的应变能。因此，当边界上位移不连续时，则有限元解就不可能收敛于精确解。可以看出，最简单的3结点三角形单元插值函数既满足完备性要求，也满足协调性要求，因此单元的解是收敛的。应当指出，对于二、三维弹性力学问题，泛函中出现导数是一阶。对于近似的位移函数的连续性要求仅是连续性，这种只要求函数自身在单元边界连续的要求很容易得到满足。而当泛函中出现导数高于一阶（如板壳，泛函中出现的导数是2阶）时，则要求试探函数在单元交界面上具有连续的一阶或高于一阶的导数，即具有或更高阶的连续性，这时构造函数比较困难。在某些情况下，可以放松对协调性的要求，只要这种单元能通过分片试验，有限元解仍然可以收敛于正确的解答。这种单元称为非协调单元。2.4.3 位移解的下限性质 以位移为基本未知量，并基于最小位能原理建立的有限元称为位移元。通过系统总位能的变分过程，可以分析位移元的近似解与精确解偏离的下限性质。 系统总位能的离散形式为 （2.4.4）由变分得到有限元求解方程 （2.3.12）将（2.3.12）式代入（2.4.4） （2.4.5）在平衡情况下，系统总位能等于负的应变能。因此，当，则。在有限元解中，由于假定的近似位移模式一般来说总是与精确解有差别，因此得到的系统总位能总会比真正的位能大。我们将有限元解的总位能、应变能、刚度矩阵和结点位移分别用表示，相应的精确解的有关量用表示。由于，则有，即 （2.4.6）对于精确解有对于近似解有 （2.4.7）将（2.4.7）式代入（2.4.6）式得到 （2.4.8）由（2.4.8）式看出，近似解应变能小于精确解应变能的原因是近似解的位移总体上要小于精确解的位移。故位移元得到的位移解总体上不大于精确解，即解具有下限性质。2.5 轴对称问题的有限元格式作为弹性力学问题有限元分析一般格式的直接应用，本节将讨论工程中经常遇到的一类实际结构问题，即它们的几何形状、约束条件以及作用的载荷都对称于某一个固定轴，我们称之为对称轴，则在载荷作用下产生的位移、应变和应力也对称于此轴。这类问题称为轴对称问题。在轴对称问题中，通常采用圆柱坐标系。以对称轴为z轴，所有的应力、应变和位移都与无关，只是r和z的函数。任一点的位移只有两个方向的分量，即沿r方向的径向位移u和沿z方向的轴向位移w。由于轴对称，方向的位移等于零。因此轴对称问题可以看成二维问题。 离散轴对称体时，采用的单元是一些圆环。这些圆环单元与rz平面正交的截面可以有不同的形状，例如3结点三角、6结点三角形、4结点四边形、8结点四边形等。单元的结点是圆周状的铰链，并且各单元在rz平面内形成网格。图2.12所示为三角形环形单元。 对轴对称问题进行计算时，只需要取出一个截面进行网格划分和分析，但注意到单元是圆环状的，所有的结点载荷都应理解为作用在单元结点所在的圆周上。2.5.1 3结点三角形环状单元的插值函数及应力应变矩阵1位移模式和插值函数取出环状单元的一个截面如图2.13所示。单元结点位移为选择位移模式为待定常数，将结点坐标代入，可以求出。将求得的广义坐标代入位移模式中，则位移函数可表示为结点位移的函数，即 （2.5.2）其中 （2.5.3）称为单元的插值函数或形函数。其中，及是坐标的函数，具体计算为 （2.5.4）上式下标轮换，可得及。 （2.5.5）称为插值函数矩阵或形函数矩阵。2单元应变和应力 （2.5.6）B称为应变矩阵，应变矩阵的分块矩阵是 （2.5.7） （2