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

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

2、3、MATLAB函數(shù)編寫(xiě),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 單元?jiǎng)偠染仃嚨男纬?y = k -k ; -k k;,function y = SpringElementStiffness(k),2020/9/7,5,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、MATLAB函數(shù)編寫(xiě),3-1 彈簧元,%Sprin

3、gAssemble 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 stiffness matrix k is assembled.,3.2 整體剛度矩陣的形成,K(i,i) = K(i,i) + k(1,1); K(i,j) = K(i,j)

4、 + 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/9/7,6,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、MATLAB函數(shù)編寫(xiě),3-1 彈簧元,%SpringElementForces This function returns the element nodal force % vector given the element stiffness matrix k % and the element nodal displacement

5、 vector u.,3.3 節(jié)點(diǎn)載荷計(jì)算,y = k * u;,function y = SpringElementForces(k,u),2020/9/7,7,結(jié)構(gòu)分析編程及軟件應(yīng)用,4、實(shí)例計(jì)算分析應(yīng)用,3-1 彈簧元,如圖所示二彈簧元結(jié)構(gòu)，假定k1=100kN/m，k2=200kN/m，P=15kN。 求：系統(tǒng)的整體剛度矩陣； 節(jié)點(diǎn)2、3的位移； 節(jié)點(diǎn)1的支反力； 每個(gè)彈簧的內(nèi)力,解：,步驟1：離散化域,2020/9/7,8,結(jié)構(gòu)分析編程及軟件應(yīng)用,4、實(shí)例計(jì)算分析應(yīng)用,3-1 彈簧元,步驟2：形成單元?jiǎng)偠染仃?k1=SpringElementStiffness(100);,k1 =

6、100 -100 -100 100,k2=SpringElementStiffness(200)；,k2 = 200 -200 -200 200,調(diào)用 function y = SpringElementStiffness(k)函數(shù),2020/9/7,9,結(jié)構(gòu)分析編程及軟件應(yīng)用,4、實(shí)例計(jì)算分析應(yīng)用,3-1 彈簧元,步驟3：集成整體剛度矩陣,調(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 = SpringAsse

7、mble(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/9/7,10,結(jié)構(gòu)分析編程及軟件應(yīng)用,4、實(shí)例計(jì)算分析應(yīng)用,3-1 彈簧元,步驟4：引入邊界條件,已知邊界條件：,2020/9/7,11,結(jié)構(gòu)分析編程及軟件應(yīng)用,5、實(shí)例計(jì)算分析應(yīng)用,3-1 彈簧元,步驟5：解方程,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/9/7,1

8、2,結(jié)構(gòu)分析編程及軟件應(yīng)用,5、實(shí)例計(jì)算分析應(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=SpringElementForces(k1,u1);,f1 = -15.0000 15.0000,u2=U(2:3); f2=SpringElementForces(k2,u2);,f2 = -15.0000 15.0000,2020/9/7,13,結(jié)構(gòu)分析編程及軟件應(yīng)用,5、實(shí)例計(jì)算分析應(yīng)用,3-1 彈簧元,k1=SpringElementStiffness(10

9、0); 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;15; 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/9/7,14,結(jié)構(gòu)分析編程及軟件應(yīng)用,1、基本方程,3-2 線性桿元,線性桿元

10、也是總體和局部坐標(biāo)一致的一維有限單元，用線性函數(shù)描述,每個(gè)線性桿元有兩個(gè)節(jié)點(diǎn)（node）,單剛矩陣為：,總剛矩陣：,結(jié)構(gòu)方程：,單元節(jié)點(diǎn)力：,2020/9/7,15,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%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 % eleme

11、nt stiffness matrix is 2 x 2.,2.1 單元?jiǎng)偠染仃嚨男纬?y = E*A/L -E*A/L ; -E*A/L E*A/L;,function y = LinearBarElementStiffness(E,A,L),3-2 線性桿元,2020/9/7,16,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%LinearBarAssemble This function assembles the element stiffness % matrix k of the linear bar with nodes i and j % into the global

12、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) = 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 線性桿

13、元,2020/9/7,17,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%LinearBarElementForces This function returns the element nodal % force vector given the element stiffness % matrix k and the element nodal % displacement vector u.,2.3 節(jié)點(diǎn)載荷計(jì)算,y = k * u;,function y = LinearBarElementForces(k,u),3-2 線性桿元,2020/9/7,18,結(jié)構(gòu)分析編程及軟件應(yīng)用,

14、2、MATLAB函數(shù)編寫(xiě),%LinearBarElementStresses This function returns the element nodal % stress vector given the element stiffness % matrix k, the element nodal displacement % vector u, and the cross-sectional area A.,2.4 節(jié)點(diǎn)應(yīng)力計(jì)算,y = k * u/A;,function y = LinearBarElementStresses(k, u, A),3-2 線性桿元,2020/9/7,1

15、9,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、實(shí)例計(jì)算分析應(yīng)用,如圖所示二線性桿元結(jié)構(gòu)，假定E=210MPa，A=0.003m2，P=10kN， 節(jié)點(diǎn)3的右位移為0.002m。 求：系統(tǒng)的整體剛度矩陣； 節(jié)點(diǎn)2的位移； 節(jié)點(diǎn)1、3的支反力； 每個(gè)桿件的應(yīng)力,解：,步驟1：離散化域,3-2 線性桿元,2020/9/7,20,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、實(shí)例計(jì)算分析應(yīng)用,步驟2：形成單元?jiǎng)偠染仃?k1=LinearBarElementStiffness(E,A,L1),k2=LinearBarElementStiffness(E,A,L2),調(diào)用 function y = LinearBarElementSti

16、ffness(E,A,L)函數(shù),3-2 線性桿元,2020/9/7,21,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、實(shí)例計(jì)算分析應(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/9/7,22,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、實(shí)例計(jì)算分析應(yīng)用,步驟4：引入邊界條件,已知邊界條件：,3-2 線性桿元,20

17、20/9/7,23,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、實(shí)例計(jì)算分析應(yīng)用,步驟5：解方程,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/9/7,24,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、實(shí)例計(jì)算分析應(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)

18、sigma1=LinearBarElementStresses(k1, u1, A),u2=U(2:3); f2= LinearBarElementForces(k2,u2) sigma2=LinearBarElementStresses(k2, u2, A),3-2 線性桿元,2020/9/7,25,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、實(shí)例計(jì)算分析應(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(

19、n,n); K = 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

20、,u2) sigma2=LinearBarElementStresses(k2, u2, A),2020/9/7,26,結(jié)構(gòu)分析編程及軟件應(yīng)用,1、基本方程,3-3 二次桿元,二次桿元也是總體和局部坐標(biāo)一致的一維有限單元，用二次方程描述,每個(gè)線性桿元有三個(gè)節(jié)點(diǎn)（node）,單剛矩陣為：,總剛矩陣：,結(jié)構(gòu)方程：,單元節(jié)點(diǎn)力：,2020/9/7,27,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%QuadraticBarElementStiffness This function returns the element % stiffness matrix for a quadratic b

21、ar % with modulus of elasticity E, % cross-sectional area A, and length L. % The size of the element stiffness % matrix is 3 x 3.,2.1 單元?jiǎng)偠染仃嚨男纬?y = E*A/(3*L)*7 1 -8 ; 1 7 -8 ; -8 -8 16;,function y = QuadraticBarElementStiffness(E,A,L),3-3 二次桿元,2020/9/7,28,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%QuadraticBarAssemb

22、le This function assembles the element stiffness % matrix k of the quadratic bar with nodes i, j % 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

23、(i,j) + k(1,2); K(i,m) = 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/9/7,29,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%Quad

24、raticBarElementForces This function returns the element nodal % force vector given the element stiffness % matrix k and the element nodal % displacement vector u.,2.3 節(jié)點(diǎn)載荷計(jì)算,y = k * u;,function y = QuadraticBarElementForces(k,u),3-3 二次桿元,2020/9/7,30,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%QuadraticBarElementStres

25、ses This function returns the element % nodal stress vector given the element % stiffness matrix k, the element nodal % displacement vector u, and the % cross-sectional area A.,2.4 節(jié)點(diǎn)應(yīng)力計(jì)算,y = k * u/A;,function y = QuadraticBarElementStresses(k, u, A),3-3 二次桿元,2020/9/7,31,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、實(shí)例計(jì)算分析應(yīng)用,如圖所示雙

26、二次桿元結(jié)構(gòu)，假定E=210MPa，A=0.003m2 求：系統(tǒng)的整體剛度矩陣； 節(jié)點(diǎn)2、3、4、5的位移； 節(jié)點(diǎn)1的支反力； 每個(gè)桿件的應(yīng)力,解：,3-3 二次桿元,2020/9/7,32,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、實(shí)例計(jì)算分析應(yīng)用,E=210E6; A=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 =QuadraticBarAss

27、emble(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= QuadraticBarElementForces(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

28、二次桿元,2020/9/7,33,結(jié)構(gòu)分析編程及軟件應(yīng)用,1、基本方程,3-4 平面桁架元,平面桁架元是既有局部坐標(biāo)又有總體坐標(biāo)二維有限元，用線性函數(shù)描述,每個(gè)平面桁架元有二個(gè)節(jié)點(diǎn)（node）,單剛矩陣為：,總剛矩陣：,結(jié)構(gòu)方程：,單元節(jié)點(diǎn)力：,2020/9/7,34,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%PlaneTrussElementLength This function returns the length of the % plane truss element whose first node has % coordinates (x1,y1) and second

29、 node has % coordinates (x2,y2).,2.1 計(jì)算單元長(zhǎng)度,y = sqrt(x2-x1)*(x2-x1) + (y2-y1)*(y2-y1);,function y = PlaneTrussElementLength(x1,y1,x2,y2),3-4 平面桁架元,2020/9/7,35,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%PlaneTrussElementStiffness This function returns the element % stiffness matrix for a plane truss % element with mo

30、dulus 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 單元?jiǎng)偠染仃嚨男纬?x = theta*pi/180; C = cos(x); 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

31、y = PlaneTrussElementStiffness(E,A,L, theta),3-4 平面桁架元,2020/9/7,36,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%PlaneTrussAssemble This function assembles the element stiffness % matrix k of the plane truss element with nodes % i and j into the global stiffness matrix K. % This function returns the global stiffness %

32、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(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(

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

34、*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/9/7,37,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%PlaneTrussElementForce This function returns the element force % given the modulus of elasticity E, the % cross-sectional a

35、rea A, the length L, % the angle theta (in degrees), and the % element nodal displacement vector u.,2.4 節(jié)點(diǎn)載荷計(jì)算,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/9/7,38,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%PlaneTrussElementStres

36、s This function returns the element stress % given the modulus of elasticity E, the % the length L, the angle theta (in % degrees), and the element nodal % displacement vector u.,2.5 節(jié)點(diǎn)應(yīng)力計(jì)算,x = theta * pi/180; C = cos(x); S = sin(x); y = E/L*-C -S C S* u;,function y = PlaneTrussElementStress(E,L,the

37、ta,u),3-4 平面桁架元,2020/9/7,39,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、實(shí)例計(jì)算分析應(yīng)用,如圖所示平面桁架結(jié)構(gòu)，假定E=210MPa，A=0.0004m2 求：系統(tǒng)的整體剛度矩陣； 節(jié)點(diǎn)2的水平位移； 節(jié)點(diǎn)3的水平豎向位移； 節(jié)點(diǎn)1、2的支反力； 每跟桿件的應(yīng)力,3-4 平面桁架元,2020/9/7,40,結(jié)構(gòu)分析編程及軟件應(yīng)用,1、基本方程,3-5 空間桁架元,空間桁架元是既有局部坐標(biāo)又有總體坐標(biāo)三維有限元，用線性函數(shù)描 述。各單元之間通過(guò)鉸接系統(tǒng)連接，只能傳遞力，而不能傳遞彎矩,每個(gè)桁架元有二個(gè)節(jié)點(diǎn)（node）,2020/9/7,41,結(jié)構(gòu)分析編程及軟件應(yīng)用,1、基本方程,3-

38、5 空間桁架元,總剛矩陣：,結(jié)構(gòu)方程：,單元節(jié)點(diǎn)力：,單剛矩陣為：,2020/9/7,42,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%SpaceTrussElementLength This function returns the length of the % space truss element whose first node has % coordinates (x1,y1,z1) and second node has % coordinates (x2,y2,z2).,2.1 計(jì)算單元長(zhǎng)度,y = sqrt(x2-x1)*(x2-x1) + (y2-y1)*(y2-y

39、1) + (z2-z1)*(z2-z1);,function y = SpaceTrussElementLength(x1,y1,z1,x2,y2,z2),3-5 空間桁架元,2020/9/7,43,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%SpaceTrussElementStiffness This function returns the element % stiffness matrix for a space truss % element with modulus of elasticity E, % cross-sectional area A, length L,

40、and % angles thetax, thetay, thetaz % (in degrees). The size of the element % stiffness matrix is 6 x 6.,2.2 單元?jiǎng)偠染仃嚨男纬?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;

41、,function y = SpaceTrussElementStiffness(E,A,L,thetax,thetay,thetaz),3-5 空間桁架元,2020/9/7,44,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%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

42、 the global 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*i-2,3*j-1) + k(

43、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);,function y =

44、SpaceTrussAssemble(K,k,i,j),3-5 空間桁架元,2020/9/7,45,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),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(3*j-1,3*j-1) + k(5

45、,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) = K(3*i,3*i-2) +

46、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,3*i) = K(3*j-2,3*i

47、) + 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/9/7,46,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%SpaceTrussElementForce This function returns the element force % given the modulus of elasticity E, the % cross-sectional area A

48、, the length L, % the angles thetax, thetay, thetaz % (in degrees), and the element nodal % displacement vector u.,2.4 節(jié)點(diǎn)載荷計(jì)算,x = thetax * pi/180; w = thetay * pi/180; v = thetaz * pi/180; Cx = cos(x); Cy = cos(w); Cz = cos(v); y = E*A/L*-Cx -Cy -Cz Cx Cy Cz*u;,function y = SpaceTrussElementForce(E,

49、A,L,thetax,thetay,thetaz,u),3-5 空間桁架元,2020/9/7,47,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%SpaceTrussElementStress This function returns the element stress % given the modulus of elasticity E, the % length L, the angles thetax, thetay, % thetaz (in degrees), and the element % nodal displacement vector u.,2.5 節(jié)點(diǎn)應(yīng)力計(jì)

50、算,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,L,thetax,thetay,thetaz,u),3-5 空間桁架元,2020/9/7,48,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、實(shí)例計(jì)算分析應(yīng)用,如圖所示空間桁架結(jié)構(gòu)，假定E=210MPa，A14=0.001m2 A24=0.002m2，A34=0.00

51、1m2，P=12kN 求：系統(tǒng)的整體剛度矩陣； 節(jié)點(diǎn)4的水平位移； 節(jié)點(diǎn)3的水平豎向位移； 節(jié)點(diǎn)1、2、3的支反力； 每跟桿件的應(yīng)力,3-5 空間桁架元,2020/9/7,49,結(jié)構(gòu)分析編程及軟件應(yīng)用,1、基本方程,3-6 梁元,梁元是總體坐標(biāo)與局部坐標(biāo)一致的二維有限元，用線性函數(shù)描 述。各單元之間通過(guò)鉸接系統(tǒng)連接，只能傳遞力，而不能傳遞彎矩,每個(gè)梁元有二個(gè)節(jié)點(diǎn)（node）,單剛矩陣為：,總剛矩陣：,結(jié)構(gòu)方程：,單元節(jié)點(diǎn)力：,2020/9/7,50,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%BeamElementStiffness This function returns the

52、element % 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單元?jiǎng)偠染仃嚨男纬?y = E*I/(L*L*L)*12 6*L -12 6*L ; 6*L 4*L*L -6*L 2*L*L ; -12 -6*L 12 -6*L ; 6*L 2*L*L -6*L 4*L*L;,function y = BeamEle

53、mentStiffness(E,I,L),3-6 梁元,2020/9/7,51,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%BeamAssemble This function assembles the element stiffness % matrix k of the beam element with nodes % i and j into the global stiffness matrix K. % This function returns the global stiffness % matrix K after the element stiffness mat

54、rix % 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(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(

55、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,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*

56、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/9/7,52,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%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é)點(diǎn)載荷計(jì)算,y

57、= k * u;,function y = BeamElementForces(k,u),3-6 梁元,2020/9/7,53,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%BeamElementShearDiagram This function plots the shear force % diagram for the beam element with nodal % force vector f and length L.,2.4 繪制剪力圖,x = 0 ; L; z = f(1) ; -f(3); hold on; title(Shear Force Diagram); p

58、lot(x,z); y1 = 0 ; 0; plot(x,y1,k),function y = BeamElementShearDiagram(f, L),3-6 梁元,2020/9/7,54,結(jié)構(gòu)分析編程及軟件應(yīng)用,2、MATLAB函數(shù)編寫(xiě),%BeamElementMomentDiagram This function plots the bending moment % diagram for the beam element with nodal % force vector f and length L.,2.4 繪制彎矩圖,x = 0 ; L; z = -f(2) ; f(4); hold on; title(Bending Moment Diagram); plot(x,z); y1 = 0 ; 0; plot(x,y1,k),function y = BeamElementMomentDiagram(f, L),3-6 梁元,2020/9/7,55,結(jié)構(gòu)分析編程及軟件應(yīng)用,3、實(shí)例計(jì)算分析應(yīng)用,3-6 梁元,