有限元绝密讲义2

1、第第2章章 平面问题的有限元方法平面问题的有限元方法2.1 弹性理论基础弹性理论基础 、基本假设：、基本假设： 连续性物质连续。相应的应力应变，位移等连续变量可以用坐标的连续函数表示； 均质各向同性物体内部各点，各方向上物理性质相同，材料常数（弹性模量，泊松比）不随坐标方向而变； 完全弹性材料服从Hooke定律； 小变形（几何假设）略去二阶小量，所有微分方程为线性的； 无初应力加载前物体内无初应力。 、基本方程、基本方程A、平衡方程：、平衡方程： (2.1)对于平面问题： (2.2)B、几何方程：、几何方程： (2.3) 00YyxXyxyyxxyx0,ijijXijjiijiiiixuxux

2、u对于平面问题： (2.4)(2.5) xvyuyvxuxyyx vuxyyxzyx00C：物理方程（应力应变关系）： (2.6)为弹性矩阵由弹性模量和泊松比构成，为常数矩阵。平面问题一般分为两类：平面应力状态 仅在平面内存在应力分量，即注意：(2.7) D D0zxzyz0z 2100010112ED：平面应变状态 仅在平面内存在应变分量，即 注意：(2.8) xyyx xyyx0z 1221000110112111ED0zxyzz实际上，只需对（2.7）式中的 、 作如下变换，即能得到（2.8）式。D、边界条件：力的边界条件（ ）：(2.9) 为边界外法线的方向余弦 位移边界条件（ ）：

3、E21EE1SuSjniiuu ijijTn 、能量原理、能量原理应变能：虚功原理：变形体中满足平衡的力系在任意满足协调条件的变形状态上作的虚功等于零。体系外力的虚功与内力的虚功相等。虚位移原理：如果弹性体处于平衡状态，当发生约束所允许的任意微小虚位移时，外力在虚位移上所作的功等于弹性体内的应力在相应虚应变上所作的功。 dU D21ddu设弹性体在体积力 ，表面力 作用下处于平衡，体内引起相应的位移 ，应力 和应变 。假定弹性体内有某约束条件所允许的虚位移 ，相应有虚应变 ，则虚位移原理可表述为：(2.10) 平面问题：设板厚t为常量，即 则，（2.11）以后证明：上式为平衡方程和力边界条件的

4、等效积分提法 f q X f ttdz0 tdxdytdsqftdxdyXfTTT ddAqfdXfTTT虚应力原理：几何方程和位移边界条件的等效积分提法应力、约束反力为虚拟量，而位移为确定量。注意：虚位移或虚应力原理中，均未涉及物理方程。也就是说，线弹性，非线弹性，弹塑性问题均适用；但是，平衡方程、几何方程都是在小变形条件下建立的，因此，二原理不能直接应用于大变形问题。更多的原理将在后续章节中阐述。2.2 分析流程分析流程对于任一平面问题，分析步骤的第一步是对结构离散化。考虑图示一平板（平面应力或平面应变不限）弹性理论求解时，取无限小微元体作为出发点，建立平板的力平衡、几何变形的微分方程形式

5、，为一无限自由度问题。而具体求解时可采用半逆解法，引入一多项式应力函数，定其系数，解出应力，积分得到位移，最后结果为一全域连续的分布函数。当采用有限元方法求解时，第一步是将平板离散成有限个小单元。A：梁结构的离散：取一段梁为一单元 单元类型：简单直线段 离散原则：几何上真实模拟原结构及其变形平板的离散：取一小面积板为一单元单元类型：由最基本的平面图形构成三角形、四边形（如正方形、长方形、梯形）而五边形、圆、扇形不宜作为单元。离散原则：几何上真实模拟原结构（无缺陷、重叠） 模拟变形状态B：节点的选择梁单元：取端点为节点，单元内任一点的几何位置可由节点坐标确定，单元内任一点的变形同样能由节点的变形

6、确定。对二节点梁单元，单元内任一点的坐标、变形可由节点的坐标、变形的线性函数插值确定。平面单元：取角点为节点，单元内任一点的几何位置、变形可由节点的坐标几变形插值确定。三角形单元节点：三个角点？矩形单元的节点：四个角点？C：节点位移梁单元中节点的位移分量对应直梁、平面梁、空间梁分别有1位移+1转动、2位移+1转动、3位移+3转动分量；平面单元的节点只有2位移：为什么？相应地，节点载荷： D：离散结果图示平板可分离成18个三角形单元、16个节点（32个位移分量）。无限自由度转换成32个自由度！ iiivu iiiYXQ 计算结果：仅16个节点的位移数值；可以将节点的位移拟合成一条曲线；当单元趋小

7、，节点增多，拟合曲线才逐渐逼近真解。 与弹性理论解的形式差别：不直观； 结果比较：不能直接反映各物理量，几何 量的影响。问题：平板的单元离散和节点位移描述都是在加载变形前实施的，各单元之间仅靠节点相连，那么加载后，能否保证单元的边界不开裂、不重叠？2.3 三角形单元三角形单元考虑图示任一三角形单元e，取角点为节点，节点编号为L、m、n(右手旋向!)，每一节点仅二个独立位移分量。单元的6个节点位移分量可按节点排序为： TnnmmlleVUVUVU梁单元：参照材料力学结果，设定位移在单元内的变化规律，再以节点位移取代待定系数，形成了由节点位移描述的位移插值函数。平面单元： 单元内任意点的变形也是坐

8、标的函数。当单元较小时，同样用节点位移的插值函数来描述单元内的位移。 一般有Lagrange、Hermite样条插值（有精度、适用范围要求）。考虑到简化微分、积分运算，多取插值函数为幂函数形式。 、选取位移模式、选取位移模式 三角形平面单元内的变形在x、y方向都是均等的，位移函数应具有相同的变化规律，即多项式函数具有相同的项，但系数可不相同！ 单元3个节点，仅6个节点位移分量，同时描述两个位移函数 ，只能假定u和v均为线性多项式，即（2.12）单元内的变形可由节点位移的插值-线性函数描述。 这里，选取最低阶的线性多项式形式，是否可以取其它形式，只要保证系数由节点位移分量取代就可以了？ ),()

9、,(yxVyxUyaxaayxVyaxaayxU654321),(),(将（2.12）改写为矩阵形式，（2.13） 简记， （2.14） 65432110000001),(),(aaaaaayxyxyxvyxu ayxHvu),( aHvu、单元节点位移、单元节点位移 与与 之关系之关系简记为（2.15） e ea654321100000011000000110000001aaaaaayxyxyxyxyxyxvuvuvunnnnmmmmllllnnmmll aAe考察 ，只要三角形的三个角点不在一条直线上，任一条边长度不显著小于其它二条边的长度，则 非奇异矩阵，其逆矩阵存在。 代入(2.14)

10、， (2.16)其中 (2.17) A A eAa1 eAHvu1 eNvu 1AHN(2.16)式即表示了单元内任意位置的位移由单元节点位移的插值函数来确定。 称为形函数矩阵，它代表和反映了位移函数与节点位移之关系，代表了单元自身的几何特性。形函数的合理性将直接影响到单元特性、求解精度、收敛性等。 N、形函数的特性、形函数的特性由（2.17）可以得到： (2.18) 式中： (对其它项， 轮换) 为三角形单元的面积引入： (2.19) nmlnmlNNNNNNN0000002ycxbaNllllnml, nnmmllyxyxyx111则：仅当右手系排序 ，能保证三角形面积恒为正值！中的分别对

11、应 第一行各元素对应的代数余子式， 分别对应第二行各元素对应的代数余子式； 第三行nml,lNlllcba, nnmmlyxyxa nmlyyb11nmlxxc11mmmcba,nnncba,21特性分析：特性分析：、将(2.18)代入(2.16)，对于单元内任一点的位移u或v，都可以同样坐标、节点位移作线性插值，形状函数均一致。、 时，代表了单元内部的位移分布函数；或者：节点 的单位位移对单元位移的贡献。在节点 ， ；在节点 ，形函数可分解为各节点的形函数，并相互独立。 iinnmmlluNuNuNuNyxU),(iinnmmllvNvNvNvNyxV),(0, 1nmluuulNyxU,l

12、l1lN0lNnm,、三个形函数之和为1，只有2个独立。（2.20）、在单元任一边上的形函数与第3个节点的坐标无关。 1nmlNNN、连续性要求连续性要求 任一结构发生变形时，结构内部的位移是一连续函数；当以有限单元离散后，位移模式应保证单元内、公共边界上的连续性。（物理现象：无缺陷、不重叠）形函数在单元内连续，则单元内的位移也连续？！相邻单元交界线上，两单元应有相同的位移？两相连的三角形单元公共边界为直线段； 变形前，边界由二公共节点唯一确定； 变形， 二节点的位移确定了边界的位移； 变形后，边界仍由二公共节点唯一确定！（设位移函数时，只取线性多项式的缘故之一） 、单元的应变应力、单元的应变

13、应力将位移函数（2.16）式直接代入平面问题的几何方程（2.5）式，可得到应变： 简写为：(2.21) ezyxNxyyxvuxyyx0000 eB式中，（2.22）记为，则,应变矩阵 的物理意义：当单元某一节点有单位位移而其它节点位移皆为零时，所引起单元内部的应变分布，亦是该点单位位移对单元应变的贡献。 nnmmllnmlnmlbcbcbccccbbbB00000021 nmlBBBB nnmmllBBB B由于 为线性函数，经微分运算后， 实际上为一常数矩阵，因此， 为常值；意指单元内部应变均匀；本三角形单元即为一种常应变单元。将(2.21)式代入(2.6)式，得到应力表达式，(2.23a

14、)(2.23b)为应力矩阵。 也可分块写成： N B eBD es s nnmmllsss的物理意义：单元内某一节点有单位位移，而其余节点位移皆为零时，所对应单元内部的应力分布；亦为该节点单位位移对单元内应力的贡献。对于平面应力问题： s iiiiiiibccbcbEs2121122),(nmli 讨论：讨论： 本三角形单元为常应变单元（亦即常应力单元）。对应力变化梯度大的区域，采用这种单元模拟时，单元之间的应力应变是跳跃的；仅当单元足够小，才能缓解单元间应力的突变，而趋于真解。 已讨论证明了位移函数在单元之间是连续分布的，但其应变应力在单元间有突变，即：只保证了函数连续，而未保证其导数连续。

15、物体内的实际应力场应该连续分布的。 这种单元简单、可靠，可以通过控制单元尺寸、疏密程度的方式，再辅以对单元节点应力作合理修正，仍能取得满意的近似结果。、单元刚度矩阵、单元刚度矩阵建立力-位移之关系？如何将单元内的变形、应力、应变与外部荷载（节点力、面力、体力等）相联系？引入虚功原理虚功原理（外力等于内力功，能量守恒）将之相关联。前述已采用节点位移插值函数描述单元内的位移、应变、应力。因此适合采用虚位移原理建立力与位移之关系单元节点位移：单元节点虚位移：（满足约束条件）单元内虚位移：(2.24) Tnnmmllevuvuvu Te eeNvuf单元内的虚应变为 ： (2.25) 单元节点力： 节

16、点力所作虚功：(2.26) 单元内应力作虚功：(2.27) eB TynxnymxmylxlePPPPPPP yllxlleTePvPuPW tdAUT eAeTTetdABDB eAeTTetdABDBU由虚位移原理： 考虑到 的任意性，与之相乘积的两部分应相等。令（2.28）则，（2.29） （2.29）即为表征节点力与节点与节点位移关系的有限元方程，单元刚度方程， 为单元刚度矩阵 。 UW Te eAeTetdABDBP ek eATetdABDBk eeekP注意：弹性矩阵 为对称矩阵，(2.28)式 中积分结果仍为对称矩阵，单元刚度矩阵为对称矩阵。在有限元方程 (2.29)式推导中，

17、没有引入或强调单元的几何性质，因而具有普遍意义，适合于所有平面单元， 表达式具有通用性。可以按节点分块描述 ek D ek ek按节点分块，(2.30) tdABDBkTAee dAtBBBDBBBeAnmlTTnTmTlnnnmnlmnmmmllmllkkkkkkkkkln tdABDBksTArrse),(nmlsr刚度元素的物理意义：当单元的 节点方向发生单位位移，而其它节点位移都保持为零的变形形态，在所有节点上需要施加的节点力。 表示 列元素中的第 个元素对应的 节点力分量。对角线元素为各节点对单元刚度的独立贡献，其余元素则为任意二节点之间的偶合效果。 ijkiijj、等效节点力、等效

18、节点力有限元方程中：节点力与节点位移之间的关系而实际外部载荷：体积力、面载、集中载三大类 外载也应转换到节点上。单元之间内力，单元边界上的分布内力也应同时转换到节点上，但公共节点上的节点力将相互抵消，毋须计算单元内力引起的节点力部分。依据能量守恒：当外载荷作的功与一组节点力作的功相等时，这组节点力就与外载荷等效。由虚功定理，所有外载荷在节点上所作虚功应与等效节点力在对应节点上所作的虚功相等。即式中， 为等效节点力。外载荷包括体积力 ，面力 ，集中力 。 eTeeTeQP)()( eQ X q F由等效节点力 ，（2.31）对于平面问题：（2.32） dXNdAqNFNQTTTe eQ tdAX

19、NdStqNFNQTTTe dXfdAqfFfQTTTeTe dXNdAqNFNTTTTe eNf逐项讨论： A、集中载集中载 作用于C点，为集中载作用点处的形函数值与集中载之积。 若集中载位于节点处，对应的形函数值为1，其等效节点力就是该点的集中载。如 作用于节点 ，则等效节点力为， F FNQTceF Fl0, 1lllNNNTyxeFFFQ0000B：表面分布力：表面分布力 注意：面积力包括了作用在单元边界上的内力，但在合成过程中内力相互抵消，因而只需考虑在弹性体边界上的表面外载所引起的等效节点力。 dAqNQTeq tdsqNtdsqNtdsqNQQQQnmlenemeleq例：考虑均

20、布侧压作用在 边界上，令 边长为 ，与x轴夹角为 ，侧压q在x，y方向的分量：nl nl llnxyylqqqsinnlyxxlqqqcos作用在单元上边界上的面积力：在单元边界上可建立局部坐标系s，根据形函数性质，沿边界的形函数可改写为：则， nllnyxxxyylqqqqlsNllsNn10mNlnxLxLllxyytqtdsqlstdsqNP2nlyLllyxxtqtdsqNP2 TnllnnllnexxyyxxyytqQ002若 为x向均布力，即有，若 为x向的三角形分布，即有， q 0qq TeqltQ01000121 q 01qlsq TeqltQ0320310021C：体积力：体

21、积力体积力 引起的等效节点力 采用（2.31）或（2.32）的益处在于，对于任意单元，计算规整化，容易实现计算机编程。 tdAXNtdAXNtdAXNQQQQnmlenemelex X、例题：、例题：考虑一直角三角形单元，节点坐标)0 , 0(), 0(),0 ,(ndmdl形函数，形函数， 三角形面积，于是，同样可得， ddyybnml011110000dyxyxannmml0010111nmlxxc221ddxycxbaNllll2dyNmdyxdNn检验形函数检验形函数 在节点 ，在节点 ，在节点 ，对任意点： 单元的形函数矩阵， (a) 应变矩阵应变矩阵， (b) )0 ,(dl0,

22、1nmlNNN), 0(dm0, 1nlmNNN)0 , 0(n0, 1mlnNNN1nmlNNN yxdyxyxdyxdN0000001 1101101010000100011dB平面应力状态下应力矩阵应力矩阵，(c) 按（2.30）式计算单元刚度矩阵刚度矩阵中的分块子矩阵，常应变状态下，与积分过程无关，直接有 (d)平面应力状态下的弹性矩阵，（2.7） 21210212101100100112dEs tBDBksTrrs 2100010112ED最后得到，（e） 对于任意等腰直角三角形的单元刚度矩阵，均与边长d无关，其刚度矩阵完全相同。 推广，对任意相似的平面单元，当材料物理性质相同时，其

23、单元刚度矩阵均相似？ 2321121212321211100212102101122Etke讨论：形状函数、应变矩阵、应力矩阵的物理意义和具体特征表述。取 ，而其余位移分量为0 ，（f） 由 ，将 (a)、(f)代入（g）由 ，将 (b)、(f)代入1lu Te000001 eNvu0dxvu eB（h）由 ，将 (c)、(f)代入（i）由 ，将 (e)、(f)代入（j）令 ，在图中，给出了节点力oodxyyx1 es0112dExyyx eeekP TeEtP1001122212EtP eP注意：显然，节点力互相平衡，单元处于平衡状态。事实上，对于任意的单元节点位移，由单元的刚度方程导出的节

24、点力都是平衡的。（方程本身就是建立在力平衡基础上！）有限元的观点：单元是在（j）式表示的节点力（或集中力）作用下，单元发生了（g）式的节点位移；相应地在单元内部引起了均匀拉伸（h）式和均匀拉应力（i）式。分析单元的边界应力边界应力：由边界条件：（k）在边界 上, 结果有： 对i 式，令则， 0ijijTnnl , 0 xn, 1yn, 0 xyyyyxxyxyxyxxqnnqnnyyq0 xqdEq21qqy在 边界上，由（k）式， 在 边界上， 由（k）式， nm1xn0yn0 xyqqxx0yqml 21xn21ynqqx21qqy2进一步合成为沿边界法线的 和沿边界切向的 ：法向应力：

25、切向应力： 上述说明：在这三组边界分布力作用下，单元内部将有均匀的应力场；反之亦然。nqtq45sin45cosyxnqqqqq221q12145cos45sinyxtqqqq121讨论： 可以验证：所有边界分布力与前述的节点力是静力平衡的，因而可以将这些节点力作为边界分布力的等效力。 从实际应力的连续分布规律看，单元之间本是以边界相连，正是边界上的分布力才能形成这种连续性。 不存在有集中的节点力，有限元中的节点力是虚构的，为单元间分布力的一种等效形式；以节点力等效代替了边界分布力。2.4 矩形单元矩形单元 矩形单元在平面问题分析中，是一种常用单元。设矩形单元 ，边长分别为2a，2b，且矩形的

26、两边分别与x，y轴平行。取4个角点为节点。 nmlk每一节点有2个位移分量 ，则单元节点位移列阵为：采用与三角形单元相似的方法来推导和分析单元的力学特性。 、局部坐标系、局部坐标系考虑到单元在平面中的位置不确定矩形尺度的变化位移、应变、应力等描述的简化积分的规则化vu, Tnnmmllkkevuvuvuvu引入自然坐标系作为矩形单元的局部坐标系，位于矩形单元中心 二坐标系之转换关系，(2.33) 或，(2.34),00, yxaxx0byy0axx0byy0单元内各点的自然坐标 ：单元的四条边分别为： 单元的四角点坐标，111111) 1, 1(k) 1, 1 ( l) 1 , 1 (m) 1

27、 , 1(n、位移模式、位移模式位移模式选取的初步考虑：位移模式在两坐标方向均应有相同变化规律； 单元变形时，其边界须由二节点的变形唯一确定。 (2.35)分析，沿二坐标方向变化规律相同；在任意一条边界上，u或v均为一线性变化，如边界 上， 87654321aaaavaaaauml 1 87654321aaaavaaaau变形前：边界 由两节点坐标确定；变 形：边界的变形是线性函数，由两节点变形确定；变形后：变形前坐标与变形量叠加，仍由两节点确定。所以，单元边界由两节点唯一确定。同时也保证了相邻两单元在公共边界上的连续性。将4个节点坐标分别代入位移模式于是， (2.36) ml aAe eAa

28、1 eAHvu1 eNvu同样可表述为分块形式：其中，(2.37)展开上式代入相应节点的自然坐标， nmlknmlkNNNNNNNNN0000000041100iNi0i0nmlki,1141kN1141lN1141mN1141nN检验形函数：满足形函数的基本特性。单元内位移函数：（2.38） jijiNji011iN eNvuf、单元应变、单元应变局部系与总体系存在坐标平动变换： 从（2.34）式， 微分算式： (2.39) 单元应变：(2.40) axx1ax1by1baababxyyx00100 eB写成分块形式，单元应力： (2.41) nmlkBBBBB 00001110014100

29、1iiiiiiiiibaababNbNaNaNbabBnmlki, D eBD eS、单元刚度、单元刚度 根据虚位移原理，可以建立单元的刚度矩阵方程。 (2.42) (2.43)按节点分块,(2.44) eeeQk tdABDBkTe nnnmnlnkmnmmmlmklmlllkknkmklkkekkkkkkkkkkkkkkkkkln tdABDBksTrrs讨论一： 矩形单元的位移为双线性模式，应变矩阵不再为常量，单元内的应变（及应力）按线性规律变化。与常应变三角形单元比较，能适应于变化梯度较大区域内的应变模拟，具有更广泛的适用性。 矩形单元的明显不足是，不能很好地模拟曲线边界，这包括与坐标

30、轴不平行的直线边界，因而应用受到局限性。 一般可将三角形单元与矩形单元混合使用。 位移模式中仅含高阶项，而无常数项、一次项，可行性？讨论二： 当选择采用某种形状的单元后，位移模式的确定自然是关键。位移函数、应变矩阵、应力矩阵、刚度矩阵、以及等效节点载荷的计算等均依赖与单元的位移模式。 如果一个单元的位移模式与真实位移分布相差很大，要得到良好的数值解是不可思议的。已经证明（在变分原理一章中详细阐述），对于一个给定的位移模式，其刚度系数的数值比精确值要大，对应的计算位移值则比真实的要小，即位移解一般偏小。 随单元网络的细分，由位移模式对应的数值解将由偏小值趋于精确解，得到真实解的下界。2.5 单元

31、的收敛准则与精度单元的收敛准则与精度 当采用一个完全的多项式作为单元内部的位移插值函数，在一个有限尺寸的单元内，可以与真实解吻合；而实际上，我们只能取有限几项的多项式，因此只能得到近似解。 需要了解的是：在什么条件下，单元尺寸趋于很小乃至零时，数值解能趋于真实解？ 如何选择单元的位移模式？其依据、判定准则是什么？ 、收敛准则：、收敛准则：A、完备性要求：、完备性要求：位移模式须包含单元的刚体位移、常应变 。位移模式除具有描述单元本身变形能力，还须具有描述其它单元形变后引发本单元刚体运动的能力；（1-4单元内有变形；5-6单元内无变形，仅有随单元4的偏转。） 当单元尺寸无限缩小时，即趋于弹性体中

32、的一点时，单元内的应变亦应趋于一常值；若单元内无常应变项，单元就不可能收敛于真实解。上述分析表明：位移模式中必须同时存在常数项、一次（线性）项。这也是选取三角形、矩形单元位移模式的主要因素之一。 B、协调性：、协调性： 单元在公共边界上位移或其导数具有连续性。 一般地，选取多项式作位移模式，在单元内部位移及其导数的连续性总是满足的；而交界面上要求位移、导数连续过渡，使得单元间在变形后不出现裂缝，重叠现象。交界界面上的位移变化应该由界面上的公共节点唯一确定。 位移在交界面不连续，其导数将不连续，引起应变奇异。平面，空间的一般问题，只需边界上位移连续；板、壳的应变由位移的二阶导数描述，则要求在交界

33、面上位移的一阶导数须连续；对于二、三维弹性体，单元交界面上的位移取决于节点的位移，即可保证单元的协调性。C、几何各向同性：、几何各向同性：位移函数从一个坐标系转换到另一个坐标系时，要保证其完备性和协调性不变。 或，单元内部相对任一方向的变形规律须一致。要求位移模式：具有完全的n次多项式；非完全多项式中，高次项具有对称性。 对于二维问题，参照pascal帕斯卡三角形取对称即可。注意： 任何一个单元的位移模式同时满足了上述三原则后，才能保证求解过程的收敛性，得到高精度的数值解。 在具体选取多项式的位移模式时，多项式中项数应等于单元节点的自由度总数，位移模式才可能唯一转换、与单元节点自由度相关联。

34、结构在某些环境中，如位移真实解为n次函数，而单元位移模式也包含了完全的n次多项式，经有限数目的单元离散后，有限元解即可得到真实解。 、精度、精度任一精确位移解可以在其邻域内作Taylor级数展开，如 在任一点 的展开式可写作： 若有限元解采用p阶完全多项式，单元尺寸为h， 。则位移模式能局部拟合精确解的前几项直到p阶。这样位移模式的误差为，相应应变、应力的误差为，收敛速度为 量级 iUyyuxxuuyxuiii,1ph),(hyxph1ph对于平面三角形三节点单元：位移模式为线性函数，即p=1,则位移U的误差为 应力应变误差为 如果单元尺寸减半，则误差将减小 相应地，应变、应力误差减小1/2收

35、敛速度亦为 量级 上述精度分析只是对单元本身而言，实际上不能对有限元解作出具体的误差估计！误差估计由综合性因素所决定。 2h 1h412122hA、类比法： 选择一存在解析解的相近类型问题，求解域尽可能与实际分析问题相近，采用形同形式的单元，网络划分密度相近。 则数值解的误差可作为实际问题求解的误差。B、推算法：当单元的位移模式满足收敛性准则， 时，有限元解应该单调收敛。 第一次网络离散的解为 ，第二次则将各单元尺寸减半后作网络离散，得到解 。而收敛速度为 ，精确解 可由下式预测， 0h1u2u rhu对于三角形单元的线性位移模式，则，2r 422221hhuuuu12431uuu rrhhu

36、uuu221注意： 这是一种离散误差，由有限单元离散结构，位移模式近似真实分布函数引起的误差。 还应注意到计算误差：计算机有限的有效位数所引起的四舍五入误差和截断误差。 如果结构刚度阵在不同方向相差过于悬殊，代数方程组有可能呈病态，导致求解误差很大，甚至可能无解。2.6 高精度三角形单元高精度三角形单元 前述的三角形三节点单元为：线性位移、常应变；与矩形单元比较：灵活。扩充三角形单元，转换为非常应变单元方式：增加位移函数的阶次，增加节点数 3节点-6节点取三个角点、三条边的中点为单元节点。单元节点位移：位移模式仍取多项式，参照帕斯卡三角形，可以取为完整的二次多项式：（2.45）待定系数代入节点

37、的坐标、位移，得到：令（2.46） ayxyxyxyxyxyxvu222210000000000001 Taaaaaaaaaaaaa121110987654321 aAe 1AHN eNvuf Tevuvu6611 式中， 及 的推导运算十分冗繁，尽管可以采用编程方式求解，但不利于应变、应力、单元刚度的表达； 与三节点三角形单元不同，不存在有 ，难以简化运算，尤其是基本矩阵的推导。 目标是建立一个通式，对任何形状的三角形单元都适用，进一步拓展三角形单元的应用范围。 参照矩形单元的分析，通过引入局部坐标、坐标变换，建立了不同尺度矩形单元的通式。 引入一种新的局部坐标系面积坐标。 1,AA N 、

38、面积坐标：、面积坐标：将三角形中任一P点与三个角点相连，形成三个子三角形。以所含大三角形边所对的角点定义为子三角形的面积。如 的面积为 ，同理有 ， 。而大三角形的面积为pmneAmAnAA定义面积坐标：且，则（2.47）P点的位置表示为： 三角形三个角点的面积坐标： 三条边的方程： AALllAALmmAALnnAAAAnmlnmlLLLP,) 1 , 0 , 0()0 , 1 , 0()0 , 0 , 1 (nmlml 0nLnm0lLln 0mL1nmlLLL注意： 三角形内与节点所对边 平行的直线上的所有点有相同的一个坐标； 三个面积坐标只有二个是独立的； 三角形的面积坐标与三角形的具

39、体形状、在总体系中的位置均无关； 面积坐标即为一种自然坐标。ml 、面积坐标与直角坐标之转换关系面积坐标与直角坐标之转换关系图示三角形角点的直角坐标： 任一点P的直角坐标： 将三角形面积 以直角坐标表示成： （2.48）llyxl,mmyxm,nnyxn,yxp,nmlAAA,nnmmlyxyxyxA11121 yxxxyyxyyxmnnmnmnm21面积坐标 ：（2.49）将(2.48)、(2.49)与(2.18)比较，面积坐标与形函数完全相同!可以将面积坐标与直角坐标关系表示成矩阵形式,（2.50）ycxbaAAALlllll21lLnmll,yxcbacbacbaALLLnnnmmmll

40、lnml121将 分别乘以后，再相加结合得到直角坐标与面积坐标关系，(2.51)nmlLLL,nmlxxx,nnmmllLxLxLxxnnmmllLyLyLyy1nmlLLLnmlnmlnmlLLLyyyxxxyx1111、面积坐标的微积分运算、面积坐标的微积分运算由(2.50)，一般复合函数求导法则： (2.52)面积坐标在三角形全面积上的积分公式， (2.53)在任一条边上如的积分公式 ，(2.54)lellelbAxNxLCAyNyL2121xLLxLLxLLxnnmmllnnmmllLbLbLbA21nnmmllLcLcLcAy21AcbacbadxdyLLLAcnbmal2!2!lb

41、abadsLLbmlal!1!三角形单元：32!2001! 0 ! 0 ! 1AAdxdyLAl62!2002! 0 ! 0 ! 22AAdxdyLAl122!2011! 0 ! 1 ! 1AAdxdyLLAml、单元形函数由面积坐标与形函数比较：直接用面积坐标表示形（插值）函数，而不必先求解待定系数与节点位移之关系，不必求解节点坐标构成的逆矩阵。对三节点单元， (2.55) 线性单元的三个形函数就是三角形单元的三个面积坐标 对六节点单元：llLN nml,节点1、2、3坐标，由单元任意点p移动到1-2边中点，既构成节点4。 在1-2边上，节点4坐标，同理有节点5、6坐标，0011，0102，

42、1003，423143212121ALL03L0 ,21,21421,21, 0521, 0 ,216各直线段方程：1-4-2：2-5-3：3-6-1：46：45：56：按节点分块形式建立形函数，对于每个节点，可选择几条直线，使得直线通过除节点本身以外其余全部节点；再利用直线方程的左部为函数因式来构造形函数。03L01L02L0122/111LL0122/122LL0122/133LL如节点1：可选二条直线，46、253，即包揽了剩余全部节点。分别取二直线方程的左部为形函数的因式，得到节点1的形函数，(2.56a)显然，对于26节点，代入 中，均有 ，仅在节点1， 。 对节点4，可选二条直线，

43、361、253，由节点4的坐标和形函数值定待定系数C(2.56b)11112LLN1N01N11N2144LLN 214LCLN 同理，(2.56c)校核：通过形函数即可建立单元位移模式：(2.57)22212LLN33312LLN3254LLN 1364LLN 13322133221161444121212LLLLLLLLLLLLNi32123212LLLLLL1 eNf2.7计算流程计算流程以平面问题为例，详细介绍有限元求解的步骤。Step1： 结构离散化A、选择正确的物理模型B、选择合理的单元、节点、位移模式C、选择高精度单元D、同一个结构，单元类型尽可能少 E、载荷突变点应选为节点F、

44、单元划分数量 G、节点编号 Step2： 单元特性分析计算针对不同单元建立其刚度矩阵，并按节点分块；同时计算单元上外载荷的等效节点载荷 Step3： 形成总刚，总载荷列阵 按总节点编号，把分块后的单元刚度子矩阵叠加存入相应的总刚中对应区域；由总刚的对称性，一般只需存入下三角元素即可；将单元的等效节点载荷按总结点编号叠加形成载荷列阵。Step4： 边界条件 引入边界约束条件，消除物体的刚体位移，消除总体刚度矩阵的奇异性。 位移条件有两类：一是给定边界位移的约束条件，排除了刚体运动的可能性；另一类是给定边界一个有限位移值，类似于位移加载。 引入已知边界条件，未知量减少，有限元方程数目减少，重新组合

45、将导致数据交换复杂化。有二种处理方式，可以维持方程维数、排序不变：局部修正法和大数法A、局部修正法、局部修正法将节点位移的指定值置入方程组中，替换相应位移变量所在序号的方程，仍保持方程组维数不变，对相关的刚度矩阵和载荷列阵作修正。如节点 在y方向有给定位移 ：令 中的元素 ，且第 行和第 列的其余元素均为零值； 令 中的第 个元素的值等于给定位移值 ， 中的第 j 元素都需减去该节点位移值 和原来 中第 行相应的第j列元素的乘积； 例： 设 ，引入上式 iiV k12,2iiki 2i2 Qi 2iV QiVi 2 k4321432144434241343332312423222114131211QQQQuuuukkkkkkkkkkkkkkkk3u43423213143214442412422211412110010000kQkQkQuuuukkkkkkkkk这样 ， 而 均为原方程的解答。这种方式保持了总刚维数不变化，求解方程组数目不减少，避免有限元方程组的重新组合、节点重新排序等烦恼。B、大数法、大数法 如节点 在y方向有给定位移