西交大有限元原理與應(yīng)用-大作業(yè)

..有限元原理及工程應(yīng)用——大作業(yè)學(xué)院：機械工程學(xué)院班級：碩4002班小組成員：李追3114001089陳草3114001080..作業(yè)題目：利用有限元方法對簡支梁問題進(jìn)行求解,梁的橫截面為矩形,其約束情況如圖1所示。已知梁的幾何尺寸和物理參數(shù)如下：〔1幾何尺寸：長度,截面尺寸；〔2物理參數(shù)：彈性模量GPa,泊松比,密度。圖1．梁及其橫截面示意圖要求：至少劃分五個節(jié)點〔四個單元；給出單元節(jié)點信息；給出單元剛度矩陣和質(zhì)量矩陣；給出總剛度矩陣和總質(zhì)量矩陣；求出梁各界固有頻率及振型〔五階；將所得結(jié)果與理論值進(jìn)行對比,驗證方法的可行性。解：由有限元知識,根據(jù)Rayleigh-Ritz法,解有限元分為四步：建立離散化、單元分析、形成總體方程、解方程,具體步驟如下：〔1建立離散化這里我們將矩形截面簡支梁等分四等分,即分為六節(jié)點的五個桿單元,如圖2所示：每個單元尺寸,這里只考慮桿在豎直平面的彎曲,每個節(jié)點只有y方向位移和繞z軸的旋轉(zhuǎn)自由度?！?單元分析構(gòu)造一組Lagrange插值基函數(shù),在本節(jié)點值為1,其他節(jié)點值為0。從Rayleigh-Ritz法可以看到,插值函數(shù)要p次可微,最高階導(dǎo)數(shù)出現(xiàn)在應(yīng)變能表達(dá)式中；同樣,我們可以這一原則適用于基函數(shù)的選擇以及形狀函數(shù),否則我們將無法正確計算應(yīng)變能當(dāng)我們使用有限元逼近方法。梁的彎曲問題,應(yīng)變能計算公式：〔1-1其中,E為彈性模量,Iz為截面慣性矩。從公式可知,位移函數(shù)必須連續(xù),并且二階導(dǎo)數(shù)平方可積。如圖3,是一維桿單元模型,每個節(jié)點兩個自由度,該單元含有四個自由度,即〔。本題中我們采用三次多項式插值函數(shù)：〔1-2因此,我們必須給出四個形函數(shù)〔位移模式。圖3一維桿單元模型構(gòu)造Hermite插值函數(shù)。選擇局部坐標(biāo)系<,>,其中l(wèi)是單元長度,轉(zhuǎn)角是撓度值的一階導(dǎo)數(shù),定義邊界條件：〔1-3因此,我們給出變形的Hermite的多項式插值函數(shù)：〔1-4其中,和分別滿足如下條件,對應(yīng)的圖形如圖4所示：〔1-5圖4一維Hermite插值多項式基于Langrage和Hermite插值多項式,寫出單元形函數(shù)〔1-6節(jié)點位移值也可以得出〔1-7同時,表達(dá)式〔1-7用矩陣表示為〔1-8其中,,用能量表達(dá)式替代表達(dá)式中的和。動能表達(dá)式：〔1-9將〔1-8帶入〔1-9,得到〔1-10從而獲得質(zhì)量矩陣：〔1-11帶入,〔1-12應(yīng)變能表達(dá)式：〔1-13剛度矩陣表達(dá)式：〔1-14帶入,可以得到〔1-15形成總體方程將每個桿單元的能量方程組裝。完整梁上的的總動能和能量的和所做的總功梁上的外力作用,所有的自由度的位移矢量可以給出：〔1-16將位移矢量轉(zhuǎn)換為全局坐標(biāo)系下的位移矢量,變換矩陣為：〔1-17〔1-18分別以矩陣形式給出動能和質(zhì)量矩陣：〔1-19〔1-20因此,總體質(zhì)量矩陣為總體剛度矩陣：求固有頻率。總應(yīng)變能<1-21><1-22>利用Lagrange方程,推導(dǎo)簡支梁自由振動方程：<1-23>這里,我們假設(shè)簡支梁做簡諧振動,則<1-24>因此,特征方程為：<1-25>其中,,為固有頻率。解方程,有限元分析結(jié)果基于上述理論,我們獲得了采用MATLAB的有限元分析程序代碼。提交邊界條件、材料特性和幾何參數(shù)到上面的方程,使用MATLAB代碼我們得到以下結(jié)果?！?離散簡支梁為5單元6節(jié)點,那么我們得到的單元質(zhì)量矩陣和剛度矩陣,如下所示：〔2總質(zhì)量矩陣和總剛度矩陣如下：<3>簡支梁的振動分析表1列出了簡支梁的五階固有振動頻率。從表中我們可以看出,有限元模擬分析方法和理論值在誤差允許范圍內(nèi)是比較吻合的。計算得離散為五單元下的簡支梁固有振動頻率：計算值〔理論值<>誤差率〔%一階固有振動頻率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='>;%輸入離散化單元的數(shù)量n%%定義形函數(shù)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;%%%%求解單元質(zhì)量矩陣、剛度矩陣以及總質(zhì)量矩陣和剛度矩陣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;%定義坐標(biāo)變換矩陣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>;%將質(zhì)量符號矩陣轉(zhuǎn)換為代數(shù)矩陣K1=matlabFunction<K0>;%將剛度符號矩陣轉(zhuǎn)換為代數(shù)矩陣%%輸入的幾何參數(shù)和物理常數(shù)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>;%輸出單元質(zhì)量矩陣disp<'ThenumericalelementstiffnessmatrixKe='>;disp<KeNumeric>;%輸出單元剛度矩陣disp<'ThenumericaltotalmassmatrixM='>;disp<MNumeric>;%輸出總質(zhì)量矩陣disp<'ThenumericaltotalstiffnessmatrixK='>;disp<KNumeric>;%輸出總剛度矩陣%%簡支梁的振動求解%%%對于簡支梁考慮約束條件：;%采用消元法：消去整體質(zhì)量矩陣的第1行,第1列；倒數(shù)第二行和倒數(shù)第二列%消去整體剛度矩陣的第1行,第1列；倒數(shù)第二行和倒數(shù)第二列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%%根據(jù)振動力學(xué)方程得到簡支梁固有頻率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'