弹性力学有限元位移法原理.docx

1、弹性力学有限元位移法原理一、课程论文：弹性力学有限元位移法原理（30分） 主要针对一维（直杆）问题，撰写一篇论文，对有限元位移法的原理作一般性虬 括和论述。要求论文论及但不限于下列内容：1）弹性力学有限元位移法的数学、 力学基础和基本思想；2）有限元法求解的原理和过程，推导所有计算列式；对 基本概念和符号进行解释和讨论；3）收敛性、收敛准则及其数学、力学意义的 讨论。弹性力学有限元位移法原理一、有限单元法的起源有限单元法的形成可以追溯到20世纪50年代甚至更早些时间， 基本思路来源于固体力学中矩阵位移法的发展和工程师对结构 相似性的直觉判断。对不同结构的杆系、不同的载荷，用矩阵 位移法求解都可

2、以得到统一的公式。在 1952-1953年期间，R- W- Clough和M- J - Turner在分析飞机三角翼振动问题时，提出了把平面应力三角形或矩形板组合起来表达机翼刚度的方法，当时被称为直接刚度法。1956年M- J - Turner，R - W- Clough， H C Martin，L - J - Topp在纽约举行的 航空学会年会上发表论文Stiffness and deflection analysis of complex structures（复杂结构的刚度和变形分 析）介绍了这种新的计算方法，从而将矩阵位移法推广推广到求 解弹性力学平面应力问题。它们把平面板壳结构划分为

3、一个个 三角形和矩形的“单元”，利用单元中近似位移函数，求得单元 节点力与结点位移关系的单元刚度矩阵。1960 年, R-W-Clough 在论文The finite element in plane stress analysis（平Finite面应力分析的有限元法 ) 中首次提出了有限单元(Eleme nt )这一术语，他也因此被称为“有限单元之父”二、 有限元法的基本思想有限元法是一种结构分析的方法，正如 O- C - Zienkiewicz所说的：“人类思维的限制在于不能通过一步运算就掌握复杂环境 和事物的行为。 因此，先把所有系统分解为它们的元件或单元， 这些元件的行为已经被充分的了

4、解，再把元件重新组装成原来 的系统来研究系统的行为” 。可以看出有限元法的基本思想是将 连续的求解区域离散为一组由有限个单元组成并按一定方式相 互连接在一起的单元组合体来加以分析。三、 有限单元法的数学基础 当有限单元法成功的应用于求解弹性力学平面问题之后，下一 步要解决的问题就是能否把这种方法应用于求解其他连续介质 问题。在寻找连续介质问题近似算法的时候，数学家们发展了 微分方程的近似解法，包括有限差分方法，变分原理和加权余 量法。四、 有限元分析的基本步骤 建立研究对象的近似模型 将研究对象分割成有限数量的单元 用标准方法对每一个单元提出一个近似解 将所有单元按标准方法组合成一个与原有系统

5、近似的系统 用数值方法求解这个近似系统计算结果处理与结构验证 五、一维杆的有限位移法分析本文以一维直杆的分析为例子，研究有限元位移法基本原理和 求解过程。 虚位移原理推到一维直杆单元的刚度方程如下图所示一维直杆，已知直杆杆长为 L，横截面积为A，材料弹性 模量为E,所受轴向分布载荷集度为q(x)。杆端位移分别记为Ui，w， 杆端力分别记为Si，Sj。q(x)Xa设局部坐标系下杆中A点的坐标为Xa，因为只有两个边界条件Ui ,Uj，因此杆轴任意一点(例如 A点)的位移可假设为式中a，b为待定常数 它们可由杆端位移条件来确定:将式代入式可得:XU = (V L)UiX Ui若引入无量纲变量:则式(

6、3)可改写成：fuilu = Nm + NjUj = Ni Nj=NUe 式中Nj = 1 - Nj =称为形函数，矩阵N称作形函数矩阵；矩阵Ue称为杆端位移矩阵或节 点位移矩阵。由式(4)可以看出，形函数具有如下性质：1、 本端为它端为零叫(0) = 1 叫(0) = 0M(1)= 0 Nj(1)= 12、 任意一点总和为1)叫()现采用虚位移原理给出该杆单元的特性公式，设杆端 i，j分别产生 虚位移Uj，由此引起的单元内任意一点的虚位移为：_ _ Tu 二 N ue = N Uj ujdudNdNidNj1-11 1 Ue-丨|Ue=|一 Uedxdx.dxdxI LL 一二【BiB2 】

7、ue = Bu式中B为应变矩阵。由此可得； = B Ue 又-= EBuEUe根据虚位移原理：对任意虚位移,外力所做的总虚功恒等于变形体所接受的总虚变形功,所以有Lq x udx 二oSTLue 0 q x N uedxT L(S q q x Ndx) ue=(S q q x NTdx)T ueL T L T T T L TW变=o Adx 二 o ueTBtEAB uedx 二 ueT( BTEABdx) ue由 W外=“变可得：(S :q x NTdx)T ueueT( BtEABdx) u(S : q x NTdx)T = u(Bt EABdx)若记L T eq x NTdx = Fe

8、o EL TBt EABdx 二 K e0 eFEe称为该杆单元等效节点载荷；K e局部坐标单元刚度矩阵。所以可得单元刚度方程:S F E 二 K eUe式中单元刚度矩阵的显式为:心宜1 TL 1可见单元刚度矩阵具有对称性。即单元刚度矩阵的每一个元素可写成LEA独 dx0 dx dx将一维直杆离散为三个单元进行分析现考虑下图所示一维直杆：长度为 L,分为三个单元，每个单元长度N J Ng Nj N4对于单元，节点位移分别为ui，U2,对应形函数为N, Nb; 由KjLdN0 dxEAdNidx得:dxKKdATj dNt 血 + 血 axdN2 dNxAE- dx + 血 dx学AE学血+ d

9、x dxdNa dN-AE1 dx + ax dxdN2 dNAE-血 + dx dxdN dN2AE1 血 + dx dxdN4 dN2AE-血 血 dx笄3dx dxdx + 血 axdN4 dNA/IE血 + dx dxdNr dN 一/IE一必 + dxAEdx dx 则 4 dNx AEdx dN* dN AE-dx dx dN2 dN殳 dxAEdx dR 4 dN, 盂ae石dNA dN2dx dx帕晰(krAE dxdN4 dN4盂ae盂dx叫+dx +dx +dx +dx +dx +dx +dx +dx +r dX1 dN-AEdx dx dxdN2 dNAE-dr dx d

10、x 学业学血dx drdg dN-AE-dx dxdN? dN2 盂AE石dN dN?盂 盂dN4 dN2AE-dx 血 血dN. A 冲,AE dx dr dxdg dN*-AE-dx dx dxdN4 dN4血. dr dxdr又对于单元、，形函数N=0;对单元,形函数N2=0;对单元,形函数N3=0;对于单元、，形函数 N=0。因此可得：dN dNx/11= ，4月一必；S =Jiii血 血f dN.dN2= I AEdx +J偽必 必 f dNA 血酹=-rAEdx + 丿血dx axf dN4 dN/44 = -rAE dx .Jq3血 血dNdNrAE dx ; Km = 0 ;

11、K4 = 0 * dbc dx誓E学血；心=/学卫e学血；隔=o q2 dx dx Jq2 ax axf dN 血V3 f dN4 dgI -AE- dx ; = -AE- dx阳必 必 /sb血 必且对于单元，Ni = i h1N2 =对于单元，N2-上h2N3 =对于单元，N3-仝h3N4 a所以式中宀,2, 3分别对应于单元，h1xh2xh3KiidN 八 dN dxEAhiKi2hi EAt dxdx2dx =0 h/EAhiEAhiKi30 , Ki4K22dN2EAdN2dxa dxhi i i hEAdx h h.dxdxi i-)EA-)dx h h.EA EAdN3EAdNd

12、xdx厂丄EA(-丄)=0 h2 h2EAh2K24K33dNo A dN33EA 3匕 dx dxh2 1 1 h3EA dx0 h2 h2 0EA EA + hz h3 dNo A dN3 ,dx 3EA 3dx门3 dx dx1 1( )EA( )dx馆 h3K34dNAEAdx13 dx dxh3 10EA(-丄)二 h3 h37EAh3K44dNEAd%x =3 dx dxh3 EA , dx 二EA所以单元整体刚度矩阵0 h32h3h11hihi心EA01 1+h1 h2h2整体单元刚度矩阵元素h2h2 h3h3Kj的物理意义为：j移时，在节点i上需施加的节点力h3节点在X方向产生

13、单位位又由单元等效节点载荷L T L _ _TF,=打 q（ x）N dx =打 q（ x ） l M NJ dx可得：每一节点等效载荷分别为：/! = / N】q 血 + / TVjq dx + / TVjq dxJ ill J Sis J 口3h= I N2q 血 + / N?q 血 + / Mq dxJ flj J Oj J O3k= I Mq 血 + / yV3q 血 + / Ng dxJ Qi J O2 v Qjf N4q(x)dx N4q(x)dx N4q(x)dx R又对于单元、，形函数N=0;对单元，形函数N2=0;对单元, 形函数N3=0；对于单元、，形函数 N=0。所以有：

14、A = / Mq 血/2 = / N 皇 q 血 + / Ar2q 血h= N：t dx+ f Mq dxJ Qa J figf4 二 N4q(x)dx RO-3以上四式可写为:h= Njq 血 +血+ / 血=於+斤Njq血+ /用q血=於+胃N2 q(x)dx R式中：N，表示单元e的形函数Nk;在本例子中，e=1、2、3,k=1、 2；而且有 X2 X X X1N； (x)=; 逍(x)=he表示单元长度。所以有9-J1 =j:hla I x 2 町 hia /云房 h23a f時一咯 h：i由单元刚度方程可得:EAhi1h1h11 丄hi h21h2h2衍（囲一叱）、 a+隔乃(啜一述

15、) . a+ r工4（诚一错）丄丄h2 h31h3h31h3u/lU23a(X2 7z 2 2、 3 3a /X2(X2 -X1 ) X2 -洛、 R( 2X1(X22 -N2)、 a /X3(X32 -X22) (3X3 X23h/ 3 2 g 2 3 )3 3 (2 2、 ,2 2. 3 3a M3X2 X2(X3X2 )、丄 a ,X4(X4X3 ) X4X3、( )十( )0 3 2 佗 2 ，3 3 22、a /X4 -X3 X3(X4 -X3 R23 3(-h3 3例：假设 A, E, L, a, 和 R 都等于 1。且 g=h2=h3=h，则有 h = 1 , X1 =30, X

16、 2 = 1 / 3, X 3 = 2 / 3, X 4 = 1,_ 2 -1 0_U20,037-1 2 一1血=01)740 -1 1 U4 0.383 由整体单元刚度方程得: 0.0 0.494U30.951_4_L333对于单元e（e二、）可写出其位移函数:应变为:殂=BN； 1 1.dxdxh h由丁二B护得每个单元的应力:b=1*48er2 = 1.37cf3= L15下图1和图2分别表示有限元解和精确解的比较:图1位移对比图2应力对比二、分析与计算（40分）1、图示两个结构和单元相似，方位相同的平面应力有限元模型，两模型的单元 厚度和材料相同。两个模型右端单元边上受均匀剪切面力。

17、对于下列 2种情况,试根据有限元法和力学有关知识来分析两个模型求解后对应节点的位移值和对应单元的应力值之间的关系：1）两个模型面力的合力相等；2）两个模型面力的 集度相等。（10分）解：建立坐标系如图所示，对（a）图，各节点坐标点号x坐标y坐标点号x坐标y坐标1400420202402050032006020单元节点信息数组可记为653346431124单元，面积 A1 =0.5 20 20 =200，由式ai =Xjym - Xmyjbi =yj-ym i, j,m （a）Ci =Xm - Xj可得= 0 & = 20b2 = -20 c2 = -20b3 = 20 c3 = 0由式12A0

18、1i,j,m (b).单元几何矩阵为00-200200 1J2 00-1 0100200-200010 -100200-20-20020 一0-1 -10i弹性矩阵大为简化,可设Cibi卩=0且为单位厚度,B 1 =400为了计算简便,0012j，可得D = E 00 01 00 0.500-2002001E?0-202010200-200010 一一 40I020-200d也00-10-100I10-1-10e= Se6e（c）,得单元的应力矩阵= tABeTDBe( d)，s1亠400t e由式 Ke =t B DB d.1Q单元的单元刚度矩阵为K二tABTS根据单元刚度矩阵的性质可得10

19、-1-101_101 120-200031-2_1-2130-10-20200_ 1_101 一二K二K二Kk(2)对（b）各节点坐标点号x坐标y坐标点号x坐标y坐标1200410102201050031006010单元节点信息数组可记为：653346431124单元，面积A二0.5 10 10=50，由式（a）可计算出bj = 0 G = 10匕2 = -10 C2=-10 b3 = 10 C3 = 0由式（c）由式（b）,单元几何矩阵为00-100100 1100-101010100-100010 一010-100100-10-1000-1-10得单元的应力矩阵00-1001001E 00

20、-20200100-100020-2005 一20.50-5-50I10-1-10S1丄100同理由单元刚度矩阵的性质可得由式（d）得单元的单元刚度矩阵-10-1-101020-200E1031-2141-2130100-202010-1-101二tABTSK(2)占北综上可知，两个模型中的单元刚度矩阵均相同， 所以它们的总刚度矩阵也相当两个模型面力的合力相等,同，即30-1-1-20010000030-2-1-1100000-103100-2-10000-1-213000-10000-2-10061-2-1-20010-10016-1-4-1-11001-20-2-16100-2-110-1

21、-1-1-416000-10000-2-10031-1000000-10013-1-2000001-20-1-130000010-1-10-203 i00它们的载荷列阵都为0 0 0 0 U5 V5 U6 V6位移列阵为- luViU2U3 V3位移约束条件为u5形成整体平衡方程 二V5二5二V6 =0，将此约束条件引入整体刚度方程，对30-1-1-200100 005 0 1030-2-1-11000 C0v1P/2-103100-2-100 00U20-1-213000-100 00V2P/2-2-10061-2-100 00U300-10016-1-400 00V3001-20-2-16

22、100 00U40110-1-1-141600 00V4000000000100 00000000000010 00000000000001 0000000000000 1一0 _ 一由上式可知在这种情况下,两种模型求解后对应节点的位移相等。有整体节点位移获取单元节点位移，所以对应的单元节点位移也相等。以单元为例，两种模型应力矩阵的关系s（b）=2s；a），又由可=six可得，模型（b）中单元的应力是模型（a）的2倍，其它单元可得到类似的结论。当两个模型面力的集度相等，可设模型（a）右端受剪力的合力为2P，模型（b）右端受剪力的合力为P,则两种模型的载荷列阵分别为F（a）j0 POPOOOOU

23、5V5U6VJ0 0 0 0 0 U5 V5 U6 V6由上述内容可得F（a）=2F（b），进而可知（a）模型中对应节点的位移是（b） 模型的2倍在单元中， 醐=2S；a），所以由2 = S】3可知模型（a）和模型（b）中单元的应力相等，其他单元可类似说明2、证明平面问题三节点三角形单元发生刚体位移（小位移平动和转动）时，单元中将不产生应力。（10分）证明：三节点三角形单元位移模式选取一次多项式:(1)u = 1 2XV = 5X 6 讨在（1）的1式中带入节点i的坐标（x，yi）得到节点i在x向的位移ui，同理可得，u 二爲 ：2人 :3yi(2)比八1 T yjUm = 1 ：2Xm ：3ym解（2）式可得到广义坐标由节点位移表示的表达式:11 =2A（aiU 二丄（bu2 Ar（Cu2A1 / （aiv2A1（bv2A（C V2 AajUj