2019/11/281第三章MATLAB有限元分析與應(yīng)用

§3-1彈簧元

§3-2線性桿元

§3-3二次桿元

§3-4平面桁架元

§3-5空間桁架元

§3-6梁元

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

1、有限元方法的步驟：

離散化域

形成單剛矩陣

集成整體剛度矩陣

引入邊界條件

求解方程

后處理

2019/11/2832、基本方程

§3-1彈簧元

彈簧元是總體和局部坐標(biāo)一致的一維有限單元

每個彈簧元有兩個節(jié)點（node）

kkkkk?????????單剛矩陣為：

22?總剛矩陣：

nn???????KUF?結(jié)構(gòu)方程：

單元節(jié)點力：

??????fku?2019/11/2843、MATLAB函數(shù)編寫

§3-1彈簧元

%SpringElementStiffnessThisfunctionreturnstheelementstiffness%matrixforaspringwithstiffnessk.%Thesizeoftheelementstiffnessmatrixis2x2.3.1單元剛度矩陣的形成

y=[k-k;-kk];functiony=SpringElementStiffness(k)2019/11/2853、MATLAB函數(shù)編寫

§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函數(shù)編寫

§3-1彈簧元

%SpringElementForcesThisfunctionreturnstheelementnodalforce%vectorgiventheelementstiffnessmatrixk%andtheelementnodaldisplacementvectoru.3.3節(jié)點載荷計算

y=k*u;functiony=SpringElementForces(k,u)2019/11/2874、實例計算分析應(yīng)用

§3-1彈簧元

如圖所示二彈簧元結(jié)構(gòu)，假定k1=100kN/m，k2=200kN/m，P=15kN。

求：系統(tǒng)的整體剛度矩陣；

節(jié)點2、3的位移；

節(jié)點1的支反力；

每個彈簧的內(nèi)力

解：

步驟1：離散化域

2019/11/2884、實例計算分析應(yīng)用

§3-1彈簧元

步驟2：形成單元剛度矩陣

k1=SpringElementStiffness(100);k1=

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

k2=

200-200

-200200調(diào)用functiony=SpringElementStiffness(k)函數(shù)

2019/11/2894、實例計算分析應(yīng)用

§3-1彈簧元

步驟3：集成整體剛度矩陣

調(diào)用functiony=SpringAssemble(K,k,i,j)函數(shù)

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、實例計算分析應(yīng)用

§3-1彈簧元

步驟4：引入邊界條件

11223310010001003002000200200UFUFUF???????????????????????????????????已知邊界條件：

1230,0,15UFF???123100100001003002000020020015FUU???????????????????????????????????2019/11/28115、實例計算分析應(yīng)用

§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、實例計算分析應(yīng)用

§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、實例計算分析應(yīng)用

§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線性桿元

線性桿元也是總體和局部坐標(biāo)一致的一維有限單元，用線性函數(shù)描述

每個線性桿元有兩個節(jié)點（node）

EAEALLkEAEALL???????????????單剛矩陣為：

22?總剛矩陣：

nn???????KUF?結(jié)構(gòu)方程：

單元節(jié)點力：

??????fku?2019/11/28152、MATLAB函數(shù)編寫

%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函數(shù)編寫

%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函數(shù)編寫

%LinearBarElementForcesThisfunctionreturnstheelementnodal%forcevectorgiventheelementstiffness%matrixkandtheelementnodal%displacementvectoru.2.3節(jié)點載荷計算

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

2019/11/28182、MATLAB函數(shù)編寫

%LinearBarElementStressesThisfunctionreturnstheelementnodal%stressvectorgiventheelementstiffness%matrixk,theelementnodaldisplacement%vectoru,andthecross-sectionalareaA.2.4節(jié)點應(yīng)力計算

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

2019/11/28193、實例計算分析應(yīng)用

如圖所示二線性桿元結(jié)構(gòu)，假定E=210MPa，A=0.003m^2，P=10kN，

節(jié)點3的右位移為0.002m。

求：系統(tǒng)的整體剛度矩陣；

節(jié)點2的位移；

節(jié)點1、3的支反力；

每個桿件的應(yīng)力

解：

步驟1：離散化域

§3-2線性桿元

2019/11/28203、實例計算分析應(yīng)用

步驟2：形成單元剛度矩陣

k1=LinearBarElementStiffness(E,A,L1)k2=LinearBarElementStiffness(E,A,L2)調(diào)用functiony=LinearBarElementStiffness(E,A,L)函數(shù)

§3-2線性桿元

2019/11/28213、實例計算分析應(yīng)用

步驟3：集成整體剛度矩陣

調(diào)用functiony=LinearBarAssemble(K,k,i,j)函數(shù)

n=3;K=zeros(n,n)K=LinearBarAssemble(K,k1,1,2)K=000000000K=LinearBarAssemble(K,k2,2,3)§3-2線性桿元

2019/11/28223、實例計算分析應(yīng)用

步驟4：引入邊界條件

已知邊界條件：

1320,0.002,10UUF????§3-2線性桿元

112233420000420000042000010500006300000630000630000UFUFUF???????????????????????????????????1234200004200000042000010500006300001006300006300000.002FUF????????????????????????????????????2019/11/28233、實例計算分析應(yīng)用

步驟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、實例計算分析應(yīng)用

步驟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、實例計算分析應(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);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二次桿元

二次桿元也是總體和局部坐標(biāo)一致的一維有限單元，用二次方程描述

每個線性桿元有三個節(jié)點（node）

71817838816EAkL???????????????單剛矩陣為：

33?總剛矩陣：

nn???????KUF?結(jié)構(gòu)方程：

單元節(jié)點力：

??????fku?31?2019/11/28272、MATLAB函數(shù)編寫

%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函數(shù)編寫

%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函數(shù)編寫

%QuadraticBarElementForcesThisfunctionreturnstheelementnodal%forcevectorgiventheelementstiffness%matrixkandtheelementnodal%displacementvectoru.2.3節(jié)點載荷計算

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

2019/11/28302、MATLAB函數(shù)編寫

%QuadraticBarElementStressesThisfunctionreturnstheelement%nodalstressvectorgiventheelement%stiffnessmatrixk,theelementnodal%displacementvectoru,andthe%cross-sectionalareaA.2.4節(jié)點應(yīng)力計算

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

2019/11/28313、實例計算分析應(yīng)用

如圖所示雙二次桿元結(jié)構(gòu)，假定E=210MPa，A=0.003m^2求：系統(tǒng)的整體剛度矩陣；

節(jié)點2、3、4、5的位移；

節(jié)點1的支反力；

每個桿件的應(yīng)力

解：

§3-3二次桿元

2019/11/28323、實例計算分析應(yīng)用

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平面桁架元

平面桁架元是既有局部坐標(biāo)又有總體坐標(biāo)二維有限元，用線性函數(shù)描述

每個平面桁架元有二個節(jié)點（node）

22222222CCSCCSCSSCSSEAkLCCSCCSCSSCSS?????????????????????單剛矩陣為：

44?總剛矩陣：

22nn???????KUF?結(jié)構(gòu)方程：

單元節(jié)點力：

??????EAfCSCSuL???41?2019/11/28342、MATLAB函數(shù)編寫

%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函數(shù)編寫

%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函數(shù)編寫

%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函數(shù)編寫

%PlaneTrussElementForceThisfunctionreturnstheelementforce%giventhemodulusofelasticityE,the%cross-sectionalareaA,thelengthL,%theangletheta(indegrees),andthe%elementnodaldisplacementvectoru.2.4節(jié)點載荷計算

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函數(shù)編寫

%PlaneTrussElementStressThisfunctionreturnstheelementstress%giventhemodulusofelasticityE,the%thelengthL,theangletheta(in%degrees),andtheelementnodal%displacementvectoru.2.5節(jié)點應(yīng)力計算

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、實例計算分析應(yīng)用

如圖所示平面桁架結(jié)構(gòu)，假定E=210MPa，A=0.0004m^2求：系統(tǒng)的整體剛度矩陣；

節(jié)點2的水平位移；

節(jié)點3的水平豎向位移；

節(jié)點1、2的支反力；

每跟桿件的應(yīng)力

§3-4平面桁架元

2019/11/28401、基本方程

§3-5空間桁架元

空間桁架元是既有局部坐標(biāo)又有總體坐標(biāo)三維有限元，用線性函數(shù)描

述。各單元之間通過鉸接系統(tǒng)連接，只能傳遞力，而不能傳遞彎矩

每個桁架元有二個節(jié)點（node）

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

§3-5空間桁架元

總剛矩陣：

33nn???????KUF?結(jié)構(gòu)方程：

單元節(jié)點力：

61?????xyzxyzEAfCCCCCCuL????????222222222222xxyxzxxyxzyxyyzyxyyzzxyzzyzyzzxxyxzxxyxzyxyyzxyyyzzxyzzxzyzzCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCEAkCCCCCCCCCCLCCCCCCCCCCCCCCCCCCCC??????????????????????????????????????66?單剛矩陣為：

2019/11/28422、MATLAB函數(shù)編寫

%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函數(shù)編寫

%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函數(shù)編寫

%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函數(shù)編寫

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函數(shù)編寫

%SpaceTrussElementForceThisfunctionreturnstheelementforce%giventhemodulusofelasticityE,the%cross-sectionalareaA,thelengthL,%theanglesthetax,thetay,thetaz%(indegrees),andtheelementnodal%displacementvectoru.2.4節(jié)點載荷計算

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函數(shù)編寫

%SpaceTrussElementStressThisfunctionreturnstheelementstress%giventhemodulusofelasticityE,the%lengthL,theanglesthetax,thetay,%thetaz(indegrees),andtheelement%nodaldisplacementvectoru.2.5節(jié)點應(yīng)力計算

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、實例計算分析應(yīng)用

如圖所示空間桁架結(jié)構(gòu)，假定E=210MPa，A14=0.001m^2

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

求：系統(tǒng)的整體剛度矩陣；

節(jié)點4的水平位移；

節(jié)點3的水平豎向位移；

節(jié)點1、2、3的支反力；

每跟桿件的應(yīng)力

§3-5空間桁架元

2019/11/28491、基本方程

§3-6梁元

梁元是總體坐標(biāo)與局部坐標(biāo)一致的二維有限元，用線性函數(shù)描

述。各單元之間通過鉸接系統(tǒng)連接，只能傳遞力，而不能傳遞彎矩

每個梁元有二個節(jié)點（node）

單剛矩陣為：

2232212612664621261266264LLLLLLEIkLLLLLLL???????????????????44?總剛矩陣：

22nn???????KUF?結(jié)構(gòu)方程：

單元節(jié)點力：

41???????fku?2019/11/28502、MATLAB函數(shù)編寫

?amElementStiffnessThisfunctionreturnstheelement%stiffnessmatrixforabeam%elementwithmodulusofelasticityE,%momentofinertiaI,andlengthL.%Thesizeoftheelementstiffness%matrixis4x4.2.1單元剛度矩陣的形成

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