有限元分析課程作業(yè)

1、有限元分析課程作業(yè)任課教師：徐亞蘭學(xué)生姓名：陳新杰學(xué) 號：班級：1304012時間：2016-01-05、問題描述及分析問題：如圖1所示，有一矩形平板，在右側(cè)受到P=10KN/m的分布力，材料常數(shù)為：彈性模量E 1 107Pa；泊松比 1/3 ； 板的厚度為t=;試按平面應(yīng)力問題利用三角形與矩形單元分 別計算各個節(jié)點位移及支座反力。P=10KN/m圖1平面矩形結(jié)構(gòu)的有限元分析分析：使用兩種方案：一、基于 3節(jié)點三角形單元的有限元建模，將矩形劃分為兩個 3節(jié)點三角形單元；二、基于4節(jié)點矩形單元的有限元建模，使用一個4節(jié)點矩形單元。利用MATLABM牛計算由各要求量，再將兩種方案的計算結(jié)果進(jìn)行比較

2、、分析、得由結(jié)論。二、有限元建模及分析1、基于3節(jié)點三角形單元的有限元建模及分析(1)結(jié)構(gòu)的離散化與編號如圖2所示，將平面矩形結(jié)構(gòu)分為兩個 3節(jié)點三角形單元。單元三個節(jié)點的編號為1, 2, 4,單元三個節(jié)點的編號為3, 4, 2,各個節(jié)點的位置坐標(biāo)為xi,yi ,i 1,2,3,4 ,各個節(jié)點的位移（分別沿X方向和y方向）為Ui，Vi ,i 1,2,3,4。圖2方案一：使用兩個 3節(jié)點三角形單元（2）各單元的剛度矩陣及剛度方程 a.單元的幾何和節(jié)點描述單元有6個節(jié)點位移自由度（DOF。將所有節(jié)點上的位移組成一個列陣，記作 q；同樣，將所有節(jié)點上的各個力也組成一個列陣，記作 F（1）,則有q（u

3、1，V1,U2,V2,U4,V4）F （Fx1，F(xiàn)y1，F(xiàn)x2，F(xiàn)y2，F(xiàn)x4，F(xiàn)y4）同理，對于單元，有q（U3,V3, U4,V4,U2,V2）F (Fx3,Fy3,Fx4,Fy4,Fx2,Fy2)b.單元的位移場描述對于單元，設(shè)位移函數(shù)u(x,y) % ax a2y(1-1)v(x,y) b0 bix b2y由節(jié)點條件，在x x,y y處，有u(xi,yi) uii 1,2,4v(xi, yi) vi(1-2)將式（1-1）代入節(jié)點條件式（1-2）中，可求生式（1-1 )中待定系數(shù)，即12A12A12AU1x1U2x2U4x4U2U4x1y1y2y4y1y2y4U1又2 U2x4 u4L

4、gv1 a2v22A1 ,、(a1U1 a2U2 a3U4) 2A1 ,.、(b1U1 b2U2 b3U4) 2A1(G3 C2U2 C3U4)2Aa3v4)(1-6)(1-3)(1-4)(1-5)1-(b!v1 b2v2 b3v4)(1-7)2A1 ,(qv1 C2v2 C3v4)(1-8)2 A在式（1-3） 式（1-8）中1.一 一 一、2(a a? a3)(1-9)1Xy11x2、21x4y4X2y2X4 V4n ； ；2y2 V、(1-10)1 x2qx2 x42 x4 4(1,2,3)上式中的符號(1,2,3)表示下標(biāo)輪換，如1 2,2 3,3 1同時更換。將單元各節(jié)點的位置坐標(biāo)x

5、1 0,y1 0,x2 1,y2 0,x4 0,y、1代入得a1 1,b11,q 1a2 0,b2 1,C2 0a30h 0,C3 1A 13將系數(shù)式(1-3)式(1-8)式代入(1-1)中，重寫位移函數(shù)，并以節(jié)點位移的形式進(jìn)行表示，有(1-11)u(x,y) N(x, y)。 N2(x,y)u2 N、(x, y)u、v(x,y) N(x,y)v1 N2(x,yM N、(x, y)v、N(x, y)令M(ai bix jy),i 1,2,3,則有形狀函數(shù)矩陣2A(1)N (x, y)(2 6)N10N20N400N10N20N41 x y 0 x 0 y 001 x y 0 x 0 y(1-1

6、2)位移函數(shù)式(1-11)寫成矩陣形式，有UiVi u(x, y) 1 x y 0 x 0 y 0U2u (x, y)(2U6)v(x, y)01 x y 0 x 0 yV2U4V4對于單元，過程同上，有形狀函數(shù)矩陣1 x y 01 x 01 y 0N("y) 01 x y 01 x 01 y(1-13)(1-14)(2)/、 u(x, y) 1+x y11(x, y)(2U6)v(x, y)01U3V301 x 01 y 0 u4x y 01 x 01 y v4位移函數(shù)U2V2(1-15)c.單元的應(yīng)變場描述/ 、(x,y)(3 1)xxyy xy對于單元，應(yīng)變函數(shù)u(x,y) v

7、(x, y)0x1 x y 0 x 0 y 0=0y 01 x y 0 x 0 yUi viU2V2U4V4uiVi-10 10 0 0= 0 -1 0 0 0 1 u2(1-16)V2-1-10 110u4V4/ 、B (x,y) (3"6)-10 10000 -1 0 001-1-10 110(1-17)其中幾何矩陣(2)(x,y)(31)u3V3u4V4u2V2對于單元，應(yīng)變方程一 0 x1 x y 01 x 01 y 00y 01 x y 01 x 01 yu3V310-10003=01000-1u4(1-18)V4110-1-10u2V2其中幾何矩陣1(2)B(x,y)0(

8、3 6)/1-10000-100-10-10(1-19)d.單元的應(yīng)力場描述xx(x,y,Z) yy (3 1) xy001""2xxyyxy=R(1)(3(1)1)(1-20)其中，彈性系數(shù)矩陣(1)D) =10012D為107-21-3(1)(3 1)=只（B）9(1)(S)9其中應(yīng)力函數(shù)矩陣s為(1)(1)(1)*D)裊=3753106 10-10-10-1-1應(yīng)力方程為=3.75106(1-21(1-22)01 =3.750-3106 -1-1-1-3-1(1-23)u1-1300-3100-1011v1- 3 -1 3 0 0 11u2- 1-310032Pav2

9、- 1-101102u4v4-3(1) =3.75 106 -1(3 1)-11(1)3 q =3.75 100 (6 1)(1-25)(1-26)(1-24)對于單元，過程同上彈性系數(shù)矩陣D為31(2) 6D =3.75 106 1 3(3 3)00應(yīng)力函數(shù)矩陣S為31-300-1(2)S=13-100-3(3 6)110-1-10應(yīng)力方程(2) 為3(2)6=3.75 106 1(3 1)1-3-1 00000-1 -1-131-300(2)6-3q=3.75 10613-1000(6 1)110-1-1-1-30u3v3u4Pav4u2v2(1-27)e. 單元的勢能表達(dá)K (1)是單元

10、剛度矩陣，即K(1)(1)B(1)TD(1)B(1)d(1)B(1)TD(1)B(1)dA t(1-28)其中薄板厚度t 0.1m o將式(1-17)、式(1-21)代入式(1-29),u1v1u2v2u4v47.53.755.6251.8751.8751.875u13.757.51.8751.8751.8755.625v1155.625K(1) 105-1.8755.625001.875u2-1.8751.87501.8751.8750v21.8751.87501.8751.8750u41.8755.6251.875005.625v4同理，得到單元的剛度陣為u3v3u4v4u2v27.53.

11、755.6251.8751.8751.875u13.757.51.8751.8751.8755.625v1255.625K(2)105-1.8755.625001.875u2-1.8751.87501.8751.8750v21.8751.87501.8751.8750u41.8755.6251.875005.625v4得到單元的剛度陣K(1) B(1)TD(1)B(1)tA 3.75 106 0.1計算得-10-10-1 -131001001000101010-10300-101-1-1100000010110將兩個單元按節(jié)點位移所對應(yīng)的位置進(jìn)行組裝，得到總體剛度矩陣為KK K(1) K(2)

12、節(jié)點力F ( FRx1 FRy1 FPx2 FPy2 FPx3 FPx3 FRx4 F Rx4 )13-0.1 10 10 (FRx1 0 10 10 FRx4 0)0.5 103(FRx1 0 10 10 FRx4 0)N系統(tǒng)的勢能U W= ；qTKq-FTq (計算結(jié)果在下面呈現(xiàn)) (4)邊界條件的處理及方程求解邊界條件為U1 V1 u4 v4 0。因此，將針對節(jié)點 2和節(jié)點 3的位移求解，節(jié)點2和節(jié)點3對應(yīng)總體剛度陣 KK中的第3 行到第6行、第3列到第6歹U,則需從KK8 8中提由，置給 k,然后生成對應(yīng)的載荷列陣 p,再采用高斯消去法進(jìn)行求解。>> k=KK(3:6,3:

13、6);>> p=500;0;500;0;>> u=kpu=將列排成了行再計算支反力。在得到整個結(jié)構(gòu)的節(jié)點位移后，由原整 體剛度方程就可以計算由對應(yīng)的支反力；先將上面得到的位 移結(jié)果與邊界條件的節(jié)點位移進(jìn)行組合，得到整體的位移列 陣U(8 1),再代回原整體剛度方程，計算由所有的節(jié)點力， 按照位置關(guān)系我由對應(yīng)的支反力。>> U=0;0;u;0;0 將列排成了行>> P=KK*UP =-500 500 0 500 0 -500 將列排成了行所以，節(jié)點1 的支反力為FRx1500N, FRy1 -176.4706N ，節(jié)點2 的支反力為FRx2500N

14、,FRy2 176.4706N。根據(jù)已求得的位移和支反力計算系統(tǒng)的勢能。>> A=*U'*KK*U-P'*UA =（5 ）結(jié)果分析上述支反力計算結(jié)果滿足靜力平衡，驗證了以上求解過程及MATLA徵法的正確性。2、基于四節(jié)點四邊形單元的有限元建模及分析（ 1）結(jié)構(gòu)的離散化與編號如圖 3 所示一個4 節(jié)點矩形單元，單元的節(jié)點位移共有8 個自由度（DOF） 。節(jié)點編號為1,2,3,4 ，各自的位置坐標(biāo)為X，yi ,i 1,2,3,4 ,各個節(jié)點的位移（分別沿x方向和y方向）為ui,vi ,i 1,2,3,4 。4e圖3方案二：使用一個4節(jié)點矩形單元(2)局部坐標(biāo)系下單元的描

15、述 a.單元的幾何和節(jié)點描述采用無量綱坐標(biāo)=x y a, b其中a 0.5,b 0.5。則單元四個節(jié)點的幾何位置為1, 11,1,1, 41111q(8 1)F1)b.單元的位移場描述將所有節(jié)點上的位移組成一個列陣，記作q；同樣，將所有節(jié)點上的各個力也組成一個列陣，記作F,則有(U1 V1 U2 V2 U3 V3 U4 V4)T (Fx1 Fy1 Fx2 Fy2 Fx3 Fy3 Fx4 Fy4)T設(shè)位移函數(shù)為u(x,y) a ax v(x, y)鳳 b)x由節(jié)點條件，在u(xi,y。 Uiiv(xi,yj via2 y a3xyb2 y b3xyx x,y y處，有1,2,3,4a0，a，a2

16、，%將位移試函數(shù)代入節(jié)點條件中，求生待定系數(shù) 和b0,b,b2,b3,再代入位移函數(shù)中，整理后得u(x,y) Ni(x, y)ui N2(x,y)U2 N3(x, y)U3 N4(x, y)u4 v(x, y) Ni(x, y)vi N?(x, yM N3(x,y)v3 N4(x, yM其中1Ni(x, y) - 1 2x 1 2y41N2(x,y)1 2x 1 2y41M(x,y) - 1 2x 1 2y41N4(x,y) 2y 0=-01 2x 1 2x 1 2y4如以無量綱坐標(biāo)來表達(dá)，可寫成1Ni -11 ,i 1,2,3,441帶入上式J等 1 1, 1 1， 21, 2 1，31,

17、31， 4 1, 4得到形狀函數(shù)矩陣1N11141N21-141N3113 41N4 1 1-14寫成矩陣形式，有UiVlH(x，y)u(x, y)v(x, y)U2N10N20N30N40v20N10N20N30N4u3N qlU4V4C.單元的應(yīng)變場描述單元應(yīng)變?yōu)閤x(31)(x，y)- Nq 且q(3 1)(3 2) 1 1 2x 1 2y 8 8 1(3 8) 8 1xy其中幾何矩陣B(x,y)為0N40xNi0N20N30N40y 0N10N20N30y x(1 2y)001 2x1 2x(1 2y)(1 2y)00(1 2x)(1 2x)(1 2y)1 2y 00(1 2x)(1

18、2x)1 2yd.單元的應(yīng)力場描述應(yīng)力表達(dá)式為DD息q(sq(3 1)(3 3)(3 1)(3 3)(3 8) (8 1)(3 8)(8 1)其中，應(yīng)力函數(shù)矩陣S DB。e.單元的勢能表達(dá)以上已將單元的三大基本變量U,用基于節(jié)點的位移列陣q來表示；將其代入單元勢能表達(dá)式中，有=lqTKq FTq ,其中K為4節(jié)點矩形單元的剛度矩陣，即kiiK A BT DBdA tk21k22symk31k32k33k41k42k43k44其中，t為薄板的厚度，t 0.1m,上式的各個字塊矩陣為krs。BTDBsr,s 1,2,3,4(2 2)(2 3) (3 3) (3 2)f.單元剛度陣及剛度方程單元剛度

19、陣在上面已經(jīng)列由。將單元的勢能對節(jié)點位移q取一階極值，可得到單元的剛度方程風(fēng)9）（F）（4）邊界條件的處理及方程求解處理方法與3節(jié)點三角形單元一致，利用上述求解程序具有的可移植性，簡化了求解過程。>> k=K(3:6,3:6);>> p=500;0;500;0;>> u=kp將列排成了行再計算支反力。同樣注意按照位置關(guān)系我由對應(yīng)的支反力。>> U=0;0;u;0;0U = *0 00 0 將列排成了行>> P=K*UP =-500 500 0 500 0 -500 將列排成了行所以，節(jié)點1的支反力為Frxi500N,FRyi -111

20、.1111N ,節(jié)點2 的支反力為FRx2500N,FRy2 111.1111N 。根據(jù)已求得的位移和支反力計算系統(tǒng)的勢能。>> A=*U'*K*U-P'*UA =(5 )結(jié)果分析基于 4 節(jié)點矩形單元計算出的勢能小于基于3 節(jié)點三角形單元計算出的結(jié)果，若將該系統(tǒng)分為更多的單元，計算精度也將提高。3. 兩種方案的比較與分析從以上計算可以看出，用三角形單元計算時，由于形函數(shù)是完全一次式，因而其應(yīng)變場在單元內(nèi)均為常數(shù)；而對于四邊形單元，其形函數(shù)帶有二次式，計算得到的應(yīng)變場和應(yīng)力場都是坐標(biāo)的一次函數(shù)，但不是完全的一次函數(shù)，對提高精度有一定作用；根據(jù)最小勢能原理，勢能越小，

21、則整體計算精度越高，比較兩種單元計算得到的系統(tǒng)勢能，可看出，在相同的節(jié)點自由度情況下，矩形單元的計算精度要比三角 形單元高。三、基于MATLAB勺編程實現(xiàn)1. 基于 3 節(jié)點三角形單元的有限元編程實現(xiàn)（ 1）程序編寫說明Triangle2D3Node_Stiffness（E,NU,t,xi,yi,xj,yj,xm ,ym,ID）該函數(shù)計算單元的剛度矩陣，輸入模量 E,泊松比NU, 厚度t ,三個節(jié)點i , j , mi平面問題性質(zhì)指示參數(shù)（1為平 面應(yīng)力，2 為平面應(yīng)變），輸出單元剛度矩陣k（ 6 6） 。Triangle2D3Node_Assemble（KK,k,i,j,m）該函數(shù)進(jìn)行單元剛

22、度矩陣的組裝，輸入單元剛度矩陣k，單元的節(jié)點編號i , j , mi輸由整體剛度矩陣 K左Triangle2D3Node_Stress（E,NU,xi,yi,xj,yj,xm,ym,u,ID）該函數(shù)計算單元的應(yīng)力，輸入彈性模量 E,泊松比NU, 厚度t ,三個節(jié)點i , j , mi平面問題性質(zhì)指示參數(shù)（1為平 面應(yīng)力，2 為平面應(yīng)變），單元的位移列陣u（6 1），輸出單元的應(yīng)力，由于它為常應(yīng)力單元，則單元的應(yīng)力分量為Sx，Sy， Sxy。（ 2）程序清單%Triangle2D3Node%begin% functionk=Triangle2D3Node_Stiffness(E,NU,t,xi,

23、yi,xj,yj,xm,ym,ID)%該函數(shù)計算單元的剛度矩陣%俞入彈T生模量E、泊松比NU和厚度t%俞入3個節(jié)點i , j , m的坐標(biāo)xi,yi,xj,yj,xm,ym%輸入平面問題性質(zhì)指示參數(shù)ID( 1 位平面應(yīng)力，2為平面應(yīng)變)%輸入單元剛度矩陣k( 6*6)%A=(xi*(yj-ym)+xj(ym-yi)+xm*(yi-yj)/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=betai 0 betaj 0 betam 0;0 gammai 0 gammaj 0 gammam;

24、gammai betai gammaj betaj gammambetam/(2*A);if ID=1D=(E/(1-NU*NU)*1 NU 0;NU 1 0;0 0 (1-NU)/2;elseif ID=2D=(E/(1+NU)/(1-2*NU)*1-NU NU 0;NU 1-NU 0;00 (1-2*NU)/2;endk=t*A*B'*D*B;end%function z=Triangle2D3Node_Assemble(KK,k,i,j,m)%該函數(shù)進(jìn)行單元剛度矩陣的組裝%輸入單元剛度矩陣k%輸入單元的節(jié)點編號i ， j ， m%輸入整體剛度矩陣KKrf%DOF(1)=2*i-1

25、;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;DOF(5)=2*m-1;DOF(6)=2*m;for n1=1:6for n2=1:6KK(DOF(n1),DOF(n2)=KK(DOF(n1),DOF(n2)+k(n1,n2);endendz=KK;% functionstress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)%該函數(shù)計算單元的應(yīng)力%俞入彈T生模量E、泊松比NU和厚度t%輸入平面問題性質(zhì)指示參數(shù)ID( 1 位平面應(yīng)力，2為平面應(yīng)變) ，單元的位移列陣u( 6*1 )%輸出單元的應(yīng)力stress

26、( 3*1 ) ， 由于它為常應(yīng)力單元，則單元的應(yīng)力分量Sx， Sy， Sxy%A=(xi*xj(ym-yi)+xm*(yi-yj)/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=batai 0 bataj 0 betam 0;0 gammai 0 gammaj 0 gammam;gammai betai gammaj betaj gammambetam/(2*A);if ID=1D=(E/(1-NU*NU)*1 NU 0;NU 1 0;0 0 (1-NU)/2;elseif ID=

27、2D=(E/(1+NU)/(1-2*NU)*1-NU NU 0;NU 1-NU 0;00 (1-2*NU)/2;endstress=D*B*u;%> > E=1E7;> > NU=1/3;> > t=;>> c=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,0,1,2)> > CC=zeros(8,8);>> CC=Triangle2D3Node_Assemble(KK,k1,1,2,4) ；>> CC=Triangle2D3Node_Assemble(KK,k1,3,4,2

28、)>> k1=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,0,1,1)>> KK=zeros(8,8);>> KK=Triangle2D3Node_Assemble(KK,k1,1,2,4)>> KK=Triangle2D3Node_Assemble(KK,k1,3,4,2)>> k=KK(3:6,3:6);>> p=500;0;500;0;>> u=kp>>U=0;0;u;0;0>>P=KK*U>> A=*U'*KK*U-P&#

29、39;*U( 3) 計算結(jié)果應(yīng)變CC =+05 *0000000000000000位移U =00節(jié)點力P =00其中，節(jié)點1 的支反力為FRx1500N, FRy1 -176.4706N ，節(jié)點2 的支反力為FRx2500N,FRy2 176.4706N。勢能A =單元剛度陣KK =+05 *00000000000000002. 基于四節(jié)點四邊形單元的有限元建模及分析 （ 1）程序編寫說明Quad2D4Node_Stiffness（E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID）該函數(shù)計算單元的剛度矩陣，輸入模量E,泊松比NU,厚度h, 4個節(jié)點i , j , m p,平面

30、問題性質(zhì)指示參數(shù)ID （ 為平面應(yīng)力，2 為平面應(yīng)變），輸出單元剛度矩陣k（ 8 8） 。Quad2D4Node_Stress（E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID）該函數(shù)計算單元的應(yīng)力，輸入彈性模量E,泊松比NU, 厚度h, 4個節(jié)點i , j , mi p,平面問題性質(zhì)指示參數(shù)ID (1 為平面應(yīng)力，2為平面應(yīng)變),單元的位移列陣u( 8 1),輸 出單元的應(yīng)力，由于它為常應(yīng)力單元，則單元的應(yīng)力分量為Sx， Sy， Sxy。( 2)程序清單%Quad2D4Node%begin%functionk=Quad2D4Node_Stiffness(E,NU,h,xi,

31、yi,xj,yj,xm,ym,xp,yp,ID)%該函數(shù)計算單元的剛度矩陣%俞入彈T生模量E、泊松比NU和厚度h%輸入4 個節(jié)點 i ， j ， m， p 的坐標(biāo)xi,yi,xj,yj,xm,ym,xp,yp%輸入平面問題性質(zhì)指示參數(shù)ID( 1 位平面應(yīng)力，2為平面應(yīng)變)%輸入單元剛度矩陣k( 8*8)%syms s t;a=(yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s)/4;b=(yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t)/4;c=(xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t)/4;d=(xi*(s-1)+

32、xj*(-1-s)+xm*(1+s)+xp*(1-s)/4;B1=a*(t-1)/4-b*(s-1)/4 0;0 c*(s-1)/4-d*(t-1)/4;c*(-1+s)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4;B2=a*(1-t)/4-b*(-1-s)/4 0;0 c*(-1-s)/4-d*(1-t)/4;c*(-1-s)/4-d*(1-t)/4 a*(1-t)/4-b*(-1-s)/4;B3=a*(t+1)/4-b*(s+1)/4 0;0 c*(s+1)/4-d*(t+1)/4;c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4;B4=a

33、*(-t-1)/4-b*(1-s)/4 0;0 c*(1-s)/4-d*(-t-1)/4;c*(1-s)/4-d*(-t-1)/4 a*(-t-1)/4-b*(1-s)/4;Bfirst=B1 B2 B3 B4;Jfirst=0 1-t t-s s-1;t-1 0 s+1 -s-t;s-t -s-1 0 t+1;1-s s+t -t-1 0;J=xi xj xm xp*yi;yj;ym;yp/8;B=Bfirst/J;if ID=1D=(E/(1-NU*NU)*1 NU 0;NU 1 0;0 0 (1-NU)/2;elseif ID=2D=(E/(1+NU)/(1-2*NU)*1-NU NU

34、0;NU 1-NU 0;0 0 (1-2*NU)/2;endBD=J*transpose(B)*D*B;r=int(int(BD,t,-1,1),s,-1,1);z=h*r;k=double(z);end% functionstress=Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID)%該函數(shù)計算單元的應(yīng)力%俞入彈T生模量E、泊松比NU和厚度h%輸入平面問題性質(zhì)指示參數(shù)ID( 1 位平面應(yīng)力，2為平面應(yīng)變)%輸入單元的位移列陣u( 8*1 )%輸出單元的應(yīng)力stress ( 3*1 ) ， 由于它為常應(yīng)力單元，則單元的應(yīng)力分量Sx， S

35、y， Sxy% sym s t;a=(yi*(s-1)+yi*(-1-s)+ym*(1+S)+yp*(1-s)/4;b=(yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t)/4c=(xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t)/4d=(xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s)/4B1=a*(t-1)/4-b*(s-1)/4 0;0 c*(s-1)/4-d*(t-1)/4;c*(-1+s)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4;B2=a*(1-t)/4-b*(-1-s)/4 0;0 c*(-1-S)/4-d*(1-t)/4;c*(-1-