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

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

2、結(jié)論。二、有限元建模及分析1、基于3節(jié)點(diǎn)三角形單元的有限元建模及分析（1）結(jié)構(gòu)的離散化與編號(hào)如圖2所示，將平面矩形結(jié)構(gòu)分為兩個(gè)3節(jié)點(diǎn)三角形單元。單元三個(gè)節(jié)點(diǎn)的編號(hào)為1，2，4，單元三個(gè)節(jié)點(diǎn)的編號(hào)為3，4，2，各個(gè)節(jié)點(diǎn)的位置坐標(biāo)為，各個(gè)節(jié)點(diǎn)的位移（分別沿方向和方向）為。Xy1432圖2 方案一：使用兩個(gè)3節(jié)點(diǎn)三角形單元（2）各單元的剛度矩陣及剛度方程 a.單元的幾何和節(jié)點(diǎn)描述單元有6個(gè)節(jié)點(diǎn)位移自由度（DOF）。將所有節(jié)點(diǎn)上的位移組成一個(gè)列陣，記作；同樣，將所有節(jié)點(diǎn)上的各個(gè)力也組成一個(gè)列陣，記作，則有 同理，對(duì)于單元，有b.單元的位移場(chǎng)描述對(duì)于單元，設(shè)位移函數(shù) (1-1)由節(jié)點(diǎn)條件，在處，有 (1

3、-2) 將式(1-1)代入節(jié)點(diǎn)條件式(1-2)中，可求出式（1-1）中待定系數(shù)，即 (1-3) (1-4) (1-5)(1-6)(1-7)(1-8)在式(1-3)式(1-8)中 (1-9) (1-10)上式中的符號(hào)(1,2,3)表示下標(biāo)輪換，如同時(shí)更換。將單元各節(jié)點(diǎn)的位置坐標(biāo)代入得 將系數(shù)式(1-3)式(1-8)式代入(1-1)中，重寫(xiě)位移函數(shù)，并以節(jié)點(diǎn)位移的形式進(jìn)行表示，有 (1-11)令，則有形狀函數(shù)矩陣 (1-12)位移函數(shù)式(1-11)寫(xiě)成矩陣形式，有 (1-13)對(duì)于單元，過(guò)程同上，有形狀函數(shù)矩陣 (1-14)位移函數(shù) (1-15)c.單元的應(yīng)變場(chǎng)描述對(duì)于單元，應(yīng)變函數(shù) (1-16)

4、其中幾何矩陣 (1-17)對(duì)于單元，應(yīng)變方程 (1-18)其中幾何矩陣 (1-19)d.單元的應(yīng)力場(chǎng)描述 (1-20)其中，彈性系數(shù)矩陣D(1)為(1-21) (1-22)其中應(yīng)力函數(shù)矩陣S(1)為 (1-23)應(yīng)力方程為 (1-24)對(duì)于單元，過(guò)程同上。彈性系數(shù)矩陣D(2)為 (1-25)應(yīng)力函數(shù)矩陣S(2)為 (1-26)應(yīng)力方程為 (1-27)e.單元的勢(shì)能表達(dá)是單元?jiǎng)偠染仃?，?(1-28)其中薄板厚度。將式(1-17)、式(1-21)代入式(1-29)，得到單元的剛度陣計(jì)算得 同理，得到單元的剛度陣為將兩個(gè)單元按節(jié)點(diǎn)位移所對(duì)應(yīng)的位置進(jìn)行組裝，得到總體剛度矩陣為節(jié)點(diǎn)力 系統(tǒng)的勢(shì)能(計(jì)算

5、結(jié)果在下面呈現(xiàn))(4)邊界條件的處理及方程求解邊界條件為。因此，將針對(duì)節(jié)點(diǎn)2和節(jié)點(diǎn)3的位移求解，節(jié)點(diǎn)2和節(jié)點(diǎn)3對(duì)應(yīng)總體剛度陣KK中的第3行到第6行、第3列到第6列，則需從KK中提出，置給k，然后生成對(duì)應(yīng)的載荷列陣p，再采用高斯消去法進(jìn)行求解。>> k=KK(3:6,3:6);>> p=500;0;500;0;>> u=kpu =0.0009 0.0001 0.0010 -0.0002 將列排成了行再計(jì)算支反力。在得到整個(gè)結(jié)構(gòu)的節(jié)點(diǎn)位移后，由原整體剛度方程就可以計(jì)算出對(duì)應(yīng)的支反力；先將上面得到的位移結(jié)果與邊界條件的節(jié)點(diǎn)位移進(jìn)行組合，得到整體的位移列陣U，再代回

6、原整體剛度方程，計(jì)算出所有的節(jié)點(diǎn)力，按照位置關(guān)系找出對(duì)應(yīng)的支反力。>> U=0;0;u;0;0U = 0 0 0.0009 0.0001 0.0010 -0.0002 0 0將列排成了行>> P=KK*UP =-500 -176.4706 500 0 500 0 -500 176.4706將列排成了行所以，節(jié)點(diǎn)1的支反力為，節(jié)點(diǎn)2的支反力為。根據(jù)已求得的位移和支反力計(jì)算系統(tǒng)的勢(shì)能。>> A=0.5*U'*KK*U-P'*UA = -0.4706(5）結(jié)果分析上述支反力計(jì)算結(jié)果滿足靜力平衡，驗(yàn)證了以上求解過(guò)程及MATLAB算法的正確性。2、基于

7、四節(jié)點(diǎn)四邊形單元的有限元建模及分析（1）結(jié)構(gòu)的離散化與編號(hào)如圖3所示一個(gè)4節(jié)點(diǎn)矩形單元，單元的節(jié)點(diǎn)位移共有8個(gè)自由度（DOF）。節(jié)點(diǎn)編號(hào)為1,2,3,4，各自的位置坐標(biāo)為，各個(gè)節(jié)點(diǎn)的位移（分別沿方向和方向）為。Xy4321圖3 方案二：使用一個(gè)4節(jié)點(diǎn)矩形單元（2）局部坐標(biāo)系下單元的描述a.單元的幾何和節(jié)點(diǎn)描述采用無(wú)量綱坐標(biāo) 其中。則單元四個(gè)節(jié)點(diǎn)的幾何位置為 將所有節(jié)點(diǎn)上的位移組成一個(gè)列陣，記作q；同樣，將所有節(jié)點(diǎn)上的各個(gè)力也組成一個(gè)列陣，記作F，則有 b.單元的位移場(chǎng)描述設(shè)位移函數(shù)為 由節(jié)點(diǎn)條件，在處，有 將位移試函數(shù)代入節(jié)點(diǎn)條件中，求出待定系數(shù)和，再代入位移函數(shù)中，整理后得 其中 如以無(wú)量綱

8、坐標(biāo)來(lái)表達(dá)，可寫(xiě)成將帶入上式得到形狀函數(shù)矩陣寫(xiě)成矩陣形式，有 c.單元的應(yīng)變場(chǎng)描述單元應(yīng)變?yōu)?其中幾何矩陣為 d.單元的應(yīng)力場(chǎng)描述應(yīng)力表達(dá)式為其中，應(yīng)力函數(shù)矩陣。e.單元的勢(shì)能表達(dá)以上已將單元的三大基本變量用基于節(jié)點(diǎn)的位移列陣來(lái)表示；將其代入單元?jiǎng)菽鼙磉_(dá)式中，有，其中為4節(jié)點(diǎn)矩形單元的剛度矩陣，即其中，t為薄板的厚度，上式的各個(gè)字塊矩陣為 f.單元?jiǎng)偠汝嚰皠偠确匠虇卧獎(jiǎng)偠汝囋谏厦嬉呀?jīng)列出。將單元的勢(shì)能對(duì)節(jié)點(diǎn)位移取一階極值，可得到單元的剛度方程(4)邊界條件的處理及方程求解 處理方法與3節(jié)點(diǎn)三角形單元一致，利用上述求解程序具有的可移植性，簡(jiǎn)化了求解過(guò)程。>> k=K(3:6,3:6)

9、;>> p=500;0;500;0;>> u=kpu =1.0e-03 *0.4815 0.1111 0.4815 -0.1111 將列排成了行再計(jì)算支反力。同樣注意按照位置關(guān)系找出對(duì)應(yīng)的支反力。>> U=0;0;u;0;0U =1.0e-03 * 0 0 0.4815 0.1111 0.4815 -0.1111 0 0將列排成了行>> P=K*UP =-500 -111.1111 500 0 500 0 -500 111.1111將列排成了行所以，節(jié)點(diǎn)1的支反力為，節(jié)點(diǎn)2的支反力為。根據(jù)已求得的位移和支反力計(jì)算系統(tǒng)的勢(shì)能。>> A=

10、0.5*U'*K*U-P'*UA = -0.2407(5）結(jié)果分析基于4節(jié)點(diǎn)矩形單元計(jì)算出的勢(shì)能小于基于3節(jié)點(diǎn)三角形單元計(jì)算出的結(jié)果，若將該系統(tǒng)分為更多的單元，計(jì)算精度也將提高。3.兩種方案的比較與分析從以上計(jì)算可以看出，用三角形單元計(jì)算時(shí)，由于形函數(shù)是完全一次式，因而其應(yīng)變場(chǎng)在單元內(nèi)均為常數(shù)；而對(duì)于四邊形單元，其形函數(shù)帶有二次式，計(jì)算得到的應(yīng)變場(chǎng)和應(yīng)力場(chǎng)都是坐標(biāo)的一次函數(shù)，但不是完全的一次函數(shù)，對(duì)提高精度有一定作用；根據(jù)最小勢(shì)能原理，勢(shì)能越小，則整體計(jì)算精度越高，比較兩種單元計(jì)算得到的系統(tǒng)勢(shì)能，可看出，在相同的節(jié)點(diǎn)自由度情況下，矩形單元的計(jì)算精度要比三角形單元高。三、基于MA

11、TLAB的編程實(shí)現(xiàn)1. 基于3節(jié)點(diǎn)三角形單元的有限元編程實(shí)現(xiàn) （1）程序編寫(xiě)說(shuō)明Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)該函數(shù)計(jì)算單元的剛度矩陣，輸入模量E，泊松比NU，厚度t，三個(gè)節(jié)點(diǎn)i，j，m，平面問(wèn)題性質(zhì)指示參數(shù)（1為平面應(yīng)力，2為平面應(yīng)變），輸出單元?jiǎng)偠染仃噆（66）。Triangle2D3Node_Assemble(KK,k,i,j,m)該函數(shù)進(jìn)行單元?jiǎng)偠染仃嚨慕M裝，輸入單元?jiǎng)偠染仃噆，單元的節(jié)點(diǎn)編號(hào)i，j，m，輸出整體剛度矩陣KK。Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm

12、,ym,u,ID)該函數(shù)計(jì)算單元的應(yīng)力，輸入彈性模量E，泊松比NU，厚度t，三個(gè)節(jié)點(diǎn)i，j，m，平面問(wèn)題性質(zhì)指示參數(shù)（1為平面應(yīng)力，2為平面應(yīng)變），單元的位移列陣u(61)，輸出單元的應(yīng)力，由于它為常應(yīng)力單元，則單元的應(yīng)力分量為Sx，Sy，Sxy。（2）程序清單%Triangle2D3Node%begin%function k=Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)%該函數(shù)計(jì)算單元的剛度矩陣%輸入彈性模量E、泊松比NU和厚度t%輸入3個(gè)節(jié)點(diǎn)i，j，m的坐標(biāo)xi,yi,xj,yj,xm,ym%輸入平面問(wèn)題性質(zhì)指示參數(shù)ID（1位

13、平面應(yīng)力，2為平面應(yīng)變）%輸入單元?jiǎng)偠染仃噆（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; gammai betai gammaj betaj gammam betam/(2*A);if ID=1 D=(E/(1-NU*NU)*1 NU 0;NU 1 0;0 0 (1-NU)/2;elseif ID

14、=2 D=(E/(1+NU)/(1-2*NU)*1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2;endk=t*A*B'*D*B;end%function z=Triangle2D3Node_Assemble(KK,k,i,j,m)%該函數(shù)進(jìn)行單元?jiǎng)偠染仃嚨慕M裝%輸入單元?jiǎng)偠染仃噆%輸入單元的節(jié)點(diǎn)編號(hào)i，j，m%輸入整體剛度矩陣KKrf%-DOF(1)=2*i-1;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:6 for n2=1:6 KK(DOF(n1),DOF(n2)=KK

15、(DOF(n1),DOF(n2)+k(n1,n2); endendz=KK;%function stress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)%該函數(shù)計(jì)算單元的應(yīng)力%輸入彈性模量E、泊松比NU和厚度t%輸入平面問(wèn)題性質(zhì)指示參數(shù)ID（1位平面應(yīng)力，2為平面應(yīng)變），單元的位移列陣u（6*1）%輸出單元的應(yīng)力stress（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=

16、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 gammam betam/(2*A);if ID=1 D=(E/(1-NU*NU)*1 NU 0;NU 1 0;0 0 (1-NU)/2;elseif ID=2 D=(E/(1+NU)/(1-2*NU)*1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2;endstress=D*B*u;%>> E=1E7;>> NU=1/3;&g

17、t;> t=0.1;>> 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)>> k1=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,0,1,1)>> KK=zeros(8,8);>> KK=Triangle2D3No

18、de_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=0.5*U'*KK*U-P'*U（3)計(jì)算結(jié)果應(yīng)變CC = 1.0e+05 * 9.3750 5.6250 -7.5000 -1.8750 0 0 -1.8750 -3.7500 5.6250 9.3750 -3.7500 -1.

19、8750 0 0 -1.8750 -7.5000 -7.5000 -3.7500 9.3750 0 -1.8750 -1.8750 0 5.6250 -1.8750 -1.8750 0 9.3750 -3.7500 -7.5000 5.6250 0 0 0 -1.8750 -3.7500 9.3750 5.6250 -7.5000 -1.8750 0 0 -1.8750 -7.5000 5.6250 9.3750 -3.7500 -1.8750 -1.8750 -1.8750 0 5.6250 -7.5000 -3.7500 9.3750 0 -3.7500 -7.5000 5.6250 0

20、-1.8750 -1.8750 0 9.3750位移U = 0 0 0.0009 0.0001 0.0010 -0.0002 0 0節(jié)點(diǎn)力P = -500.0000 -176.4706 500.0000 0 500.0000 0 -500.0000 176.4706其中，節(jié)點(diǎn)1的支反力為，節(jié)點(diǎn)2的支反力為。勢(shì)能A = -0.4706單元?jiǎng)偠汝嘖K = 1.0e+05 * 7.5000 3.7500 -5.6250 -1.8750 0 0 -1.8750 -1.8750 3.7500 7.5000 -1.8750 -1.8750 0 0 -1.8750 -5.6250 -5.6250 -1.87

21、50 7.5000 0 -1.8750 -1.8750 0 3.7500 -1.8750 -1.8750 0 7.5000 -1.8750 -5.6250 3.7500 0 0 0 -1.8750 -1.8750 7.5000 3.7500 -5.6250 -1.8750 0 0 -1.8750 -5.6250 3.7500 7.5000 -1.8750 -1.8750 -1.8750 -1.8750 0 3.7500 -5.6250 -1.8750 7.5000 0 -1.8750 -5.6250 3.7500 0 -1.8750 -1.8750 0 7.50002. 基于四節(jié)點(diǎn)四邊形單元的

22、有限元建模及分析（1）程序編寫(xiě)說(shuō)明Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID)該函數(shù)計(jì)算單元的剛度矩陣，輸入模量E，泊松比NU，厚度h，4個(gè)節(jié)點(diǎn)i，j，m，p，平面問(wèn)題性質(zhì)指示參數(shù)ID（1為平面應(yīng)力，2為平面應(yīng)變），輸出單元?jiǎng)偠染仃噆（88）。Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID)該函數(shù)計(jì)算單元的應(yīng)力，輸入彈性模量E，泊松比NU，厚度h，4個(gè)節(jié)點(diǎn)i，j，m，p，平面問(wèn)題性質(zhì)指示參數(shù)ID（1為平面應(yīng)力，2為平面應(yīng)變），單元的位移列陣u(1)，輸出單元的應(yīng)力，由于

23、它為常應(yīng)力單元，則單元的應(yīng)力分量為Sx，Sy，Sxy。（2）程序清單%Quad2D4Node%begin%function k=Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID)%該函數(shù)計(jì)算單元的剛度矩陣%輸入彈性模量E、泊松比NU和厚度h%輸入4個(gè)節(jié)點(diǎn)i，j，m，p的坐標(biāo)xi,yi,xj,yj,xm,ym,xp,yp%輸入平面問(wèn)題性質(zhì)指示參數(shù)ID（1位平面應(yīng)力，2為平面應(yīng)變）%輸入單元?jiǎng)偠染仃噆（8*8)%-syms s t;a=(yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s)/4;b=(yi*(t-1)

24、+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)+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

25、*(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*(-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

26、;if ID=1 D=(E/(1-NU*NU)*1 NU 0;NU 1 0;0 0 (1-NU)/2;elseif ID=2 D=(E/(1+NU)/(1-2*NU)*1-NU NU 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%function stress=Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID)%該函數(shù)計(jì)算單元的應(yīng)力%輸入彈性模量E、泊松比NU和厚度h%輸入平面問(wèn)題

27、性質(zhì)指示參數(shù)ID（1位平面應(yīng)力，2為平面應(yīng)變）%輸入單元的位移列陣u（8*1）%輸出單元的應(yīng)力stress（3*1），由于它為常應(yīng)力單元，則單元的應(yīng)力分量Sx，Sy，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*

28、(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*(-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

29、 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=1 D=(E/(1-NU*NU)*1 NU 0;NU 1 0;0 0 (1-NU)/2;elseif ID=2 D=(E/(1+NU)/(1-2*NU)*1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2;endstr1=D*B*u;str2=subs(str1,s,t,0,

30、0);stress=double(str2);end>>=Quad2D4Node_Stiffness(E,NU,h,0,0,1,0,1,1,0,1，)>> K=Quad2D4Node_Stiffness(E,NU,h,0,0,1,0,1,1,0,1,1)>> k=K(3:6,3:6);>> p=500;0;500;0;>> u=kp>> U=0;0;u;0;0>> P=K*U>> A=0.5*U'*K*U-P'*U>>stress=Quad2D4Node_Stress(E

31、,NU,0,0,1,0,1,1,0,1,10(-3)*0;0;0.4815;0.1111;0.4815;-0.1111;0;0,1)（3)計(jì)算結(jié)果應(yīng)變C = 1.0e+06 * 1.2500 0.5625 -0.8750 0.1875 -0.6250 -0.5625 0.2500 -0.1875 0.5625 1.2500 -0.1875 0.2500 -0.5625 -0.6250 0.1875 -0.8750 -0.8750 -0.1875 1.2500 -0.5625 0.2500 0.1875 -0.6250 0.5625 0.1875 0.2500 -0.5625 1.2500 -0.1875 -0.8750 0.5625 -0.6250 -0.6250 -0.5625 0.2500 -0.1875 1.2500 0.5625 -0.8750 0.1875 -0.5625 -0.6250 0.1875 -0.87