[数学]第三章-MATLAB有限元分析与应用课件.ppt

《[数学]第三章-MATLAB有限元分析与应用课件.ppt》由会员分享，可在线阅读，更多相关《[数学]第三章-MATLAB有限元分析与应用课件.ppt（55页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第三章第三章 MATLABMATLAB有限元分析与应用有限元分析与应用3-1 3-1 弹簧元弹簧元3-2 3-2 线性杆元线性杆元3-3 3-3 二次杆元二次杆元3-4 3-4 平面桁架元平面桁架元3-5 3-5 空间桁架元空间桁架元3-6 3-6 梁元梁元2/10/202313-1 弹簧元弹簧元 1、有限元方法的步骤：离散化域形成单刚矩阵集成整体刚度矩阵引入边界条件求解方程后处理2/10/202322、基本方程3-1 弹簧元弹簧元 弹簧元是总体和局部坐标一致的一维有限单元每个弹簧元有两个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：2/10/202333、MATLAB函数编写3

2、-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)2/10/202343、MATLAB函数编写3-1 弹簧元弹簧元%SpringAssemble This function assembles the

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

4、function y=SpringAssemble(K,k,i,j)2/10/202353、MATLAB函数编写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 节点载荷计算y=k*u;function y=SpringElementForces(k,u)2/10/202364、实例计算分析应用3-1 弹簧元弹

5、簧元 如图所示二弹簧元结构，假定k1=100kN/m，k2=200kN/m，P=15kN。求：系统的整体刚度矩阵；节点2、3的位移；节点1的支反力；每个弹簧的内力解：步骤1：离散化域2/10/202374、实例计算分析应用3-1 弹簧元弹簧元 步骤2：形成单元刚度矩阵k1=SpringElementStiffness(100);k1=100 -100 -100 100k2=SpringElementStiffness(200)；k2=200 -200 -200 200调用 function y=SpringElementStiffness(k)函数2/10/202384、实例计算分析应用3-1

6、 弹簧元弹簧元 步骤3：集成整体刚度矩阵调用 function y=SpringAssemble(K,k,i,j)函数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 0 -200 2002/10/202394、实例计算分析应用3-1 弹簧元弹簧元 步骤4：引入边界条件已知边界条件：2/10/2023105、实例计算分析应用3-1 弹簧元弹簧元 步骤5：解方程

7、U=zeros(2,1);F=0;15;K=K(2:3,2:3);U=KFU=inv(K)*FK(1,:)=;K(:,1)=;U=0.1500 0.22502/10/2023115、实例计算分析应用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.00002/10/2023125、

8、实例计算分析应用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;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)2/10/2023131、基本方

9、程3-2 线性杆元线性杆元 线性杆元也是总体和局部坐标一致的一维有限单元，用线性函数描述每个线性杆元有两个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：2/10/2023142、MATLAB函数编写%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

10、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 线性杆元线性杆元 2/10/2023152、MATLAB函数编写%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 retu

11、rns 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 线性杆元线性杆元 2/10/2023162、MATLAB函数编写%LinearBarElementForces This function retu

12、rns the element nodal%force vector given the element stiffness%matrix k and the element nodal%displacement vector u.2.3 节点载荷计算y=k*u;function y=LinearBarElementForces(k,u)3-2 线性杆元线性杆元 2/10/2023172、MATLAB函数编写%LinearBarElementStresses This function returns the element nodal%stress vector given the elem

13、ent stiffness%matrix k,the element nodal displacement%vector u,and the cross-sectional area A.2.4 节点应力计算y=k*u/A;function y=LinearBarElementStresses(k,u,A)3-2 线性杆元线性杆元 2/10/2023183、实例计算分析应用如图所示二线性杆元结构，假定E=210MPa，A=0.003m2，P=10kN，节点3的右位移为0.002m。求：系统的整体刚度矩阵；节点2的位移；节点1、3的支反力；每个杆件的应力解：步骤1：离散化域3-2 线性杆元线性杆

14、元 2/10/2023193、实例计算分析应用步骤2：形成单元刚度矩阵k1=LinearBarElementStiffness(E,A,L1)k2=LinearBarElementStiffness(E,A,L2)调用 function y=LinearBarElementStiffness(E,A,L)函数3-2 线性杆元线性杆元 2/10/2023203、实例计算分析应用步骤3：集成整体刚度矩阵调用 function y=LinearBarAssemble(K,k,i,j)函数n=3;K=zeros(n,n)K=LinearBarAssemble(K,k1,1,2)K=0 0 0 0 0

15、0 0 0 0K=LinearBarAssemble(K,k2,2,3)3-2 线性杆元线性杆元 2/10/2023213、实例计算分析应用步骤4：引入边界条件已知边界条件：3-2 线性杆元线性杆元 2/10/2023223、实例计算分析应用步骤5：解方程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 线性杆元线性杆元 2/10/2023233、实例计算分析应用步骤6：后处理U=0;U;0.002U=0 0.0012 0.0002F=K*UF=-500.0000 -10.000

16、0 510.0000u1=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)3-2 线性杆元线性杆元 2/10/2023243、实例计算分析应用E=210E6;A=0.003;L1=1.5;L2=1;k1=LinearBarElementStiffness(E,A,L1);k2=LinearBarElementStiffnes

17、s(E,A,L2);n=3;K=zeros(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*Uu1=U(1:2);f1=LinearBarElementForces(k1,u1)sigma1=LinearBarElementStresses(k1,u1,A)u2=U(2:3);f2=LinearBarElementForces(k2

18、,u2)sigma2=LinearBarElementStresses(k2,u2,A)2/10/2023251、基本方程3-3 二次杆元二次杆元 二次杆元也是总体和局部坐标一致的一维有限单元，用二次方程描述每个线性杆元有三个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：2/10/2023262、MATLAB函数编写%QuadraticBarElementStiffness This function returns the element%stiffness matrix for a quadratic bar%with modulus of elasticity E,%cro

19、ss-sectional area A,and length L.%The size of the element stiffness%matrix 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 二次杆元二次杆元 2/10/2023272、MATLAB函数编写%QuadraticBarAssemble This function assembles the element stiffness%matrix k of the quadr

20、atic 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(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=Qua

21、draticBarAssemble(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;2/10/2023282、MATLAB函数编写%QuadraticBarElementForces This function returns the element nodal%force vector given the element stiffness%matrix k and the element nodal%displacem

22、ent vector u.2.3 节点载荷计算y=k*u;function y=QuadraticBarElementForces(k,u)3-3 二次杆元二次杆元 2/10/2023292、MATLAB函数编写%QuadraticBarElementStresses 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

23、 节点应力计算y=k*u/A;function y=QuadraticBarElementStresses(k,u,A)3-3 二次杆元二次杆元 2/10/2023303、实例计算分析应用如图所示双二次杆元结构，假定E=210MPa，A=0.003m2求：系统的整体刚度矩阵；节点2、3、4、5的位移；节点1的支反力；每个杆件的应力解：3-3 二次杆元二次杆元 2/10/2023313、实例计算分析应用E=210E6;A=0.003;L=2;k1=QuadraticBarElementStiffness(E,A,L);k2=QuadraticBarElementStiffness(E,A,L);

24、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(k1,u1,A);u2=U(3);U(5);U(4);f2=QuadraticBarElementForces(k2,u2

25、);sigma2=QuadraticBarElementStresses(k2,u2,A);3-3 二次杆元二次杆元 2/10/2023321、基本方程3-4 平面桁架元平面桁架元 平面桁架元是既有局部坐标又有总体坐标二维有限元，用线性函数描述每个平面桁架元有二个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：2/10/2023332、MATLAB函数编写%PlaneTrussElementLength This function returns the length of the%plane truss element whose first node has%coordinat

26、es(x1,y1)and second 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 平面桁架元平面桁架元 2/10/2023342、MATLAB函数编写%PlaneTrussElementStiffness This function returns the element%stiffness matrix for a plane truss%element with modulus o

27、f 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);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

28、 平面桁架元平面桁架元 2/10/2023352、MATLAB函数编写%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%matrix K after the element stiffness matrix k is assembled.2.3 整体

29、刚度矩阵的形成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(2,3);K(2*i,2*j)=K(2*i,2*j)+k(2,4);function y=PlaneTrussAssemble

30、(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*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 平面桁架元平面桁架元 2/10/202

31、3362、MATLAB函数编写%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.2.4 节点载荷计算x=theta*pi/180;C=cos(x);S=sin(x);y=E*A/L*-C-S C S*u;function y=P

32、laneTrussElementForce(E,A,L,theta,u)3-4 平面桁架元平面桁架元 2/10/2023372、MATLAB函数编写%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 节点应力计算x=theta*pi/180;C=cos(x);S=s

33、in(x);y=E/L*-C-S C S*u;function y=PlaneTrussElementStress(E,L,theta,u)3-4 平面桁架元平面桁架元 2/10/2023383、实例计算分析应用如图所示平面桁架结构，假定E=210MPa，A=0.0004m2求：系统的整体刚度矩阵；节点2的水平位移；节点3的水平竖向位移；节点1、2的支反力；每跟杆件的应力3-4 平面桁架元平面桁架元 2/10/2023391、基本方程3-5 空间桁架元空间桁架元 空间桁架元是既有局部坐标又有总体坐标三维有限元，用线性函数描 述。各单元之间通过铰接系统连接，只能传递力，而不能传递弯矩 每个桁架元

34、有二个节点（node）2/10/2023401、基本方程3-5 空间桁架元空间桁架元 总刚矩阵：结构方程：单元节点力：单刚矩阵为：2/10/2023412、MATLAB函数编写%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 计算单元长度y=sqrt(x2-x1)*(x2-x1)+(y2-y1)*(y2

35、-y1)+(z2-z1)*(z2-z1);function y=SpaceTrussElementLength(x1,y1,z1,x2,y2,z2)3-5 空间桁架元空间桁架元 2/10/2023422、MATLAB函数编写%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,the

36、tay,thetaz%(in degrees).The size of the element%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=SpaceTrussElementStiffness(E,A,L,thetax,thetay,thet

37、az)3-5 空间桁架元空间桁架元 2/10/2023432、MATLAB函数编写%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.

38、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(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)

39、;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 空间桁架元空间桁架元 2/10/2023442、MATLAB函数编写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

40、*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(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)=

41、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-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

42、-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);2/10/2023452、MATLAB函数编写%SpaceTrussElementForce This function returns the element force%given the modulus of elasticity E,the%cross-sectional ar

43、ea A,the length L,%the angles thetax,thetay,thetaz%(in degrees),and the element nodal%displacement vector u.2.4 节点载荷计算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,A,L,thetax,thetay,thetaz,u)3-5 空间桁架元空间桁

44、架元 2/10/2023462、MATLAB函数编写%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 节点应力计算x=thetax*pi/180;w=thetay*pi/180;v=thetaz*pi/180;Cx=cos(x);Cy=c

45、os(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 空间桁架元空间桁架元 2/10/2023473、实例计算分析应用如图所示空间桁架结构，假定E=210MPa，A14=0.001m2 A24=0.002m2，A34=0.001m2，P=12kN求：系统的整体刚度矩阵；节点4的水平位移；节点3的水平竖向位移；节点1、2、3的支反力；每跟杆件的应力3-5 空间桁架元空间桁架元 2/10/2023481、基本方程3-6 梁元梁元 梁元是

46、总体坐标与局部坐标一致的二维有限元，用线性函数描 述。各单元之间通过铰接系统连接，只能传递力，而不能传递弯矩 每个梁元有二个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：2/10/2023492、MATLAB函数编写%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 stiff

47、ness%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*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 梁元梁元 2/10/2023502、MATLAB函数编写%BeamAssemble This function assembles the element stiffness%matrix k of the beam element with nodes%i and j into t

48、he global stiffness matrix K.%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(1,4);K(2*i,2*i-1)=K(

49、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,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-

50、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;2/10/2023512、MATLAB函数编写%BeamElementForces This function returns the element nodal force%vector given the element stiffness matrix k%and the element nodal displacement vector u.2.4 节点