第四章 平面问题有限单元法

第四章

平面问题的有限单元法弹性力学的平面应力问题基本条件（1）等厚度的薄板；（2）体力作用于体内，平行于板的中面，沿板厚不变；（3）面力作用于板边，平行于板的中面，沿板厚不变；（4）约束作用于板边，平行于板的中面，沿板厚不变。弹性力学的平面应力问题坐标系：由于两板面上无面力和约束作用：由于薄板很薄，应力是连续变化的，又无z向外力，可认为：简化为平面应力问题，仅剩：其值与z无关弹性力学的平面应变问题基本条件（1）很长的常截面柱；（2）体力作用于体内，平行于横截面，沿柱体长度方向不变；（3）面力作用于柱面，平行于横截面，沿柱体长度方向不变；（4）约束作用于柱面，平行于横截面，沿柱体长度方向不变。弹性力学的平面应变问题坐标系：由于截面、外力、约束沿z向不变，外力、约束平行xy面，柱体非常长：故任何z面（截面）均为对称面。简化为平面应变问题：其值与z无关以平面应力问题为例介绍

平面问题的有限单元法平面应力单元类型平面应力单元类型简介3节点三角形单元4节点4边形单元8节点4边形曲边单元节点位移分量每节点2个位移分量（自由度）x方向的位移u，y方向的位移v单元位移分量（4节点）jik三角形单元单元ekijl单元e四边形单元123456788节点单元单元e平面应力单元网格划分应力梯度变化比较大的地方，网格应密一些有应力集中的地方，网格应密一些单元边界长度不要相差过大单元各边夹角不要太大集中载荷处要设置节点结构不同材料交界面处要设置节点并作为单元边界结构厚度突变处要设置节点并作为单元边界分布载荷突变处要设置节点施加位移约束处要设置节点注意单元间的连接平面应力单元网格划分设置节点设置节点材料A材料B界面这样不行病态单元a-边长差别太大b-边长差别太大c-边夹角太大abc单元节点信息节点信息节点号xyz100021003200421051106010700.50810.50920.50单元拓扑信息单元号节点i节点j节点k节点l材料编号其它常数112871278561358941439881582331以最经典的三角形单元为例单元位移函数（位移模式）单元位移模式概念单元内任一点的位移要用节点上的位移值近似表达出来，这就需要假定一个近似函数来表示单元内的位移分布，所选择的近似函数就称为单元位移函数或单元位移模式。对于弹性力学平面问题，一般选择多项式(polynomial)来作为单元内的位移解或插值函数或位移模式。一维单元二维单元多项式的项数越多，结果就越精确，但取多少项由单元形式决定。三角单元的位移函数节点上只有六个位移分量，所以单元内部位移函数的待定参数不能超过这个数目。可假设单元内部位移为x、y的线性函数：参数ai由位移边界条件确定。三角单元的位移函数节点i节点j节点k于是：三角单元的位移函数如果令则：根据线性代数的知识，可知：三角单元的位移函数T*为T的伴随矩阵其中：三角单元的位移函数把求得的系数代入位移函数公式：得到：三角单元的位移函数表达为矩阵形式：这里：Ni,Nj,Nk是坐标的函数，它们反映了单元的位移形态，故称为三角单元的形态函数（或形函数）三角单元的位移函数形函数具有明确的几何意义：如图所示三角单元IJK，P为三角单元内任意一点，其坐标为(x,y)P点在三角单元各角点上产生的形函数分别是Ni,Nj,Nk同理：三角单元的位移函数位移函数运用示例：已知各节点位移为：求P点位移P点的位移可由节点位移近似表示为三角单元的位移函数于是：三角单元的位移函数形函数的本质计算点(x,y)的位移u(x,y)、v(x,y)可用单元内各节点的位移值ui,vi的加权之和来近似表示，其中，各节点位移加权系数为关于计算点(x,y)的函数，即为形函数三角形单元形函数的性质1、单元节点产生的形函数值为1或02、形函数之和等于1位移函数应满足的条件应满足：单元内位移模式必须是连续的，公共边上位移必须协调位移模式必须反映单元的刚体位移位移模式必须反映单元的常应变可以证明三节点三角形单元是收敛的完备单元和协调单元三条准则：1、位移模式必须包含单元的刚体位移2、位移模式必须能包含单元的常应变3、位移模式在单元内要连续、并使相邻单元间的位移必须协调满足条件1、2的单元为完备单元满足条件3的单元为协调单元应变的离散过程应变的离散过程根据弹性力学中的几何关系，单元内任一点(x,y)的应变表达式为矩阵形式应变的离散过程应变的离散过程单元内任一点(x,y)的位移(u,v)可以采用节点位移近似表示：将其代入应变表达式，则应变的离散过程应变的离散过程为书写方便，应变分量矩阵可用分块矩阵表示简写为：B也称为“应变矩阵”应变的离散过程应变的离散过程由于形函数所以刚才的矩阵事实上可表示为B矩阵中的所有元素已经由三角形单元的节点坐标确定。应变在单元内为常数，所以又称为常应变单元。应力的离散过程应力的离散过程根据广义虎克定律，对于平面应力问题：矩阵形式物理方程应力的离散过程应力的离散过程如果令：物理方程简写为D又称为“弹性矩阵”将前面应变的表达式代入，则虚位移与虚应变我们已经知道了应变与位移的关系那么很自然的如果发生了虚位移则会发生虚应变虚功原理建立控制方程外力虚功等于内力虚功。外力虚功内力虚功其具体可计算为：虚功原理建立控制方程外力虚功等于内力虚功。结果：考虑到节点虚位移的任意性：上式即为有限元控制方程。此处K称为“刚度矩阵”刚度矩阵如果将求解域划分为多个单元，则即总体刚度矩阵（总刚）单元刚度矩阵（单刚）单元刚度矩阵三节点等厚三角形单元中B和D的分量均为常量，则单元刚度矩阵可以表示为其具体形式为单元刚度矩阵对于平面应力问题，其具体可计算如下：单元刚度矩阵的物理意义把前面获得的有限元控制方程展开：那么事实上就是当节点j产生单位位移时，在节点i上需要施加的节点力。单元刚度矩阵的物理意义更具体一点：当节点i在垂直方向产生单位位移时，在节点i上需要施加的垂直节点力

当节点j在水平方向产生单位位移时，在节点j上需要施加的水平节点力当节点j在垂直方向产生单位位移时，在节点i上需要施加的水平节点力当节点i在水平方向产生单位位移时，在节点j上需要施加的垂直节点力单元刚度矩阵的性质性质1：对称性单元刚度矩阵的性质性质2：对角线上元素恒为正单元刚度矩阵的性质性质3：此矩阵为奇异矩阵意义：没有对节点施加位移约束，所以单元产生任何的刚性位移都是可以的，由力得不到位移的唯一解。性质4：此矩阵的各行元素之和为零，由于对称性，各列元素之和也为零。整体刚度矩阵的形成单元刚度矩阵形成后，要将单元组成一个整体结构，即整体分析，基本方法是刚度集成法，即整体刚度矩阵是单元刚度矩阵的集成。整体刚度矩阵的集成是按对号入座的方式叠加的。用下面的三角形薄板作为示例：共计4个单元，单元节点编号为：整体刚度矩阵的形成各个单元的刚度矩阵为：整体刚度矩阵的形成设单元节点总数为N，每个节点的自由度数为NDOF。(对于一维情况，NDOF=1；对于二维情况，NDOF=2；三维，NDOF=3)。如果是二维问题，则总自由度数为2N个，相应的整体刚度矩阵大小为2N×2N阶方阵。整体刚度矩阵的意义与性质Kij表示j自由度发生单位位移，其他位移为零时，第i个自由度上必须施加的节点力。总体刚度矩阵中的元素具有如下性质：（1）主对角元素Kii

>0（2）总体刚度矩阵K是对称的奇异矩阵（3）总体刚度矩阵K是带状稀疏矩阵。

整体刚度矩阵的存储由于整体刚度矩阵具有对称性、稀疏性和非零元素带状分布的特点，所以没有必要将全部的整体刚度矩阵进行存储。（1）利用对称性：只保存整体刚度矩阵上(下)三角的带宽内元素即可；（2）利用稀疏性：在用分块表示的整体刚度矩阵中，与相关节点对应的分块才能具有非零的元素，其他位置上的分块矩阵的元素为零（3）利用带状分布：整体刚度矩阵的非零元素分布在以对角线为中心的带状区域中，每行具有元素的元素的个数叫做“半带宽”，用D表示。

整体刚度矩阵的存储各行的半带宽D怎么计算：整体刚度矩阵的存储可用一维数组A来存储半带宽内的元素，而不必储存所有元素。本例中：总带宽则可以采用如下方式存储：整体刚度矩阵的存储最大半带宽是多少？相邻节点的编码最大差值+1）×NDOFDmax=(10-6+1)×2=10设整体刚度矩阵K是一个n×n的矩阵,其最大半带宽为D，那么利用带状矩阵的特点和对称性，只需要存储以D为固定宽度的元素，这种存储方法称为二维等带宽存储。整体刚度矩阵的存储元素Krs在整刚矩阵K中的行列编码记为r，s，在二维等带宽矩阵K*中的行列编码为r*，s*边界条件的处理非节点载荷的移植：集中力：最好在集中力处设置节点分布面力：分布体力：边界条件的处理事实上，由于整体刚度矩阵的奇异性，仍然是没有办法求解。原因在于位移边界条件没有引入。方法一：划零置一法若已知边界条件：则边界条件的处理事实上，由于整体刚度矩阵的奇异性，仍然是没有办法求解。原因在于位移边界条件没有引入。方法一：乘大数法若已知边界条件：则M是很大的