基于MATLAB的有限元法分析平面應(yīng)力應(yīng)變問(wèn)題劉剛0001

1、姓名：劉剛學(xué)號(hào)：15平面應(yīng)力應(yīng)變分析有限元法Abstruct:本文通過(guò)對(duì)平面應(yīng)力/應(yīng)變問(wèn)題的簡(jiǎn)要理論闡述，使讀者對(duì)要分析的問(wèn) 題有大致的印象，然后結(jié)合兩個(gè)實(shí)例，通過(guò) MATLAB軟件的計(jì)算，將有限元分 析平面應(yīng)力/應(yīng)變問(wèn)題的過(guò)程形象的展示給讀者，讓人一目了然，快速了解有限 元解決這類(lèi)問(wèn)題的方法和步驟！一 基本理論有限元法的基本思路和基本原則以結(jié)構(gòu)力學(xué)中的位移法為基礎(chǔ)，把復(fù)雜的結(jié)構(gòu)或連續(xù) 體看成有限個(gè)單元的組合，各單元彼此在節(jié)點(diǎn)出連接而組成整體。把連續(xù)體分成有限 個(gè)單元和節(jié)點(diǎn)，稱(chēng)為離散化。先對(duì)單元進(jìn)行特性分析，然后根據(jù)節(jié)點(diǎn)處的平衡和協(xié)調(diào) 條件建立方程，綜合后做整體分析。這樣一分一合，先離散再綜合

2、的過(guò)程，就是把復(fù) 雜結(jié)構(gòu)或連續(xù)體的計(jì)算問(wèn)題轉(zhuǎn)化簡(jiǎn)單單元分析與綜合問(wèn)題。因此，一般的有限揭發(fā)包 括三個(gè)主要步驟：離散化 單元分析整體分析。二.用到的函數(shù)1. Li nearTria ngleEleme ntStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p)2. Li nearBarAssemble(K k I f)3. L in earBarEleme ntForces(k u)4. L i nearBarEleme ntStresses(k u A)5. L in earTria ngleEleme ntArea(E NU t)三實(shí)例例1考慮如圖所示的受均布載荷作用的薄

3、平板結(jié)構(gòu)。將平板離散化成兩個(gè)線(xiàn) 性三角元，假定 E=200GPa, v=0.3，t=0.025m,w=3000kN/m.1. 離散化2. 寫(xiě)出單元?jiǎng)偠染仃嚭瘮?shù)，得到兩個(gè)單元?jiǎng)偠染仃囃ㄟ^(guò) matlab 的 Li nearTria ngleEleme ntStiffnessk1和k2，每個(gè)矩陣都是6 6的。>> E=210e6210000000>>k1=Li nearTria ngleEleme ntStiffness(E,NU,t,0,0,0.5,0.25,0,0.25,1)k1 =1.0e+006 *Columns 1 through 52.019200 -1.0096

4、-2.019205.7692-0.865400.86540 -0.86541.44230-1.44231.0096000.50481.0096-2.0192 0.8654 -1.4423 1.0096 3.46151.0096 -5.7692 0.8654 -0.5048 -1.8750Colu mn 61.0096-5.76920.8654-0.5048-1.87506.2740>> NU=0.3NU =0.3000>> t=0.0250.0250>> k2=LinearTriangleElementStiffness(E,NU,t,0,0,0.5,0,0

5、.5,0.25,1)k2 =1.0e+006 *Columns 1 through 51.44230-1.44230.8654000.50481.0096-0.5048-1.0096-1.44231.00963.4615-1.8750-2.01920.8654-0.5048-1.87506.27401.00960-1.0096-2.01921.00962.0192-0.865400.8654-5.76920Column 6-0.865400.8654-5.769205.76923. 集成整體剛度矩陣8*8 零矩陣K =00000000000000000000000000000000000000

6、00000000000000000000000000>> K=LinearTriangleAssemble(K,k1,1,3,4)K =1.0e+006 *Columns 1 through 52.0192000005.769200-0.865400000000000-0.8654001.4423-1.00960000-2.01920.865400-1.44231.0096-5.7692000.86541.0e+007 *0.40380.10100-0.576900-0.10100-0.20190.08650-0.0505-0.2019-0.18750.10100.62740 01

7、.1538 -0.0865-0.0865 0.14420 00 -0.1442-0.57690.0865-0.57690.086500-0.1010 -0.2019000.05050.1010-0.0505000-0.14420.10100.4904-0.1875-0.14420.08650-0.2019-0.57690.08650.08650-0.05050-0.1875-0.14420.67790.10100.10100.3462-0.0505-0.1875Columns 6 through 8-1.0096-2.01921.009600.8654-5.76920000000-1.4423

8、0.86540.50481.0096-0.50481.00963.4615-1.8750-0.5048-1.87506.2740>> K=LinearTriangleAssemble(K,k1,1,2,3)4.引入邊界條件 .用上步得到的整體剛度矩陣，可以得到該結(jié)構(gòu)的方程組如下形式3.46150 -1.4423 0.86540 -1.8750 -2.0192 1.0096 'U 1X% 10 6.2740 1.0096 - 0.5048 -1.87500 0.8654 -5.7692SyF1y-1.4423 1.0096 3.4615 -1.8750 -2.01920.86

9、5400U 2XF2X0.8654 -0.5048 -1.87506.2740 1.0096 -5.769200U2yF2y0 -1.8750 -2.01921.0096 3.46150 -1.4423 0.8654U 3XF3X-1.875000.8654 -5.76920 6.2740 1.0096 -0.5048U3yF3y-2.01920.865400 -1.4423 1.00963.4615 -1.8750U 4XF4X1.0096 -5.769200 0.8654 - 0.5048 -1.8750 6.2740 一U4y -IF4y本題的邊界條件:1063yF?x = 9.375,

10、 F?y = 0, F3X = 9.375, F將邊界條件帶入，得到:3.46150 -1.4423 0.86540 -1.8750-2.01920 6.2740 1.0096 - 0.5048-1.87500 0.8654-1.4423 1.00963.4615 -1.8750 -2.01920.86546100 -1.8750 -2.01921.00963.4615， 0-1.44230.8654-1.87500 0.8654-5.769206.27401.0096-0.5048-2.01920.865400 -1.44231.00963.4615-1.87501.0096-5.76920

11、0 0.8654 -0.5048-1.87506.27400.8654 -0.5048 -1.87506.2740 1.0096-5.76920001耳0F1yU2X9.375U2y0U3X9.375U3y00F4X0F4y1.0096-5.76920K的第3-6行的第3-6列作為子矩陣3.4615- 1.87506-1.87506.274010-2.01921.0096-0.8654-5.7692Matlab命令-2.01921.00963.46155.解方程分解上述方程組，提取總體剛度矩陣>> k=K(3:6,3:6)0.8654-5.76926.2740U 2XU2yU 3X

12、9. 3759. 375k =1.0e+006 *3.4615-1.8750-2.01920.8654-1.87506.27401.0096-5.7692-2.01921.00963.461500.8654-5.76920 6.2740>> f=9.375;0;9.375;09.3750 09.37500>> u=kfu =1.0e-005 *0.71110.11150.65310.0045現(xiàn)在可以清楚的看出，節(jié)點(diǎn) 2的水平位移和垂直位移分別是0.7111m和0.1115m。節(jié)點(diǎn)3的水平位移和垂直位移分別是0.6531m 和0.0045m。6. 后處理用matlab命令

13、求出節(jié)點(diǎn)1和節(jié)點(diǎn)4的支反力以及每個(gè)單元的應(yīng)力。首先建立總體節(jié)點(diǎn)位移矢量 U，U=0;0;u;0;0U =1.0e-005 *000.71110.11150.65310.004500>> F=K*U-9.3750-5.62959.37500.00009.37500.0000-9.37505.6295由以上知，節(jié)點(diǎn)1的水平反力和垂直反力分別是9.375k n(指向左邊)和5.6295kn (作用力方向向下)，節(jié)點(diǎn) 4的水平反力和垂直反力分別是9.375kn (指向左邊)和 5.6295kn(作用力方向向下) . 滿(mǎn)足力平衡條件。接著，建立單元節(jié)點(diǎn)位移矢量比和u2，然后調(diào)用matlab命

14、令LinearTriangleElementStresses計(jì)算單元應(yīng)力 sigma1 和 sigma2>> u1=U(1);U(2);U(5);U(6);U(7);U(8)u1 =1.0e-005 *000.65310.004500>> u2=U(1);U(2);U(3);U(4);U(5);U(6)u2 =1.0e-005 *000.71110.11150.65310.0045>>sigma1=LinearTriangleElementStresses(E,NU,0.025,0,0,0.5,0.25,0,0.25,1,u1)sigma1 =1.0e+00

15、3 *3.01440.90430.0072>>sigma2=Li nearTria ngleEleme ntStresses(E,NU,0.025,0,0,0.5,0,0.5,0.25,1,u2)sigma2 =1.0e+003 *2.9856-0.0036-0.0072由以上可知，單元1的應(yīng)力二x =3.0144MPa（拉應(yīng)力），二y =0.9043MPa （拉應(yīng)力），xy = 0.0072MPa（正值）。單元2的應(yīng)力是二 x =2.9856MPa（拉應(yīng)力）：=0.0036MPa（壓應(yīng)力）岑=0.0072MPa（負(fù)值）。顯然，在x方向的應(yīng)力（拉應(yīng)力）接近于正確的值3MPa （拉應(yīng)

16、力）。接著調(diào)用Lin earTria ngleEleme ntStresses函數(shù)計(jì)算每個(gè)單元的主應(yīng)力和主應(yīng)力方向角。>> s1= Li nearTria ngleEleme ntPStresses（sigmal）si =1.0e+003 *3.01440.90430.0002>> s2= Lin earTria ngleEleme ntPStresses(sigma2)s2 =2.98561.0e+3 *-0.0036-0.00016 =3.0144MPa（拉應(yīng)力），二0.9043MPa（拉應(yīng)力）主應(yīng)力方向角 =0.2二1 = 2.9856MPa（拉應(yīng)力），二0.00

17、36MPa（壓應(yīng)力），二p - -0.1例2考慮如圖3.1所示的由均勻分布載荷和集中載荷作用的薄平板結(jié)構(gòu)。將平板離散化成12個(gè)線(xiàn)性三角單元，如圖 4所示。假定E=210GPa，v=03，t=0025m，w=100kN/m 和 P=125kN。1.離散化2. 寫(xiě)出單元?jiǎng)偠染仃?gt;> E=201e6;>> NU=0.3;>> t=0.025;>> k1= Li nearTria ngleEleme ntStiffness(E,NU,t,0,0.5,0.125,0.375,0.25,0.5,1);>> k2= Lin earTria ngle

18、Eleme ntStiffness(E,NU,t,0,0.5,0,0.25,0.125,0.375,1);>> k3= Li nearTria ngleEleme ntStiffness(E,NU,t,0.125,0.375,0.25,0.25,0.25,0.5,1);>> k4= Li nearTria ngleEleme ntStiffness(E,NU,t,0.125,0.375,0,0.25,0.25,0.25,1);>> k5= Li nearTria ngleEleme ntStiffness(E,NU,t,0,0.25,0.125,0.125,

19、0.25,0.25,1);>> k6= Li nearTria ngleEleme ntStiffness(E,NU,t,0,0.25,0,0,0.125,0.125,1);>> k7= Li nearTria ngleEleme ntStiffness(E,NU,t,0.25,0.25,0.125,0.125,0.25,0,1);>> k8= Lin earTria ngleEleme ntStiffness(E,NU,t,0.125,0.125,0,0,0.25,0,1);>> k9= Lin earTria ngleEleme ntStif

20、fness(E,NU,t,025,0.25,0.25,0,0.375,0.125,1);>> k10= Li nearTria ngleEleme ntStiffness(E,NU,t,0.25,0.25,0.375,0.125,0.5,0.25,1);>> k11= Lin earTria ngleEleme ntStiffness(E,NU,t,0.25,0,0.5,0,0.375,0.125,1);>> k12= Li nearTria ngleEleme ntStiffness(E,NU,t,0.375,0.125,0.5,0,0.5,0.25,1)

21、k1 =1.8637-0.8973-0.96630.8283-0.89730.0690-0.89731.86370.9663-2.7610-0.06900.89731.0e+6 *-0.9663 0.9663 1.9327 0 -0.9663 -0.96630.8283 -2.7610 0 5.5220 -0.8283 -2.7610-0.8973 -0.0690 -0.9663 -0.8283 1.8637 0.89730.0690 0.8973 -0.9663 -2.7610 0.8973 1.86373. 集成整體剛度矩陣：>>K=zero(22,22);>>K=

22、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,

23、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) 運(yùn)行得1.0e+008 *Columns 1 through 70.0389-0.0187-0.00940.0007-0.03

24、890.01870.0094-0.01870.0389-0.00070.00940.0187-0.03890.0007-0.0094-0.00070.03890.0187-0.0389-0.018700.00070.00940.01870.0389-0.0187-0.03890-0.03890.0187-0.0389-0.01870.15580-0.03890.0187-0.0389-0.0187-0.038900.1558-0.01870.00940.000700-0.0389-0.01870.0779-0.0007-0.009400-0.0187-0.03890000.0094-0.000

25、7-0.03890.0187-0.0187000.0007-0.00940.0187-0.03890000000-0.03890000000.01870000000.0094000000-0.000700000000000000000000000000000000000000000000000000000000Columns 8 through 14-0.0007000000-0.009400000000.00940.000700000-0.0007-0.00940000-0.0187-0.03890.01870000-0.03890.0187-0.038900000-0.01870-0.03

26、890.01870.0094-0.00070.077900.01870.0187-0.03890.0007-0.009400.0972-0.0093-0.0389-0.0187000.0187-0.00930.0972-0.0187-0.0389000.0187-0.0389-0.01870.15580-0.0389-0.0187-0.0389-0.0187-0.038900.1558-0.0187-0.03890.000700-0.0389-0.01870.03890.0187-0.009400-0.0187-0.03890.01870.03890-0.00090.0095-0.03890.

27、0187-0.00940.000700.0095-0.03840.0187-0.0389-0.00070.009400.0004-0.000200000-0.00020.000400000-0.0094-0.0007000000.00070.0094000000000000000000Columns 15 through 2100000000000000000000000000000000000000000000000000000000-0.00090.00950.0004-0.0002 -0.00940.000700.0095-0.0384-0.00020.0004 -0.00070.009

28、40-0.03890.0187000000.0187-0.038900000-0.0094-0.0007000000.00070.009400000-1.94080.00951.9994-0.037700-0.00940.0095-5.6533-0.03775.7119000.00071.9994-0.0377-1.92190.0379-0.0389-0.0187-0.0389-0.03775.71190.0379-5.6344-0.0187-0.03890.018700-0.0389-0.01870.03890.01870.009400-0.0187-0.03890.01870.0389-0

29、.0007-0.00940.0007-0.03890.01870.0094-0.00070.0389-0.00070.00940.0187-0.03890.0007-0.0094-0.0187Column 2200000000000000-0.00070.00940.0187-0.03890.0007-0.0094-0.01870.0389 0.0007-0.0094-0.01870.03894. 引入邊界條件：U1x= U 1y= U 4x= U4y=U 7x=U7y=0F2x= F2y= F3x= F3y=F6x=F6y=F8x= F8y= F9x= F 9y =F10x=F 10y= F

30、11x= F11y= 0 F5x= 0,F5y= -12.55. 解方程： >>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);01.0e+008 *Columns 1 through 80.0094 -0.0007-0.01870.0389-0.00070.00940.0187-0.0389-0.0094-0.00070.03890.0187-0.0389-0.01870.00070.00

31、940.01870.0389-0.0187-0.0389-0.03890.0187-0.0389-0.01870.155800.0187-0.0389-0.0187-0.038900.15580.00940.000700-0.0389-0.0187-0.0007-0.009400-0.0187-0.0389000.0094-0.0007-0.03890.0187000.0007-0.00940.0187-0.03890000000000000000000000000000000000000000000000000000000000000.018700000.0007 -0.03890.0389

32、 -0.0187 -0.00940.0007 -0.0094-0.0389 -0.0187-0.0187 -0.03890.0779-0.0187-0.03890.07790.01870.01870.0187 -0.03890.00940.0007-0.0007-0.0094Columns 9 through 16000000000000000.0094 0.00070000000000000000-0.0007-0.0094000000-0.03890.01870000000.0187-0.0389000000-0.01870-0.03890.01870.0094-0.00070000.01

33、870.0187-0.03890.0007-0.0094000.0972-0.0093-0.0389-0.018700-0.00090.0095-0.00930.0972-0.0187-0.0389000.0095-0.0384-0.0389-0.01870.15580-0.0389-0.0187-0.03890.0187-0.0187-0.038900.1558-0.0187-0.03890.0187-0.038900-0.0389-0.01870.03890.0187-0.0094-0.000700-0.0187-0.03890.01870.03890.00070.0094-0.00090.0095-0.03890.0187-0.00940.0007-1.94080.00950.0095-0.03840.0187-0.0389-0.00070.00940.0095-5.65330.0004-0.000200001.9994-0.0377-0.00020.00040000-0.03775.7119-0.0094-0.0