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

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

2、： fku2022-3-1643、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)偠染仃嚨男纬蓎 = k -k ; -k k;function y = SpringElementStiffness(k)2022-3-1653、MATLAB函數(shù)編寫(xiě)3-1 彈簧元彈簧元 %Spr

3、ingAssemble 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) + k

4、(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)2022-3-1663、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 vector u.3.3 節(jié)點(diǎn)載荷計(jì)算y

5、 = k * u;function y = SpringElementForces(k,u)2022-3-1674、實(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：離散化域2022-3-1684、實(shí)例計(jì)算分析應(yīng)用3-1 彈簧元彈簧元 步驟2：形成單元?jiǎng)偠染仃噆1=SpringElementStiffness(100);k1 = 100 -100 -100 100k2=SpringElementStiffness(200)；k2 =

6、200 -200 -200 200調(diào)用 function y = SpringElementStiffness(k)函數(shù)2022-3-1694、實(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 0K = SpringAssemble(K,k2,2,3)K = 100 -100 0 -100 100 0 0 0 0K = 100 -100 0 -100 300 -200

7、 0 -200 2002022-3-16104、實(shí)例計(jì)算分析應(yīng)用3-1 彈簧元彈簧元 步驟4：引入邊界條件11223310010001003002000200200UFUFUF已知邊界條件：1230,0,15UFF123100100001003002000020020015FUU2022-3-16115、實(shí)例計(jì)算分析應(yīng)用3-1 彈簧元彈簧元 步驟5：解方程23300200020020015UUU=zeros(2,1);F=0;15;K = K(2:3,2:3);U=KFU=inv(K)*FK(1,:)=;K(:,1)=;U = 0.1500 0.22502022-3-16125、實(shí)例計(jì)算分析應(yīng)

8、用2-1 彈簧元彈簧元 步驟6：后處理U=0;UU = 0 0.1500 0.2250F=K*UF = -15.0000 0.0000 15.0000u1=U(1:2);f1=SpringElementForces(k1,u1);f1 = -15.0000 15.0000u2=U(2:3);f2=SpringElementForces(k2,u2);f2 = -15.0000 15.00002022-3-16135、實(shí)例計(jì)算分析應(yīng)用3-1 彈簧元彈簧元 k1=SpringElementStiffness(100);k2=SpringElementStiffness(200);n=3;K=zer

9、os(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=KFU=0;U;F=K*U;u1=U(1:2);f1=SpringElementForces(k1,u1)u2=U(2:3);f2=SpringElementForces(k2,u2)2022-3-16141、基本方程3-2 線性桿元線性桿元 線性桿元也是總體和局部坐標(biāo)一致的一維有限單元，用線性函數(shù)描述每個(gè)線性桿元有兩個(gè)節(jié)點(diǎn)（node）EAEALLkEAEALL 單剛矩陣為：2 2總剛矩陣

10、：nn KUF結(jié)構(gòu)方程：?jiǎn)卧?jié)點(diǎn)力： fku2022-3-16152、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 % element stiffness matrix is 2 x 2.2.1 單元?jiǎng)偠染仃嚨男纬蓎 = E*A/L -E*A/L ; -E

11、*A/L E*A/L;function y = LinearBarElementStiffness(E,A,L)3-2 線性桿元線性桿元 2022-3-16162、MATLAB函數(shù)編寫(xiě)%LinearBarAssemble This function assembles the element stiffness% matrix k of the linear bar with nodes i and j % into the global stiffness matrix K.% This function returns the global stiffness % matrix K aft

12、er 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 線性桿元線性桿元 2022-3-16172、MATLAB函數(shù)編寫(xiě)%LinearBarElementForces This function returns the element

13、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 線性桿元線性桿元 2022-3-16182、MATLAB函數(shù)編寫(xiě)%LinearBarElementStresses This function returns the element nodal % stress vector given the element

14、 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 線性桿元線性桿元 2022-3-16193、實(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：

15、離散化域3-2 線性桿元線性桿元 2022-3-16203、實(shí)例計(jì)算分析應(yīng)用步驟2：形成單元?jiǎng)偠染仃噆1=LinearBarElementStiffness(E,A,L1)k2=LinearBarElementStiffness(E,A,L2)調(diào)用 function y = LinearBarElementStiffness(E,A,L)函數(shù)3-2 線性桿元線性桿元 2022-3-16213、實(shí)例計(jì)算分析應(yīng)用步驟3：集成整體剛度矩陣調(diào)用 function y = LinearBarAssemble(K,k,i,j)函數(shù)n=3; K = zeros(n,n)K = LinearBarAssemb

16、le (K,k1,1,2)K = 0 0 0 0 0 0 0 0 0K = LinearBarAssemble (K,k2,2,3)3-2 線性桿元線性桿元 2022-3-16223、實(shí)例計(jì)算分析應(yīng)用步驟4：引入邊界條件已知邊界條件：1320,0.002,10UUF 3-2 線性桿元線性桿元 112233420000420000042000010500006300000630000630000UFUFUF1234200004200000042000010500006300001006300006300000.002FUF 2022-3-16233、實(shí)例計(jì)算分析應(yīng)用步驟5：解方程 2105000

17、6300000.00210U U=zeros(1,1);U3=0.002F=-10;K = K(2,2) 105000K0 = K(2,3); -630000U=K(F-K0*U3)U =0.00123-2 線性桿元線性桿元 2022-3-16243、實(shí)例計(jì)算分析應(yīng)用步驟6：后處理U=0;U;0.002U = 0 0.0012 0.0002F=K*UF = -500.0000 -10.0000 510.0000u1=U(1:2);f1= LinearBarElementForces(k1,u1)sigma1=LinearBarElementStresses(k1, u1, A)u2=U(2:3

18、);f2= LinearBarElementForces(k2,u2)sigma2=LinearBarElementStresses(k2, u2, A)3-2 線性桿元線性桿元 2022-3-16253、實(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(n,n);K = LinearBarAssemble (K,k1,1,2);K = LinearBarAssemble (K,k2,2,3)

19、;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*Uu1=U(1:2);f1= LinearBarElementForces(k1,u1)sigma1=LinearBarElementStresses(k1, u1, A)u2=U(2:3);f2= LinearBarElementForces(k2,u2)sigma2=LinearBarElementStresses(k2, u2, A)2022-3-16261、基本方程3-3 二次桿元二次桿元 二次桿元也是

20、總體和局部坐標(biāo)一致的一維有限單元，用二次方程描述每個(gè)線性桿元有三個(gè)節(jié)點(diǎn)（node）71817838816EAkL單剛矩陣為：3 3總剛矩陣：nn KUF結(jié)構(gòu)方程：?jiǎn)卧?jié)點(diǎn)力： fku3 12022-3-16272、MATLAB函數(shù)編寫(xiě)%QuadraticBarElementStiffness This function returns the element % stiffness matrix for a quadratic bar % with modulus of elasticity E, % cross-sectional area A, and length L. % The si

21、ze of the element stiffness % matrix is 3 x 3.2.1 單元?jiǎng)偠染仃嚨男纬蓎 = E*A/(3*L)*7 1 -8 ; 1 7 -8 ; -8 -8 16;function y = QuadraticBarElementStiffness(E,A,L)3-3 二次桿元二次桿元 2022-3-16282、MATLAB函數(shù)編寫(xiě)%QuadraticBarAssemble This function assembles the element stiffness% matrix k of the quadratic bar with nodes i, j %

22、 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) = 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 =

23、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;2022-3-16292、MATLAB函數(shù)編寫(xiě)%QuadraticBarElementForces This function returns the element nodal % force vector given the element stiffness % matrix k and th

24、e element nodal % displacement vector u.2.3 節(jié)點(diǎn)載荷計(jì)算y = k * u;function y = QuadraticBarElementForces(k,u)3-3 二次桿元二次桿元 2022-3-16302、MATLAB函數(shù)編寫(xiě)%QuadraticBarElementStresses This function returns the element % nodal stress vector given the element % stiffness matrix k, the element nodal % displacement vec

25、tor 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 二次桿元二次桿元 2022-3-16313、實(shí)例計(jì)算分析應(yīng)用如圖所示雙二次桿元結(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 二次桿元二次桿元 2022-3-16323、實(shí)例計(jì)算分析應(yīng)用E=210E6;A=0.003;L=2;k1= QuadraticBarElementS

26、tiffness(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= QuadraticBarElementForces(k1,u1);sigma1=QuadraticBarElementStresses

27、(k1, u1, A);u2=U(3);U(5);U(4);f2=QuadraticBarElementForces(k2,u2);sigma2=QuadraticBarElementStresses(k2, u2, A);3-3 二次桿元二次桿元 2022-3-16331、基本方程3-4 平面桁架元平面桁架元 平面桁架元是既有局部坐標(biāo)又有總體坐標(biāo)二維有限元，用線性函數(shù)描述每個(gè)平面桁架元有二個(gè)節(jié)點(diǎn)（node）22222222CCSCCSCSSCSSEAkLCCSCCSCSSCSS單剛矩陣為：4 4總剛矩陣：22nn KUF結(jié)構(gòu)方程：?jiǎn)卧?jié)點(diǎn)力： EAfCSCSuL4 12022-3-16342

28、、MATLAB函數(shù)編寫(xiě)%PlaneTrussElementLength This function returns the length of the% plane truss element whose first node has % coordinates (x1,y1) and second 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 平面桁架元平面桁架元 2

29、022-3-16352、MATLAB函數(shù)編寫(xiě)%PlaneTrussElementStiffness This function returns the element % stiffness matrix 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 單元?jiǎng)偠染仃嚨男纬蓌 = th

30、eta*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 y = PlaneTrussElementStiffness(E,A,L, theta)3-4 平面桁架元平面桁架元 2022-3-16362、MATLAB函數(shù)編寫(xiě)%PlaneTrussAssemble This function assembles the element stiffness% matrix k of the plane

31、 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(1,2);K(2*i-1,2*j-1) = K(2*i-1,2*j-1) +

32、 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) + k(3,1);K(2*j-1,2*i) = K(2*j-1,2*i) + k(3,2

33、);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 平面桁架元平面桁架元 2022-3-16372、MATLAB函數(shù)編寫(xiě)%PlaneTrussElementForce This functio

34、n 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.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,

35、L,theta,u)3-4 平面桁架元平面桁架元 2022-3-16382、MATLAB函數(shù)編寫(xiě)%PlaneTrussElementStress 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

36、 = E/L*-C -S C S* u;function y = PlaneTrussElementStress(E,L,theta,u)3-4 平面桁架元平面桁架元 2022-3-16393、實(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 平面桁架元平面桁架元 2022-3-16401、基本方程3-5 空間桁架元空間桁架元 空間桁架元是既有局部坐標(biāo)又有總體坐標(biāo)三維有限元，用線性函數(shù)描 述。各單元之間通過(guò)鉸接系統(tǒng)連接，只能傳遞力，而不能傳遞彎矩 每個(gè)

37、桁架元有二個(gè)節(jié)點(diǎn)（node）cos,cos,cosxxyyzzCCC2022-3-16411、基本方程3-5 空間桁架元空間桁架元 總剛矩陣：33nn KUF結(jié)構(gòu)方程：?jiǎn)卧?jié)點(diǎn)力：6 1 xyzxyzEAfCCCCCCuL222222222222xxyxzxxyxzyxyyzyxyyzzxyzzyzyzzxxyxzxxyxzyxyyzxyyyzzxyzzxzyzzCC CC CCC CC CC CCC CC CCC CC CC CCC CC CCEAkCC CC CCC CC CLC CCC CC CCC CC CC CCC CC CC6 6單剛矩陣為：2022-3-16422、MATLAB函

38、數(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-y1) + (z2-z1)*(z2-z1);function y = SpaceTrussElementLength(x1,y1,z1,x2,

39、y2,z2)3-5 空間桁架元空間桁架元 2022-3-16432、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, and% angles thetax, thetay, thetaz % (in degrees). The size of the element % stif

40、fness matrix is 6 x 6.2.2 單元?jiǎng)偠染仃嚨男纬蓌 = 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 = SpaceTrussElementStiffness(E,A,L,thetax,thetay,thetaz)3-5 空間桁架元空間桁架元 2022-3-16

41、442、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 the global stiffness % matrix K after the element stiffness matrix % k is assembled.2.3 整體剛度矩陣的形成K(3*i-2

42、,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(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

43、-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 =SpaceTrussAssemble(K,k,i,j)3-5 空間桁架元空間桁架元 2022-3-16452、MATLAB函數(shù)編寫(xiě)2.3 整體剛度矩陣的形成3-5 空間桁架元空間桁架元 K(3*j-1,3*i-2) = K(3*j-1,3

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

45、(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) + 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

46、-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) + 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);2022-3-16462、MAT

47、LAB函數(shù)編寫(xiě)%SpaceTrussElementForce This function returns the element force% given the modulus of elasticity E, the % cross-sectional area A, 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

48、= 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,A,L,thetax,thetay,thetaz,u)3-5 空間桁架元空間桁架元 2022-3-16472、MATLAB函數(shù)編寫(xiě)%SpaceTrussElementStress This function returns the element stress% given the modulus of elasticity E, the % le

49、ngth L, the angles thetax, thetay, % thetaz (in degrees), and the element % nodal displacement vector u.2.5 節(jié)點(diǎn)應(yīng)力計(jì)算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,thet

50、az,u)3-5 空間桁架元空間桁架元 2022-3-16483、實(shí)例計(jì)算分析應(yīng)用如圖所示空間桁架結(jié)構(gòu)，假定E=210MPa，A14=0.001m2 A24=0.002m2，A34=0.001m2，P=12kN求：系統(tǒng)的整體剛度矩陣； 節(jié)點(diǎn)4的水平位移； 節(jié)點(diǎn)3的水平豎向位移； 節(jié)點(diǎn)1、2、3的支反力； 每跟桿件的應(yīng)力3-5 空間桁架元空間桁架元 2022-3-16491、基本方程3-6 梁元梁元 梁元是總體坐標(biāo)與局部坐標(biāo)一致的二維有限元，用線性函數(shù)描 述。各單元之間通過(guò)鉸接系統(tǒng)連接，只能傳遞力，而不能傳遞彎矩 每個(gè)梁元有二個(gè)節(jié)點(diǎn)（node）單剛矩陣為：22322126126646212612

51、66264LLLLLLEIkLLLLLLL4 4總剛矩陣：22nn KUF結(jié)構(gòu)方程：?jiǎn)卧?jié)點(diǎn)力：4 1 fku2022-3-16502、MATLAB函數(shù)編寫(xiě)%BeamElementStiffness This function returns the 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)偠染仃嚨男纬蓎 = 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 = BeamElementStiffness(E,I,L)3-6 梁元梁元 2022-3-16512、MATLAB函數(shù)編寫(xiě)%BeamAssemble This function assembles the element stiffness% matrix k of the b