弹性力学与有限元分析-第四章 平面问题有限元分析及程序设计
弹性力学与有限元完整版
• 第一章 弹性力学基本方程
1.1 绪论 1.2 弹性力学的基本假定 1.3 几个基本概念 1.4 弹性力学基本方程
• 第二章 弹性力学平面问题
2.1 平面应力问题 2.2 平面应变问题 2.3 平面问题的基本方程
• 第三章 弹性力学问题求解方法简述
• 第一章 弹性力学基本方程
1.1 绪论 1.2 弹性力学的基本假定 1.3 几个基本概念 1.4 弹性力学基本方程
3、平面应力问题应力、应变
• 应力分量
x、 y、 xy
• 应变分量
z 0 yz = zx 0 x、 y、 xy
x
{} y
xy
x
y
xy
2.2 平面应变问题
1 平面应变问题的概念
– 弹性体是具有很长的纵向轴的柱形物体，横截面大 小和形状沿轴线长度不变；作用外力与纵向轴垂直， 并且沿长度不变；柱体的两端受固定约束。
应力分量——6个
x、 y、 z、 xy、 yz、 zx
应变分量——6个
x、 y、 z、 xy、yz、 zx
位移分量——3个
u、v、w
合计 15
• 第二章 弹性力学平面问题
2.1 平面应力问题 2.2 平面应变问题 2.3 平面问题的基本方程
2.1 平面应力问题
1、平面应力问题的概念
平面应力问题讨论的弹性 体为薄板。薄壁厚度远小于 结构另外两个方向的尺度。 薄板的中面为平面，其所受 外力，包括体力均平行于中 面O-xy面内，并沿厚度方向 z不变。而且薄板的两个表 面不受外力作用。
2、平面应变问题的位移
• 沿纵向轴的位移恒等于零； • 由于无限长，所以任一个横截面都是一样的，与z
轴无关。
有限元经典PPT第4章
Pii Kiiui
Ki1u1 Ki2u2 Kiiui K u i,i1 i1
ui
n
Kiiui Kiiui
Kiju j
4.1.2 平面应力问题有限元的基本思想和瑞雷－里兹法
v3 f3y
3
u3
f3x
f1y v1 u1
1 f1x
v2 f2y u2
2 f2x
给定一个三角形单元和作用在角点上 的六个力，要求得六个角点的位移。 或者是要求三角形角点发生指定的位 移，在三角形三个角点如何加力？
很显然，问题的精确解很困难。采用 瑞雷－里兹法求近似式解
e号单元的三个节点I,j,k的力对应的 力的平衡方程是第2i-1,2i;2j-1,2j;2k1,2k个平衡方程
e号单元的三个节点I,j,k的位移是第 2i-1,2i;2j-1,2j;2k-1,2k个未知数
弹性模量：E 横截面积：A
1
1 L
2
2L
3
局部系单元刚度阵：
k
1
EA L
1 -1
-1
1
2 集成总刚：
0 1
解得：
ux uy
L EA
3.8284L
EA
i
j
第一类位移条件：
Ki1u1 Ki2u2 Kiiui Ki1ui1
ui 0
令： Kij 0 i j
m
vi 0
Kii 1
um 0
Pi 0
ui 0
第二类位移条件：um um
大数
充大数法： Kii Kii
第一步：求转换矩阵
k2
EA 1 2L -1
-1
1
P
cos 0
T sin
有限元程序求解弹性力学平面问题
计算力学课程设计报告有限元程序求解弹性力学平面问题专业：班级：姓名：学号：指导教师：有限元程序求解弹性力学平面问题设计目的:1、学习有限元程序求解弹性力学平面问题的方法；2、学习有限元程序编写技巧；3、加深对有限元方法的理解；4、锻炼处理复杂弹性力学问题的能力。
题一：例3.9设深梁承受均布荷载，如下图（a）所示。
假定E=1，泊松比=0.17μ，不计容重，厚度t=1m，为平面应力问题。
因对称去半边结构进行计算，结构支承、单元划分、节点编号如图（b）所示。
试画出y=0及y=6m截面的竖向位移图，x=3m截面的xσ应力分布图。
50kN100kN100kN50kNy1、有限元Fortran源程序如下：COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),*S(3,6),TKZ(200,20),EKE(6,6),P(200)CALL DATAIF(IND.EQ.0)GOTO 10EO=EO/(1.0-UN*UN)UN=UN/(1.0-UN)10 CALL TOTSTICALL LOADCALL SUPPORCALL SOLVEQCALL STRESSSTOPENDSUBROUTINE DATACOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3), *S(3,6),TKZ(200,20),EKE(6,6),P(200)PRINT*,"INPUT:NJ,NE,NZ,NDD,NPJ,IND"READ(5,*)NJ,NE,NZ,NDD,NPJ,INDNJ2=NJ*2NPJ1=NPJ+1PRINT*,"INPUT:EO,UN,GAMA,TE"READ(5,*)EO,UN,GAMA,TEPRINT*,"INPUT:JM"READ(5,*)((JM(I,J),J=1,3),I=1,NE)PRINT*,"INPUT:CJZ"READ(5,*)((CJZ(I,J),J=1,2),I=1,NJ)PRINT*,"INPUT:NZC"READ(5,*)(NZC(I),I=1,NZ)PRINT*,"INPUT:PJ"READ(5,*)((PJ(I,J),J=1,2),I=1,NPJ1)OPEN(100,FILE='1.TXT')WRITE(100,10)(I,(CJZ(I,J),J=1,2),I=1,NJ)10 FORMAT(4X,2HNO,6X,1HX,6X,1HY/(I6,2X,F7.2,F7.2))RETURNENDSUBROUTINE ELEST(MEO,IASK)COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3), *S(3,6),TKZ(200,20),EKE(6,6),P(200)IE=JM(MEO,1)JE=JM(MEO,2)ME=JM(MEO,3)CM=CJZ(JE,1)-CJZ(IE,1)BM=CJZ(IE,2)-CJZ(JE,2)CJ=CJZ(IE,1)-CJZ(ME,1)BJ=CJZ(ME,2)-CJZ(IE,2)AE=(BJ*CM-BM*CJ)/2.0IF(IASK.LE.1) GOTO 50DO 10 I=1,3DO 10 J=1,6B(I,J)=0.010 CONTINUEB(1,1)=-BJ-BMB(1,3)=BJB(1,5)=BMB(2,2)=-CJ-CMB(2,4)=CJB(2,6)=CMB(3,1)=B(2,2)B(3,2)=B(1,1)B(3,3)=B(2,4)B(3,4)=B(1,3)B(3,5)=B(2,6)B(3,6)=B(1,5)DO 20 I=1,3DO 20 J=1,6B(I,J)=B(I,J)/(2.0*AE)20 CONTINUED(1,1)=EO/(1.0-UN*UN)D(1,2)=EO*UN/(1.0-UN*UN)D(2,1)=D(1,2)D(2,2)=D(1,1)D(1,3)=0.0D(2,3)=0.0D(3,1)=0.0D(3,2)=0.0D(3,3)=EO/(2.0*(1.0+UN))DO 30 I=1,3DO 30 J=1,6S(I,J)=0.0DO 30 K=1,3S(I,J)=S(I,J)+D(I,K)*B(K,J)30 CONTINUEIF(IASK.LE.2) GOTO 50DO 40 I=1,6DO 40 J=1,6DO 40 K=1,3EKE(I,J)=EKE(I,J)+S(K,I)*B(K,J)*AE*TE40 CONTINUE50 CONTINUERETURNENDSUBROUTINE TOTSTICOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3), *S(3,6),TKZ(200,20),EKE(6,6),P(200)DO 20 I=1,NJ2DO 20 J=1,NDDTKZ(I,J)=0.020 CONTINUEDO 30 MEO=1,NECALL ELEST(MEO,3)DO 30 I=1,3DO 30 II=1,2LH=2*(I-1)+IILDH=2*(JM(MEO,I)-1)+IIDO 30 J=1,3DO 30 JJ=1,2L=2*(J-1)+JJLZ=2*(JM(MEO,J)-1)+JJLD=LZ-LDH+1IF(LD.LE.0) GOTO 30TKZ(LDH,LD)=TKZ(LDH,LD)+EKE(LH,L)30 CONTINUERETURNENDSUBROUTINE LOADCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3), *S(3,6),TKZ(200,20),EKE(6,6),P(200)DO 10 I=1,NJ2P(I)=0.010 CONTINUEIF(NPJ.EQ.0) GOTO 30DO 20 I=1,NPJI1=I+1J=IFIX(PJ(I1,2))20 CONTINUE30 IF(GAMA.LE.0.0) GOTO 50DO 40 MEO=1,NECALL ELEST(MEO,1)PE=-GAMA*AE*TE/3.0IE=JM(MEO,1)JE=JM(MEO,2)ME=JM(MEO,3)P(2*IE)=P(2*IE)+PEP(2*JE)=P(2*JE)+PEP(2*ME)=P(2*ME)+PE40 CONTINUE50 CONTINUERETURNENDSUBROUTINE SUPPORCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3), *S(3,6),TKZ(200,20),EKE(6,6),P(200)DO 60 I=1,NZMZ=NZC(I)TKZ(MZ,1)=1.0DO 10 J=2,NDDTKZ(MZ,J)=0.010 CONTINUEIF(MZ-NDD)20,20,3020 JO=MZGOTO 4030 JO=NDD40 DO 50 J=2,JOJ1=MZ-JTKZ(J1+1,J)=0.050 CONTINUEP(MZ)=0.060 CONTINUERETURNENDSUBROUTINE SOLVEQCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3), *S(3,6),TKZ(200,20),EKE(6,6),P(200)DO 50 K=1,NJ1IF(NJ2-K-NDD+1)10,10,2010 IM=NJ2GOTO 3020 IM=K+NDD-130 K1=K+1DO 50 I=K1,IML=I-K+1C=TKZ(K,L)/TKZ(K,1)LD1=NDD-L+1DO 40 J=1,LD1M=J+I-KTKZ(I,J)=TKZ(I,J)-C*TKZ(K,M)40 CONTINUEP(I)=P(I)-C*P(K)50 CONTINUEP(NJ2)=P(NJ2)/TKZ(NJ2,1)DO 100 I1=1,NJ1I=NJ2-I1IF(NDD-NJ2+I-1)60,60,7060 JO=NDDGOTO 8070 JO=NJ2-I+180 DO 90 J=2,JOLH=J+I-1P(I)=P(I)-TKZ(I,J)*P(LH)90 CONTINUEP(I)=P(I)/TKZ(I,1)100 CONTINUEOPEN(101,FILE='2.TXT')WRITE(101,110)(I,P(2*I-1),P(2*I),I=1,NJ)110 FORMAT(2X,3HJD=,5X,2HU=,12X,2HV=/(I4,3X,F11.3,2X,F11.3)) RETURNENDSUBROUTINE STRESSCOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AE COMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),*S(3,6),TKZ(200,20),EKE(6,6),P(200)DIMENSION WY(6),YL(3)DO 60 MEO=1,NECALL ELEST(MEO,2)DO 10 I=1,3LH=2*(I-1)+JLDH=2*(JM(MEO,I)-1)+JWY(LH)=P(LDH)10 CONTINUEDO 20 I=1,3YL(I)=0.0DO 20 J=1,6YL(I)=YL(I)+S(I,J)*WY(J)20 CONTINUESIGX=YL(1)SIGY=YL(2)TOXY=YL(3)PYL=(SIGX+SIGY)/2.0SIG=(SIGX-SIGY)**2/4.0+TOXY*TOXYRYL=SQRT(SIG)SIG1=PYL+RYLSIG2=PYL-RYLIF(SIGY.EQ.SIG2) GOTO 30CETA1=TOXY/(SIGY-SIG2)CETA=90.0-57.29578*ATAN(CETA1)GOTO 4030 CETA=0.0OPEN(102,FILE='3.TXT')40 WRITE(102,50)MEO,SIGX,SIGY,TOXY,SIG1,SIG2,CETA50 FORMAT(4X,2HE=,I3/2X,3HSX=,F11.3,3X,3HSY=,F11.3,3X,4HTAU=,*F11.3/2X,3HS1=,F11.3,3X,3HS2=,F11.3,3X,4HCET=,F11.3)60 CONTINUERETURNEND2、输入数据：28,36,9,10,4,0 1,0.17,0,1 1,5,22,5,62,6,33,6,73,7,44,7,85,9,66,9,106,10,7 7,10,11 7,11,8 8,11,12 9,13,10 10,13,14 10,14,11 11,14,15 11,15,12 12,15,16 13,17,14 14,17,1814,18,15 15,18,19 15,19,16 16,19,20 17,21,18 18,21,22 18,22,19 19,22,23 19,23,20 20,23,24 21,25,22 22,25,26 22,26,23 23,26,27 23,27,24 24,27,28 0,61,62,63,60,51,52,53,50,4 1,42,43,40,31,32,33,30,21,22,23,20,11,12,13,10,01,02,03,07,15,23,31,39,47,49,50,55 0,0-50000,2-100000,4-100000,6-50000,83、输出结果：⑴、节点坐标：NO X Y1 0.00 6.002 1.00 6.003 2.00 6.004 3.00 6.005 0.00 5.006 1.00 5.007 2.00 5.008 3.00 5.009 0.00 4.0010 1.00 4.0011 2.00 4.0012 3.00 4.0013 0.00 3.0014 1.00 3.0015 2.00 3.0016 3.00 3.0017 0.00 2.0018 1.00 2.0019 2.00 2.0020 3.00 2.0021 0.00 1.0022 1.00 1.0023 2.00 1.0024 3.00 1.0025 0.00 0.0026 1.00 0.0027 2.00 0.0028 3.00 0.00⑵、节点位移（m）：JD= U= V=1 -29767.031 -1173919.8752 -14003.248 -1174020.7503 -3753.289 -1179520.1254 0.000 -1181721.8755 -26382.504 -1072683.5006 -10746.987 -1073616.8757 -2064.578 -1082362.6258 0.000 -1085875.3759 -13536.982 -964012.00010 3372.830 -970057.00011 7268.439 -989270.93812 0.000 -998403.81213 7816.650 -835385.12514 27176.352 -861715.688 15 22063.291 -905727.93816 0.000 -927167.25017 29514.604 -665604.25018 53419.770 -747342.18819 34876.914 -839808.62520 0.000 -881221.00021 29580.316 -416289.65622 52945.020 -632602.68823 17504.215 -803767.50024 0.000 -859483.75025 0.000 0.00026 -120103.062 -583506.75027 -76202.484 -787348.93828 0.000 -829172.562⑶、单元应力（Pa）：E= 1SX= -1489.449 SY=-101489.578 TAU= -1489.493 S1= -1467.270 S2=-101511.766 CET= 179.147 E= 2SX= -1475.786 SY=-100654.781 TAU= -1790.448 S1= -1443.477 S2=-100687.094 CET= 178.966 E= 3SX= -7021.617 SY=-101597.578 TAU= -3741.747 S1= -6873.812 S2=-101745.383 CET= 177.738 E= 4SX= -8067.512 SY= -98528.883 TAU= -4459.167 S1= -7848.234 S2= -98748.156 CET= 177.185E= 5SX= -13143.324 SY= -99391.883 TAU= -1662.572 S1= -13111.285 S2= -99423.922 CET= 178.896E= 6SX= -14652.791 SY= -98337.594 TAU= -1501.185 S1= -14625.871 S2= -98364.516 CET= 178.973E= 7SX= -2923.120 SY=-109168.492 TAU= -5888.416 S1= -2597.762 S2=-109493.844 CET= 176.837 E= 8SX= -716.068 SY=-103681.609 TAU= -8617.434 S1= 0.160 S2=-104397.844 CET= 175.249 E= 9SX= -9188.308 SY=-105121.906 TAU= -9771.604 S1= -8203.109 S2=-106107.102 CET= 174.243E= 10SX= -12285.007 SY= -95180.242 TAU= -12199.550S1= -10526.906 S2= -96938.344 CET= 171.799 E= 11SX= -14170.538 SY= -95500.742 TAU= -5489.623 S1= -13801.672 S2= -95869.609 CET= 176.156 E= 12SX= -22797.447 SY= -91347.070 TAU= -3902.925 S1= -22575.945 S2= -91568.570 CET= 176.752 E= 13SX= -5104.274 SY=-129494.578 TAU= -11708.809 S1= -4011.730 S2=-130587.125 CET= 174.669 E= 14SX= 969.693 SY=-108176.422 TAU= -21424.820 S1= 5024.629 S2=-112231.359 CET= 169.283 E= 15SX= -14954.604 SY=-110883.578 TAU= -18383.518 S1= -11552.312 S2=-114285.867 CET= 169.515 E= 16SX= -19890.213 SY= -86924.367 TAU= -25131.225S1= -11514.891 S2= -95299.688 CET= 161.569 E= 17SX= -22109.713 SY= -87301.680 TAU= -10225.519S1= -20543.445 S2= -88867.945 CET= 171.292 E= 18SX= -35190.500 SY= -77218.953 TAU= -9162.105 S1= -33280.023 S2= -79129.430 CET= 168.221 E= 19SX= -9785.861 SY=-171444.500 TAU= -20525.008 S1= -7220.609 S2=-174009.750 CET= 172.876 E= 20SX= 4594.449 SY=-113592.445 TAU= -46145.883 S1= 20477.516 S2=-129475.516 CET= 161.007 E= 21SX= -25287.357 SY=-118672.375 TAU= -30023.787 S1= -16467.543 S2=-127492.188 CET= 163.629 E= 22SX= -30634.465 SY= -71127.195 TAU= -44991.484 S1= -1543.734 S2=-100217.922 CET= 147.114 E= 23SX= -34259.684 SY= -71743.445 TAU= -14638.012 S1= -29220.699 S2= -76782.422 CET= 161.005 E= 24SX= -43958.176 SY= -53419.129 TAU= -17697.580S1= -30369.762 S2= -67007.547 CET= 142.483 E= 25SX= -19028.232 SY=-252549.391 TAU= -34958.820 S1= -13907.102 S2=-257670.531 CET= 171.666 E= 26SX= 3973.836 SY=-114063.977 TAU= -92238.578 S1= 54459.203 S2=-164549.344 CET= 151.307 E= 27SX= -39180.902 SY=-121400.266 TAU= -39312.672 S1= -23409.199 S2=-137171.969 CET= 158.140 E= 28SX= -42804.859 SY= -43317.969 TAU= -65723.141S1= 22662.227 S2=-108785.055 CET= 135.112 E= 29SX= -42224.188 SY= -43219.219 TAU= -10273.362S1= -32436.301 S2= -53007.105 CET= 136.386 E= 30SX= -21830.449 SY= -25448.408 TAU= -23810.377S1= 239.566 S2= -47518.426 CET= 137.172 E= 31SX= -48815.301 SY=-424588.281 TAU= -79800.312 S1= -32570.906 S2=-440832.688 CET= 168.494 E= 32SX=-132272.047 SY= -71582.141 TAU=-175409.688 S1= 76088.000 S2=-279942.188 CET= 130.093 E= 33SX= -45090.219 SY= -56761.309 TAU= 804.799 S1= -45034.984 S2= -56816.547 CET= 3.926 E= 34SX= 42332.844 SY= -9222.024 TAU= -47066.441 S1= 70218.484 S2= -37107.668 CET= 149.354 E= 35SX= -20899.361 SY= -19971.461 TAU= 16235.217 S1= -4193.565 S2= -36677.254 CET= 45.818 E= 36SX= 73164.031 SY= -17873.285 TAU= -17873.346 S1= 76547.367 S2= -21256.619 CET= 169.2814、数据处理：由以上运行结果，并结合对称特点可得y=0及 y=6面上节点竖向位移值如下表：根据上表位移值，可作出y=0及 y=6截面的竖向位移图如下：整理x=3m 处各点的x σ值时，图（b ）中的b 、c 、d 、e 、f 、g 六点按二单元平均法整理；a 、h 两点为边界点，在a 点附近应力x σ变化显著，用四点插值公式计算a 点的应力；而h 点用三点插值公式计算应力。
有限元分析——平面问题
Re=
NT
s
Pstds
江西五十铃发动机有限公司
技术中心 12 /33
4、整体分析 整体刚度矩阵 整体刚度矩阵组装的基本步骤：
先求出各个单元的单元刚度矩阵； 将单元刚度矩阵中的每个子块放在整体刚度矩阵中的对应位置上，得到单 元的扩大刚度矩阵； 将全部单元的扩大矩阵相加得到整体刚度矩阵。
不失一般性，仅考虑模型中有四个单元，如图所示，四个单元的整体节点位 移列阵为
τZX z= + t/2 =0
因板很薄，载荷又不沿厚度变化，应力沿板 的厚度方向是连续分布的，可以认为，在整
Z
个板内各点都有
σZ=0 τYZ=0 τZX=0
O
tX
图1 平面应力问题
根据剪应力的互等性、物理方程，可得描述平面应力问题的八个独立的基本变量 为
江西五十铃发动机有限公司
技术中心 4 /33
σ=[σX σY τXY]T ε=[εX εY γXY]T
x2 y2 ɑ1= x 3 y 3
1 y2 b1=- 1 y 3
1 c1= 1
x2 x3
(1,2,3)
上式表示下标轮换，即1 2，2 3，3 1同时更换。
江西五十铃发动机有限公司
技术中心 9 /33
重写位移函数，并以节点位移的形式进行表达，有
uv((xx,,yy))N（x,y）qe
其中形函数矩阵为
Y
江西五十铃发动机有限公司
图2 平面应变问题
技术中心 5 /33
根据几何方程、物理方程可得，描述平面应变问题的独立变量也是八个，且与 平面应力问题的一样。只是弹性矩阵变为
1
D=
E1
1 1 2 1
1
第4章 平面问题的有限元法-4收敛准则
8
9 10 11 12 13 14
2
4
6
8 10 12 14
(a)
(b)
图4-13
四. 单元节点i、j、m的次序 在前面章节中，我们曾指出，为了在计算中保证单元的 面积 不会出现负值，节点i、j、m的编号次序必须是逆时 针方向。事实上，节点i、j、m的编号次序是可以任意安排 的，只要在计算刚度矩阵的各元素时，对取绝对值，即可 得到正确的计算结果。在实际计算时，应该注意所选有限元 分析软件的使用要求。 五. 边界条件的处理及整体刚度矩阵的修正 在前面讨论整体刚度矩阵时，已经提到，整体刚度矩阵 的奇异性可以提高考虑边界约束条件来排除弹性体的刚体位 移，以达到求解的目的。
B =2(d+1)
若采取带宽压缩存储，则整体刚度矩阵的存储量N 最 多为N =2nB = 4n(d+1) 其中：d为相邻节点的最大差值，n为节点总数。 例如在图4-13中，(a)与(b)的单元划分相同，且节点 总数都等于14，但两者的节点编号方式却完全不同。(a) 是按长边进行编号， d =7， N =488；而(b)是按短边进行 编号，d =2，N =168。显然(b)的编号方式可比(a)的编号 方式节省280个存储单元。
为了保证解答的收敛性，要求位移模式必须满足以下三 个条件，即 ⑴ 位移模式必须包含单元的刚体位移。也就是说，当 节点位移是由某个刚体位移所引起时，弹性体内将不会产生 应变。所以，位移模式不但要具有描述单元本身形变的能力 ，而且还要具有描述由于其它单元形变而通过节点位移引起 单元刚体位移的能力。 例如，三角形三节点单元位移模式中，常数项1、4 就 是用于提供刚体位移的。 ⑵ 位移模式必须能包含单元的常应变。每个单元的应变 一般都是包含着两个部分：一部分是与该单元中各点的坐标 位置有关的应变（即所谓各点的变应变）；另一部分是与位 置坐标无关的应变（即所谓的常应变）。从物理意义上看，
有限元讲义弹性力学平面问题有限单元法(四结点四边形等参元,八结点曲线四边形等参元,问题补充)分析
2.6 四结点四边形单元(The four-node quadrilateral element)前面介绍了四结点的矩形单元其位移函数:xy y x U 4321αααα+++=xy y x V8765αααα+++=为双线性函数，应力，应变在单元内呈线性变化，比常应力三角形单元精度高。
但它对边界要求严格。
本节介绍的四结点四边形等参元，它不但具有较高的精度，而且其网格划分也不受边界的影响。
对任意四边形单元（图见下面）若仍直接采用前面矩形单元的位移函数，在边界上它便不再是线性的（因边界不与x,y 轴一致），这样会使得相邻两单元在公共边界上的位移可能会出现不连续现象（非协调元），而使收敛性受到影响。
可以验证，利用坐标变换就能解决这个问题，即可以通过坐标变换将整体坐标中的四边形（图a ）变换成在局部坐标系中与四边形方向无关的边长为2的正方形。
正方形四个结点i,j,m,p 按反时钟顺序对应四边形的四个结点i j m p 。
正方形的 1-=η 和 1=η 二条边界，分别对应四边形的i ，j 边界和p,m 边界；ξ=-1和ξ=+1分别对应四边形的i ，p 边界和j ，m 边界。
如果用二组直线等分四边形的四个边界线段，使四边形绘成一个非正交网格，那么该非正交网格在正方形上对应着一个等距离的规则网格（见图a, b ）。
当然, 局部坐标上的A 点与整体坐标的A 点对应。
一、四结点四边形等参单元的形函数及坐标变换由于可以将整体坐标下的四边形单元变换成局部坐标下的正方形单元，对于这种正方形单元，自然仍取形函数为： ξηαηαξαα2321+++=U ξηαηαξαα8765+++=V引入边界条件，即可得位移函数：∑=ijmpi i U N Ui ijmpi V N V ∑==写成矩阵形式：{}{}[]{}ee p i p i ed N d N N N N V U f =⎥⎦⎤⎢⎣⎡=⎭⎬⎫⎩⎨⎧=000 式中形函数: ()()()ηηξξηξi i i N ++=1141, ()p m j i ,,, 按照等参元的定义,我们将坐标变换式亦取为: p p m m j j i i i ijmpi x N x N x N x N x N x +++==∑p p m m j j i i i ijmpi y N y N y N y N y N y +++==∑ ()162-- 式中形函数N 与位移函数中的完全一致。
现代设计方法4-1弹性力学平面问题的基本方程
些概念和方程，作为弹性力学有限单元法
的预备知识。
弹性力学—区别与联系—材料力学
1、研究的内容：基本上没有什么区别。
弹性力学也是研究弹性体在外力作用下的平衡和运动，
以及由此产生的应力和变形。
2、研究的对象：有相同也有区别。
材料力学基本上只研究杆、梁、柱、轴等杆状构件，即 长度远大于宽度和厚度的构件。弹性力学虽然也研究杆 状构件，但还研究材料力学无法研究的板与壳及其它实 体结构，即两个尺寸远大于第三个尺寸，或三个尺寸相 当的构件。
y
x
某一个截面上的外法线方向是沿坐标轴的正方向，这个截面称 正面，面上的应力沿正向为正，负方向为负。相反如果某截 面上的外法线是沿坐标轴的负方向，截面为负面，面上的应 力以沿坐标轴负向为正，正向为负。 空间问题有九个应力分量：三个正应力&六个剪应力三个独立
剪应力：txy = tyx 线应变：ex
s De
e x e e y g xy
1 E D 2 1 0
1 0
0 0 1 2
[D]平面应力问题的弹性矩阵，对称与E，有关
(3)平面应变问题的物理方程
角应变：gxy
tyz = tzy tzx = txz
空间应力状态有6个独立应力分量，对应6个应变分量：
ey
ez
gyz gzx
完全弹性,各向同性物体应变与应力关系(胡克定律导出) 胡克定律：在单向应力 状态下，处于弹性阶段
1 e x E s x (s y s z ) 1 e y s y (s x s z ) E 1 e z s z (s x s y ) E 1 1 1 g xy G t xy , g yz G t yz , g zx G t zx
有限元分析第四章
19
4）形函数的性质
形函数是有限单元法中的一个重要函数，它具 有以下性质： 性质1 形函数Ni在节点i上的值等于1，在其它节点 上的值等于0。对于本单元，有
20
Ni ( xi , yi ) 1 Ni ( x j , y j ) 0 Ni ( xm , ym ) 0
（i、j、m）
利用 N i 1 (ai bi x ci y )和ai、bi、ci公式证明 2A
对于一个具体问题进行分析，不管采用什么样的单元， 分析过程与思路是一样的，所不同的只是各种单元的位移模 式和单元刚度矩阵不一样，其他的包括整体刚度矩阵的组装 过程都完全一样，所以我们仅仅对矩形单元位移模式的求取 和单元刚度矩阵的求解加以介绍。
4.7 收敛准则
可以证明，对于一个给定的位移模式，其刚度系统的数 值要比精确值大。所以，在给定载荷的作用下，有限元计算 模型的变形要比实际结构的变形小。因而，当单元网格分得 越来越细时，位移的近似解将由下方收敛于精确解，即得到 真实解的下界。 为了保证解答的收敛性，要求选取的位移模式必须满足 以下三个条件： 1）位移模式必须包含单元的刚体位移 也就是说，当节点位移是某个刚体位移所引起时，弹 性体内将不会产生应变。所以位移模式不但要具有描述单元 本身形变的能力，而且还要具有描述由其他变形而通过节点 位移引起单元刚体位移的能力。例如，三角形三节点位移模 式中，常数项就是用于提供刚体位移的。
Ni(x、y)
1 i(xi，yi) x xi
x xi N i ( x, y ) 1 x j xi
N m ( x, y ) 0
证
N
y j (xj，yj)
m (xm，ym)
xj
x
N i ( x, y )
有限元分析第4章 平面问题有限单元法1
6
P
3
4 5
4
2
位移协调条件：各单元共享节点的位移相等 节点平衡条件：各节点单元内力与节点外力构成平衡力系
最终数学模型： K Q
基本概念
单元（element） 节点 (node)
回顾
单元节点位移 (node displacement)
单元节点内力 (node force)
单元刚度矩阵 (element stiffness matrix)
e
bx u by v
d
S
e p
px u py v dS
代入
u v
N
e
{} [B]{ }e
{ } [S]{ }e
得
内力虚功＝
e x x y y xy xy d
T d
cj
y)v j
(am
bmx
cm y)vm ]
二、平面问题三角形单元分析
三角形单元形函数
形函数
u x,
y
1 2A
[(ai
bi x
ci
y)ui
(a j
bj x
cj
y)u j
(am
bm x
cm
y)um ]
v x,
y
1 2A
[(ai
bi x
ci
y)vi
(a j
插值系数的确定：待定系数法
ui a1 a2 xi a3 yi u j a1 a2 x j a3 y j um a1 a2 xm a3 ym
第4章 平面问题的有限元法-1离散化
e
T i
T j
T T m
u
i
vi
T
uj
vj
um
vm
T
(4-7)
其中的子矩阵
i ui
vi
（i，j，m 轮换） (a)
式中 ui、vi 是节点i在x轴和y轴方向的位移。
在有限单元法中，虽然是用离散化模型来代替原来 的连续体，但每一个单元体仍是一个弹性体，所以在其 内部依然是符合弹性力学基本假设的，弹性力学的基本 方程在每个单元内部同样适用。 从弹性力学平面问题的解析解法中可知，如果弹性 体内的位移分量函数已知，则应变分量和应力分量也就 确定了。但是，如果只知道弹性体中某几个点的位移分 量的值，那么就不能直接求得应变分量和应力分量。因 此，在进行有限元分析时，必须先假定一个位移模式。 由于在弹性体内，各点的位移变化情况非常复杂，很难 在整个弹性体内选取一个恰当的位移函数来表示位移的 复杂变化，但是如果将整个区域分割成许多小单元，那 么在每个单元的局部范围内就可以采用比较简单的函数 来近似地表示单元的真实位移，将各单元的位移式连接
结构离散化后，要用单元内节点的位移通过插值(?)来 获得单元内各点的位移。在有限元法中，通常都是假定 单元的位移模式是多项式，一般来说，单元位移多项式 的项数应与单元的自由度数相等。它的阶数至少包含常 数项和一次项。至于高次项要选取多少项，则应视单元 的类型而定。 （4-1） f N e
(c)
由(c)式左边的三个方程可以求得
1
1 uj 2 um ui xi xj xm 1 y j ,2 1 uj 2 ym 1 um yi 1 ui 1 y j , 1 1 xj 2 ym 1 xm yi 1 xi ui uj um
弹性力学与有限元分析.ppt
上式建立了单元中任意一点的位移与节点位移的关系，
即通过单元节点位移 e 插值求出单元中任一点位移
f (x, y)，把位移函数的这种描述形式称为插值函数形
式。 形函数具有以下两个性质： 1、形函数 N i在节点 i处的值为1，而在其余两个节点 处的值为0。
2、在单元中任意一点，3个形函数之和为1，即：
差太大，即单元划分中不应出现过大的钝角或过 小的锐角，否则，计算误差较大。 在应力较大和应力集中的区域，单元应划分细一 些，以提高精度。 如果边界上有集中力作用，则该点应被划分为点。
单元的大小和数目应根据精度要求来确定，在保证
精度的前提下，力求采用较少的单元。
当物体的厚度有突变或物体由不同材料组成时，不 要把厚度不同或材料不同的区域划分在统一单元。
x y xy
且它们只是
x, y 的函数，与 z 无关。工程实际中，炮
筒、桥梁支座的柱形辊轴等都可简化为平面应变问题。
所以无论是平面应力问题还是平面应变问题，都只 需研究3个应力分量 x ,y ,xy，3个应变分量 x , y , xy
2个位移分量 U和 V。
四、单元划分
单元划分是有限元分析的基本前提，也是有限元 法解题的重要步骤。常用的单元类型有： 杆单元 平面单元 轴对称单元
空间单元 对平面问题，一般采用三角形单元，此时单元划
分应注意以下问题：
任一三角形单元的顶点必须同时也是其相邻三角
形单元的顶点，而不能是其内点。
三角形单元的3条边长（或3个顶角）之间不应相
x y xy
x y xy
8 第四章 用常应变三角形单元解力学平面问题 (2)解析
um xm ym
1 um ym
1 xm um
其中
1 xi yi
2 1 x j y j
1 xm ym
(c) (d) (1)
从解析几何可知，式中的 就是三角形i、j、m的面积。
为保证求得的面积为正值，节点i、j、m的编排次序必须是逆 时针方向，如图1所示。
7. 由单元的节点位移列阵计算单元应力
解出整体结构的节点位移列阵 后，再根据单元节点的 编号找出对应于单元的位移列阵 e，将 e代入（3-3）式就
可求出各单元的应力分量值。
8. 计算结果输出
求解出整体结构的位移和应力后，可有选择 地整理输出某些关键点的位移值和应力值，特别 要输出结构的 变形图、应力图、应变图、结构仿 真变形过程动画图及整体结构的弯矩、剪力图等 等。
平面问题可用三角元，四边元等。
例如：
3. 选择单元的位移模式
结构离散化后，要用单元内结点的位移通过插值来获得 单元内各点的位移。在有限元法中，通常都是假定单元的位 移模式是多项式，一般来说，单元位移多项式的项数应与单 元的自由度数相等。它的阶数至少包含常数项和一次项。至 于高次项要选取多少项，则应视单元的类型而定。
有限元法的实质是：把有无限个自由度的连续体， 理想化为只有有限个自由度的单元集合体，使问题简化 为适合于数值解法的结构型问题。
二、经典解与有限元解的区别：
微分 经 典 解 法 —— （解析法）
数目增到∞ 大小趋于 0
建立一个描述连续体 性质的偏微分方程
有限单元 离散化 集合
总体分析解
有限元法——连续体——单元——代替原连续体
式中：
Re ke e
（3-4）
——单元刚度矩阵
ke BT DBdxdydz
弹性力学平面问题的有限元法
用于描述四节点四边形单元内任意一点的位移和 应力状态。
刚度矩阵
由四节点四边形单元的形状函数和弹性力学基本 公式构建，用于描述单元的刚度特性。
平面六面体八节点单元
六面体八节点单元
是一种三维有限元单元， 具有六个面和八个节点。
形状函数
用于描述六面体八节点 单元内任意一点的位移 和应力状态。
刚度矩阵
对复杂问题的处理能力有限
对于一些高度非线性或耦合问题，有限元法可能难以获得准确解，需要采用其他数值方法 或实验手段。
对高维问题的处理难度较大
随着问题维度的增加，有限元法的计算量和内存消耗会急剧增加，限制了其在高维问题中 的应用。
未来发展方向与挑战
高效算法设计
研究更高效的有限元算法，提高计算速度和精度，降低计算成本。
载荷向量的确定
根据边界条件和外力分布，确定每个节点的载荷 向量。
3
系统刚度矩阵与总载荷向量
将各个单元的刚度矩阵和载荷向量组合起来，形 成系统刚度矩阵和总载荷向量。
求解线性方程组
线性方程组的求解
利用数值方法（如Gauss消去法、迭代法等）求解由 系统刚度矩阵和总载荷向量构成的线性方程组。
解的收敛性与稳定性
02 弹性力学基本方程
应力和应变的关系
01
02
03
胡克定律
在弹性范围内，应力与应 变之间存在线性关系，即 应力与应变成正比。
应变分量
描述物体变形的量，包括 线应变和角应变。
应力分量
描述物体内部受力情况的 量，包括正应力和剪切应 力。
平衡方程
静力平衡
物体在无外力作用下保持静止状态， 即合力为零。
弹性力学平面问题的有限元法
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第四章 平面问题有限元分析及程序设计
§4.1 平面问题单元离散 §4.2 平面问题单元位移模式 §4.3 平面问题单元分析 §4.4 平面问题整体分析 §4.5 平面问题有限元程序设计
有限元网格划分的基本原则
• 网格数目 • 网格疏密 • 单元阶次 • 网格质量 • 网格分界面和分界点 • 位移协调性 • 网格布局 • 结点和单元编号 • 网格自动剖分
f
y
面力
f
f y
xy
xy
基本量和方程的矩阵表示
位移
d
u
v
物理方程 简写为
x y
xy
E
1 2
1
0
1
0
0 0
x y
1
xy
2
D
§4.2 单元位移模式
几何方程：
ux
v y
xvuyT
只要知道了单元的位移函数，就可由几何方程求出应变，再由物理 方程就可求出应力。
（1）位移模式必须能够反映单元的刚体位移； （2）位移模式必须能够反映单元的常应变；
必要条件
（3）位移模式尽可能反映位移的连续性；
u12x3y12x5 23y5 23y v4 5x6y46y5 23x5 23x
u0 1
v0 4
5 3
2
刚体平动
刚体转动
充分条件
u
v
u0 v0
y x
作业： P141 6-1
u12x3y N iuiNjujN m um
其中， N i 、N j 、N m 是系数，是 x、 y 的线性函数；
可以求得：
N i a i b ix ciy2A (i, j, m )
其中：
ai
xj xm
yj ym
1 bi - 1
yj ym
1 ci 1
xj xm
1 A11
2
xi xj
yi yj
1 xm ym
i 结点
j 结点
12xi3yiui 12xj3yjuj
45xi6yivi 45xj6yjvj
m 结点
12xm3ymum 45xm 6ymvm
§4.2 单元位移模式
• 写成矩阵形式
1 xi 0 0
yi 0 0 0 1 xi
0 1 ui
yi
2
vi
1 0
xj 0
因此， u、v 在坐标空间应该为一平面。
um
vm
m
m
ui
vi
i
uj
i
vj
j
j
§4.2 单元位移模式
位移写成向量形式：
d u v N Niiu vii N Njju vjj N Nm mv um m
0 N iu u i i 0 N iv v i i N 0 j u u j j N 0 jv v j j N 0 m u u m m N 0 m v v m m
ui
vi
Ni 0
0 Ni
Nj 0
0 Nj
Nm 0
0 Nm
u v
j j
um
vm
dNe
N 称为形函数矩阵。
§4.2 单元位移模式
有限元分析中，应力转换矩阵、刚度矩阵都是依赖于位移模式建立 起来的，因此，位移模式必须能够反映弹性体的真实位移形态，才 能得到正确的解答。
位移模式需要满足的条件：
1 2
ij
ANidxdy
1 3
A
m 1
1/3 i
1/2
j
(i, j, m)
§4.2 单元位移模式
位移函数：
u 1 2 x 3 y N iu i N ju j N m u m v 4 5 x 6 y N iv i N jv j N m v m
由于 N i 、N j 、N m 是坐标 (x 、y ) 的线性函数， 因此， u、v 也是 (ui,uj,um)、(vi,vj,vm) 的线性函数。
网格数量20万 最小网格尺度150m 最大网格尺度3500m
§4.1 平面问题单元离散
几个重要概念：
平面问题单元： 三角形单元 平面应力： 三角形板 平面应变： 三棱柱
平面问题结点： 绞结点 平面问题约束： 绞支座、链杆 平面问题荷载： 结点荷载和非结点荷载
基本量和方程的矩阵表示
体积力
f
fx
j (xj, yj) uj
vm um
m (xm, ym)
x
假定位移模式如下所示：
ufu(x,y) vfv(x,y)
三个结点的位移也必定满足 位移场函数，因此有：
ui fu(xi,yi) vi fv(xi,yi) uj fu(xj,yj) vj fv(xj,yj)
umfu(xm,ym) vmfv(xm,ym)
有限单元法：未知量是结点的位移分量
e u i v i u j v j u mv m T
那么单元内任意一点的位移跟结点位移有什么关系呢？
因此说，只要知道了位移场的分布，即可解决上述问题。
§4.2 单元位移模式
位移模式：单元位移场分布形式
建立一个坐标系，如下图所示：
y
vj
vi ui i (xi, yi)
yj 0
00 1 xj
0 yj
3 4
=
u v
j j
1
xm
ym
0
0
0
5
u
m
0 0 0 1 xm ym 6 vm
§4.2 单元位移模式
位移模式的选取
因此， 1 、 2 、 3 是 u i 、u j 、u m 的线性函数； 同样， 4 、 5 、 5 是 v i 、v j 、v m 的线性函数； 代入位移场函数，则 u是 u i 、u j 、u m 的线性函数，即：
§4.2 单元位移模式
位移模式的选取
位移函数的选取是任意的，所选取的位移函数越接近于真实情况，所
求得的形变和内力结果就越准确。
最简单的位移场函数是线性函数，即：
u1 2x3y
v4 5x6y
其中， 1 、 2 、 3 、 1 、 2 、 3 是系数，由边界条件求得。
边界条件：在三个结点也应满足位移场函数；
同理，可求得 N j 、N m ，且下标可轮换 ；
注意：i, j, m 必须是逆时针 排列，否则面 积为负。
同理可得：
v 4 5 x 6 y N iv i N jv j N m v m
§4.2 单元位移模式
上式也可以写成：
1x y
1 xj yj
1 Ni 1
xm xi
ym yi
(i, j, m)
1 xj yj
1 xm ym
N i 、N j 、 N m 是坐标 (x 、y ) 的 线性函数；
N i 、N j 、 N m 表明了单元的位移 形态（位移在单元的变化规律）
—称为形态函数，简称形函数
形函数的性质：
Nii 1 Ni j 0 Nim0
Ni oij
1 2
Nio
1 3
N ids
i j