姓名:刘刚学号:15平面应力应变分析有限元法Abstruct:本文通过对平面应力/应变问题的简要理论阐述,使读者对要分析的问题有大致的印象,然后结合两个实例,通过MATLAB软件的计算,将有限元分析平面应力/应变问题的过程形象的展示给读者,让人一目了然,快速了解有限元解决这类问题的方法和步骤!一.基本理论有限元法的基本思路和基本原则以结构力学中的位移法为基础,把复杂的结构或连续体看成有限个单元的组合,各单元彼此在节点出连接而组成整体。
把连续体分成有限个单元和节点,称为离散化。
先对单元进行特性分析,然后根据节点处的平衡和协调条件建立方程,综合后做整体分析。
这样一分一合,先离散再综合的过程,就是把复杂结构或连续体的计算问题转化简单单元分析与综合问题。
因此,一般的有限揭发包括三个主要步骤:离散化单元分析整体分析。
二.用到的函数1. LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p)(K k I f)(k u)(k u A)(E NU t)三.实例例1.考虑如图所示的受均布载荷作用的薄平板结构。
将平板离散化成两个线性三角元,假定E=200GPa,v=,t=0.025m,w=3000kN/m.1.离散化2.写出单元刚度矩阵通过matlab 的LinearTriangleElementStiffness 函数,得到两个单元刚度矩阵1k 和2k ,每个矩阵都是6 6的。
>> E=210e6 E =>> k1=LinearTriangleElementStiffness(E,NU,t,0,0,,,0,,1) k1 = +006 *Columns 1 through 50 0 0 0 0 0 0 0Column 6 >> NU= NU = >> t= t =>> k2=LinearTriangleElementStiffness(E,NU,t,0,0,,0,,,1)k2 =+006 *Columns 1 through 50 00 0Column 63.集成整体刚度矩阵 8*8零矩阵K =0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> K=LinearTriangleAssemble(K,k1,1,3,4)K =+006 *Columns 1 through 50 0 0 00 0 00 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0Columns 6 through 80 0 0 0 0 0 0>> K=LinearTriangleAssemble(K,k1,1,2,3) K =+007 *0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 04.引入边界条件.用上一步得到的整体刚度矩阵,可以得到该结构的方程组如下形式⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡4y 4X 3y 3X 2y 2X 1y 1X 4y 4X 3y 3X 2y 2X 1y 1X 6F F F F F F F F U U U U U U U U 6.2740 1.8750- 0.5048- 0.8654 0 0 5.7692- 1.0096 1.8750- 3.4615 1.0096 1.4423- 0 0 0.8654 2.0192- 0.5048- 1.0096 6.2740 0 5.7692- 0.8654 01.8750- 0.8654 1.4423- 0 3.4615 1.00962.0192- 1.8750- 0 0 0 5.7692- 1.0096 6.2740 1.8750- 0.5048- 0.8654 0 0 0.8654 2.0192- 1.8750-3.4615 1.0096 1.4423- 5.7692- 0.8654 0 1.8750- 0.5048- 1.0096 6.2740 0 1.0096 2.0192- 1.8750- 0 0.8654 1.4423- 03.461510 本题的边界条件:04411====y x y x U U U U0,375.9,0,375.93322====y x y x F F F F将边界条件带入,得到:⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡4y 4X 1y 1X 3y 3X 2y 2X 6F F 09.37509.375F F 0 0 U U U U 0 0 6.2740 1.8750- 0.5048- 0.8654 0 0 5.7692- 1.0096 1.8750- 3.4615 1.0096 1.4423- 0 0 0.8654 2.0192- 0.5048- 1.0096 6.2740 0 5.7692- 0.8654 01.8750- 0.8654 1.4423- 0 3.4615 1.00962.0192- 1.8750- 0 0 0 5.7692- 1.0096 6.2740 1.8750- 0.5048- 0.8654 0 0 0.8654 2.0192- 1.8750-3.4615 1.0096 1.4423- 5.7692- 0.8654 0 1.8750- 0.5048- 1.0096 6.2740 0 1.0096 2.0192- 1.8750- 0 0.8654 1.4423- 03.4615105.解方程分解上述方程组,提取总体刚度矩阵K 的第3-6行的第3-6列作为子矩阵⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡ 09.37509.375 U U U U 6.2740 0 5.7692- 0.8654 0 3.4615 1.0096 2.0192- 5.7692- 1.0096 6.2740 1.8750- 0.8654 2.0192- 1.8750-3.4615103y 3X 2y 2X 6 Matlab 命令>> k=K(3:6,3:6) k =+006 *0 0>> f=[;0;;0] f = 0 0 >> u=k\f u =*现在可以清楚的看出,节点2的水平位移和垂直位移分别是0.7111m 和0.1115m 。
节点3的水平位移和垂直位移分别是0.6531m 和0.0045m 。
6.后处理用matlab 命令求出节点1和节点4的支反力以及每个单元的应力。
首先建立总体节点位移矢量U ,U=[0;0;u;0;0] U = *0 0 0 0 >> F=K*U F =由以上知,节点1的水平反力和垂直反力分别是(指向左边)和(作用力方向向下),节点4的水平反力和垂直反力分别是(指向左边)和(作用力方向向下).满足力平衡条件。
接着,建立单元节点位移矢量21u 和u ,然后调用matlab 命令LinearTriangleElementStresses 计算单元应力sigma1和sigma2>> u1=[U(1);U(2);U(5);U(6);U(7);U(8)]u1 =*>> u2=[U(1);U(2);U(3);U(4);U(5);U(6)]u2 =*>> sigma1=LinearTriangleElementStresses(E,NU,,0,0,,,0,,1,u1) sigma1 =+003 *>> sigma2=LinearTriangleElementStresses(E,NU,,0,0,,0,,,1,u2) sigma2 =+003 *由以上可知,单元1的应力)(0144.3拉应力MPa x =σ,)(9043.0y 拉应力MPa =σ,)(0072.0y 正值MPa x =τ 。
单元2的应力是)(9856.2拉应力MPa x =σ)(0036.0y 压应力MPa =σ)(0072.0y 负值MPa x =τ。
显然,在x方向的应力(拉应力)接近于正确的值3MPa (拉应力)。
接着调用LinearTriangleElementStresses 函数计算每个单元的主应力和主应力方向角。
>> s1= LinearTriangleElementPStresses(sigma1) s1 = +003 *>> s2= LinearTriangleElementPStresses(sigma2) s2 = +003 *)(0144.31拉应力MPa =σ,)(9043.02拉应力MPa =σ主应力方向角ο2.0=p θ )(M 9856.21拉应力Pa =σ,)(0036.02压应力MPa =σ,ο1.0-=p θ例 2.考虑如图所示的由均匀分布载荷和集中载荷作用的薄平板结构。
将平板离散化成12个线性三角单元,如图4所示。
假定E=210GPa,v=,t=0.025m,w=100kN/m和P=。
1.离散化2.写出单元刚度矩阵>> E=201e6;>> NU=;>> t=;>> k1= LinearTriangleElementStiffness(E,NU,t,0,,,,,,1); >> k2= LinearTriangleElementStiffness(E,NU,t,0,,0,,,,1); >> k3= LinearTriangleElementStiffness(E,NU,t,,,,,,,1); >> k4= LinearTriangleElementStiffness(E,NU,t,,,0,,,,1); >> k5= LinearTriangleElementStiffness(E,NU,t,0,,,,,,1); >> k6= LinearTriangleElementStiffness(E,NU,t,0,,0,0,,,1); >> k7= LinearTriangleElementStiffness(E,NU,t,,,,,,0,1); >> k8= LinearTriangleElementStiffness(E,NU,t,,,0,0,,0,1); >> k9= LinearTriangleElementStiffness(E,NU,t,025,,,0,,,1); >> k10= LinearTriangleElementStiffness(E,NU,t,,,,,,,1); >> k11= LinearTriangleElementStiffness(E,NU,t,,0,,0,,,1); >> k12= LinearTriangleElementStiffness(E,NU,t,,,,0,,,1)k1 =+006 *3.集成整体刚度矩阵:>>K=zero(22,22);>>K=LinearTriangleAssemble(K,k1,1,3,2);>>K=LinearTriangleAssemble(K,k2,1,4,3);>>K=LinearTriangleAssemble(K,k3,3,5,2);>>K=LinearTriangleAssemble(K,k4,3,4,5);>>K=LinearTriangleAssemble(K,k5,4,6,5);>>K=LinearTriangleAssemble(K,k6,4,7,6);>>K=LinearTriangleAssemble(K,k7,5,6,8);>>K=LinearTriangleAssemble(K,k8,6,7,8);>>K=LinearTriangleAssemble(K,k9,5,8,9);>>K=LinearTriangleAssemble(K,k10,5,9,10);>>K=LinearTriangleAssemble(K,k11,8,11,9);>>K=LinearTriangleAssemble(K,k12,9,11,10)运行得+008 *Columns 1 through 70 00 0 00 00 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0Columns 8 through 140 0 0 0 0 00 0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 00 0 0 00 00 0 00 00 00 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0Columns 15 through 210 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 00 00 00 0Column 224.引入边界条件:U1x= U1y= U4x= U4y=U7x=U7y=0F2x= F2y= F3x= F3y=F6x=F6y=F8x= F8y= F9x= F9y=F10x=F10y= F11x= F11y= 0F5x= 0,F5y=5.解方程:>>k=[K(3:6,3:6),K(3:6,9:12),K(3:6,15:22);K(9:12,3:6),K(9:12,9:12),K(9:12,15:22) ;K(15:22,3:6),K(15:22,9:12) ,K(15:22,15:22)];K =+008 *Columns 1 through 80 00 00 0 00 0 00 0 00 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0Columns 9 through 160 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 00 0 00 00 00 00 00 0 0 00 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0Columns 17 through 220 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 00 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 00 0>>f=[0;0;0;0;0;;0;0;0;0;0;0;0;0;0;0];>>u=k\fu=*[ ]T6.后处理>>U=[0;0;u(1:4);0;0;u(5:8);0;0;u(9:16)];>>F=K*UF=[ 0 0 ]T三.小结通过这次练习,对有限元结合MATLAB软件解决平面应力/应变问题的方法和过程逐渐清晰,特别是在应用MATLAB软件的过程中学到了很多东西。