有限单元法课件

有限单元法和ANSYS程序(FiniteElementMethodandANSYSProgram)§

5.1变分法(VariationMethods)简介

变分法的定义:主要是研究泛函(functional)及其极值的求解方法.所谓泛函,就是以函数作为自变量的函数,即泛函就是函数的函数.弹性力学变分法所研究的泛函,主要是指能量(energy),能量是应力应变的函数,而应力应变又是坐标、力、材料性质等的函数，所以能量是泛函。因此，弹性力学中的变分法又称为能量法(energymethod).一、形变势能

按照材料力学(themechanicsofmaterials),当弹性体仅受一个方向力作用,例如x方向，其stresscomponent为

x,相应的

straincomponent为

x,那末单位体积的应变能（

strainenergy）

is

x

x/2.即U1又称为变形势能密度或单位比能（specificenergy）．同样，当单元体只受剪stresscomponents

xy

作用时，相应的剪straincomponent

xy时,单位体积的应变能is

xy

xy/2.

设弹性受全部六个应力分量的作用，对于平面问题,,而如果是平面应力问题,σz=0,如果是平面应变问题，εz=0,所以无论对任何平面问题，单位体积的应变能都为：(a)(b)因此，物体总的应变能为：利用物理方程可得到用应变表示的单位体积的应变能为(c)(e)(e)式分别对三个应变分量求导,得到上式说明，弹性体每单位体积的应变能，对于任一形变分量的微分就等于相应的应力分量，将几何方程代入方程(e),得用位移分量表示的单位比能的形式：(f)(5.1)又由方程(c),得在以上式子中，将Ｅ变为，将就得平面应变问题的相应方程．二、外力势能（potentialenergyofexternalforce)

若弹性体受外力作用，平面区域Ａ内的体力分量为Ｘ，Ｙ，边界sc上的面力分量为，则外力所做的功称为外力功(workofexternalforce)：(5.2)由于外力做了功，消耗了外力势能，因此，在发生实际位移时，弹性体的外力势能是三、变分方程(variation

equation)

１、虚位移(virtualdisplacements)的概念设有平面问题中的任一单位厚度的弹性体，在一定外力作用下处于平衡状态．命u,v为实际的位移分量，它们满足平衡方程和边界条件．现假想这些位移分量发生了边界条件所容许的微小位移，即所谓的虚位移或位移变分，则有

2、位移变分方程(equationofdisplacements

variation)

由于位移的变分，引起外力功的变分δＷ和外力势能的变分δＶ为由于位移的变分，引起应变的变分为从而引起形变势能的变分为根据能量守恒定理，形变势能的增加应该等于外力势能的减少，也就等于外力所做的功，故有因而得这就是位移变分方程,又称为位移变分Lagrange’s方程３、虚功方程(virtualworkequation)

(5.3)将(5.1)代进上式,得将上式代入位移变分方程(5.3)中,得(5.4)(5.4)式就叫做虚功方程(virtualworkequation).它表示：如果在虚位移发生之前，弹性体处于平衡状态，那么，在虚位移过程中，外力在虚位移上所做的功就等于应力在虚应变上所做的虚功．§5.2基本量和基本方程的矩阵表示(MatrixRepresentationsofBasicQuantitiesandEquations)

在有限单元法中,为了简洁清晰地表示各个基本量以及它们之间的关系,也为了编写程序以便利用电子计算机进行计算,广泛地采用矩阵表示和矩阵计算.1体力(bodyforces)的矩阵表示:设物体受体力P的作用,它的两个分量为X和Y,则有:

２

面力(faceforces)的矩阵表示:设物体受面力的作用,它的两个分量为

则有:(5.5)(5.6)３

应力(stress)的矩阵表示:４

应变(strian)的矩阵表示:５

位移(displacement)的矩阵表示:(5.7)(5.8)(5.9)６几何方程(geometricalequations)的矩阵表示:７

物理方程(physicalequations)的矩阵表示:平面应力情况下简写的形式：(5.10)(5.11)(5.12)这里［Ｄ］称为弹性矩阵(elasticmatrix)对于平面应变问题，可将上式中的Ｅ变为，即得平面应变问题的弹性矩阵（８）虚功方程(virtualworkequation)的矩阵表示虚功方程(5.4)可以表示为：(5.13)(5.1４)这里，在finiteelementmethod中,作用于弹性体的分布力（distributedforces）

一般用作用在某些点等价的集中力（equivalentconcentratedforces）来代替，如

Fig5.1所示．假定作用在

ipoint的集中力是

FxiandFyi

，假定作用在

jpoint的集中力

是

FxjandFyj

respectively,etc.这些分量和相应的虚位移分别用以下列阵表示：Fig.5.1(5.17)

(5.1８)

于是各外力在虚位移的虚功为：

于是由（５.14),得这就是在集中力作用下的虚功方程的矩阵表示形式(5.19)§5.3有限元导言

(IntroductionofFiniteElementMethod)

有限单元法出现于40年代，被应用于飞机结构分析，有限元这个术语是1956年Turner首先使用的。一般来说，任何能用微分方程描述的物理现象，都能够通过变分原理建立的有限元方法来模拟，所以有限元的应用非常广泛。不但可以解决像荷载—位移问题，还可以解决渗流、热传导等问题。有限元法的基本思想是用分段逼近(piecewiseapproximations)的方法,如利用多项式(polinomials)来代替连续函数(continuousfunctions),即是将连续体无限个节点变成有限个足够多的节点，通过建立节点力与节点位移之间的关系，求解节点位移，又通过给出一定的假设建立单元内任一点位移与节点位移之间的关系，计算出应变，运用应力—应变关系，从应变就可以计算出应力。Fig.5.2

在Fig.5.1(a)中，表示一个承受初始应力px,py,pxy

的一个地下洞室(undergroundopening)；在Fig.5.1(b)中，表示将地下洞室周围区域划分为很多三角形单元(triangularelements)；在Fig.5.1(c)中,表示节点(nodes)为I,j,m

的一个三角形单元．有限单元法的基本解题步骤为：1.划分单元；2.建立位移模型，即建立单元内任一点位移与节点位移之间的关系

设三角形单元三个节点的位移分别为：(ui,vi),(uj,vj),(um,vm),则三角形单元任何一点的位移u与v为：这样，这里［Ｎ］称为插值函数矩阵(amatrixofinterpolationfunctions,{Ue}称为节点位移列阵(acolumnvectorlistingthenodalDisplacements)或(5.20)3.推导单元的刚度矩阵方程。既建立节点力与节点位移之间的关系．如果知道任一点的位移，则其应变为：或(5.21)将(5.20)式代入(5.21)式，得：(5.22)这里，(5.23)[B]称为应变与节点位移关系矩阵(therelationshipmatrixbetweenstrainsandnodaldisplacements.而应力由下式确定：(5.2４)[Ｄ]称为弹性矩阵（theelasticitymatrix）.在有限单元法中，作用在单元所有的力都用等价的节点力矢量(equivalentnodalforcevector)来表示:根据虚功原理:在虚位移过程中，外力在虚位移上所做的功就等于应力在虚应变上所做的虚功．而外力的虚功为:(a)应力在虚应变上所做的虚功为：(b)根据虚功原理，这两个功应相等，故有：

称为单元的刚度矩阵（thestiffnessmatrixofaelement）.

我们就得到单元的刚度方程(thestiffnessequationofaelement):

如果记：

4.建立有限单元法的基本方程：将所有单元组合起来，就得到有限单元法的总体刚度方程

其中：[K]—总体刚度矩阵，

—位移列阵，各个节点的位移，

—荷载列阵，将荷载化为节点力列阵5.求解节点位移(考虑边界条件)。6.通过位移求解应力和应变。§5.4有限元网格(FiniteElementNetwork)1.常用单元(1)线性单元(linearelements)：三角形、矩形或其他四边形。形函数是线性函数，即单元内任一点的坐标可用单元节点坐标的线性函数来表示。(2)抛物形单元(parabolaelements)：除了角上有节点（主外节点），边缘上也有节点（副外节点）的单元。有时在内部也有节点（内节点），直边或曲边均可。单元内任一点的坐标，可用一个抛物线内插法来求得。121231234图5.3线形单元及节点位置2.单元划分要注意的几个问题

(1)相邻两个单元的节点要与节点重合（外节点与外节点、内节点与内节点），不能与无节点边重合。

(2)单元不必是相同尺寸，应力有突变的地方，单元划分应较小。

(3)任何一个单元必须只能在一种材料区，即它不能跨越两种材料区。

(4)同一单元的各个边长，一般不要相差太大。5216436312456783124578125463图5.4抛物线形单元及节点位置§5.5位移插值函数DisplacementInterpolationFunction(形函数,插值函数）(ShapeFunction)1.三节点三角单元的位移插值函数

所谓位移插值函数是指表是节点位移与任一点位移之间的关系的函数,即:从上一节我们知道,三节点三角单元的位移插值函数矩阵为:现在讨论如何确定[N]的各个分量,我们知道,组成[N]各分量的插值函数,当它的变量等于某个节点的值时,它应该等于该节点的位移值,这就要求:这里,等这里[I]和[0]表示单位矩阵和零矩阵.同时,一个点位移的两个分量可通过同样的方式插值而得,因此很明显有:最简单的位移插值函数为线性插值函数,即假定为(5.28)

这些系数等是待定系数,通过解下列联立方程(simultaneousequations)确定:

从以上六个方程解出6个待定系数这里一定是的线性函数,即(5.29)因而可将(5.28)写成节点位移函数的形式：而一定是的线性函数,即这里，

(i,j,m)(5.30)而(5.3１a)(i,j,m)(5.31c)(5.3１b)A=areaofthetriangularelement

2.四节点四边形单元插值函数对于四节点四边形单元，单元内任意一点的位移与节点位移的关系，也是通过形函数建立的，以下以一个四节点的正方形单元为例，说明形函数的含义。取局部坐标如下图所示，四个节点的坐标分别是(1,1),

(1,-1),

(-1,-1),

(-1,1).（1，1）（-1，1）（-1，-1）（1，-1）1234（0，0）ξη图5.5四边形单元映射为面积2×2个单位的正方形单元假定四个节点分别有位移，则单元内任一点的位移可用下式表示：即：（1）则是分别与节点1，2，3，4相关的形函数，形函数Ni（i=1,2,3,4）是坐标点的函数，即

形函数有如下的性质，当把第i个节点坐标代入时：如将节点1的坐标（1，1）代入，那么就有形函数怎么得来的：（1）从观察和直觉中得到例如：统一的式子可以写为：

（2）用数学过程来导出形函数的表达式对于矩形单元，可假设位移是坐标的双线性函数。

式中，是待定系数。将4个节点位移值代入上述方程得:

这里8个方程8个未知数，可以求得，并将其代入上述假定式（2），与（1）式相比较，就可以求得值。

3.坐标变换(CoordinateTransformation)

在进行每个单元的计算时往往采用局部坐标，而每个单元的局部坐标与整体坐标往往是不一致的，故需要进行坐标变换。为什么要对单元采用局部坐标，这里因为在局部坐标下很多积分计算比较容易，例如函数在正方形区域上积分，可写为：比在一个不规则的原四方形域内的积分容易的多。

从整体坐标到局部坐标的变换，可以看成是一种映射，例如从整体坐标下任意形状的四节点单元映射成局部坐标下的两个单位边长的正方形单元。我们可以通过形函数来映射，即可假定单元内点的坐标满足：式中，

分别为整体坐标下的四节点单元四个节点坐标，

代表形函数，它是局部坐标的函数，这就建立了整体坐标与的变换式，上式可写为：

或：以下，我们介绍等参数单元的概念：对于位移，我们已有

如果对于单元内任一点的坐标也存在如下关系：即位移的插值公式与坐标的插值公式具有相同的形状（位移插值函数=坐标插值函数），则这样的单元就叫做等参数单元。我们一般都采用等参数单元.由数学我们知道，对于无限小区域dA有，

雅可比矩阵的元素很容易通过下式求得：式中

叫做雅可比矩阵的行列式。这样整体坐标下的积分就可以变成局部坐标下的积分：在求刚度矩阵中，我们还需要求形函数对笛卡儿系坐标的导数，如求：写成矩阵形式：⒋应变位移关系(Therelationshipbetweenthestrainanddisplacement)①.平面应力或平面应变：②.轴对称情况：r为径向，y为对称轴§

5.6等效节点力（EquivalentNodalForces)

在有限单元法中，所有的力，无论是载荷重力、温度应力、初应力还是由初应变、孔隙压力、离心力、惯性力等引起的外力都要转换到作用于节点上的等效力，以下我们来具体分析。5.6.1集中力引起的等效节点力设单元边上作用集中力R,则利用力矩平衡条件可求得:

图5.6集中力的等效节点力(5.32)5.6.2面力的等效节点力(Theequivalentnodalforcescausingbydistributedforces)

由作用在单元边界或表面上的面力引起的等效节点力(这是由虚功方程推得的,证明从略)

如果知道表面正应力σ和剪应力τ,则其等效节点力为:θ为单元的荷载边界与总体坐标x轴的夹角。式中,T是转换矩阵，(5.33)5.6.3单元体力引起的等效节点力（Theequivalentnodalforcescausingbybodyforces

）

根据虚功原理,可求出单元体力引起的等效节点力.设体力为则等效节点力为:(5.34)

5.6.4初应变引起的等效节点力（Theequivalentnodalforce仿

causingbyinitialstrains

）初应变通常当作零来处理，但当我们试图模拟洞室开挖、岩石加固、灌浆等过程，在这些过程中的各个阶段，其前一阶段产生的应变就被当作初应变，初应变也可能是外力以外的其他因素引起的，如蠕变、渗流。

根据虚功原理可以推导出初应变引起的等效节点力

若某种原因，初应变产生变化，则采用增量公式来计算由产生初应变的变化引起的节点力增量。(5.35)5.6.5初应力（或初应力增量）引起的等效节点力（节点力增量）（Theequivalentnodalforcescausingbyinitialstresses

）根据虚功原理可以推导出初应力引起的等效节点力或写成增量的形式：

(5.36)5.6.6孔隙压力引起的等效节点力5.6.7温度应力在受到约束的岩体中，温度的改变会引起应力和应变。如果岩土的热膨胀系数为，温度的改变产生的线应变为单元由于温度改变而引起的等效节点力以类似于初应变方式处理。

与初应力的处理方法一样。

这里

为孔隙压力的增量，，，

§5.7单元刚度矩阵(StiffnessMatrixofanElement)

假如设插值函数为线性,将

Eq.(5.30)代入

Eq.(5.23),我们就得到三节点三角形单元应变与节点位移的关系矩阵[B](relationshipmatrixbetweenstrainsandnodaldisplacementsforthethreenodestriangularelement),(5.37)这里A为三角形的面积,bi,bj

等由(5.3１)确定．将(5.13)和(5.37)代如方程(5.25)，就得(5.38)这里

[krs]isamatrix(r=i,j,k;s=i,j,k),对于

planestress问题:

(5.39)对于planestrainproblems,弹性矩阵(elasticitymatrix)变成:(5.40)所以在planestrain情况下,三角单元的刚度矩阵(stiffnessmatrix)变为:

(5.41)

§5.8有限单元法的总体刚度矩阵

(

GlobalStiffnessEquationofFiniteElementMethod)有限单元法包含着对每个单元都要计算它的刚度矩阵(stiffnessmatrix[ke])和等效节点力vector{Fe}.对于单元刚度矩阵，简单的单元可以直接积分解出；对于高阶单元，还需要数值积分。将这些刚度“叠加”在一起，就可以得到总体结构的刚度矩阵。相加时是把对应点的对应项相加。设整个结构的节点数为