第三章-MATLAB有限元分析与应用课件

第三章MATLAB有限元分析与应用§3-1弹簧元结构分析编程及软件应用§3-2线性杆元§3-3二次杆元§3-4平面桁架元§3-5空间桁架元§3-6梁元7/31/20231第三章MATLAB有限元分析与应用§3-1弹簧元结构分§3-1弹簧元结构分析编程及软件应用1、有限元方法的步骤：离散化域形成单刚矩阵集成整体刚度矩阵引入边界条件求解方程后处理7/31/20232§3-1弹簧元结构分析编程及软件应用1、有限元结构分析编程及软件应用2、基本方程§3-1弹簧元弹簧元是总体和局部坐标一致的一维有限单元每个弹簧元有两个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：7/31/20233结构分析编程及软件应用2、基本方程§3-1弹簧元结构分析编程及软件应用3、MATLAB函数编写§3-1弹簧元%SpringElementStiffnessThisfunctionreturnstheelementstiffness%matrixforaspringwithstiffnessk.%Thesizeoftheelementstiffnessmatrixis2x2.3.1单元刚度矩阵的形成y=[k-k;-kk];functiony=SpringElementStiffness(k)7/31/20234结构分析编程及软件应用3、MATLAB函数编写§3-1结构分析编程及软件应用3、MATLAB函数编写§3-1弹簧元%SpringAssembleThisfunctionassemblestheelementstiffness%matrixkofthespringwithnodesiandjintothe%globalstiffnessmatrixK.%ThisfunctionreturnstheglobalstiffnessmatrixK%aftertheelementstiffnessmatrixkisassembled.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;functiony=SpringAssemble(K,k,i,j)7/31/20235结构分析编程及软件应用3、MATLAB函数编写§3-1结构分析编程及软件应用3、MATLAB函数编写§3-1弹簧元%SpringElementForcesThisfunctionreturnstheelementnodalforce%vectorgiventheelementstiffnessmatrixk%andtheelementnodaldisplacementvectoru.3.3节点载荷计算y=k*u;functiony=SpringElementForces(k,u)7/31/20236结构分析编程及软件应用3、MATLAB函数编写§3-1结构分析编程及软件应用4、实例计算分析应用§3-1弹簧元如图所示二弹簧元结构，假定k1=100kN/m，k2=200kN/m，P=15kN。求：系统的整体刚度矩阵；节点2、3的位移；节点1的支反力；每个弹簧的内力解：步骤1：离散化域7/31/20237结构分析编程及软件应用4、实例计算分析应用§3-1结构分析编程及软件应用4、实例计算分析应用§3-1弹簧元步骤2：形成单元刚度矩阵k1=SpringElementStiffness(100);k1=

100-100-100100k2=SpringElementStiffness(200)；k2=

200-200 -200200调用functiony=SpringElementStiffness(k)函数7/31/20238结构分析编程及软件应用4、实例计算分析应用§3-1结构分析编程及软件应用4、实例计算分析应用§3-1弹簧元步骤3：集成整体刚度矩阵调用functiony=SpringAssemble(K,k,i,j)函数n=3;K=zeros(n,n);K=SpringAssemble(K,k1,1,2)K=000000000K=SpringAssemble(K,k2,2,3)K=100-1000-1001000000K=100-1000-100300-2000-2002007/31/20239结构分析编程及软件应用4、实例计算分析应用§3-1结构分析编程及软件应用4、实例计算分析应用§3-1弹簧元步骤4：引入边界条件已知边界条件：7/31/202310结构分析编程及软件应用4、实例计算分析应用§3-1结构分析编程及软件应用5、实例计算分析应用§3-1弹簧元步骤5：解方程U=zeros(2,1);F=[0;15];K=K(2:3,2:3);U=K\FU=inv(K)*FK(1,:)=[];K(:,1)=[];U=0.15000.22507/31/202311结构分析编程及软件应用5、实例计算分析应用§3-1结构分析编程及软件应用5、实例计算分析应用§2-1弹簧元步骤6：后处理U=[0;U]U=00.15000.2250F=K*UF=-15.00000.000015.0000u1=U(1:2);f1=SpringElementForces(k1,u1);f1=-15.000015.0000u2=U(2:3);f2=SpringElementForces(k2,u2);f2=-15.000015.00007/31/202312结构分析编程及软件应用5、实例计算分析应用§2-1结构分析编程及软件应用5、实例计算分析应用§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=K\FU=[0;U];F=K*U;u1=U(1:2);f1=SpringElementForces(k1,u1)u2=U(2:3);f2=SpringElementForces(k2,u2)7/31/202313结构分析编程及软件应用5、实例计算分析应用§3-1结构分析编程及软件应用1、基本方程§3-2线性杆元线性杆元也是总体和局部坐标一致的一维有限单元，用线性函数描述每个线性杆元有两个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：7/31/202314结构分析编程及软件应用1、基本方程§3-2线性杆元结构分析编程及软件应用2、MATLAB函数编写%LinearBarElementStiffnessThisfunctionreturnstheelement%stiffnessmatrixforalinearbarwith%modulusofelasticityE,cross-sectional%areaA,andlengthL.Thesizeofthe%elementstiffnessmatrixis2x2.2.1单元刚度矩阵的形成y=[E*A/L-E*A/L;-E*A/LE*A/L];functiony=LinearBarElementStiffness(E,A,L)§3-2线性杆元7/31/202315结构分析编程及软件应用2、MATLAB函数编写%Linear结构分析编程及软件应用2、MATLAB函数编写%LinearBarAssembleThisfunctionassemblestheelementstiffness%matrixkofthelinearbarwithnodesiandj%intotheglobalstiffnessmatrixK.%Thisfunctionreturnstheglobalstiffness%matrixKaftertheelementstiffnessmatrix%kisassembled.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;functiony=LinearBarAssemble(K,k,i,j)§3-2线性杆元7/31/202316结构分析编程及软件应用2、MATLAB函数编写%Linear结构分析编程及软件应用2、MATLAB函数编写%LinearBarElementForcesThisfunctionreturnstheelementnodal%forcevectorgiventheelementstiffness%matrixkandtheelementnodal%displacementvectoru.2.3节点载荷计算y=k*u;functiony=LinearBarElementForces(k,u)§3-2线性杆元7/31/202317结构分析编程及软件应用2、MATLAB函数编写%Linear结构分析编程及软件应用2、MATLAB函数编写%LinearBarElementStressesThisfunctionreturnstheelementnodal%stressvectorgiventheelementstiffness%matrixk,theelementnodaldisplacement%vectoru,andthecross-sectionalareaA.2.4节点应力计算y=k*u/A;functiony=LinearBarElementStresses(k,u,A)§3-2线性杆元7/31/202318结构分析编程及软件应用2、MATLAB函数编写%Linear结构分析编程及软件应用3、实例计算分析应用如图所示二线性杆元结构，假定E=210MPa，A=0.003m^2，P=10kN，节点3的右位移为0.002m。求：系统的整体刚度矩阵；节点2的位移；节点1、3的支反力；每个杆件的应力解：步骤1：离散化域§3-2线性杆元7/31/202319结构分析编程及软件应用3、实例计算分析应用如图所示二线性杆元结构分析编程及软件应用3、实例计算分析应用步骤2：形成单元刚度矩阵k1=LinearBarElementStiffness(E,A,L1)k2=LinearBarElementStiffness(E,A,L2)调用functiony=LinearBarElementStiffness(E,A,L)函数§3-2线性杆元7/31/202320结构分析编程及软件应用3、实例计算分析应用步骤2：形成单元刚结构分析编程及软件应用3、实例计算分析应用步骤3：集成整体刚度矩阵调用functiony=LinearBarAssemble(K,k,i,j)函数n=3;K=zeros(n,n)K=LinearBarAssemble(K,k1,1,2)K=000000000K=LinearBarAssemble(K,k2,2,3)§3-2线性杆元7/31/202321结构分析编程及软件应用3、实例计算分析应用步骤3：集成整体刚结构分析编程及软件应用3、实例计算分析应用步骤4：引入边界条件已知边界条件：§3-2线性杆元7/31/202322结构分析编程及软件应用3、实例计算分析应用步骤4：引入边界条结构分析编程及软件应用3、实例计算分析应用步骤5：解方程U=zeros(1,1);U3=0.002F=[-10];K=K(2,2)

105000K0=K(2,3);-630000U=K\(F-K0*U3)U=0.0012§3-2线性杆元7/31/202323结构分析编程及软件应用3、实例计算分析应用步骤5：解方程U=结构分析编程及软件应用3、实例计算分析应用步骤6：后处理U=[0;U;0.002]U=00.00120.0002F=K*UF=-500.0000-10.0000510.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线性杆元7/31/202324结构分析编程及软件应用3、实例计算分析应用步骤6：后处理U=结构分析编程及软件应用3、实例计算分析应用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);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/202325结构分析编程及软件应用3、实例计算分析应用E=210E6;§结构分析编程及软件应用1、基本方程§3-3二次杆元二次杆元也是总体和局部坐标一致的一维有限单元，用二次方程描述每个线性杆元有三个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：7/31/202326结构分析编程及软件应用1、基本方程§3-3二次杆元结构分析编程及软件应用2、MATLAB函数编写%QuadraticBarElementStiffnessThisfunctionreturnstheelement%stiffnessmatrixforaquadraticbar%withmodulusofelasticityE,%cross-sectionalareaA,andlengthL.%Thesizeoftheelementstiffness%matrixis3x3.2.1单元刚度矩阵的形成y=E*A/(3*L)*[71-8;17-8;-8-816];functiony=QuadraticBarElementStiffness(E,A,L)§3-3二次杆元7/31/202327结构分析编程及软件应用2、MATLAB函数编写%Quadra结构分析编程及软件应用2、MATLAB函数编写%QuadraticBarAssembleThisfunctionassemblestheelementstiffness%matrixkofthequadraticbarwithnodesi,j%andmintotheglobalstiffnessmatrixK.%Thisfunctionreturnstheglobalstiffness%matrixKaftertheelementstiffnessmatrix%kisassembled.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);functiony=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;7/31/202328结构分析编程及软件应用2、MATLAB函数编写%Quadra结构分析编程及软件应用2、MATLAB函数编写%QuadraticBarElementForcesThisfunctionreturnstheelementnodal%forcevectorgiventheelementstiffness%matrixkandtheelementnodal%displacementvectoru.2.3节点载荷计算y=k*u;functiony=QuadraticBarElementForces(k,u)§3-3二次杆元7/31/202329结构分析编程及软件应用2、MATLAB函数编写%Quadra结构分析编程及软件应用2、MATLAB函数编写%QuadraticBarElementStressesThisfunctionreturnstheelement%nodalstressvectorgiventheelement%stiffnessmatrixk,theelementnodal%displacementvectoru,andthe%cross-sectionalareaA.2.4节点应力计算y=k*u/A;functiony=QuadraticBarElementStresses(k,u,A)§3-3二次杆元7/31/202330结构分析编程及软件应用2、MATLAB函数编写%Quadra结构分析编程及软件应用3、实例计算分析应用如图所示双二次杆元结构，假定E=210MPa，A=0.003m^2求：系统的整体刚度矩阵；节点2、3、4、5的位移；节点1的支反力；每个杆件的应力解：§3-3二次杆元7/31/202331结构分析编程及软件应用3、实例计算分析应用如图所示双二次杆元结构分析编程及软件应用3、实例计算分析应用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=QuadraticBarAssemble(K,k2,3,5,4);U=zeros(4,1);F=[5;-10;-7;10];KK=K;K=K(2:n,2:n);U=K\F;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二次杆元7/31/202332结构分析编程及软件应用3、实例计算分析应用E=210E6;K结构分析编程及软件应用1、基本方程§3-4平面桁架元平面桁架元是既有局部坐标又有总体坐标二维有限元，用线性函数描述每个平面桁架元有二个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：7/31/202333结构分析编程及软件应用1、基本方程§3-4平面桁架结构分析编程及软件应用2、MATLAB函数编写%PlaneTrussElementLengthThisfunctionreturnsthelengthofthe%planetrusselementwhosefirstnodehas%coordinates(x1,y1)andsecondnodehas%coordinates(x2,y2).2.1计算单元长度y=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));functiony=PlaneTrussElementLength(x1,y1,x2,y2)§3-4平面桁架元7/31/202334结构分析编程及软件应用2、MATLAB函数编写%PlaneT结构分析编程及软件应用2、MATLAB函数编写%PlaneTrussElementStiffnessThisfunctionreturnstheelement%stiffnessmatrixforaplanetruss%elementwithmodulusofelasticityE,%cross-sectionalareaA,lengthL,and%angletheta(indegrees).%Thesizeoftheelementstiffness%matrixis4x4.2.2单元刚度矩阵的形成x=theta*pi/180;C=cos(x);S=sin(x);y=E*A/L*[C*CC*S-C*C-C*S;C*SS*S-C*S-S*S;-C*C-C*SC*CC*S;-C*S-S*SC*SS*S];functiony=PlaneTrussElementStiffness(E,A,L,theta)§3-4平面桁架元7/31/202335结构分析编程及软件应用2、MATLAB函数编写%PlaneT结构分析编程及软件应用2、MATLAB函数编写%PlaneTrussAssembleThisfunctionassemblestheelementstiffness%matrixkoftheplanetrusselementwithnodes%iandjintotheglobalstiffnessmatrixK.%Thisfunctionreturnstheglobalstiffness%matrixKaftertheelementstiffnessmatrixkisassembled.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(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);functiony=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*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/202336结构分析编程及软件应用2、MATLAB函数编写%PlaneT结构分析编程及软件应用2、MATLAB函数编写%PlaneTrussElementForceThisfunctionreturnstheelementforce%giventhemodulusofelasticityE,the%cross-sectionalareaA,thelengthL,%theangletheta(indegrees),andthe%elementnodaldisplacementvectoru.2.4节点载荷计算x=theta*pi/180;C=cos(x);S=sin(x);y=E*A/L*[-C-SCS]*u;functiony=PlaneTrussElementForce(E,A,L,theta,u)§3-4平面桁架元7/31/202337结构分析编程及软件应用2、MATLAB函数编写%PlaneT结构分析编程及软件应用2、MATLAB函数编写%PlaneTrussElementStressThisfunctionreturnstheelementstress%giventhemodulusofelasticityE,the%thelengthL,theangletheta(in%degrees),andtheelementnodal%displacementvectoru.2.5节点应力计算x=theta*pi/180;C=cos(x);S=sin(x);y=E/L*[-C-SCS]*u;functiony=PlaneTrussElementStress(E,L,theta,u)§3-4平面桁架元7/31/202338结构分析编程及软件应用2、MATLAB函数编写%PlaneT结构分析编程及软件应用3、实例计算分析应用如图所示平面桁架结构，假定E=210MPa，A=0.0004m^2求：系统的整体刚度矩阵；节点2的水平位移；节点3的水平竖向位移；节点1、2的支反力；每跟杆件的应力§3-4平面桁架元7/31/202339结构分析编程及软件应用3、实例计算分析应用如图所示平面桁架结结构分析编程及软件应用1、基本方程§3-5空间桁架元空间桁架元是既有局部坐标又有总体坐标三维有限元，用线性函数描述。各单元之间通过铰接系统连接，只能传递力，而不能传递弯矩每个桁架元有二个节点（node）7/31/202340结构分析编程及软件应用1、基本方程§3-5空间桁架结构分析编程及软件应用1、基本方程§3-5空间桁架元总刚矩阵：结构方程：单元节点力：单刚矩阵为：7/31/202341结构分析编程及软件应用1、基本方程§3-5空间桁架结构分析编程及软件应用2、MATLAB函数编写%SpaceTrussElementLengthThisfunctionreturnsthelengthofthe%spacetrusselementwhosefirstnodehas%coordinates(x1,y1,z1)andsecondnodehas%coordinates(x2,y2,z2).2.1计算单元长度y=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));functiony=SpaceTrussElementLength(x1,y1,z1,x2,y2,z2)§3-5空间桁架元7/31/202342结构分析编程及软件应用2、MATLAB函数编写%SpaceT结构分析编程及软件应用2、MATLAB函数编写%SpaceTrussElementStiffnessThisfunctionreturnstheelement%stiffnessmatrixforaspacetruss%elementwithmodulusofelasticityE,%cross-sectionalareaA,lengthL,and%anglesthetax,thetay,thetaz%(indegrees).Thesizeoftheelement%stiffnessmatrixis6x6.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*CxCx*CyCx*Cz;Cy*CxCy*CyCy*Cz;Cz*CxCz*CyCz*Cz];y=E*A/L*[w-w;-ww];functiony=SpaceTrussElementStiffness(E,A,L,thetax,thetay,thetaz)§3-5空间桁架元7/31/202343结构分析编程及软件应用2、MATLAB函数编写%SpaceT结构分析编程及软件应用2、MATLAB函数编写%SpaceTrussAssembleThisfunctionassemblestheelementstiffness%matrixkofthespacetrusselementwithnodes%iandjintotheglobalstiffnessmatrixK.%Thisfunctionreturnstheglobalstiffness%matrixKaftertheelementstiffnessmatrix%kisassembled.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);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);functiony=SpaceTrussAssemble(K,k,i,j)§3-5空间桁架元7/31/202344结构分析编程及软件应用2、MATLAB函数编写%SpaceT结构分析编程及软件应用2、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*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)=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-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/202345结构分析编程及软件应用2、MATLAB函数编写2.3整体刚结构分析编程及软件应用2、MATLAB函数编写%SpaceTrussElementForceThisfunctionreturnstheelementforce%giventhemodulusofelasticityE,the%cross-sectionalareaA,thelengthL,%theanglesthetax,thetay,thetaz%(indegrees),andtheelementnodal%displacementvectoru.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-CzCxCyCz]*u;functiony=SpaceTrussElementForce(E,A,L,thetax,thetay,thetaz,u)§3-5空间桁架元7/31/202346结构分析编程及软件应用2、MATLAB函数编写%SpaceT结构分析编程及软件应用2、MATLAB函数编写%SpaceTrussElementStressThisfunctionreturnstheelementstress%giventhemodulusofelasticityE,the%lengthL,theanglesthetax,thetay,%thetaz(indegrees),andtheelement%nodaldisplacementvectoru.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-CzCxCyCz]*u;functiony=SpaceTrussElementStress(E,L,thetax,thetay,thetaz,u)§3-5空间桁架元7/31/202347结构分析编程及软件应用2、MATLAB函数编写%SpaceT结构分析编程及软件应用3、实例计算分析应用如图所示空间桁架结构，假定E=210MPa，A14=0.001m^2

A24=0.002m^2，A34=0.001m^2，P=12kN求：系统的整体刚度矩阵；节点4的水平位移；节点3的水平竖向位移；节点1、2、3的支反力；每跟杆件的应力§3-5空间桁架元7/31/202348结构分析编程及软件应用3、实例计算分析应用如图所示空间桁架结结构分析编程及软件应用1、基本方程§3-6梁元梁元是总体坐标与局部坐标一致的二维有限元，用线性函数描述。各单元之间通过铰接系统连接，只能传递力，而不能传递弯矩每个梁元有二个节点（node）单刚矩阵为：总刚矩阵：结构方程：单元节点力：7/31/202349结构分析编程及软件应用1、基本方程§3-6梁元结构分析编程及软件应用2、MATLAB函数编写%BeamElementStiffnessThisfunctionreturnstheelement%stiffnessmatrixforabeam%elementwithmodulusofelasticityE,%momentofinertiaI,andlengthL.%Thesizeoftheelementstiffness%matrixis4x4.2.1单元刚度矩阵的形成y=E*I/(L*L*L)*[126*L-126*L;6*L4*L*L-6*L2*L*L;-12-6*L12-6*L;6*L2*L*L-6*L4*L*L];functiony=BeamElementStiffness(E,I,L)§3-6梁元7/31/202350结构分析编程及软件应用2、MATLAB函数编写%BeamEl结构分析编程及软件应用2、MATLAB函数编写%BeamAssembleThisfunctionassemblest