MATLAB有限元分析与应用

2020/3/191第三章MATLAB有限元分析与应用§3-1弹簧元§3-2线性杆元§3-3二次杆元§3-4平面桁架元§3-5空间桁架元§3-6梁元2020/3/192§3-1弹簧元1、有限元方法的步骤：离散化域形成单刚矩阵集成整体刚度矩阵引入边界条件求解方程后处理2020/3/1932、基本方程§3-1弹簧元弹簧元是总体和局部坐标一致的一维有限单元每个弹簧元有两个节点（node）kkkkk单刚矩阵为：22总刚矩阵：nnKUF结构方程：单元节点力：fku2020/3/1943、MATLAB函数编写§3-1弹簧元%SpringElementStiffnessThisfunctionreturnstheelementstiffness%matrixforaspringwithstiffnessk.%Thesizeoftheelementstiffnessmatrixis2x2.3.1单元刚度矩阵的形成y=[k-k;-kk];functiony=SpringElementStiffness(k)2020/3/1953、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)2020/3/1963、MATLAB函数编写§3-1弹簧元%SpringElementForcesThisfunctionreturnstheelementnodalforce%vectorgiventheelementstiffnessmatrixk%andtheelementnodaldisplacementvectoru.3.3节点载荷计算y=k*u;functiony=SpringElementForces(k,u)2020/3/1974、实例计算分析应用§3-1弹簧元如图所示二弹簧元结构，假定k1=100kN/m，k2=200kN/m，P=15kN。求：系统的整体刚度矩阵；节点2、3的位移；节点1的支反力；每个弹簧的内力解：步骤1：离散化域2020/3/1984、实例计算分析应用§3-1弹簧元步骤2：形成单元刚度矩阵k1=SpringElementStiffness(100);k1=100-100-100100k2=SpringElementStiffness(200)；k2=200-200-200200调用functiony=SpringElementStiffness(k)函数2020/3/1994、实例计算分析应用§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-2002002020/3/19104、实例计算分析应用§3-1弹簧元步骤4：引入边界条件11223310010001003002000200200UFUFUF已知边界条件：1230,0,15UFF123100100001003002000020020015FUU2020/3/19115、实例计算分析应用§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.22502020/3/19125、实例计算分析应用§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.00002020/3/19135、实例计算分析应用§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)2020/3/19141、基本方程§3-2线性杆元线性杆元也是总体和局部坐标一致的一维有限单元，用线性函数描述每个线性杆元有两个节点（node）EAEALLkEAEALL单刚矩阵为：22总刚矩阵：nnKUF结构方程：单元节点力：fku2020/3/19152、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线性杆元2020/3/19162、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线性杆元2020/3/19172、MATLAB函数编写%LinearBarElementForcesThisfunctionreturnstheelementnodal%forcevectorgiventheelementstiffness%matrixkandtheelementnodal%displacementvectoru.2.3节点载荷计算y=k*u;functiony=LinearBarElementForces(k,u)§3-2线性杆元2020/3/19182、MATLAB函数编写%LinearBarElementStressesThisfunctionreturnstheelementnodal%stressvectorgiventheelementstiffness%matrixk,theelementnodaldisplacement%vectoru,andthecross-sectionalareaA.2.4节点应力计算y=k*u/A;functiony=LinearBarElementStresses(k,u,A)§3-2线性杆元2020/3/19193、实例计算分析应用如图所示二线性杆元结构，假定E=210MPa，A=0.003m^2，P=10kN，节点3的右位移为0.002m。求：系统的整体刚度矩阵；节点2的位移；节点1、3的支反力；每个杆件的应力解：步骤1：离散化域§3-2线性杆元2020/3/19203、实例计算分析应用步骤2：形成单元刚度矩阵k1=LinearBarElementStiffness(E,A,L1)k2=LinearBarElementStiffness(E,A,L2)调用functiony=LinearBarElementStiffness(E,A,L)函数§3-2线性杆元2020/3/19213、实例计算分析应用步骤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线性杆元2020/3/19223、实例计算分析应用步骤4：引入边界条件已知边界条件：1320,0.002,10UUF§3-2线性杆元112233420000420000042000010500006300000630000630000UFUFUF1234200004200000042000010500006300001006300006300000.002FUF2020/3/19233、实例计算分析应用步骤5：解方程21050006300000.00210U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线性杆元2020/3/19243、实例计算分析应用步骤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线性杆元2020/3/19253、实例计算分析应用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