FEMch有限元法求解平面问题改前课件

FEMch有限元法求解平面问题改前FEMch有限元法求解平面问题改前FEMch有限元法求解平面问题改前有限元法根本思路:离散化构造单元内位移函数；单元位移模式单元分析；划分网格，将连续体划分为有限数量的单元。单元内位移节点位移单元刚度矩阵单元节点力节点位移变分法思想整体分析；总体刚度矩阵节点位移外载荷静力平衡求解；节点位移单元位移模式几何方程物理方程差分法思想FEMch有限元法求解平面问题改前FEMch有限元法求解平面1有限元法根本思路:

离散化

构造单元内

位移函数；单元位移模式

单元分析；划分网格，将连续体划分为有限数量的单元。单元内位移节点位移单元刚度矩阵单元节点力节点位移变分法思想

整体分析；总体刚度矩阵节点位移外载荷静力平衡

求解；节点位移单元位移模式几何方程物理方程差分法思想有限元法根本思路:离散化构造单元内单元位移模式单第一节根本物理量和方程的

矩阵表示第一节根本物理量和方程的

矩阵表示3第一节根本物理量和方程的矩阵表示1根本物理量的矩阵表示外力：应力：应变：位移：虚应变：虚位移：第一节根本物理量和方程的矩阵表示1根本物理量的矩阵表示平面应变问题弹性矩阵：第一节根本物理量和方程的矩阵表示2根本方程几何方程：物理方程：符号矩阵：平面应力问题弹性矩阵：平面应变问第一节根本物理量和方程的矩阵表示2根本方程几第一节根本物理量和方程的矩阵表示位移变分方程虚功方程极小势能原理应力边界条件平衡微分方程等价!节点力:节点对单元的作用力。对单元来说，节点力是作用在单元上的外力。节点位移:节点产生的位移。当单元进入虚位移状态：节点虚位移：单元内应力：单元内位移：单元虚应变：第一节根本物理量和方程的矩阵表示位移变分方程应力边界条件第一节根本物理量和方程的矩阵表示对单元，外力虚功等于内力虚功。虚功方程矩阵表示外力虚功：节点力在节点虚位移上所做的功内力虚功：单元内应力在虚应变上所做的功第一节根本物理量和方程的矩阵表示对单元，外力虚功等于内力第二节构造离散化第二节构造离散化8第二节构造离散化

深梁（离散化结构）将连续体变换为离散化构造：将连续体划分为有限多个、有限大小的单元，并使这些单元仅在一些节点处连接，构成所谓“离散化构造〞。单元节点1单元要素节点：单元与单元之间的连接点〔i,j,m〕。节点位移：节点产生的位移。ij

m节点力：通过节点传递的内力。节点载荷：作用在节点上的载荷〔外力〕。单元位移：单元内位移分布〔u(x,y),v(x,y)〕第二节构造离散化深梁（离散化结构）将第二节构造离散化2单元类型一维单元：如杆单元，梁单元

二维单元：如三角形单元，四边形单元。

三维单元：如四面体单元，六面体单元，棱柱单元。

3连续体离散化模型

单元间仅通过节点连接，没有其它联系；位移，载荷仅通过节点传递；单元内依然是连续体，位移是坐标的连接函数。第二节构造离散化2单元类型一维单元：如杆单元，梁单元第三节单元位移模式解的收敛性第三节单元位移模式解的收敛性11第三节单元位移模式解的收敛性1单元位移函数

将连续体离散化后，各单元间通过节点相连，但每个单元仍然是连续的，即单元内位移是坐标的连续函数：物理意义表示了单元内位移

u,v的分布形式，称单元位移模式。数学意义构造各单元节点位移间的位移插值函数。2形函数设位移函数：变分法在单元上取位移试函数。第三节单元位移模式解的收敛性1单元位移函数第三节单元位移模式解的收敛性假设边界：求解三元一次方程组2）为的代数余子式。

其中：1〕A为三角形ijm的面积，(i,j,m按逆时针编号)第三节单元位移模式解的收敛性假设边界：求解三元一次第三节单元位移模式解的收敛性这里：规律：i,j,m

安逆时针替换。一般写成：第三节单元位移模式解的收敛性这里：规律：i,j,第三节单元位移模式解的收敛性用矩阵形式表示：这里：形函数形函数矩阵那么：几何意义：反映单元位移形态。对三角形单元，在图形上是一个平面。第三节单元位移模式解的收敛性用矩阵形式表示：这里：第三节单元位移模式解的收敛性

位移函数表示了位移u,v

在单元内的分布形式，因此又称位移模式。沟通了单元中离散点的位移和单元内位移之间的关系。3形函数的性质i点：j点：m点：第三节单元位移模式解的收敛性位移函数表示第三节单元位移模式解的收敛性性质1:在节点r处在其它节点处物理意义：不同单元在同一节点处位移相等，即节点处位移与形函数无关，反映了单元在节点处的连续性。总结性质2:总结第三节单元位移模式解的收敛性性质1:在节点r处在其第三节单元位移模式解的收敛性令：各节点位移相等时，单元内位移为常数，且等于节点位移。物理意义：各节点形函数之和为1，反映了单元的刚体位移。是x,y

的线性函数。jmijim由性质1：1m11第三节单元位移模式解的收敛性令：各节点位移相等时，第三节单元位移模式解的收敛性4位移函数的收敛性(1)位移函数收敛性的概念位移模式建立以后，便可逐步求应变、应力、结点力等一系列工作。所以位移模式是有限元法的工作根底。…逐步将单元细分，可以得到同一问题近似解的一个序列：准确解？第三节单元位移模式解的收敛性4位移函数的收敛性(第三节单元位移模式解的收敛性当时，有限元解答收敛于精确解有限元解的收敛性位移函数的收敛性位移函数收敛于正确解(2)位移函数收敛于正确解的条件必须反映单元的刚体位移。单元位移自身变形引起位移刚体位移=+由其它单元变形引起

必须反映单元的常量应变。单元应变常量应变变量应变=+与位置无关与位置坐标有关a,b

两个条件为位移函数收敛的必要条件。第三节单元位移模式解的收敛性当第三节单元位移模式解的收敛性

位移模式应尽可能反映位移的连续性。单元内连续〔位移函数取坐标的连续函数〕；单元公共节点上位移相等；单元公共边界上位移相等。保证离散化后的构造仍为连续弹性体。c条件为位移收敛性的充分条件。满足a,b条件的单元完备单元还满足c条件的单元协条单元第三节单元位移模式解的收敛性位移模式应尽可能反复习：复习：复习性质1:在节点r处在其它节点处性质2:位移函数收敛于正确解的条件:必须反映单元的刚体位移。

必须反映单元的常量应变。

位移模式应尽可能反映位移的连续性。jmijim复习性质1:在节点r处在其它节点处性质2:位移函数收敛于正确例题：例：判断三角形单元位移函数是否满足位移函数的收敛性？解：〔1〕令：因此，u,v

包含了单元的刚体位移。〔2〕因此，应变分量包含了常应变。例题：例：判断三角形单元位移函数例题：〔3〕单元内，u,v

显然是连续函数。对单元①：对单元②：所以，不同单元在公共节点处位移相等，保证了节点处的连续性。u,v

在任一单元内都是坐标的线性函数，在公共边界ij

上当然也是线性的。

经过两点的直线只有一条，所以公共边界ij

上任一点位移相等，保证位移函数在公共边界上是连续的。例题：〔3〕单元内，u,v显然是连续函数。对单元①：对例题：例题：例题：例题：第四节单元分析单元刚度矩阵第四节单元分析单元刚度矩阵28第四节单元分析单元刚度矩阵1.单元的应变和应力(1)单元的应变单元的位移模式反映了单元内任一点的应变与节点位移之间的关系第四节单元分析单元刚度矩阵1.单元的应变和应第四节单元分析单元刚度矩阵应变矩阵应变矩阵子矩阵〔i,j,m〕对三角形单元，均仅与节点坐标有关为常量矩阵

为常量三角形单元为常应变单元(2)单元的应力反映了单元内任一点的应力与节点位移之间的关系第四节单元分析单元刚度矩阵应变矩阵应变矩阵子矩第四节单元分析单元刚度矩阵应力矩阵应力矩阵子矩阵〔i,j,m〕说明：

为常量，单元应力为常量。三角形单元为常应力单元；

对三角形单元，不同单元之间应力、应变有突变。位移函数误差量级：应力应变误差数量级：即使位移精度满足了要求，应力精度不一定满足要求措施：a,细化网格；b,采用高次单元权衡计算精度和计算量之间的关系第四节单元分析单元刚度矩阵应力矩阵应力矩阵子矩第四节单元分析单元刚度矩阵2.单元刚度矩阵建立节点力与节点位移之间的关系。(1)单元刚度矩阵的推导

以单元为研究对象，假想将单元与节点i切开，节点作用于单元上的力，称为节点力。令节点处产生了一组虚位移：相应单元虚位移：相应单元虚应变：

由虚功方程：外力（节点力）在虚位移（节点虚位移）上的虚功，等于应力在虚应变上的虚功，即：第四节单元分析单元刚度矩阵2.单元刚度矩阵建第四节单元分析单元刚度矩阵（与坐标x,y

无关）对任意虚位移上式都必须满足，则：令：则：单元刚度矩阵

有限元实际运算中大量的计算工作量是形成单元刚度矩阵单元刚度方程第四节单元分析单元刚度矩阵（第四节单元分析单元刚度矩阵3.单元刚度矩阵的计算ijm

ijm每个分块子矩阵位置与节点号码（逆时针）顺序相对应平面应力单元刚度矩阵计算公式平面应变问题的刚度矩阵第四节单元分析单元刚度矩阵3.单元刚度矩阵的第四节单元分析单元刚度矩阵4.单元刚度矩阵的物理意义当节点j

处产生单位位移，而其余点完全被约束时，节点i

处所需施加的节点力。当节点j

处产生单位水平位移时，节点i

处所需施加的水平节点力。即时

时

时

时对子矩阵，如：第四节单元分析单元刚度矩阵4.单元刚度矩阵的第四节单元分析单元刚度矩阵4.单元刚度矩阵的性质

对称性：

奇异性：从物理意义证明：直接由单刚计算公式证明：令：则：单元内无应力应变产生以第1行为例：第四节单元分析单元刚度矩阵4.单元刚度矩阵的第四节单元分析单元刚度矩阵上式对任一都成立，则一定有：为奇异矩阵

三角形单元不变性：

与无关。第四节单元分析单元刚度矩阵上式对任一第四节单元分析单元刚度矩阵

与无关。

与无关。

代表单元刚体平移。单元刚度矩阵与单元的刚体平移无关另外可以证明：单元刚度矩阵与单元的刚体转动也无关。影响单元刚度矩阵的因素：单元形状单元所取位移函数材料属性：第四节单元分析单元刚度矩阵例题：例题：应用变分原理推导单元刚度方程(1)单元应变能由于上式中：代入应变能公式：仅与节点坐标有关：应用变分原理推导单元刚度方程(1)单元应变能由于上式中：应用变分原理推导单元刚度方程令：,则：(2)单元外力势能(3)单元总势能根据极小势能原理：应用变分原理推导单元刚度方程令：第五节载荷移置第五节载荷移置42第五节载荷移置1.静力等效原那么刚体静力等效原那么只从运动效应来考虑，得出移置荷载不是唯一的解；变形体的静力等效原那么考虑了变形效应，在一定的位移模式下，其结果是唯一的，且也满足了前者条件的。所以在FEM中，采用变形体的静力等效原那么。刚体静力等效原那么:变形体静力等效原那么:使原荷载与移置荷载的主矢量以及对同一点的主矩也一样。在任意的虚位移上，使原荷载与移置荷载在虚位移上所做的虚功相等。单元原载荷在虚位移上做的虚功=移置后节点载荷在相应虚位移上做的虚功第五节载荷移置1.静力等效原那么刚体第五节载荷移置2.移置方法(1)集中力的移置集中力按形函数在各节点分配第五节载荷移置2.移置方法(1)集中力的移置集中力按第五节载荷移置(2)体力的移置取微单元，那么微单元上的体力：可看作集中力。例1：重力移置。令密度为,则:第五节载荷移置(2)体力的移置取微单元，那么微单元上的第五节载荷移置(3)面力的移置可看作集中力。取微单元，面积为。则微单元上的力：则:例2：三角形单元一边受沿x方向均布载荷。q第五节载荷移置(3)面力的移置可看作集中力。取微单元，第五节载荷移置q例3：三角形单元一边受沿x方向线性分布载荷，在i

点为q，在j点为0，求移置后的节点载荷。令ij边长为as第五节载荷移置q例3：三角形单元一边受沿x方向线性分布第六节整体分析总体刚度矩阵第六节整体分析总体刚度矩阵48〔1〕节点协调原那么：节点处保证协调连接不同单元在同一节点处位移相等〔2〕节点平衡原那么：①②第六节整体分析总体刚度矩阵节点平衡方程〔1〕节点协调原那么：节点处保证协调连接2.总体刚度矩阵形成节点载荷：节点1：①①③③节点2：①①②②节点3：①①②②④③③④节点4：③③④④节点5：①③①③①②①②③①②④③①②④③④③④②④②④第六节整体分析总体刚度矩阵②②④④节点平衡原那么2.总体刚度矩阵形成节点载荷：节点1：①①③③节点2：①

将上面10个方程用矩阵表示：①③①③①②①②③①②④③①②④③④③④②④②④①①①①①①②②②②②②③③③③③③④④④④④④单元①单元②单元④单元③12345节点第六节整体分析总体刚度矩阵节点平衡方程组将上面10个方程用矩阵表示：①③①③①②①②③①②④③①②用分块矩阵表示：①①①②②②③③③④④④单元①单元②单元④单元③12345节点用零（矩阵）升阶使节点力矩阵与节点载荷矩阵同阶第六节整体分析总体刚度矩阵用分块矩阵表示：①①①②②②③③③④④④单元①单元②单元④单列出各个单元单元刚度方程：单元①①①①单元②②②②①①①①①①①①①②②②②②②②②②第六节整体分析总体刚度矩阵单元位移列阵用整体位移列阵代替单元刚度矩阵用零升阶与总体刚度矩阵同阶，又称单元奉献矩阵。①②列出各个单元单元刚度方程：单元①①①①单元②②②②①①①①①单元③③③③单元④④④④③③③③③③③③③④④④④④④④④④(1)

扩阶后的单元刚度方程只不过增加了两个0=0的方程，因此与原方程等价。第六节整体分析总体刚度矩阵〔3〕用构造的整体位移矩阵代替单元节点位移矩阵-----节点协调原那么(2)不影响虚功及应变能的计算。③④单元③③③③单元④④④④③③③③③③③③③④④④④④④④④④将以上四式代入节点平衡方程组：①②③④①①①②②②③③③④④④节点平衡方程组令：①②③④总体刚度矩阵整体平衡方程,沟通节点载荷与节点位移关系第六节整体分析总体刚度矩阵将以上四式代入节点平衡方程组：①②③④①①①②②②③③③④①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④①②③④总体刚度矩阵的形成：所有升阶后的单元刚度矩阵〔单元奉献矩阵〕的叠加！第六节整体分析总体刚度矩阵①①①①①①①①①②②②②②②②②②③③③③③③③③③④④④④④④④④④①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③3.总体刚度矩阵的计算规律

令：为总体刚度矩阵的任一子块:当时：为主对角线上子块，由环绕节点r或s的单元在单刚矩阵相应节点上子矩阵的叠加如：①+②+③+④

(2)当且rs为相邻单元的公共边时：为相邻两单元在单刚矩阵相应子块矩阵的叠加。如：①+③

(3)当且rs仅为某一单元的一边时：为该单元单刚矩阵相应子块矩阵。如：①

(4)当且rs为互不相关的两节点时：如：①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④第六节整体分析总体刚度矩阵①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④3.总体刚度矩阵的计算规律令：3.总体刚度矩阵的性质

对称性①+③①+③①③①③①+③①③①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④存贮一半元素加上对角线上的元素第六节整体分析总体刚度矩阵3.总体刚度矩阵的性质对称性①+③①+③①③①③①

奇异性①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④以第一行为例,考察第一行之和：①+③③①①+③第一行子块矩阵有：①③单元刚度矩阵第一行元素之和①单元刚度矩阵第一行元素之和③反映了弹性体的刚体位移①①①①①①①①①③③③③③③③③③第六节整体分析总体刚度矩阵奇异性①+③③①①+③①①+②①+②②①+③①+②①+

稀疏性当且r、s为互不相关的两节点时：那么单元越多，不相关的节点越多，等于零的元素越多，称总刚矩阵的稀疏性例：总刚为24×24阶矩阵，非0元素仅占37.5%稀疏矩阵一般一个节点的相关结点不会超过九个，如果网格中有200个节点，那么一行中非零子块的个数与该行的子块总数相比不大于9/200，即非零子块在5%以下，零子块在95%以上。高稀疏矩阵第六节整体分析总体刚度矩阵稀疏性当且r、s为互不相关的两

带状性00

如果编码合适，的非零元素分布在以对角线为中心的带状区域内。B半带宽B包括对角元素在内的半个带状区域内，每行〔或列〕包含的最多元素个数B=〔相邻节点号码的最大差值+1〕×2存贮包括主对角线在内的半个带状区域元素第六节整体分析总体刚度矩阵带状性00如果编码合适，的非零元4.四节点矩形单元三角形单元：对曲线边界适应性好，

计算精度低位移函数的项数的多少要受到单元型式的限制，项数应与单元的自由度数一致。i,j,m,p四个点的坐标值位移函数求得代回位移函数代入位移模式第六节整体分析总体刚度矩阵4.四节点矩形单元三角形单元：对曲线边界适应性好，式中形函数为：其中：第六节整体分析总体刚度矩阵式中形函数为：其中：第六节整体分析总体刚度矩阵矩形单元位移函数的收敛性讨论：1〕刚体位移〔同三角形单元〕2〕常量应变：第六节整体分析总体刚度矩阵矩形单元位移函数的收敛性讨论：1〕刚体位移〔同三角形单元〕23〕位移连续性：a.单元内连续b.节点连续在节点r处在其它节点处显然u,v是x,y的连续函数。第六节整体分析总体刚度矩阵3〕位移连续性：a.单元内连续b.节点连续在节点r处在其不同单元在同一节点处位移相等。c.边界连续以ij边为例：y=-b代入：所以，u,v在边界上为一直线，经过两点有且只有一条直线，边界连续。第六节整体分析总体刚度矩阵不同单元在同一节点处位移相等。c.边界连续以ij边为例：第七节边界条件处理

计算成果整理第七节边界条件处理

计算成果整671.进展边界处理的原因总刚矩阵的奇异性：无逆矩阵无法求解的奇异性是由于其中包含了弹性体的刚体位移考虑边界约束，消除刚体位移边界处理

给定位移边界2.边界的约束情况对称构造：第七节边界条件处理计算成果整理几何形状，边界约束及载荷均关于x轴或y轴对称。1.进展边界处理的原因总刚矩阵的奇异性：无逆矩阵无法求解2.进展边界处理的方法(1〕降阶法：将总体平衡方程中位移为0的相应行和列划去第七节边界条件处理计算成果整理1.降阶：2n×2n降为(2n-r)×(2n-r)2.奇异矩阵变为非奇异矩阵3.打乱[K][FL]的储存顺序,不利于编程2.进展边界处理的方法(1〕降阶法：将总体平衡方程中位移(2)对角元素置1法第七节边界条件处理计算成果整理将位移为0的所在行和列置0，对角元素置1如：那么：(2)对角元素置1法第七节边界条件处理计算成(3)对角元素乘大数法令：第七节边界条件处理计算成果整理如：那么：(3)对角元素乘大数法令：第七节边界条件处理3.计算成果整理

在有限元法中，位移的精度较高，其误差量级是，即与单元尺度的二次幂成正比。应力的误差量级是，即与单元的大小成正比。计算成果(1)位移(2)应力〔应变〕对于结点位移的成果，可以直接采用。

三结点三角形单元的应力的成果，不但应力的精度较低，而且还产生了所谓应力的波动性。假设单元①的应力成果为，其中为真解，为误差。由于在结点都列出了平衡方程并令其满足，从而使相邻的②单元的应力趋近于。应力的波动性第七节边界条件处理计算成果整理3.计算成果整理在有限元法中，位移的精度应力成果的整理方法：〔1〕绕结点平均法：把环绕某一结点的各单元中的常量应力加以平均，用来表征该结点处的应力；〔2〕两相邻单元平均法：把两个相邻单元中的常量应力加以平均，用来表征公共边中点处的应力。提高应力的精度方法：〔1〕加密网格，减少单元的尺寸，以提高应力的精度。〔2〕采用较多结点的单元，并使位移模式中包含一些高幂次的项，从而提高位移和应力的精度。第七节边界条件处理计算成果整理应力成果的整理方法：〔1〕绕结点平均法：把环绕某一结点的各单由具体构造及给定条件，绘计算简图---尺寸、外载荷、约束等；有限元法的求解步骤：选定坐标系，划分单元，单元与节点编号〔整体编号，局部编号〕；

计算单刚矩阵：单元面积A，,，形成单刚；

组集总刚；

载荷移置；

边界处理；

求解线性方程组。

计算成果整理。第七节边界条件处理计算成果整理由具体构造及给定条件，绘计算简图---尺寸、外载荷、约束等第八节有限元法求解实例第八节有限元法求解实例解:

(1)对称问题处理：几何形状，边界约束，和载荷均关于x轴和y轴对称。单元号总体局部节点单元①单元②ijmijm123341

如图为厚度为t的矩形薄板。两端受均布拉力1N/m，板长2m，宽1m,材料常数为E，，在不计自重的情况下，试用有限元法求板内应力分布。31=m第八节有限元法求解实例〔2〕划分单元，单元及节点编号i,j,m逆时针方向编码，保证单元面积为正解:(1)对称问题处理：几何形状，边界约束，单元号总体局(3)求单刚矩阵：以单元①为例：第八节有限元法求解实例(3)求单刚矩阵：以单元①为例：第八节有限元法求解实①②单元②可通过单元①的逆时针旋转而得到：第八节有限元法求解实例①②单元②可通过单元①的逆时针旋转而得到：第八节有〔4〕组集总刚矩阵ijmijm①②①①②对称②以子块为例：第八节有限元法求解实例单元号总体局部节点单元①单元②ijmijm123341〔4〕组集总刚矩阵ijmijm①②①①②对称②以子块(5)边界条件处理：

分别对应于节点位移矩阵的第1、2、4、7行，划去总刚矩阵的第1、2、4、7行与列，得到：对称(6)载荷移置：对称对称第八节有限元法求解实例(5)边界条件处理：分别对应于节点位移矩阵的第(7)求解线性方程组：①①①②②②第八节有限元法求解实例(7)求解线性方程组：①①①②②②第八节有限元法求解实第九节有限元法程序设计简介有限元法程序设计:划分单元输入参数:〔坐标，编码，材料常数，载荷，约束〕计算：形成单刚子程序形成总刚子程序形成载荷子程序边界处理子程序解线性方程组子程序计算结果处理〔整理，分析，显示〕前处理〔Preprocessor〕求解〔Solution〕后处理(Postprocessor)第九节有限元法程序设计简介有限元法程序设计:划分单元前第九节有限元法程序设计简介第九节有限元法程序设计简介输入数据信息形成总刚K形成载荷R计算应力输出应力求解总刚方程，输出位移开场完毕第九节有限元法程序设计简介程序流程图程序框图有限元法程序设计:输入数据信息形成总刚K形成载荷R计算应力求解总刚方程，输出，INPUT——输入结点坐标、单元信息和材料参数；MR——形成结点自由度序号矩阵；FORMMA——形成指标矩阵MA〔N〕并调用其他功能子程序，相当于主控程序；DIV——取出单元的3个结点号码和该单元的材料号并计算单元的等；等；MGK——形成整体刚度矩阵并存放在SK〔NH〕中；LOAD——形成整体结点荷载列阵F；OUTPUT——输出结点位移或结点荷载；TREAT——由于有非零位移，对K和F进展处理；DECOMP——整体刚度矩阵的

分解运算；FOBA——前代、回代求出未知

结点位移

；ERFAC——计算约束结点的

支座反力；KRS——计算单元劲度矩阵中

的子块Krs。第九节有限元法程序设计简介，INPUT——输入结点坐标、单元信息和材料参数；MR——形CFINITEELEMENTPROGRAMFORTWODIMENSIONALCTRIANGLEELEMENTCDIMENSIONK(800000),COOR(2,3000),AE(4,11),*MEL(5,2000),MA(6000)CHARACTER*32datCOMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NCWRITE(*,300)300FORMAT(///''*':****'/'+PLEASEINPUTFILENAMEOFDATA')READ(*,*)dataOPEN(4,FILE=data,STATUS='OLD')OPEN(7,FILE='OUT',STATUS='UNKNOWN')READ(4,*)NP,NE,NM,NR,NI,NL,NG,ND,NCCWRITE(*,400)NP,NE,NM,NR,NI,NL,NG,ND,NCCWRITE(7,400)NP,NE,NM,NR,NI,NL,NG,ND,NCCALLINPUT(JR,COOR,MEL,AE)CALLCBAND(MA,JR,MEL)CALLSK0(SK,R,COOR,MEL,MA,JR,AE) CALLLOAD(COOR,MEL,R,JR,AE)IF(ND.GT.0)CALLTREAT(SK,MA,JR,R)CALLDECOMP(SK,MA)CALLFOBA(SK,MA,R)WRITE(*,650)WRITE(7,650)CALLOUTPUT(JR,R)WRITE(*,700)WRITE(7,700)CALLCES(COOR,MEL,JR,R,AE)IF(NC.GT.0)callERFAC(COOR,MEL,JR,R,AE)400FORMAT(/2X,'NP=',I3,2X,'NE=',I3,2X,'NM='*,I3,2X,'NR=',I3,2X,'NI='I3,2X,'NL=',I3,2X,*'NG=',I3,2X,'ND=',I3,2X,'NC=',I3)500FORMAT(1X,'TOTALDEGREESOFFREEDOMN=',*I4,1X,'TOTAL-STORAGE','NH=',I5,1X,*'MAX-SEMI-BANDWIDTHMX=',I3)550FORMAT(/20X,'TOTALSTORAGEIS*GREATERTHAN50000')600FORMAT(30X,'NODALFORCES'/8X,'NODE',*11X,'X-COMP.',14X,'Y-COMP.')650FORMAT(/30X,'NODALDISPLACEMENTS'/8X,*'NODE',13X,'X-COMP.',12X,'Y-COMP.')700FORMAT(/30X,'ELEMENTSTRESSES'/5X,*'ELEMENT',5X,'X-STRESS',3X,'Y-STRESS',*2X,'XY-STRESS',1X,'MAX-STRESS',1X,*'MIN-STRESS',6X,'ANGLE'/)STOPENDC*********************************************SUBROUTINEKRS(BR,BS,CR,CS)COMMON/CB/EO,VO,W,T,A,H11,H12,H21,H22*,ME(3),BI(3),CI(3)ET=EO*T/(1.0-VO*VO)/A/4.0V=(1.0-VO)/2.0H11=ET*(BR*BS+V*CR*CS)H12=ET*(VO*BR*CS+V*BS*CR)H21=ET*(VO*CR*BS+V*BR*CS)H22=ET*(CR*CS+V*BR*BS)RETURNEND第九节有限元法程序设计简介CFINITEELEMENTPROGRAMFOR第九节有限元法程序设计简介第九节有限元法程序设计简介作业：本章作业,P.141:6-1,6-2,6-5,6-7(a图)，

6-9可直接用的单刚(注意单元编码要一致)作业：本章作业,P.141:可直接用的单刚课堂测试：1。单元刚度矩阵中子块的物理意义是什么？2。在有限元法中，构造满足哪些条件可以用轴对称问题进展求解？课堂测试：谢谢欣赏谢谢欣赏90FEMch有限元法求解平面问题改前FEMch有限元法求解平面问题改前FEMch有限元法求解平面问题改前有限元法根本思路:离散化构造单元内位移函数；单元位移模式单元分析；划分网格，将连续体划分为有限数量的单元。单元内位移节点位移单元刚度矩阵单元节点力节点位移变分法思想整体分析；总体刚度矩阵节点位移外载荷静力平衡求解；节点位移单元位移模式几何方程物理方程差分法思想FEMch有限元法求解平面问题改前FEMch有限元法求解平面91有限元法根本思路:

离散化

构造单元内

位移函数；单元位移模式

单元分析；划分网格，将连续体划分为有限数量的单元。单元内位移节点位移单元刚度矩阵单元节点力节点位移变分法思想

整体分析；总体刚度矩阵节点位移外载荷静力平衡

求解；节点位移单元位移模式几何方程物理方程差分法思想有限元法根本思路:离散化构造单元内单元位移模式单第一节根本物理量和方程的

矩阵表示第一节根本物理量和方程的

矩阵表示93第一节根本物理量和方程的矩阵表示1根本物理量的矩阵表示外力：应力：应变：位移：虚应变：虚位移：第一节根本物理量和方程的矩阵表示1根本物理量的矩阵表示平面应变问题弹性矩阵：第一节根本物理量和方程的矩阵表示2根本方程几何方程：物理方程：符号矩阵：平面应力问题弹性矩阵：平面应变问第一节根本物理量和方程的矩阵表示2根本方程几第一节根本物理量和方程的矩阵表示位移变分方程虚功方程极小势能原理应力边界条件平衡微分方程等价!节点力:节点对单元的作用力。对单元来说，节点力是作用在单元上的外力。节点位移:节点产生的位移。当单元进入虚位移状态：节点虚位移：单元内应力：单元内位移：单元虚应变：第一节根本物理量和方程的矩阵表示位移变分方程应力边界条件第一节根本物理量和方程的矩阵表示对单元，外力虚功等于内力虚功。虚功方程矩阵表示外力虚功：节点力在节点虚位移上所做的功内力虚功：单元内应力在虚应变上所做的功第一节根本物理量和方程的矩阵表示对单元，外力虚功等于内力第二节构造离散化第二节构造离散化98第二节构造离散化

深梁（离散化结构）将连续体变换为离散化构造：将连续体划分为有限多个、有限大小的单元，并使这些单元仅在一些节点处连接，构成所谓“离散化构造〞。单元节点1单元要素节点：单元与单元之间的连接点〔i,j,m〕。节点位移：节点产生的位移。ij

m节点力：通过节点传递的内力。节点载荷：作用在节点上的载荷〔外力〕。单元位移：单元内位移分布〔u(x,y),v(x,y)〕第二节构造离散化深梁（离散化结构）将第二节构造离散化2单元类型一维单元：如杆单元，梁单元

二维单元：如三角形单元，四边形单元。

三维单元：如四面体单元，六面体单元，棱柱单元。

3连续体离散化模型

单元间仅通过节点连接，没有其它联系；位移，载荷仅通过节点传递；单元内依然是连续体，位移是坐标的连接函数。第二节构造离散化2单元类型一维单元：如杆单元，梁单元第三节单元位移模式解的收敛性第三节单元位移模式解的收敛性101第三节单元位移模式解的收敛性1单元位移函数

将连续体离散化后，各单元间通过节点相连，但每个单元仍然是连续的，即单元内位移是坐标的连续函数：物理意义表示了单元内位移

u,v的分布形式，称单元位移模式。数学意义构造各单元节点位移间的位移插值函数。2形函数设位移函数：变分法在单元上取位移试函数。第三节单元位移模式解的收敛性1单元位移函数第三节单元位移模式解的收敛性假设边界：求解三元一次方程组2）为的代数余子式。

其中：1〕A为三角形ijm的面积，(i,j,m按逆时针编号)第三节单元位移模式解的收敛性假设边界：求解三元一次第三节单元位移模式解的收敛性这里：规律：i,j,m

安逆时针替换。一般写成：第三节单元位移模式解的收敛性这里：规律：i,j,第三节单元位移模式解的收敛性用矩阵形式表示：这里：形函数形函数矩阵那么：几何意义：反映单元位移形态。对三角形单元，在图形上是一个平面。第三节单元位移模式解的收敛性用矩阵形式表示：这里：第三节单元位移模式解的收敛性

位移函数表示了位移u,v

在单元内的分布形式，因此又称位移模式。沟通了单元中离散点的位移和单元内位移之间的关系。3形函数的性质i点：j点：m点：第三节单元位移模式解的收敛性位移函数表示第三节单元位移模式解的收敛性性质1:在节点r处在其它节点处物理意义：不同单元在同一节点处位移相等，即节点处位移与形函数无关，反映了单元在节点处的连续性。总结性质2:总结第三节单元位移模式解的收敛性性质1:在节点r处在其第三节单元位移模式解的收敛性令：各节点位移相等时，单元内位移为常数，且等于节点位移。物理意义：各节点形函数之和为1，反映了单元的刚体位移。是x,y

的线性函数。jmijim由性质1：1m11第三节单元位移模式解的收敛性令：各节点位移相等时，第三节单元位移模式解的收敛性4位移函数的收敛性(1)位移函数收敛性的概念位移模式建立以后，便可逐步求应变、应力、结点力等一系列工作。所以位移模式是有限元法的工作根底。…逐步将单元细分，可以得到同一问题近似解的一个序列：准确解？第三节单元位移模式解的收敛性4位移函数的收敛性(第三节单元位移模式解的收敛性当时，有限元解答收敛于精确解有限元解的收敛性位移函数的收敛性位移函数收敛于正确解(2)位移函数收敛于正确解的条件必须反映单元的刚体位移。单元位移自身变形引起位移刚体位移=+由其它单元变形引起

必须反映单元的常量应变。单元应变常量应变变量应变=+与位置无关与位置坐标有关a,b

两个条件为位移函数收敛的必要条件。第三节单元位移模式解的收敛性当第三节单元位移模式解的收敛性

位移模式应尽可能反映位移的连续性。单元内连续〔位移函数取坐标的连续函数〕；单元公共节点上位移相等；单元公共边界上位移相等。保证离散化后的构造仍为连续弹性体。c条件为位移收敛性的充分条件。满足a,b条件的单元完备单元还满足c条件的单元协条单元第三节单元位移模式解的收敛性位移模式应尽可能反复习：复习：复习性质1:在节点r处在其它节点处性质2:位移函数收敛于正确解的条件:必须反映单元的刚体位移。

必须反映单元的常量应变。

位移模式应尽可能反映位移的连续性。jmijim复习性质1:在节点r处在其它节点处性质2:位移函数收敛于正确例题：例：判断三角形单元位移函数是否满足位移函数的收敛性？解：〔1〕令：因此，u,v

包含了单元的刚体位移。〔2〕因此，应变分量包含了常应变。例题：例：判断三角形单元位移函数例题：〔3〕单元内，u,v

显然是连续函数。对单元①：对单元②：所以，不同单元在公共节点处位移相等，保证了节点处的连续性。u,v

在任一单元内都是坐标的线性函数，在公共边界ij

上当然也是线性的。

经过两点的直线只有一条，所以公共边界ij

上任一点位移相等，保证位移函数在公共边界上是连续的。例题：〔3〕单元内，u,v显然是连续函数。对单元①：对例题：例题：例题：例题：第四节单元分析单元刚度矩阵第四节单元分析单元刚度矩阵118第四节单元分析单元刚度矩阵1.单元的应变和应力(1)单元的应变单元的位移模式反映了单元内任一点的应变与节点位移之间的关系第四节单元分析单元刚度矩阵1.单元的应变和应第四节单元分析单元刚度矩阵应变矩阵应变矩阵子矩阵〔i,j,m〕对三角形单元，均仅与节点坐标有关为常量矩阵

为常量三角形单元为常应变单元(2)单元的应力反映了单元内任一点的应力与节点位移之间的关系第四节单元分析单元刚度矩阵应变矩阵应变矩阵子矩第四节单元分析单元刚度矩阵应力矩阵应力矩阵子矩阵〔i,j,m〕说明：

为常量，单元应力为常量。三角形单元为常应力单元；

对三角形单元，不同单元之间应力、应变有突变。位移函数误差量级：应力应变误差数量级：即使位移精度满足了要求，应力精度不一定满足要求措施：a,细化网格；b,采用高次单元权衡计算精度和计算量之间的关系第四节单元分析单元刚度矩阵应力矩阵应力矩阵子矩第四节单元分析单元刚度矩阵2.单元刚度矩阵建立节点力与节点位移之间的关系。(1)单元刚度矩阵的推导

以单元为研究对象，假想将单元与节点i切开，节点作用于单元上的力，称为节点力。令节点处产生了一组虚位移：相应单元虚位移：相应单元虚应变：

由虚功方程：外力（节点力）在虚位移（节点虚位移）上的虚功，等于应力在虚应变上的虚功，即：第四节单元分析单元刚度矩阵2.单元刚度矩阵建第四节单元分析单元刚度矩阵（与坐标x,y

无关）对任意虚位移上式都必须满足，则：令：则：单元刚度矩阵

有限元实际运算中大量的计算工作量是形成单元刚度矩阵单元刚度方程第四节单元分析单元刚度矩阵（第四节单元分析单元刚度矩阵3.单元刚度矩阵的计算ijm

ijm每个分块子矩阵位置与节点号码（逆时针）顺序相对应平面应力单元刚度矩阵计算公式平面应变问题的刚度矩阵第四节单元分析单元刚度矩阵3.单元刚度矩阵的第四节单元分析单元刚度矩阵4.单元刚度矩阵的物理意义当节点j

处产生单位位移，而其余点完全被约束时，节点i

处所需施加的节点力。当节点j

处产生单位水平位移时，节点i

处所需施加的水平节点力。即时

时

时

时对子矩阵，如：第四节单元分析单元刚度矩阵4.单元刚度矩阵的第四节单元分析单元刚度矩阵4.单元刚度矩阵的性质

对称性：

奇异性：从物理意义证明：直接由单刚计算公式证明：令：则：单元内无应力应变产生以第1行为例：第四节单元分析单元刚度矩阵4.单元刚度矩阵的第四节单元分析单元刚度矩阵上式对任一都成立，则一定有：为奇异矩阵

三角形单元不变性：

与无关。第四节单元分析单元刚度矩阵上式对任一第四节单元分析单元刚度矩阵

与无关。

与无关。

代表单元刚体平移。单元刚度矩阵与单元的刚体平移无关另外可以证明：单元刚度矩阵与单元的刚体转动也无关。影响单元刚度矩阵的因素：单元形状单元所取位移函数材料属性：第四节单元分析单元刚度矩阵例题：例题：应用变分原理推导单元刚度方程(1)单元应变能由于上式中：代入应变能公式：仅与节点坐标有关：应用变分原理推导单元刚度方程(1)单元应变能由于上式中：应用变分原理推导单元刚度方程令：,则：(2)单元外力势能(3)单元总势能根据极小势能原理：应用变分原理推导单元刚度方程令：第五节载荷移置第五节载荷移置132第五节载荷移置1.静力等效原那么刚体静力等效原那么只从运动效应来考虑，得出移置荷载不是唯一的解；变形体的静力等效原那么考虑了变形效应，在一定的位移模式下，其结果是唯一的，且也满足了前者条件的。所以在FEM中，采用变形体的静力等效原那么。刚体静力等效原那么:变形体静力等效原那么:使原荷载与移置荷载的主矢量以及对同一点的主矩也一样。在任意的虚位移上，使原荷载与移置荷载在虚位移上所做的虚功相等。单元原载荷在虚位移上做的虚功=移置后节点载荷在相应虚位移上做的虚功第五节载荷移置1.静力等效原那么刚体第五节载荷移置2.移置方法(1)集中力的移置集中力按形函数在各节点分配第五节载荷移置2.移置方法(1)集中力的移置集中力按第五节载荷移置(2)体力的移置取微单元，那么微单元上的体力：可看作集中力。例1：重力移置。令密度为,则:第五节载荷移置(2)体力的移置取微单元，那么微单元上的第五节载荷移置(3)面力的移置可看作集中力。取微单元，面积为。则微单元上的力：则:例2：三角形单元一边受沿x方向均布载荷。q第五节载荷移置(3)面力的移置可看作集中力。取微单元，第五节载荷移置q例3：三角形单元一边受沿x方向线性分布载荷，在i

点为q，在j点为0，求移置后的节点载荷。令ij边长为as第五节载荷移置q例3：三角形单元一边受沿x方向线性分布第六节整体分析总体刚度矩阵第六节整体分析总体刚度矩阵138〔1〕节点协调原那么：节点处保证协调连接不同单元在同一节点处位移相等〔2〕节点平衡原那么：①②第六节整体分析总体刚度矩阵节点平衡方程〔1〕节点协调原那么：节点处保证协调连接2.总体刚度矩阵形成节点载荷：节点1：①①③③节点2：①①②②节点3：①①②②④③③④节点4：③③④④节点5：①③①③①②①②③①②④③①②④③④③④②④②④第六节整体分析总体刚度矩阵②②④④节点平衡原那么2.总体刚度矩阵形成节点载荷：节点1：①①③③节点2：①

将上面10个方程用矩阵表示：①③①③①②①②③①②④③①②④③④③④②④②④①①①①①①②②②②②②③③③③③③④④④④④④单元①单元②单元④单元③12345节点第六节整体分析总体刚度矩阵节点平衡方程组将上面10个方程用矩阵表示：①③①③①②①②③①②④③①②用分块矩阵表示：①①①②②②③③③④④④单元①单元②单元④单元③12345节点用零（矩阵）升阶使节点力矩阵与节点载荷矩阵同阶第六节整体分析总体刚度矩阵用分块矩阵表示：①①①②②②③③③④④④单元①单元②单元④单列出各个单元单元刚度方程：单元①①①①单元②②②②①①①①①①①①①②②②②②②②②②第六节整体分析总体刚度矩阵单元位移列阵用整体位移列阵代替单元刚度矩阵用零升阶与总体刚度矩阵同阶，又称单元奉献矩阵。①②列出各个单元单元刚度方程：单元①①①①单元②②②②①①①①①单元③③③③单元④④④④③③③③③③③③③④④④④④④④④④(1)

扩阶后的单元刚度方程只不过增加了两个0=0的方程，因此与原方程等价。第六节整体分析总体刚度矩阵〔3〕用构造的整体位移矩阵代替单元节点位移矩阵-----节点协调原那么(2)不影响虚功及应变能的计算。③④单元③③③③单元④④④④③③③③③③③③③④④④④④④④④④将以上四式代入节点平衡方程组：①②③④①①①②②②③③③④④④节点平衡方程组令：①②③④总体刚度矩阵整体平衡方程,沟通节点载荷与节点位移关系第六节整体分析总体刚度矩阵将以上四式代入节点平衡方程组：①②③④①①①②②②③③③④①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④①②③④总体刚度矩阵的形成：所有升阶后的单元刚度矩阵〔单元奉献矩阵〕的叠加！第六节整体分析总体刚度矩阵①①①①①①①①①②②②②②②②②②③③③③③③③③③④④④④④④④④④①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③3.总体刚度矩阵的计算规律

令：为总体刚度矩阵的任一子块:当时：为主对角线上子块，由环绕节点r或s的单元在单刚矩阵相应节点上子矩阵的叠加如：①+②+③+④

(2)当且rs为相邻单元的公共边时：为相邻两单元在单刚矩阵相应子块矩阵的叠加。如：①+③

(3)当且rs仅为某一单元的一边时：为该单元单刚矩阵相应子块矩阵。如：①

(4)当且rs为互不相关的两节点时：如：①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④第六节整体分析总体刚度矩阵①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④3.总体刚度矩阵的计算规律令：3.总体刚度矩阵的性质

对称性①+③①+③①③①③①+③①③①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④存贮一半元素加上对角线上的元素第六节整体分析总体刚度矩阵3.总体刚度矩阵的性质对称性①+③①+③①③①③①

奇异性①+③③①①+③①①+②①+②②①+③①+②①+②+③+④③+④②+④③③+④③+④④②②+④④②+④以第一行为例,考察第一行之和：①+③③①①+③第一行子块矩阵有：①③单元刚度矩阵第一行元素之和①单元刚度矩阵第一行元素之和③反映了弹性体的刚体位移①①①①①①①①①③③③③③③③③③第六节整体分析总