基于MATLAB懸臂梁擾度分析

范玲玲，樂騰勝，童智能，項炳泉，張金鐘

（1.南昌市京東學校，江西 南昌 330029；2.安徽省建筑科學研究設計院，安徽 合肥 230001；3.江西科技師范大學建筑工程學院，江西 南昌 330013）

1 引言

數學建模就是構造數學模型的過程，即用數學的語言-公式、符號、圖表等刻畫和描述一個實際問題，然后經過數學的處理-計算、迭代等得到定量的結果，以供人們做分析、預報、決策和控制[1]。隨著計算機廣泛使用與科學技術迅速發展，科學計算已是科學研究、工程設計中的一個重要的手段。合理的利用MATLAB數學軟件來輔助數學演算和繪圖，已成為與理論分析、科學實驗并駕齊驅的科學研究方法[2～5]。懸臂梁在各種荷載及不同支撐下的擾度分析是工民建中常見的問題，主要涉及到了常微分方程初值問題的數值求解。對于一般的初值問題，通常采用改進的Euler公式，就能保證二階精度[6-7]。如果方程右端項f（x，y）足夠光滑，計算精度較高，經典四階RK方法是不錯的方法，這一點在分析懸臂梁的擾度問題中得到了很好的應用[8-9]。

本文針對懸臂梁的擾度隨著梁的長度的變化而變化，在MATLAB中運用經典四階RK公式計算，作出懸臂梁擾度圖（z-y圖），確定出懸臂梁中擾度最大的位置以及求出相對應的最大擾度近似解，并與撓度微分方程計算得出的精確解相對比，驗證了運用MATLAB軟件近似計算懸臂梁擾度的可行性，促進了懸臂梁擾度計算的研究與分析。

2 經典四階龍格-庫塔算法的介紹

經典龍格-庫塔（RK）方法是一種在工程上應用廣泛的高精度單步算法。由于此算法精度高，采取措施對誤差進行抑制，所以其實現原理也較復雜。同Euler等算法一樣，該算法也是構建在數學支持的基礎之上的。對于一階精度的歐拉公式有：

當用點Xn處的率近似值K1與右端點Xn+1處的斜率K2的算術平均值作為平均斜率K*的近似值，那么就會得到二階精度的改進歐拉公式：

依次類推，如果在區間[Xn，Xn+1]內多預估幾個點上的斜率值K1、K2、……Km，并用他們的加權平均數作為平均斜率K*的近似值，顯然能構造出具有很高精度的高階計算公式。經數學推導、求解，可以得出四階龍格－庫塔公式，也就是在工程中應用廣泛的經典四階龍格－庫塔算法[10]：

3 問題描述

如圖1所示的一段被嵌入墻內的懸臂梁，在固定點A，位移y和斜度dz皆為零。設梁是一根均勻細桿，其長度為L。梁關于垂直方向y的擾度滿足微分方程：

其中：I是梁的橫截面關于其主軸的慣性力矩；E是彈性模量；常數ρ是梁的線密度；g是重力加速度。

選取參數L=2m，線密度ρ=10kg/m，慣性力矩與彈性模量的乘積 IE=2400kg·m3/s2[11]。

根據以上條件，確定出懸臂梁中擾度最大的位置以及求出相對應的最大擾度，作出懸臂梁擾度圖（z-y圖）。

圖1 懸臂梁示意圖

4 MATLAB近似求解

將區間[0,2]進行100等分，采用經典四階RK公式計算。

用RK求解微分方程的MATLAB程序[12]如下：function RK4=RK4(a,b,m)

運行MATLAB程序得出懸臂梁擾度圖（z-y圖）如圖2所示。

圖2 懸臂梁繞度圖

再運行程序得出近似解：最大擾度在端點B出現，即ymax1（2）=16.48cm。

5 解析法精確求解

將已知參數L=2m，線密度ρ=10kg/m，慣性力矩與彈性模量的乘積 IE=2400kg·m3/s2代入（11）式，化簡求解得：

對（15）式求最值得出：當z=2時，即在懸臂梁最右端，撓度最大，ymax2（2）=16.67cm。

6 結論