基于MATLAB懸臂梁擾度分析

范玲玲，樂騰勝，童智能，項(xiàng)炳泉，張金鐘

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

1 引言

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

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

2 經(jīng)典四階龍格-庫塔算法的介紹

經(jīng)典龍格-庫塔（RK）方法是一種在工程上應(yīng)用廣泛的高精度單步算法。由于此算法精度高，采取措施對(duì)誤差進(jìn)行抑制，所以其實(shí)現(xiàn)原理也較復(fù)雜。同Euler等算法一樣，該算法也是構(gòu)建在數(shù)學(xué)支持的基礎(chǔ)之上的。對(duì)于一階精度的歐拉公式有：

當(dāng)用點(diǎn)Xn處的率近似值K1與右端點(diǎn)Xn+1處的斜率K2的算術(shù)平均值作為平均斜率K*的近似值，那么就會(huì)得到二階精度的改進(jìn)歐拉公式：

依次類推，如果在區(qū)間[Xn，Xn+1]內(nèi)多預(yù)估幾個(gè)點(diǎn)上的斜率值K1、K2、……Km，并用他們的加權(quán)平均數(shù)作為平均斜率K*的近似值，顯然能構(gòu)造出具有很高精度的高階計(jì)算公式。經(jīng)數(shù)學(xué)推導(dǎo)、求解，可以得出四階龍格－庫塔公式，也就是在工程中應(yīng)用廣泛的經(jīng)典四階龍格－庫塔算法[10]：

3 問題描述

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

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

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

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

圖1 懸臂梁示意圖

4 MATLAB近似求解

將區(qū)間[0,2]進(jìn)行100等分，采用經(jīng)典四階RK公式計(jì)算。

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

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

圖2 懸臂梁繞度圖

再運(yùn)行程序得出近似解：最大擾度在端點(diǎn)B出現(xiàn)，即ymax1（2）=16.48cm。

5 解析法精確求解

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

對(duì)（15）式求最值得出：當(dāng)z=2時(shí)，即在懸臂梁最右端，撓度最大，ymax2（2）=16.67cm。

6 結(jié)論