MATLAB有限元分析與應(yīng)用.ppt

《MATLAB有限元分析與應(yīng)用.ppt》由會員分享，可在線閱讀，更多相關(guān)《MATLAB有限元分析與應(yīng)用.ppt（55頁珍藏版）》請在裝配圖網(wǎng)上搜索。

1、2020/11/11 1 第三章 MATLAB有限元分析與應(yīng)用 3-1 彈簧元 3-2 線性桿元 3-3 二次桿元 3-4 平面桁架元 3-5 空間桁架元 3-6 梁元 2020/11/11 2 3-1 彈簧元 1、有限元方法的步驟： 離散化域 形成單剛矩陣 集成整體剛度矩陣 引入邊界條件 求解方程 后處理 2020/11/11 3 2、基本方程 3-1 彈簧元 彈簧元是總體和局部坐標(biāo)一致的一維有限單元 每個彈簧元有兩個節(jié)點（ node） kkk kk 單剛矩陣為： 22 總剛矩陣： nn K U F結(jié)構(gòu)方程：

2、 單元節(jié)點力： f k u 2020/11/11 4 3、 MATLAB函數(shù)編寫 3-1 彈簧元 %SpringElementStiffness This function returns the element stiffness %matrix for a spring with stiffness k. %The size of the element stiffness matrix is 2 x 2. 3.1 單元剛度矩陣的形成 y = k -k ; -k k; function y = SpringElementStiffness(k)

3、 2020/11/11 5 3、 MATLAB函數(shù)編寫 3-1 彈簧元 %SpringAssemble This function assembles the element stiffness % matrix k of the spring with nodes i and j into the % global stiffness matrix K. % This function returns the global stiffness matrix K % after the element stiffn

4、ess matrix k is assembled. 3.2 整體剛度矩陣的形成 K(i,i) = K(i,i) + k(1,1); K(i,j) = K(i,j) + k(1,2); K(j,i) = K(j,i) + k(2,1); K(j,j) = K(j,j) + k(2,2); y = K; function y = SpringAssemble(K,k,i,j) 2020/11/11 6 3、 MATLAB函數(shù)編寫 3-1 彈簧元 %SpringElementForces This function returns the element nodal force %

5、 vector given the element stiffness matrix k % and the element nodal displacement vector u. 3.3 節(jié)點載荷計算 y = k * u; function y = SpringElementForces(k,u) 2020/11/11 7 4、實例計算分析應(yīng)用 3-1 彈簧元 如圖所示二彈簧元結(jié)構(gòu)，假定 k1=100kN/m， k2=200kN/m， P=15kN。 求：系統(tǒng)的整體剛度矩陣； 節(jié)點 2、 3的位移； 節(jié)點 1的支

6、反力； 每個彈簧的內(nèi)力 解： 步驟 1：離散化域 2020/11/11 8 4、實例計算分析應(yīng)用 3-1 彈簧元 步驟 2：形成單元剛度矩陣 k1=SpringElementStiffness(100); k1 = 100 -100 -100 100 k2=SpringElementStiffness(200)； k2 = 200 -200 -200 200 調(diào)用 function y = SpringElementStiffness(k)函數(shù) 2020/11/11 9 4、實例計算分析應(yīng)用 3-1 彈簧元 步驟 3：集成整體

7、剛度矩陣 調(diào)用 function y = SpringAssemble(K,k,i,j)函數(shù) n=3; K = zeros(n,n); K = SpringAssemble(K,k1,1,2) K = 0 0 0 0 0 0 0 0 0 K = SpringAssemble(K,k2,2,3) K = 100 -100 0 -100 100 0 0 0 0 K = 100 -100 0 -100 300 -200 0 -200 200 2020/11/

8、11 10 4、實例計算分析應(yīng)用 3-1 彈簧元 步驟 4：引入邊界條件 11 22 33 1 0 0 1 0 0 0 1 0 0 3 0 0 2 0 0 0 2 0 0 2 0 0 UF UF UF 已知邊界條件： 1 2 30 , 0 , 1 5U F F 1 2 3 1 0 0 1 0 0 0 0 1 0 0 3 0 0 2 0 0 0 0 2 0 0 2 0 0 1 5 F U U 2020/11/11 11 5、實例計算分析應(yīng)用 3-1 彈簧元 步驟 5

9、：解方程 2 3 300 200 0 200 200 15 U U U=zeros(2,1); F=0;15; K = K(2:3,2:3); U=KF U=inv(K)*F K(1,:)=; K(:,1)=; U = 0.1500 0.2250 2020/11/11 12 5、實例計算分析應(yīng)用 2-1 彈簧元 步驟 6：后處理 U=0;U U = 0 0.1500 0.2250 F=K*U F = -15.0000 0.0000 15.0000 u1=U(1:2); f1=SpringElementForc

10、es(k1,u1); f1 = -15.0000 15.0000 u2=U(2:3); f2=SpringElementForces(k2,u2); f2 = -15.0000 15.0000 2020/11/11 13 5、實例計算分析應(yīng)用 3-1 彈簧元 k1=SpringElementStiffness(100); k2=SpringElementStiffness(200); n=3; K=zeros(n,n); K=SpringAssemble(K,k1,1,2); K=SpringAssemble(K,k2,2,3); U=zeros(2,1); F=0;1

11、5; K = K(2:3,2:3); KK=K; U=KF U=0;U; F=K*U; u1=U(1:2); f1=SpringElementForces(k1,u1) u2=U(2:3); f2=SpringElementForces(k2,u2) 2020/11/11 14 1、基本方程 3-2 線性桿元 線性桿元也是總體和局部坐標(biāo)一致的一維有限單元，用線性函數(shù)描述 每個線性桿元有兩個節(jié)點（ node） E A E A LLk E A E A LL 單剛矩陣為： 22 總剛矩陣： nn K U F結(jié)構(gòu)方程： 單元節(jié)點力： f k u 2020

12、/11/11 15 2、 MATLAB函數(shù)編寫 %LinearBarElementStiffness This function returns the element % stiffness matrix for a linear bar with % modulus of elasticity E, cross-sectional % area A, and length L. The size of the %

13、 element stiffness matrix is 2 x 2. 2.1 單元剛度矩陣的形成 y = E*A/L -E*A/L ; -E*A/L E*A/L; function y = LinearBarElementStiffness(E,A,L) 3-2 線性桿元 2020/11/11 16 2、 MATLAB函數(shù)編寫 %LinearBarAssemble This function assembles the element stiffness % matrix k of the linear bar with node

14、s i and j % into the global stiffness matrix K. % This function returns the global stiffness % matrix K after the element stiffness matrix % k is assembled. 2.2 整體剛度矩陣的形成 K(i,i) = K(i,i) + k(1,1); K(i,j) = K(i,j) + k(1,2); K(j,i

15、) = K(j,i) + k(2,1); K(j,j) = K(j,j) + k(2,2); y = K; function y =LinearBarAssemble(K,k,i,j) 3-2 線性桿元 2020/11/11 17 2、 MATLAB函數(shù)編寫 %LinearBarElementForces This function returns the element nodal % force vector given the element stiffness % matrix k an

16、d the element nodal % displacement vector u. 2.3 節(jié)點載荷計算 y = k * u; function y = LinearBarElementForces(k,u) 3-2 線性桿元 2020/11/11 18 2、 MATLAB函數(shù)編寫 %LinearBarElementStresses This function returns the element nodal % stress vector given the element stif

17、fness % matrix k, the element nodal displacement % vector u, and the cross-sectional area A. 2.4 節(jié)點應(yīng)力計算 y = k * u/A; function y = LinearBarElementStresses(k, u, A) 3-2 線性桿元 2020/11/11 19 3、實例計算分析應(yīng)用 如圖所示二線性桿元結(jié)構(gòu)，假定 E=210MPa， A=0.003m2， P=10kN，

18、 節(jié)點 3的右位移為 0.002m。 求：系統(tǒng)的整體剛度矩陣； 節(jié)點 2的位移； 節(jié)點 1、 3的支反力； 每個桿件的應(yīng)力 解： 步驟 1：離散化域 3-2 線性桿元 2020/11/11 20 3、實例計算分析應(yīng)用 步驟 2：形成單元剛度矩陣 k1=LinearBarElementStiffness(E,A,L1) k2=LinearBarElementStiffness(E,A,L2) 調(diào)用 function y = LinearBarElementStiffness(E,A,L)函數(shù) 3-2 線性桿元 2020/11/11

19、21 3、實例計算分析應(yīng)用 步驟 3：集成整體剛度矩陣 調(diào)用 function y = LinearBarAssemble(K,k,i,j)函數(shù) n=3; K = zeros(n,n) K = LinearBarAssemble (K,k1,1,2) K = 0 0 0 0 0 0 0 0 0 K = LinearBarAssemble (K,k2,2,3) 3-2 線性桿元 2020/11/11 22 3、實例計算分析應(yīng)用 步驟 4：引入邊界條件 已知邊界條件： 1 3 20 , 0 . 0 0 2 , 1 0U U F

20、 3-2 線性桿元 11 22 33 420 000 420 000 0 420 000 105 000 0 630 000 0 630 000 630 000 UF UF UF 1 2 3 420 000 420 000 0 0 420 000 105 000 0 630 000 10 0 630 000 630 000 0.00 2 F U F 2020/11/11 23 3、實例計算分析應(yīng)用 步驟 5：解方程 21 0 5 0 0 0 6 3 0 0 0 0

21、 0 .0 0 2 1 0U U=zeros(1,1); U3=0.002 F=-10; K = K(2,2) 105000 K0 = K(2,3); -630000 U=K(F-K0*U3) U =0.0012 3-2 線性桿元 2020/11/11 24 3、實例計算分析應(yīng)用 步驟 6：后處理 U=0;U;0.002 U = 0 0.0012 0.0002 F=K*U F = -500.0000 -10.0000 510.0000 u1=U(1:2); f1= LinearBarElementForces(k1,u1) sigma1

22、=LinearBarElementStresses(k1, u1, A) u2=U(2:3); f2= LinearBarElementForces(k2,u2) sigma2=LinearBarElementStresses(k2, u2, A) 3-2 線性桿元 2020/11/11 25 3、實例計算分析應(yīng)用 E=210E6; A=0.003; L1=1.5; L2=1; k1= LinearBarElementStiffness(E,A,L1); k2= LinearBarElementStiffness(E,A,L2); n=3; K = zeros(n,n); K =

23、 LinearBarAssemble (K,k1,1,2); K = LinearBarAssemble (K,k2,2,3); U=zeros(1,1); U3=0.002; F=-10; 3-2 線性桿元 KK=K; K=K(2,2); K0=K(2,3); U=K(F-K0*U3); U=0;U;U3; F=KK*U u1=U(1:2); f1= LinearBarElementForces(k1,u1) sigma1=LinearBarElementStresses(k1, u1, A) u2=U(2:3); f2= LinearBarElementForces(k2,u2)

24、 sigma2=LinearBarElementStresses(k2, u2, A) 2020/11/11 26 1、基本方程 3-3 二次桿元 二次桿元也是總體和局部坐標(biāo)一致的一維有限單元，用二次方程描述 每個線性桿元有三個節(jié)點（ node） 7 1 8 1 7 8 3 8 8 16 EAk L 單剛矩陣為： 33 總剛矩陣： nn K U F結(jié)構(gòu)方程： 單元節(jié)點力： f k u 31 2020/11/11 27 2、 MATLAB函數(shù)編寫 %QuadraticBarElementStiffness This function returns

25、 the element % stiffness matrix for a quadratic bar % with modulus of elasticity E, % cross-sectional area A, and length L. % The size of the element stiffness % matrix

26、is 3 x 3. 2.1 單元剛度矩陣的形成 y = E*A/(3*L)*7 1 -8 ; 1 7 -8 ; -8 -8 16; function y = QuadraticBarElementStiffness(E,A,L) 3-3 二次桿元 2020/11/11 28 2、 MATLAB函數(shù)編寫 %QuadraticBarAssemble This function assembles the element stiffness % matrix k of the quadratic bar with nodes i, j %

27、 and m into the global stiffness matrix K. % This function returns the global stiffness % matrix K after the element stiffness matrix % k is assembled. 2.2 整體剛度矩陣的形成 K(i,i) = K(i,i) + k(1,1); K(i,j) = K(i,j) + k(1,2); K(i,m

28、) = K(i,m) + k(1,3); K(j,i) = K(j,i) + k(2,1); K(j,j) = K(j,j) + k(2,2); function y =QuadraticBarAssemble(K,k,i,j,m) 3-3 二次桿元 K(j,m) = K(j,m) + k(2,3); K(m,i) = K(m,i) + k(3,1); K(m,j) = K(m,j) + k(3,2); K(m,m) = K(m,m) + k(3,3); y = K; 2020/11/11 29 2、 MATLAB函數(shù)編寫 %QuadraticBarElementForces T

29、his function returns the element nodal % force vector given the element stiffness % matrix k and the element nodal % displacement vector u. 2.3 節(jié)點載荷計算 y = k * u; function y = QuadraticBarElementForces(k,u) 3-3 二次桿元 2020

30、/11/11 30 2、 MATLAB函數(shù)編寫 %QuadraticBarElementStresses This function returns the element % nodal stress vector given the element % stiffness matrix k, the element nodal % displacement vector u, and the %

31、 cross-sectional area A. 2.4 節(jié)點應(yīng)力計算 y = k * u/A; function y = QuadraticBarElementStresses(k, u, A) 3-3 二次桿元 2020/11/11 31 3、實例計算分析應(yīng)用 如圖所示雙二次桿元結(jié)構(gòu)，假定 E=210MPa， A=0.003m2 求：系統(tǒng)的整體剛度矩陣； 節(jié)點 2、 3、 4、 5的位移； 節(jié)點 1的支反力； 每個桿件的應(yīng)力 解： 3-3 二次桿元 2020/11/11 32 3、實例計算分析應(yīng)用 E=210E6; A=

32、0.003; L=2; k1= QuadraticBarElementStiffness(E,A,L); k2= QuadraticBarElementStiffness(E,A,L); n=5; K = zeros(n,n); K =QuadraticBarAssemble(K,k1,1,3,2); K =QuadraticBarAssemble(K,k2,3,5,4); U=zeros(4,1); F=5;-10;-7;10; KK=K; K=K(2:n,2:n); U=KF; U=0;U; F=KK*U; u1=U(1);U(3);U(2); f1= QuadraticBarEleme

33、ntForces(k1,u1); sigma1=QuadraticBarElementStresses(k1, u1, A); u2=U(3);U(5);U(4); f2=QuadraticBarElementForces(k2,u2); sigma2=QuadraticBarElementStresses(k2, u2, A); 3-3 二次桿元 2020/11/11 33 1、基本方程 3-4 平面桁架元 平面桁架元是既有局部坐標(biāo)又有總體坐標(biāo)二維有限元，用線性函數(shù)描述 每個平面桁架元有二個節(jié)點（ node） 22 22 22 22 C CS C CS CS S C

34、S SEA k L C CS C CS CS S CS S 單剛矩陣為： 44 總剛矩陣： 22nn K U F結(jié)構(gòu)方程： 單元節(jié)點力： EAf C S C S uL 41 2020/11/11 34 2、 MATLAB函數(shù)編寫 %PlaneTrussElementLength This function returns the length of the % plane truss element whose first node has % coordinates (x1,y1) and se

35、cond node has % coordinates (x2,y2). 2.1 計算單元長度 y = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)); function y = PlaneTrussElementLength(x1,y1,x2,y2) 3-4 平面桁架元 2020/11/11 35 2、 MATLAB函數(shù)編寫 %PlaneTrussElementStiffness This function returns the element % stiffness matrix

36、for a plane truss % element with modulus of elasticity E, % cross-sectional area A, length L, and % angle theta (in degrees). % The size of the element stiffness % matrix is 4 x 4. 2.2 單元剛度矩陣的形成 x = theta*pi/180; C = cos(x);

37、S = sin(x); y = E*A/L*C*C C*S -C*C -C*S ; C*S S*S -C*S -S*S ; -C*C -C*S C*C C*S ; -C*S -S*S C*S S*S; function y = PlaneTrussElementStiffness(E,A,L, theta) 3-4 平面桁架元 2020/11/11 36 2、 MATLAB函數(shù)編寫 %PlaneTrussAssemble This function assembles the element stiffness % matrix k of the plan

38、e truss element with nodes % i and j into the global stiffness matrix K. % This function returns the global stiffness % matrix K after the element stiffness matrix k is assembled. 2.3 整體剛度矩陣的形成 K(2*i-1,2*i-1) = K(2*i-1,2*i-1) + k(1,1); K(2*i-1,2*i) = K(2*i-1,2*i) + k(

39、1,2); K(2*i-1,2*j-1) = K(2*i-1,2*j-1) + k(1,3); K(2*i-1,2*j) = K(2*i-1,2*j) + k(1,4); K(2*i,2*i-1) = K(2*i,2*i-1) + k(2,1); K(2*i,2*i) = K(2*i,2*i) + k(2,2); K(2*i,2*j-1) = K(2*i,2*j-1) + k(2,3); K(2*i,2*j) = K(2*i,2*j) + k(2,4); function y =PlaneTrussAssemble(K,k,i,j) K(2*j-1,2*i-1) = K(2*j-1,2*i-1

40、) + k(3,1); K(2*j-1,2*i) = K(2*j-1,2*i) + k(3,2); K(2*j-1,2*j-1) = K(2*j-1,2*j-1) + k(3,3); K(2*j-1,2*j) = K(2*j-1,2*j) + k(3,4); K(2*j,2*i-1) = K(2*j,2*i-1) + k(4,1); K(2*j,2*i) = K(2*j,2*i) + k(4,2); K(2*j,2*j-1) = K(2*j,2*j-1) + k(4,3); K(2*j,2*j) = K(2*j,2*j) + k(4,4); y = K; 3-4 平面桁架元 2020

41、/11/11 37 2、 MATLAB函數(shù)編寫 %PlaneTrussElementForce This function returns the element force % given the modulus of elasticity E, the % cross-sectional area A, the length L, % the angle theta (in degrees), and the % element nodal displacement vector u.

42、 2.4 節(jié)點載荷計算 x = theta * pi/180; C = cos(x); S = sin(x); y = E*A/L*-C -S C S* u; function y = PlaneTrussElementForce(E,A,L,theta,u) 3-4 平面桁架元 2020/11/11 38 2、 MATLAB函數(shù)編寫 %PlaneTrussElementStress This function returns the element stress % given the modulus of elasticity E, the %

43、 the length L, the angle theta (in % degrees), and the element nodal % displacement vector u. 2.5 節(jié)點應(yīng)力計算 x = theta * pi/180; C = cos(x); S = sin(x); y = E/L*-C -S C S* u; function y = PlaneTrussElementStress(E,L,theta,u) 3-4 平面桁架元 2020/11/11 39 3、實例計算分析應(yīng)用 如

44、圖所示平面桁架結(jié)構(gòu)，假定 E=210MPa， A=0.0004m2 求：系統(tǒng)的整體剛度矩陣； 節(jié)點 2的水平位移； 節(jié)點 3的水平豎向位移； 節(jié)點 1、 2的支反力； 每跟桿件的應(yīng)力 3-4 平面桁架元 2020/11/11 40 1、基本方程 3-5 空間桁架元 空間桁架元是既有局部坐標(biāo)又有總體坐標(biāo)三維有限元，用線性函數(shù)描 述。 各單元之間通過鉸接系統(tǒng)連接，只能傳遞力，而不能傳遞彎矩 每個桁架元有二個節(jié)點（ node） c o s , c o s , c o sx x y y z zC C C 2020/1

45、1/11 41 1、基本方程 3-5 空間桁架元 總剛矩陣： 33nn K U F結(jié)構(gòu)方程： 單元節(jié)點力： 61 x y z x y zEAf C C C C C C uL 22 22 22 22 22 22 x x y x z x x y x z y x y y z y x y y z z x y z z y z y z z x x y x z x x y x z y x y y z x y y y z z x y z z x z y z z C C C C C C C C C C C C C C C C C C C C C C C C C C C C

46、 C CEA k C C C C C C C C C CL C C C C C C C C C C C C C C C C C C C C 66 單剛矩陣為： 2020/11/11 42 2、 MATLAB函數(shù)編寫 %SpaceTrussElementLength This function returns the length of the % space truss element whose first node has % coordinates (x1,y1,z1) and second

47、node has % coordinates (x2,y2,z2). 2.1 計算單元長度 y = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1)); function y = SpaceTrussElementLength(x1,y1,z1,x2,y2,z2) 3-5 空間桁架元 2020/11/11 43 2、 MATLAB函數(shù)編寫 %SpaceTrussElementStiffness This function returns the element %

48、 stiffness matrix for a space truss % element with modulus of elasticity E, % cross-sectional area A, length L, and % angles thetax, thetay, thetaz % (in degrees). The size of the eleme

49、nt % stiffness matrix is 6 x 6. 2.2 單元剛度矩陣的形成 x = thetax*pi/180; u = thetay*pi/180; v = thetaz*pi/180; Cx = cos(x); Cy = cos(u); Cz = cos(v); w = Cx*Cx Cx*Cy Cx*Cz ; Cy*Cx Cy*Cy Cy*Cz ; Cz*Cx Cz*Cy Cz*Cz; y = E*A/L*w -w ; -w w; function y = SpaceTrussElementStiff

50、ness(E,A,L,thetax,thetay,thetaz) 3-5 空間桁架元 2020/11/11 44 2、 MATLAB函數(shù)編寫 %SpaceTrussAssemble This function assembles the element stiffness % matrix k of the space truss element with nodes % i and j into the global stiffness matrix K. % This function returns the glo

51、bal stiffness % matrix K after the element stiffness matrix % k is assembled. 2.3 整體剛度矩陣的形成 K(3*i-2,3*i-2) = K(3*i-2,3*i-2) + k(1,1); K(3*i-2,3*i-1) = K(3*i-2,3*i-1) + k(1,2); K(3*i-2,3*i) = K(3*i-2,3*i) + k(1,3); K(3*i-2,3*j-2) = K(3*i-2,3*j-2) + k(1,4); K(3*i-2,3*j-1) = K(3*

52、i-2,3*j-1) + k(1,5); K(3*i-2,3*j) = K(3*i-2,3*j) + k(1,6); K(3*i-1,3*i-2) = K(3*i-1,3*i-2) + k(2,1); K(3*i-1,3*i-1) = K(3*i-1,3*i-1) + k(2,2); K(3*i-1,3*i) = K(3*i-1,3*i) + k(2,3); K(3*i-1,3*j-2) = K(3*i-1,3*j-2) + k(2,4); K(3*i-1,3*j-1) = K(3*i-1,3*j-1) + k(2,5); K(3*i-1,3*j) = K(3*i-1,3*j) + k(2,6

53、); function y =SpaceTrussAssemble(K,k,i,j) 3-5 空間桁架元 2020/11/11 45 2、 MATLAB函數(shù)編寫 2.3 整體剛度矩陣的形成 3-5 空間桁架元 K(3*j-1,3*i-2) = K(3*j-1,3*i-2) + k(5,1); K(3*j-1,3*i-1) = K(3*j-1,3*i-1) + k(5,2); K(3*j-1,3*i) = K(3*j-1,3*i) + k(5,3); K(3*j-1,3*j-2) = K(3*j-1,3*j-2) + k(5,4); K(3*j-1,3*j-1) = K(

54、3*j-1,3*j-1) + k(5,5); K(3*j-1,3*j) = K(3*j-1,3*j) + k(5,6); K(3*j,3*i-2) = K(3*j,3*i-2) + k(6,1); K(3*j,3*i-1) = K(3*j,3*i-1) + k(6,2); K(3*j,3*i) = K(3*j,3*i) + k(6,3); K(3*j,3*j-2) = K(3*j,3*j-2) + k(6,4); K(3*j,3*j-1) = K(3*j,3*j-1) + k(6,5); K(3*j,3*j) = K(3*j,3*j) + k(6,6); y = K; K(3*i,3*i-2)

55、 = K(3*i,3*i-2) + k(3,1); K(3*i,3*i-1) = K(3*i,3*i-1) + k(3,2); K(3*i,3*i) = K(3*i,3*i) + k(3,3); K(3*i,3*j-2) = K(3*i,3*j-2) + k(3,4); K(3*i,3*j-1) = K(3*i,3*j-1) + k(3,5); K(3*i,3*j) = K(3*i,3*j) + k(3,6); K(3*j-2,3*i-2) = K(3*j-2,3*i-2) + k(4,1); K(3*j-2,3*i-1) = K(3*j-2,3*i-1) + k(4,2); K(3*j-2,

56、3*i) = K(3*j-2,3*i) + k(4,3); K(3*j-2,3*j-2) = K(3*j-2,3*j-2) + k(4,4); K(3*j-2,3*j-1) = K(3*j-2,3*j-1) + k(4,5); K(3*j-2,3*j) = K(3*j-2,3*j) + k(4,6); 2020/11/11 46 2、 MATLAB函數(shù)編寫 %SpaceTrussElementForce This function returns the element force % given the modulus of elasticity E, the

57、% cross-sectional area A, the length L, % the angles thetax, thetay, thetaz % (in degrees), and the element nodal % displacement vector u. 2.4 節(jié)點載荷計算 x = thetax * pi/180; w = thetay * pi/180; v = thetaz * pi/180; Cx = cos(x); Cy = cos(w); Cz = cos(v

58、); y = E*A/L*-Cx -Cy -Cz Cx Cy Cz*u; function y = SpaceTrussElementForce(E,A,L,thetax,thetay,thetaz,u) 3-5 空間桁架元 2020/11/11 47 2、 MATLAB函數(shù)編寫 %SpaceTrussElementStress This function returns the element stress % given the modulus of elasticity E, the % length L, the angl

59、es thetax, thetay, % thetaz (in degrees), and the element % nodal displacement vector u. 2.5 節(jié)點應(yīng)力計算 x = thetax * pi/180; w = thetay * pi/180; v = thetaz * pi/180; Cx = cos(x); Cy = cos(w); Cz = cos(v); y = E/L*-Cx -Cy -Cz Cx Cy Cz*u; function y = SpaceTrussElementStress(E,

60、L,thetax,thetay,thetaz,u) 3-5 空間桁架元 2020/11/11 48 3、實例計算分析應(yīng)用 如圖所示空間桁架結(jié)構(gòu)，假定 E=210MPa， A14=0.001m2 A24=0.002m2， A34=0.001m2， P=12kN 求：系統(tǒng)的整體剛度矩陣； 節(jié)點 4的水平位移； 節(jié)點 3的水平豎向位移； 節(jié)點 1、 2、 3的支反力； 每跟桿件的應(yīng)力 3-5 空間桁架元 2020/11/11 49 1、基本方程 3-6 梁元 梁元是總體坐標(biāo)與局部坐標(biāo)一致的二維有限元，用線性函數(shù)描

61、 述。 各單元之間通過鉸接系統(tǒng)連接，只能傳遞力，而不能傳遞彎矩 每個梁元有二個節(jié)點（ node） 單剛矩陣為： 22 3 22 12 6 12 6 6 4 6 2 12 6 12 6 6 2 6 4 LL L L L LEI k LLL L L L L 44 總剛矩陣： 22nn K U F結(jié)構(gòu)方程： 單元節(jié)點力： 41 f k u 2020/11/11 50 2、 MATLAB函數(shù)編寫 %BeamElementStiffness This function returns the element %

62、 stiffness matrix for a beam % element with modulus of elasticity E, % moment of inertia I, and length L. % The size of the element stiffness % matrix is 4 x 4. 2.1單元剛度矩陣的形成 y = E*I/(L*L*L)*12 6*L -12 6*L ; 6*L 4*L*L -6

63、*L 2*L*L ; -12 -6*L 12 -6*L ; 6*L 2*L*L -6*L 4*L*L; function y = BeamElementStiffness(E,I,L) 3-6 梁元 2020/11/11 51 2、 MATLAB函數(shù)編寫 %BeamAssemble This function assembles the element stiffness % matrix k of the beam element with nodes % i and j into the global stiffness matrix K. %

64、 This function returns the global stiffness % matrix K after the element stiffness matrix % k is assembled. 2.2 整體剛度矩陣的形成 K(2*i-1,2*i-1) = K(2*i-1,2*i-1) + k(1,1); K(2*i-1,2*i) = K(2*i-1,2*i) + k(1,2); K(2*i-1,2*j-1) = K(2*i-1,2*j-1) + k(1,3); K(2*i-1,2*j) = K(2*i-1,2*j) + k(

65、1,4); K(2*i,2*i-1) = K(2*i,2*i-1) + k(2,1); K(2*i,2*i) = K(2*i,2*i) + k(2,2); K(2*i,2*j-1) = K(2*i,2*j-1) + k(2,3); K(2*i,2*j) = K(2*i,2*j) + k(2,4); function y =BeamAssemble(K,k,i,j) 3-6 梁元 K(2*j-1,2*i-1) = K(2*j-1,2*i-1) + k(3,1); K(2*j-1,2*i) = K(2*j-1,2*i) + k(3,2); K(2*j-1,2*j-1) = K(2*j-1

66、,2*j-1) + k(3,3); K(2*j-1,2*j) = K(2*j-1,2*j) + k(3,4); K(2*j,2*i-1) = K(2*j,2*i-1) + k(4,1); K(2*j,2*i) = K(2*j,2*i) + k(4,2); K(2*j,2*j-1) = K(2*j,2*j-1) + k(4,3); K(2*j,2*j) = K(2*j,2*j) + k(4,4); y = K; 2020/11/11 52 2、 MATLAB函數(shù)編寫 %BeamElementForces This function returns the element nodal force % vector given the element stiffness matrix k % and the element nodal displacement vector u. 2.4 節(jié)點載荷計算 y = k * u; function y = BeamElementForces(k,u) 3-6 梁元 2020/11/11 53 2