西交大有限元原理与应用-大作业

..有限元原理及工程应用——大作业学院：机械工程学院班级：硕4002班小组成员：李追3114001089陈草3114001080..作业题目：利用有限元方法对简支梁问题进行求解,梁的横截面为矩形,其约束情况如图1所示。已知梁的几何尺寸和物理参数如下：〔1几何尺寸：长度,截面尺寸；〔2物理参数：弹性模量GPa,泊松比,密度。图1．梁及其横截面示意图要求：至少划分五个节点〔四个单元；给出单元节点信息；给出单元刚度矩阵和质量矩阵；给出总刚度矩阵和总质量矩阵；求出梁各界固有频率及振型〔五阶；将所得结果与理论值进行对比,验证方法的可行性。解：由有限元知识,根据Rayleigh-Ritz法,解有限元分为四步：建立离散化、单元分析、形成总体方程、解方程,具体步骤如下：〔1建立离散化这里我们将矩形截面简支梁等分四等分,即分为六节点的五个杆单元,如图2所示：每个单元尺寸,这里只考虑杆在竖直平面的弯曲,每个节点只有y方向位移和绕z轴的旋转自由度。〔2单元分析构造一组Lagrange插值基函数,在本节点值为1,其他节点值为0。从Rayleigh-Ritz法可以看到,插值函数要p次可微,最高阶导数出现在应变能表达式中；同样,我们可以这一原则适用于基函数的选择以及形状函数,否则我们将无法正确计算应变能当我们使用有限元逼近方法。梁的弯曲问题,应变能计算公式：〔1-1其中,E为弹性模量,Iz为截面惯性矩。从公式可知,位移函数必须连续,并且二阶导数平方可积。如图3,是一维杆单元模型,每个节点两个自由度,该单元含有四个自由度,即〔。本题中我们采用三次多项式插值函数：〔1-2因此,我们必须给出四个形函数〔位移模式。图3一维杆单元模型构造Hermite插值函数。选择局部坐标系<,>,其中l是单元长度,转角是挠度值的一阶导数,定义边界条件：〔1-3因此,我们给出变形的Hermite的多项式插值函数：〔1-4其中,和分别满足如下条件,对应的图形如图4所示：〔1-5图4一维Hermite插值多项式基于Langrage和Hermite插值多项式,写出单元形函数〔1-6节点位移值也可以得出〔1-7同时,表达式〔1-7用矩阵表示为〔1-8其中,,用能量表达式替代表达式中的和。动能表达式：〔1-9将〔1-8带入〔1-9,得到〔1-10从而获得质量矩阵：〔1-11带入,〔1-12应变能表达式：〔1-13刚度矩阵表达式：〔1-14带入,可以得到〔1-15形成总体方程将每个杆单元的能量方程组装。完整梁上的的总动能和能量的和所做的总功梁上的外力作用,所有的自由度的位移矢量可以给出：〔1-16将位移矢量转换为全局坐标系下的位移矢量,变换矩阵为：〔1-17〔1-18分别以矩阵形式给出动能和质量矩阵：〔1-19〔1-20因此,总体质量矩阵为总体刚度矩阵：求固有频率。总应变能<1-21><1-22>利用Lagrange方程,推导简支梁自由振动方程：<1-23>这里,我们假设简支梁做简谐振动,则<1-24>因此,特征方程为：<1-25>其中,,为固有频率。解方程,有限元分析结果基于上述理论,我们获得了采用MATLAB的有限元分析程序代码。提交边界条件、材料特性和几何参数到上面的方程,使用MATLAB代码我们得到以下结果。〔1离散简支梁为5单元6节点,那么我们得到的单元质量矩阵和刚度矩阵,如下所示：〔2总质量矩阵和总刚度矩阵如下：<3>简支梁的振动分析表1列出了简支梁的五阶固有振动频率。从表中我们可以看出,有限元模拟分析方法和理论值在误差允许范围内是比较吻合的。计算得离散为五单元下的简支梁固有振动频率：计算值〔理论值<>误差率〔%一阶固有振动频率0.01810.01810.0107二阶固有振动频率0.07250.07270.1657三阶固有振动频率0.16320.16450.7942四阶固有振动频率0.29010.29682.3037五阶固有振动频率0.45330.503210.9918离散为30单元的简支梁的振动模式如下图所示。一阶振动图像二阶振动图像三阶振动图像四阶振动图像五阶振动图像六阶振动图像采用Matlab编写的程序代码%%利用有限元方法求解简支梁的振动问题%%clc;clearall;symsxlrhobtEA=b*t;I=b*t^3/12;n=input<'Pleaseinputthenumberofdiscreteelementsn='>;%输入离散化单元的数量n%%定义形函数N1=1-3*<x/l>^2+2*<x/l>^3;N2=<x/l-2*<x/l>^2+<x/l>^3>*l;N3=3*<x/l>^2-2*<x/l>^3;N4=<<x/l>^3-<x/l>^2>*l;%%%%求解单元质量矩阵、刚度矩阵以及总质量矩阵和刚度矩阵N=[N1,N2,N3,N4];Me0=int<N'*N,x,0,l>;Me=rho*A*Me0Ke0=int<diff<N.',x,2>*diff<N,x,2>,x,0,l>Ke=<E*I>*Ke0;M0=zeros<2*<n+1>,2*<n+1>>;K0=zeros<2*<n+1>,2*<n+1>>;fori=1:1:nae=zeros<4,2*<n+1>>;forj=1:1:4ae<j,2*i+j-2>=1;%定义坐标变换矩阵endM0=M0+ae.'*Me*ae;K0=K0+ae.'*Ke*ae;enddisp<'TheelementmassmatrixMe='>;disp<Me>;disp<'TheelementstiffnessmatrixKe='>;disp<Ke>;disp<'ThetotalmassmatrixM='>;disp<M0>;disp<'ThetotalstiffnessmatrixK='>;disp<K0>;%%Me1=matlabFunction<Me>;Ke1=matlabFunction<Ke>;M1=matlabFunction<M0>;%将质量符号矩阵转换为代数矩阵K1=matlabFunction<K0>;%将刚度符号矩阵转换为代数矩阵%%输入的几何参数和物理常数L=0.4;%梁的长度,mb=0.02;%梁的宽度,mt=0.002;%梁的厚度,mE=0.7*10^11;%梁的弹性模量,GParho=2700;%梁的密度,kg/m^3l=L/n;MeNumeric=Me1<b,l,t,rho>;KeNumeric=Ke1<E,b,l,t>;MNumeric=M1<b,l,t,rho>;KNumeric=K1<E,b,l,t>;disp<'ThenumericalelementmassmatrixMe='>;disp<MeNumeric>;%输出单元质量矩阵disp<'ThenumericalelementstiffnessmatrixKe='>;disp<KeNumeric>;%输出单元刚度矩阵disp<'ThenumericaltotalmassmatrixM='>;disp<MNumeric>;%输出总质量矩阵disp<'ThenumericaltotalstiffnessmatrixK='>;disp<KNumeric>;%输出总刚度矩阵%%简支梁的振动求解%%%对于简支梁考虑约束条件：;%采用消元法：消去整体质量矩阵的第1行,第1列；倒数第二行和倒数第二列%消去整体刚度矩阵的第1行,第1列；倒数第二行和倒数第二列Mssb=MNumeric<[2:2*n,2*<n+1>],[2:2*n,2*<n+1>]>;Kssb=KNumeric<[2:2*n,2*<n+1>],[2:2*n,2*<n+1>]>;[Xssb,lamdassb]=eig<Kssb,Mssb>;omegassb=sort<sqrt<diag<lamdassb>>>;XXssb=zeros<n,n>;XXssb<1,:>=0;XXssb<n,:>=0;fori=2:n-1forj=1:nXXssb<i,j>=Xssb<2*<i-1>,j>;endend%%根据振动力学方程得到简支梁固有频率betalssb=[];omega_realssb=[];errorssb=[];fori=1:10betalssb<i>=pi*i;omega_realssb<i>=<betalssb<i>/L>^2*sqrt<E*<b*t^3/12>/<rho*b*t>>;enddisp<'Trueoutput='>;disp<omega_realssb>;%输出理论固有频率disp<'CalculatedOutput='>;disp<omegassb>;%输出计算固有频率值fori=1:10errorssb<i>=<abs<omega_realssb<i>-omegassb<i>>/omega_realssb<i>>*100;%输出误差率enddisp<'Outputpercentageerror='>;disp<errorssb>;%%%画出前五阶振型r=1:n;fori=1:nfigure<'Color',[111]>;plot<r,XXssb<:,i>>;xlabel<'\itx\rm/\itl','FontName','TimesNewRoman','FontSize',20>;ylabel<'u\rm<\itx\rm>','FontName','TimesNewRoman'