FEM在求解连续梁振动问题中的应用

ApplicationofFEMinsolvingthecontinuousbeamvibrationproblem)Ifwewrite,thisisasetoflinearhomogenousequationsintheunknowns,,…,.Bysolvingthefollowingeigenequation,wewillhavenrootsalsocalledtheeigenvalue,,…,.Notingthat,thenaturalfrequencyisgivenby.TheresultsoftheFEManalysisBasedonthetheorymentionedabove,weobtainaFEMcodeaccordingtofigure5usingMATLAB.Submittingtheboundarycondition,thematerialpropertiesandgeometryparametersintotheaboveequations,weobtainthefollowingresultsusingtheMATLABcode.Discretethebeaminto5elements,thenwegettheelementmassmatrixandstiffnessmatrixasfollows FigureSEQFigure\*ARABIC5theprogramchartThetotalmassmatrixandthetotalstiffnessmatrixThecantileverbeamTable1listthetop5inherentfrequenciesofthecantileverbeam.Aswecanseeintable1,theFEMresultsagreedwellwiththetheoreticalvalue.Inaddition,theerrorachievedbydividingthebeaminto5elementsislargerthantheerrorobtainedfromthemodelwith50discreteelements.Thevibrationmode(basedonthedatafromthemodelwith50discreteelements)correspondingtothetop6naturalfrequenciesarepresentedfigure6.TableSEQTable\*ARABIC1thenaturalfrequencyofcantileverbeamOrderTheoreticalvalueFEMresults（5elements）FEMresults（50elements）Error(%)Error(%)132.912032.91610.01245332.91570.011101602206.271206.3820.0538736206.2790.003888393577.473579.6620.379044577.5870.019798941131.851145.111.172051131.840.00044828151871.011900.601.581491871.020.000462699ThefirstordervibrationmodeThesecondordervibrationmodeThethirdordervibrationmodeTheforthordervibrationmodeThefifthordervibrationmodeThesixthordervibrationmodeFigureSEQFigure\*ARABIC6vibrationmodeofthecantileverbeamThesimplysupportedbeamTable2listthetop5inherentfrequenciesofthesimplysupportedbeam.Aswecanseeintable2,theFEMresultsagreedwellwiththetheoreticalvalue.Additionally,theerrorachievedbydividingthebeaminto5elementsislargerthantheerrorobtainedfromthemodelwith50discreteelements.Thevibrationmode(basedonthedatafromthemodelwith50discreteelements)correspondingtothetop6naturalfrequenciesarepresentedfigure7.TableSEQTable\*ARABIC2thenaturalfrequencyofsimplysupportedbeamOrderTheoreticalvalueFEMresults（5elements）FEMresults（50elements）Error(%)Error(%)192.405592.39560.010702592.39560.01110162370.195369.5830.165673369.5830.003888393838.165831.5610.794154831.5620.019798941512.391478.332.303681478.330.00044828152563.792309.8910.99182309.910.000462699ThefirstordervibrationmodeThesecondordervibrationmodeThethirdordervibrationmodeTheforthordervibrationmodeThefifthordervibrationmodeThesixthordervibrationmodeFigureSEQFigure\*ARABIC7vibrationmodeofthesimplysupportedbeamAppendix:MATLABcode%%利用有限元方法求解悬臂梁和简支梁的振动问题%%clc;clearall;symsxlrhobtEA=b*t;I=b*t^3/12;n=input('Pleaseinputthenumberofdiscreteelementsn=');%inputthetotalelementsofthediscreteelement%%difinetheshapefunctionN1=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;%%%%solvetheelementmassmatrixandstiffnessmatix,totalmassmatrixandstiffnessmatrixN=[N1,N2,N3,N4];Me0=int(N'*N,x,0,l);%obtaintheelementmassmatrixMe=rho*A*Me0Ke0=int(diff(N.',x,2)*diff(N,x,2),x,0,l)%obtaintheelementstiffnessmatrixKe=(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;%difinethecoordinatetransformationmatrixendM0=M0+ae.'*Me*ae;%solvethetotalmassstiffnessmatrixK0=K0+ae.'*Ke*ae;%solvethetotalstiffnessmatixenddisp('TheelementmassmatrixMe=');disp(Me);%showtheelementmassmatrixdisp('TheelementstiffnessmatrixKe=');disp(Ke);%showtheelementtiffnessmatrixdisp('ThetotalmassmatrixM=');disp(M0);%showthetotalmassmatrixdisp('ThetotalstiffnessmatrixK=');disp(K0);%showthetotaltiffnessmatrix%%Me1=matlabFunction(Me);Ke1=matlabFunction(Ke);M1=matlabFunction(M0);%将质量符号矩阵转换为代数矩阵K1=matlabFunction(K0);%将刚度符号矩阵转换为代数矩阵%%inputthegeometryparametersandphysicalconstantsL=0.4;%thetotallengthofthebeam,mb=0.02;%thewidthofthebeam,mt=0.001;%thethicknessofthebeam,mE=2.1*10^11;%theelasticmodulus,GParho=7800;%thedensityofthebeam,kg/m^3l=L/n;MNumeric=M1(b,l,t,rho);KNumeric=K1(E,b,l,t);disp('ThenumericaltotalmassmatrixM=');disp(MNumeric);%showthetotalmassmatrixdisp('ThenumericaltotalstiffnessmatrixK=');disp(KNumeric);%showthetotaltiffnessmatrix%%1.悬臂梁的振动求解%对于悬臂梁考虑约束条件节点1：挠度为零，转角为零；%采用消元法：消去整体质量矩阵的前两行，前两列；消去整体刚度矩阵的前两行，前两列；M=MNumeric([3:2*n+2],[3:2*n+2]);K=KNumeric([3:2*n+2],[3:2*n+2]);[X,lamda]=eig(K,M);omega=sort(sqrt(diag(lamda)));XX=zeros(n,n);XX(1,:)=0;fori=2:nforj=1:nXX(i,j)=X(2*(i-1)+1,j);endend%%画出前五阶振型r=1:n;fori=1:6figure('Color',[111]);plot(r,XX(:,i));xlabel('\itx\rm/\itl','FontName','TimesNewRoman','FontSize',36);ylabel('u\rm(\itx\rm)','FontName','TimesNewRoman','FontSize',36);set(gca,'Linewidth',2);%settingupthelinewidthofthecoordinateset(gca,'FontName','TimesNewRoman','FontSize',36);set(get(gca,'Children'),'linewidth',3);%settingupthelinewidthofthecurveend%%根据振动力学方程得到悬臂梁固有频率betal1=1.875;betal2=4.694;betal=[];omega_real=[];error=[];fori=1:20betal(1)=betal1;betal(2)=betal2;ifi>=3betal(i)=pi*(i-0.5);endomega_real(i)=(betal(i)/L)^2*sqrt(E*(b*t^3/12)/(rho*b*t));endifn>=20m=20;elsem=n;endfori=1:merror(i)=(abs(omega_real(i)-omega(i))/omega_real(i))*100;%outputthepercentageerrorend%%简支梁的振动求解%%%对于简支梁考虑约束条件：第一个节点的挠度为零，最后一个节点的挠度为零；%采用消元法：消去整体质量矩阵的第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%%根据振动力学方程得到简