第2章-弹性力学基础及有限元法的基本原理1

第2章有限元法的基本原理

—平面问题有限元法（1）力学应用引言（2）力学回顾力学回顾力学回顾第一节弹性力学基础知识1.1弹性力学中的物理量1、载荷定义：作用在弹性体上的力（力矩），又称外力。载荷可分为：体力、面力、集中力1）体力定义：分布于整个弹性体体积内的外力。如：重力可以分解为三个坐标系上的分量，用向量表示为：2）面力定义：作用于弹性体表面上的外力，如：流体压力可以分解为三个坐标系上的分量，用向量表示为：3）集中力定义：集中在某一点上的外力，如：牵引力可以分解为三个坐标系上的分量，用向量表示为：2、应力（注意下标）定义：弹性体内某一点作用于某个截面单位面积上的内力，反映了内力在截面上的分布密度。微元体表面上的应力：一个正应力（拉压）两个切应力（剪切）切应力互等定律：符号标注规则：第一个角标是所在面垂直于哪个轴；

第二个是指向哪个轴弹性体内某一点的应力状态由六个应力所决定应力向量可以表示为：3、应变（对应于应力）定义：微元体体发生变形后，单位长度的变形量。对应于应力，应变向量可以表示为：4、位移定义：弹性体内质点位置的变化。位移向量可以表示为：1.2弹性力学的基本方程主要是描述应力、应变、位移及外力间的相互关系1、平衡方程（应力间的关系）[L]微分算子2、几何方程（应变与位移的关系）3、物理方程（应力与应变之间的关系）杨氏弹性模量泊松比弹性矩阵E杨氏弹性模量1D问题的弹性模量泊松比是指材料在单向受拉或受压时，横向正应变与轴向正应变的绝对值的比值，也叫横向变形系数，它是反映材料横向变形的弹性常数。可以这样记忆：空气的泊松比为0，45#钢0.3，水的泊松比为0.5，中间的可以推出。若在弹性范围内加载，横向应变εx与纵向应变εy之间存在下列关系：εx=-νεy式中ν为材料的一个弹性常数，称为泊松比。泊松比是量纲为一的量。未知数应力6个+应变6个+位移3个=15个方程个数平衡方程3个+几何方程6个+物理方程6个=15个原则上可以根据15个方程求出15个未知物理量但实际求解时先求出一部分再通过方程求解剩下的。目前有限元法主要采用的是位移法，以三个位移分量为基本未知量。位移-应变-应力，应力和外力平衡1.3虚位移原理（1）虚功与虚应变能弹性体在外力作用下变形，外力对弹性体做功，所做的功以应变能的形式储存于弹性体中。弹性体单位体积的应变能为：虚位移定义：在约束条件允许的范围内弹性体可能发生的任意微小位移。虚位移与时间及外载荷无关实际位移是在外载荷作用下可能的虚位移弹性体在平衡状态下发生虚位移1）外力所做的虚功为：2）应力在虚应变上所做的虚功，也就是存储在弹性体内的虚应变能为：（2）虚位移原理表述：如果在在虚位移发生之前弹性体是平衡的，那么在虚位移发生时外力在虚位移上所做的功就等于弹性体的虚应变能。即：当外力的形式是多样的时，外力的虚功等于：1.4平面问题定义严格地讲，任何结构都是空间的。对于某些特殊情况，空间问题可以转化为平面问题。（1）平面应力问题满足条件：1）几何条件厚度尺寸远远小于截面尺寸；2）载荷条件载荷平行于板平面且沿厚度方向均匀分布，而板平面不受任何外力作用。此时，应力应变分量变为：几何方程物理方程简化为8个方程（2）平面应变问题满足条件1）几何条件结构呈等截面的细长形；2）载荷条件载荷垂直于厚度方向(平行横截面)且沿厚度均匀分布，两个端面不受力。此时，应力应变分量变为：物理方程简化为8个方程1.5弹性力学的参量及方程汇总1）参量位移：应力：应变：2）方程平衡方程：几何方程：物理方程：其中：

第2章作业一：手画微分体应力分量图。回顾弹性力学的参量及方程汇总1）参量位移：应力：应变：2）方程平衡方程：几何方程：物理方程：3）虚位移、虚功原理4）平面应力问题；平面应力问题第二节平面问题的有限元法2.1平面问题的有限元分析步骤（平面应力问题）1、结构离散离散：将一个连续的弹性体分割为一定形状和数量的单元的组合

单元也称为网格连续体→有限个单元的组合体可用于离散的单元：三角形单元；矩形单元；不规则四边形单元。DOF节点的自由度：节点所具有的位移分量的数量。一个单元所有节点的自由度总和称为单元自由度。（1）单元参数只能通过节点传递到相邻单元（2）单元和节点必须统一编号2.2单元分析（位移、应力、应变）任务：形成单元刚度矩阵，建立单元特性方程因此必须建立坐标系，如下图：1）位移函数分片插值→假设一种函数来表示单元位移分布一般选取多项式（简单而且易求导）对于三角形单节点单元（DOF=6）根据克莱姆法则（Cramer，瑞士数学家），可求得：[说明1][说明2]位移函数应该满足以下几个条件（1）包括常数项（保证刚体位移）（2）包括一次项（保证常应变）（3）保证位移的连续性（性质3保证）（4）各项几何同性（x,y应该是可以互换的）满足上述三个条件的目的是满足有限元的收敛性（1）和（2）是收敛的必要条件完备性条件（3）是收敛的充分条件协调条件注意：非协调单元的解不一定不收敛2)单元应力和应变将位移表达式（2-3）和（2-4）代入几何方程得：则应变向量可以表示为：

[D]、[B]均为常数矩阵，因此三角形三节点单元为常应力单元。3）单元刚度矩阵设作用在单元节点上的单元节点力列阵为：而节点发生的虚位移为：则节点力在虚位移上做的虚功为：将（2-5）展开得：单元刚度矩阵的物理意义：在一个节点处产生单位位移而其他点为零时，在该节点上需要的外力的大小。单元某个元素的影响单元刚度矩阵的特性（1）对称性：（弹性力学互等定理）（2）奇异性：（刚体位移）不等于零，无力下单元无位移，与弹性力学理论相悖。2.3总刚度矩阵的集成通过单元特性方程并不能求出单元节点位移。因为包含单元间的作用力。因此，必须将每个单元的特性方程相加消除内力的影响。这就是总刚度矩阵集成的目的。1）总刚集成原理在整个结构中，一个节点为几个单元共有。在第i个节点处的平衡方程为：例：总刚的形成过程1234562）总刚度矩阵集成过程（1）扩阶过程（可由转换矩阵完成）（2）叠加过程：总刚矩阵的特点：（1）对称性节省存储容量（2）稀疏性可能存在大量零元素（3）带状性半带宽与节点的编号有关（4）奇异性保证刚体位移2.4载荷移置移置可能在局部产生误差，但不会影响整个结构的力学特性。1）集中力的移置（虚功等效）2）面力的移置3）体力的移置

2.5约束处理1）边界位移为零2）边界位移为已知量2.6求解线性方程组

2.7计算其它物理量2.8计算结果处理2.9结果显示、打印、分析%-----------------------------------------------------------functionk=Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)%该函数计算单元的刚度矩阵%输入弹性模量E，泊松比NU，厚度t%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)%输出单元刚度矩阵k(6X6)%---------------------------------------------------------------A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=[betai0betaj0betam0;0gammai0gammaj0gammam;gammaibetaigammajbetajgammambetam]/(2*A);ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];endk=t*A*B'*D*B;（1）Triangle2D3Node_Stiffness第三节平面3节点三角形单元及MATLAB编程%-----------------------------------------------------------functionz=Triangle2D3Node_Assembly(KK,k,i,j,m)%该函数进行单元刚度矩阵的组装%输入单元刚度矩阵k%输入单元的节点编号I、j、m%输出整体刚度矩阵KK%---------------------------------------------------------------DOF(1)=2*i-1;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;DOF(5)=2*m-1;DOF(6)=2*m;forn1=1:6forn2=1:6KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);endendz=KK;（2）Triangle2D3Node_Assembly%-----------------------------------------------------------functionstress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)%该函数计算单元的应力%输入弹性模量E，泊松比NU，厚度t%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)，单元的位移列阵u(6X1)%输出单元的应力stress(3X1)，由于它为常应力单元，则单元的应力分量为Sx,Sy,Sz%---------------------------------------------------------------A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=[betai0betaj0betam0;0gammai0gammaj0gammam;gammaibetaigammajbetajgammambetam]/(2*A);ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];endstress=D*B*u;（