MATLAB有限元分析与应用

2019/11/281第三章MATLAB有限元分析与应用

§3-1弹簧元

§3-2线性杆元

§3-3二次杆元

§3-4平面桁架元

§3-5空间桁架元

§3-6梁元

2019/11/282§3-1弹簧元

1、有限元方法的步骤：

离散化域

形成单刚矩阵

集成整体刚度矩阵

引入边界条件

求解方程

后处理

2019/11/2832、基本方程

§3-1弹簧元

弹簧元是总体和局部坐标一致的一维有限单元

每个弹簧元有两个节点（node）

kkkkk?????????单刚矩阵为：

22?总刚矩阵：

nn???????KUF?结构方程：

单元节点力：

??????fku?2019/11/2843、MATLAB函数编写

§3-1弹簧元

%SpringElementStiffnessThisfunctionreturnstheelementstiffness%matrixforaspringwithstiffnessk.%Thesizeoftheelementstiffnessmatrixis2x2.3.1单元刚度矩阵的形成

y=[k-k;-kk];functiony=SpringElementStiffness(k)2019/11/2853、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)2019/11/2863、MATLAB函数编写

§3-1弹簧元

%SpringElementForcesThisfunctionreturnstheelementnodalforce%vectorgiventheelementstiffnessmatrixk%andtheelementnodaldisplacementvectoru.3.3节点载荷计算

y=k*u;functiony=SpringElementForces(k,u)2019/11/2874、实例计算分析应用

§3-1弹簧元

如图所示二弹簧元结构，假定k1=100kN/m，k2=200kN/m，P=15kN。

求：系统的整体刚度矩阵；

节点2、3的位移；

节点1的支反力；

每个弹簧的内力

解：

步骤1：离散化域

2019/11/2884、实例计算分析应用

§3-1弹簧元

步骤2：形成单元刚度矩阵

k1=SpringElementStiffness(100);k1=

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

k2=

200-200

-200200调用functiony=SpringElementStiffness(k)函数

2019/11/2894、实例计算分析应用

§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-2002002019/11/28104、实例计算分析应用

§3-1弹簧元

步骤4：引入边界条件

11223310010001003002000200200UFUFUF???????????????????????????????????已知边界条件：

1230,0,15UFF???123100100001003002000020020015FUU???????????????????????????????????2019/11/28115、实例计算分析应用

§3-1弹簧元

步骤5：解方程

23300200020020015UU?????????????????????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.22502019/11/28125、实例计算分析应用

§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.00002019/11/28135、实例计算分析应用

§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)2019/11/28141、基本方程

§3-2线性杆元

线性杆元也是总体和局部坐标一致的一维有限单元，用线性函数描述

每个线性杆元有两个节点（node）

EAEALLkEAEALL???????????????单刚矩阵为：

22?总刚矩阵：

nn???????KUF?结构方程：

单元节点力：

??????fku?2019/11/28152、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线性杆元

2019/11/28162、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线性杆元

2019/11/28172、MATLAB函数编写

%LinearBarElementForcesThisfunctionreturnstheelementnodal%forcevectorgiventheelementstiffness%matrixkandtheelementnodal%displacementvectoru.2.3节点载荷计算

y=k*u;functiony=LinearBarElementForces(k,u)§3-2线性杆元

2019/11/28182、MATLAB函数编写

%LinearBarElementStressesThisfunctionreturnstheelementnodal%stressvectorgiventheelementstiffness%matrixk,theelementnodaldisplacement%vectoru,andthecross-sectionalareaA.2.4节点应力计算

y=k*u/A;functiony=LinearBarElementStresses(k,u,A)§3-2线性杆元

2019/11/28193、实例计算分析应用

如图所示二线性杆元结构，假定E=210MPa，A=0.003m^2，P=10kN，

节点3的右位移为0.002m。

求：系统的整体刚度矩阵；

节点2的位移；

节点1、3的支反力；

每个杆件的应力

解：

步骤1：离散化域

§3-2线性杆元

2019/11/28203、实例计算分析应用

步骤2：形成单元刚度矩阵

k1=LinearBarElementStiffness(E,A,L1)k2=LinearBarElementStiffness(E,A,L2)调用functiony=LinearBarElementStiffness(E,A,L)函数

§3-2线性杆元

2019/11/28213、实例计算分析应用

步骤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线性杆元

2019/11/28223、实例计算分析应用

步骤4：引入边界条件

已知边界条件：

1320,0.002,10UUF????§3-2线性杆元

112233420000420000042000010500006300000630000630000UFUFUF???????????????????????????????????1234200004200000042000010500006300001006300006300000.002FUF????????????????????????????????????2019/11/28233、实例计算分析应用

步骤5：解方程

??????????21050006300000.00210U????U=zeros(1,1);U3=0.002F=[-10];K=K(2,2)

105000K0=K(2,3);-630000

U=K\(F-K0*U3)U=0.0012§3-2线性杆元

2019/11/28243、实例计算分析应用

步骤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线性杆元

2019/11/28253、实例计算分析应用

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)2019/11/28261、基本方程

§3-3二次杆元

二次杆元也是总体和局部坐标一致的一维有限单元，用二次方程描述

每个线性杆元有三个节点（node）

71817838816EAkL???????????????单刚矩阵为：

33?总刚矩阵：

nn???????KUF?结构方程：

单元节点力：

??????fku?31?2019/11/28272、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二次杆元

2019/11/28282、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;2019/11/28292、MATLAB函数编写

%QuadraticBarElementForcesThisfunctionreturnstheelementnodal%forcevectorgiventheelementstiffness%matrixkandtheelementnodal%displacementvectoru.2.3节点载荷计算

y=k*u;functiony=QuadraticBarElementForces(k,u)§3-3二次杆元

2019/11/28302、MATLAB函数编写

%QuadraticBarElementStressesThisfunctionreturnstheelement%nodalstressvectorgiventheelement%stiffnessmatrixk,theelementnodal%displacementvectoru,andthe%cross-sectionalareaA.2.4节点应力计算

y=k*u/A;functiony=QuadraticBarElementStresses(k,u,A)§3-3二次杆元

2019/11/28313、实例计算分析应用

如图所示双二次杆元结构，假定E=210MPa，A=0.003m^2求：系统的整体刚度矩阵；

节点2、3、4、5的位移；

节点1的支反力；

每个杆件的应力

解：

§3-3二次杆元

2019/11/28323、实例计算分析应用

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二次杆元

2019/11/28331、基本方程

§3-4平面桁架元

平面桁架元是既有局部坐标又有总体坐标二维有限元，用线性函数描述

每个平面桁架元有二个节点（node）

22222222CCSCCSCSSCSSEAkLCCSCCSCSSCSS?????????????????????单刚矩阵为：

44?总刚矩阵：

22nn???????KUF?结构方程：

单元节点力：

??????EAfCSCSuL???41?2019/11/28342、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平面桁架元

2019/11/28352、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平面桁架元

2019/11/28362、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平面桁架元

2019/11/28372、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平面桁架元

2019/11/28382、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平面桁架元

2019/11/28393、实例计算分析应用

如图所示平面桁架结构，假定E=210MPa，A=0.0004m^2求：系统的整体刚度矩阵；

节点2的水平位移；

节点3的水平竖向位移；

节点1、2的支反力；

每跟杆件的应力

§3-4平面桁架元

2019/11/28401、基本方程

§3-5空间桁架元

空间桁架元是既有局部坐标又有总体坐标三维有限元，用线性函数描

述。各单元之间通过铰接系统连接，只能传递力，而不能传递弯矩

每个桁架元有二个节点（node）

cos,cos,cosxxyyzzCCC??????2019/11/28411、基本方程

§3-5空间桁架元

总刚矩阵：

33nn???????KUF?结构方程：

单元节点力：

61?????xyzxyzEAfCCCCCCuL????????222222222222xxyxzxxyxzyxyyzyxyyzzxyzzyzyzzxxyxzxxyxzyxyyzxyyyzzxyzzxzyzzCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCEAkCCCCCCCCCCLCCCCCCCCCCCCCCCCCCCC??????????????????????????????????????66?单刚矩阵为：

2019/11/28422、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空间桁架元

2019/11/28432、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空间桁架元

2019/11/28442、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空间桁架元

2019/11/28452、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);2019/11/28462、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空间桁架元

2019/11/28472、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空间桁架元

2019/11/28483、实例计算分析应用

如图所示空间桁架结构，假定E=210MPa，A14=0.001m^2

A24=0.002m^2，A34=0.001m^2，P=12kN

求：系统的整体刚度矩阵；

节点4的水平位移；

节点3的水平竖向位移；

节点1、2、3的支反力；

每跟杆件的应力

§3-5空间桁架元

2019/11/28491、基本方程

§3-6梁元

梁元是总体坐标与局部坐标一致的二维有限元，用线性函数描

述。各单元之间通过铰接系统连接，只能传递力，而不能传递弯矩

每个梁元有二个节点（node）

单刚矩阵为：

2232212612664621261266264LLLLLLEIkLLLLLLL???????????????????44?总刚矩阵：

22nn???????KUF?结构方程：

单元节点力：

41???????fku?2019/11/28502、MATLAB函数编写

?amElementStiffnessThisfunctionreturnstheelement%stiffnessmatrixforabeam%elementwithmodulusofelasticityE,%momentofinertiaI,andlengthL.%Thesizeoftheelementstiffness%matrixis4x4.2.1单元刚度矩阵的形成

y=E*I/(L*L*L)*[126*L-12