《第三章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 梁元梁元7/31/202413-1 弹簧元弹簧元 1、有限元方法的步骤：离散化域形成单刚矩阵集成整体刚度矩阵引入边界条件求解方程后处理7/31/202422、基本方程3-1 弹簧元弹簧元 弹簧元是总体和局部坐标一致的一维有限单元每个弹簧元有两个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：7/31/202433、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)7/31/202443、MATLAB函数编写3-1 弹簧元弹簧元 %SpringAssemble This function as
3、sembles 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(1,2);K(j,i) = K(j,i) + k(2,
4、1);K(j,j) = K(j,j) + k(2,2);y = K;function y = SpringAssemble(K,k,i,j)7/31/202453、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 = Spring
5、ElementForces(k,u)7/31/202464、实例计算分析应用3-1 弹簧元弹簧元 如图所示二弹簧元结构，假定k1=100kN/m，k2=200kN/m，P=15kN。求：系统的整体刚度矩阵； 节点2、3的位移； 节点1的支反力； 每个弹簧的内力解：步骤1：离散化域7/31/202474、实例计算分析应用3-1 弹簧元弹簧元 步骤2：形成单元刚度矩阵k1=SpringElementStiffness(100);k1 = 100 -100 -100 100k2=SpringElementStiffness(200)；k2 = 200 -200 -200 200调用 function
6、 y = SpringElementStiffness(k)函数7/31/202484、实例计算分析应用3-1 弹簧元弹簧元 步骤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 2007/31/202494、实例计算分
7、析应用3-1 弹簧元弹簧元 步骤4：引入边界条件已知边界条件：7/31/2024105、实例计算分析应用3-1 弹簧元弹簧元 步骤5：解方程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.22507/31/2024115、实例计算分析应用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
8、 15.0000u2=U(2:3);f2=SpringElementForces(k2,u2);f2 = -15.0000 15.00007/31/2024125、实例计算分析应用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);
9、f1=SpringElementForces(k1,u1)u2=U(2:3);f2=SpringElementForces(k2,u2)7/31/2024131、基本方程3-2 线性杆元线性杆元 线性杆元也是总体和局部坐标一致的一维有限单元，用线性函数描述每个线性杆元有两个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：7/31/2024142、MATLAB函数编写%LinearBarElementStiffness This function returns the element % stiffness matrix for a linear bar with % modulu
10、s of elasticity E, cross-sectional % area A, and length L. The size of the % 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 线性杆元线性杆元 7/31/2024152、MATLAB函数编写%LinearBarAssemble This function assembles the element stiffnes
11、s% 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 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) +
12、k(2,2);y = K;function y =LinearBarAssemble(K,k,i,j)3-2 线性杆元线性杆元 7/31/2024162、MATLAB函数编写%LinearBarElementForces This function returns 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 = LinearBarElementF
13、orces(k,u)3-2 线性杆元线性杆元 7/31/2024172、MATLAB函数编写%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 节点应力计算y = k * u/A;function y = LinearBarElementStresse
14、s(k, u, A)3-2 线性杆元线性杆元 7/31/2024183、实例计算分析应用如图所示二线性杆元结构，假定E=210MPa，A=0.003m2，P=10kN， 节点3的右位移为0.002m。求：系统的整体刚度矩阵； 节点2的位移； 节点1、3的支反力； 每个杆件的应力解：步骤1：离散化域3-2 线性杆元线性杆元 7/31/2024193、实例计算分析应用步骤2：形成单元刚度矩阵k1=LinearBarElementStiffness(E,A,L1)k2=LinearBarElementStiffness(E,A,L2)调用 function y = LinearBarElementS
15、tiffness(E,A,L)函数3-2 线性杆元线性杆元 7/31/2024203、实例计算分析应用步骤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 0 0 0 0K = LinearBarAssemble (K,k2,2,3)3-2 线性杆元线性杆元 7/31/2024213、实例计算分析应用步骤4：引入边界条件已知边界条件：3-2 线性杆元线性杆元 7/31/2024223、实例计算分析应用步骤5：
16、解方程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 线性杆元线性杆元 7/31/2024233、实例计算分析应用步骤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);f2= Linear
17、BarElementForces(k2,u2)sigma2=LinearBarElementStresses(k2, u2, A)3-2 线性杆元线性杆元 7/31/2024243、实例计算分析应用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);U=zeros(1,1
18、);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)7/31/2024251、基本方程3-3 二次杆元二次杆元 二次杆元也是总体和局部坐标一致的一维
19、有限单元，用二次方程描述每个线性杆元有三个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：7/31/2024262、MATLAB函数编写%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 size of the element stiffness % matrix is 3 x
20、 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 二次杆元二次杆元 7/31/2024272、MATLAB函数编写%QuadraticBarAssemble This function assembles the element stiffness% matrix k of the quadratic bar with nodes i, j % and m into the global stiffness matrix K.%
21、 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 =QuadraticBarAssemble(K,k,i,j,m)3-3 二次杆元二次杆元
22、 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;7/31/2024282、MATLAB函数编写%QuadraticBarElementForces This function returns the element nodal % force vector given the element stiffness % matrix k and the element nodal % displacement vector u.2.3
23、 节点载荷计算y = k * u;function y = QuadraticBarElementForces(k,u)3-3 二次杆元二次杆元 7/31/2024292、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
24、 节点应力计算y = k * u/A;function y = QuadraticBarElementStresses(k, u, A)3-3 二次杆元二次杆元 7/31/2024303、实例计算分析应用如图所示双二次杆元结构，假定E=210MPa，A=0.003m2求：系统的整体刚度矩阵； 节点2、3、4、5的位移； 节点1的支反力； 每个杆件的应力解：3-3 二次杆元二次杆元 7/31/2024313、实例计算分析应用E=210E6;A=0.003;L=2;k1= QuadraticBarElementStiffness(E,A,L);k2= QuadraticBarElementStif
25、fness(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(k1, u1, A);u2=U(3);U(5);U(4);f2=QuadraticB
26、arElementForces(k2,u2);sigma2=QuadraticBarElementStresses(k2, u2, A);3-3 二次杆元二次杆元 7/31/2024321、基本方程3-4 平面桁架元平面桁架元 平面桁架元是既有局部坐标又有总体坐标二维有限元，用线性函数描述每个平面桁架元有二个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：7/31/2024332、MATLAB函数编写%PlaneTrussElementLength This function returns the length of the% plane truss element whose
27、first node has % coordinates (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 平面桁架元平面桁架元 7/31/2024342、MATLAB函数编写%PlaneTrussElementStiffness This function returns the element % stiffness matrix
28、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);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 ;
29、-C*S -S*S C*S S*S;function y = PlaneTrussElementStiffness(E,A,L, theta)3-4 平面桁架元平面桁架元 7/31/2024352、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 globa
30、l 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) + 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*
31、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);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*
32、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 平面桁架元平面桁架元 7/31/2024362、MATLAB函数编写%PlaneTrussElementForce This function returns the element force% given the modulus of elasticity E, the % cross-sectional area A, the lengt
33、h 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 = PlaneTrussElementForce(E,A,L,theta,u)3-4 平面桁架元平面桁架元 7/31/2024372、MATLAB函数编写%PlaneTrussElementStress This function returns the elem
34、ent 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 = sin(x);y = E/L*-C -S C S* u;function y = PlaneTrussElementStress(E,L,theta,u)3-4 平面桁架元平面桁架元 7/31/2024383、实例计算分析
35、应用如图所示平面桁架结构，假定E=210MPa，A=0.0004m2求：系统的整体刚度矩阵； 节点2的水平位移； 节点3的水平竖向位移； 节点1、2的支反力； 每跟杆件的应力3-4 平面桁架元平面桁架元 7/31/2024391、基本方程3-5 空间桁架元空间桁架元 空间桁架元是既有局部坐标又有总体坐标三维有限元，用线性函数描 述。各单元之间通过铰接系统连接，只能传递力，而不能传递弯矩 每个桁架元有二个节点（node）7/31/2024401、基本方程3-5 空间桁架元空间桁架元 总刚矩阵：结构方程：单元节点力：单刚矩阵为：7/31/2024412、MATLAB函数编写%SpaceTrussE
36、lementLength 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-y1) + (z2-z1)*(z2-z1);function y = SpaceTrussElementLength(x1,y1,z1,x2,y2,z2)3-5 空间桁架元
37、空间桁架元 7/31/2024422、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, thetay, thetaz % (in degrees). The size of the element % stiffness matrix is
38、 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,thetaz)3-5 空间桁架元空间桁架元 7/31/2024432、MATLAB函数编写%
39、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,3*i-2) = K(3*i
40、-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);K(
41、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 空间桁架元空间桁架元 7/31/2024442、MATLAB函数编写2.3 整体刚度矩阵的形成3-5 空间桁架元空间桁架元 K(3*j-1,3*i-2) = K(3*j-1,3*i-2) + k(5,1);
42、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(3*j,3*i) + k(6
43、,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-1) + k(3,5);K(
44、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);7/31/2024452、MATLAB函数编写%SpaceTr
45、ussElementForce 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 节点载荷计算x = thetax * pi/180;w = thetay * pi/180;v = thetaz * pi/1
46、80;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 空间桁架元空间桁架元 7/31/2024462、MATLAB函数编写%SpaceTrussElementStress This function returns the element stress% given the modulus of elasticity E, the % length L, the ang
47、les 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 = 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 空间桁架元空
48、间桁架元 7/31/2024473、实例计算分析应用如图所示空间桁架结构，假定E=210MPa，A14=0.001m2 A24=0.002m2，A34=0.001m2，P=12kN求：系统的整体刚度矩阵； 节点4的水平位移； 节点3的水平竖向位移； 节点1、2、3的支反力； 每跟杆件的应力3-5 空间桁架元空间桁架元 7/31/2024481、基本方程3-6 梁元梁元 梁元是总体坐标与局部坐标一致的二维有限元，用线性函数描 述。各单元之间通过铰接系统连接，只能传递力，而不能传递弯矩 每个梁元有二个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：7/31/2024492、MATLAB
49、函数编写%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单元刚度矩阵的形成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
50、-6*L ; 6*L 2*L*L -6*L 4*L*L;function y = BeamElementStiffness(E,I,L)3-6 梁元梁元 7/31/2024502、MATLAB函数编写%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
51、 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(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*
52、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-1) = K(2*j,2*i-1) + k(4,1);K(2
53、*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;7/31/2024512、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 节点载荷计算y = k
54、 * u;function y = BeamElementForces(k,u)3-6 梁元梁元 7/31/2024522、MATLAB函数编写%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);plot(x,z);y1 = 0 ; 0;plot
55、(x,y1,k)function y = BeamElementShearDiagram(f, L)3-6 梁元梁元 7/31/2024532、MATLAB函数编写%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 梁元梁元 7/31/2024543、实例计算分析应用3-6 梁元梁元 7/31/202455
