范玲玲,樂騰勝,童智能,項炳泉,張金鐘
(1.南昌市京東學校,江西 南昌 330029;2.安徽省建筑科學研究設計院,安徽 合肥 230001;3.江西科技師范大學建筑工程學院,江西 南昌 330013)
數(shù)學建模就是構(gòu)造數(shù)學模型的過程,即用數(shù)學的語言-公式、符號、圖表等刻畫和描述一個實際問題,然后經(jīng)過數(shù)學的處理-計算、迭代等得到定量的結(jié)果,以供人們做分析、預報、決策和控制[1]。隨著計算機廣泛使用與科學技術迅速發(fā)展,科學計算已是科學研究、工程設計中的一個重要的手段。合理的利用MATLAB數(shù)學軟件來輔助數(shù)學演算和繪圖,已成為與理論分析、科學實驗并駕齊驅(qū)的科學研究方法[2~5]。懸臂梁在各種荷載及不同支撐下的擾度分析是工民建中常見的問題,主要涉及到了常微分方程初值問題的數(shù)值求解。對于一般的初值問題,通常采用改進的Euler公式,就能保證二階精度[6-7]。如果方程右端項f(x,y)足夠光滑,計算精度較高,經(jīng)典四階RK方法是不錯的方法,這一點在分析懸臂梁的擾度問題中得到了很好的應用[8-9]。
本文針對懸臂梁的擾度隨著梁的長度的變化而變化,在MATLAB中運用經(jīng)典四階RK公式計算,作出懸臂梁擾度圖(z-y圖),確定出懸臂梁中擾度最大的位置以及求出相對應的最大擾度近似解,并與撓度微分方程計算得出的精確解相對比,驗證了運用MATLAB軟件近似計算懸臂梁擾度的可行性,促進了懸臂梁擾度計算的研究與分析。
經(jīng)典龍格-庫塔(RK)方法是一種在工程上應用廣泛的高精度單步算法。由于此算法精度高,采取措施對誤差進行抑制,所以其實現(xiàn)原理也較復雜。同Euler等算法一樣,該算法也是構(gòu)建在數(shù)學支持的基礎之上的。對于一階精度的歐拉公式有:
當用點Xn處的率近似值K1與右端點Xn+1處的斜率K2的算術平均值作為平均斜率K*的近似值,那么就會得到二階精度的改進歐拉公式:
依次類推,如果在區(qū)間[Xn,Xn+1]內(nèi)多預估幾個點上的斜率值K1、K2、……Km,并用他們的加權(quán)平均數(shù)作為平均斜率K*的近似值,顯然能構(gòu)造出具有很高精度的高階計算公式。經(jīng)數(shù)學推導、求解,可以得出四階龍格-庫塔公式,也就是在工程中應用廣泛的經(jīng)典四階龍格-庫塔算法[10]:
如圖1所示的一段被嵌入墻內(nèi)的懸臂梁,在固定點A,位移y和斜度dz皆為零。設梁是一根均勻細桿,其長度為L。梁關于垂直方向y的擾度滿足微分方程:
其中:I是梁的橫截面關于其主軸的慣性力矩;E是彈性模量;常數(shù)ρ是梁的線密度;g是重力加速度。
選取參數(shù)L=2m,線密度ρ=10kg/m,慣性力矩與彈性模量的乘積 IE=2400kg·m3/s2[11]。
根據(jù)以上條件,確定出懸臂梁中擾度最大的位置以及求出相對應的最大擾度,作出懸臂梁擾度圖(z-y圖)。
圖1 懸臂梁示意圖
將區(qū)間[0,2]進行100等分,采用經(jīng)典四階RK公式計算。
用RK求解微分方程的MATLAB程序[12]如下:function RK4=RK4(a,b,m)
運行MATLAB程序得出懸臂梁擾度圖(z-y圖)如圖2所示。
圖2 懸臂梁繞度圖
再運行程序得出近似解:最大擾度在端點B出現(xiàn),即ymax1(2)=16.48cm。
將已知參數(shù)L=2m,線密度ρ=10kg/m,慣性力矩與彈性模量的乘積 IE=2400kg·m3/s2代入(11)式,化簡求解得:
對(15)式求最值得出:當z=2時,即在懸臂梁最右端,撓度最大,ymax2(2)=16.67cm。