Matlab 求解理論力學(xué)問題系列 (一)剛體系統(tǒng)及桁架受力問題

高云峰

(清華大學(xué)航天航空學(xué)院, 北京100084)

如果在理論力學(xué)教學(xué)中引入Matlab，根據(jù)經(jīng)驗(yàn)，只需要三次課，就可以讓學(xué)生掌握代數(shù)方程和微分方程的數(shù)值求解、符號推導(dǎo)、動畫演示等，讓學(xué)生對理論力學(xué)問題的理解有飛躍式的提升；而教學(xué)中某些解題技巧性的內(nèi)容則可以壓縮，保持總學(xué)時不變。具體來說：

(1)在靜力學(xué)中，以往對于復(fù)雜系統(tǒng)的受力分析通常要適當(dāng)取分離體，有時需要高度的技巧[1]；同時由于傳統(tǒng)計(jì)算能力的限制，往往只要求解出某些部件的受力；如果采用Matlab 處理，可以采用統(tǒng)一的處理方式，把系統(tǒng)全部拆開，快速求出所有部件的受力，對系統(tǒng)的整體和各部件受力有更全面的了解。

(2) 在運(yùn)動學(xué)中，以往分析系統(tǒng)運(yùn)動時，強(qiáng)調(diào)求特定時刻或特定位置某點(diǎn)或剛體的速度和加速度，而系統(tǒng)的整體運(yùn)動特點(diǎn)、某些點(diǎn)的運(yùn)動軌跡有時難以想象；而采用Matlab 處理，可以求出系統(tǒng)任意點(diǎn)或剛體在任意時刻的速度和加速度等運(yùn)動量，特別是其畫圖和動畫演示功能，可以快速直觀地顯示系統(tǒng)的整個運(yùn)動過程、給出任意點(diǎn)的運(yùn)動軌跡。

(3)在動力學(xué)中，以往絕大部分問題都只能列寫動力學(xué)方程，通常沒有解析解，傳統(tǒng)數(shù)學(xué)分析的方法也用不上，系統(tǒng)豐富復(fù)雜的動力學(xué)現(xiàn)象很難從方程中看出；而采用Matlab 處理，可以獲得系統(tǒng)整個運(yùn)動過程中的受力、速度和加速度等量，還可以快速直觀地演示系統(tǒng)的運(yùn)動過程。

考慮到目前理論力學(xué)教學(xué)中對于數(shù)值計(jì)算、符號推導(dǎo)很少介紹，為此專門準(zhǔn)備系列理論力學(xué)教學(xué)文章，每篇介紹1～2 個典型的理論力學(xué)問題及如何利用Matlab 進(jìn)行處理。系列文章具體計(jì)劃分為如下專題：

(1)靜力學(xué)專題1 篇：剛體系統(tǒng)及桁架的受力問題(著重介紹Matlab 中代數(shù)方程的數(shù)值求解和符號求解)；

(2) 運(yùn)動學(xué)專題1 篇：典型機(jī)構(gòu)的運(yùn)動分析(著重介紹Matlab 中非線性方程組的求解、動畫顯示，如何對運(yùn)動方程求導(dǎo)數(shù))；

(3)動力學(xué)專題2 篇：單擺和橢圓擺的運(yùn)動和周期(著重介紹Matlab 中微分方程的數(shù)值求解、計(jì)算可靠性、根據(jù)數(shù)據(jù)的快速傅里葉變換求周期)、乒乓球滾動問題(著重介紹Matlab 中分段積分的處理方法，以及與分段對應(yīng)的積分中斷點(diǎn)問題)；

(4) 綜合運(yùn)用專題1 篇：數(shù)據(jù)轉(zhuǎn)換問題(著重介紹在不同坐標(biāo)系中看到結(jié)果，包括運(yùn)動和動力學(xué)問題)。

通過這幾篇文章，可以讓學(xué)生們了解、熟悉Matlab，大大提高解決問題的能力。

下面首先從靜力學(xué)開始。根據(jù)教學(xué)經(jīng)驗(yàn)，學(xué)生在靜力學(xué)中對于剛體系統(tǒng)和桁架的受力分析感到相對比較困難，通常要適當(dāng)拆開，否則解不出來。而Matlab 可以采用統(tǒng)一的方法求解，降低了解題的技巧，但是得到的解答更全面。

1 Matlab 中代數(shù)方程的求解

靜力學(xué)問題的求解一般可以化為代數(shù)方程的求解，代數(shù)方程一般可以寫為

其中A是n×n階的矩陣，由系統(tǒng)的位置、尺寸等參數(shù)構(gòu)成，X是n×1 階的列陣，由系統(tǒng)中待求解的未知數(shù)構(gòu)成，B是n×1 階的列陣，由系統(tǒng)中已知載荷、尺寸等參數(shù)構(gòu)成。具體內(nèi)容和形式見下面案例中的具體表達(dá)式。

列寫力和力矩的平衡方程是理論力學(xué)教學(xué)中的重點(diǎn)。一旦有了平衡方程就可以獲得式(1) 中的A矩陣和B列陣，而Matlab 處理矩陣運(yùn)算特別方便，其求解格式為

其中inv 是Matlab 中矩陣求逆的函數(shù)[2]，運(yùn)行后就能直接解出系統(tǒng)中所有的未知力。

案例 1：圖示桁架系統(tǒng)中 (圖 1)，ABC是正三角形，邊長為 1 m，DEF也是正三角形，且∠ACD= ∠BAE= ∠CBF=15?，水平力P=10 N，垂直力Q=20 N，求 1，2，3 桿的內(nèi)力[3]。

圖1 桁架系統(tǒng)

從理論力學(xué)教學(xué)的角度，希望學(xué)生采用特殊截面法，把 1，2，3 桿截斷，把三角形DEF“挖出來”(圖2)，把DEF看作剛體，三個未知數(shù)正好可以求解 (求解過程略)。但是這個特殊截面對很多同學(xué)而言有一定的難度，不容易想到。而且每個桁架問題都可能有特殊性，求解時需要經(jīng)驗(yàn)和技巧。

而采用節(jié)點(diǎn)法就不需要什么技巧：將所有的桿件都編號(圖3)，全部拆開，設(shè)桿件受拉為正，對各節(jié)點(diǎn)列寫平衡方程(為節(jié)省篇幅只以外部A節(jié)點(diǎn)和內(nèi)部E節(jié)點(diǎn)為例，見圖4 和圖 5)。

圖2 特殊截面法

圖3 所有桿件編號

圖4 A 節(jié)點(diǎn)受力圖

圖5 E 節(jié)點(diǎn)受力圖

對每個節(jié)點(diǎn)，根據(jù)水平和豎直方向力的平衡方程，分別有

其他節(jié)點(diǎn)的平衡方程類似，最后合在一起，寫為AX=B的形式，有

在Matlab 中運(yùn)行X=inv(A)*B，就得到所有桿件以及A和B鉸鏈處的力，具體桁架問題求解程序源代碼見圖6。源代碼中clc 是清除屏幕；Matlab 中表示矩陣很方便，例如[1, 2, 3] 是3×1 的行陣，而[1; 2; 3]是1×3 的列陣；zeros(12,12)是生成12×12的零矩陣，里面元素全是0；A(i,j)表示A矩陣中第i行第j列的元素；在屏幕上顯示的格式為 disp(可以顯示特定的文字)，所以用num2str 命令把具體數(shù)值轉(zhuǎn)換為符號。如果想更簡單些，X=inv(A)*B 后面不寫分號“；”就能直接顯示結(jié)果(如顯示為3.4509，而不是S1=?3.4509 N)。

其解答的截圖見圖7。

圖6 桁架問題源代碼

圖7 全部解答的截圖

因此使用Matlab 求解靜力學(xué)問題，關(guān)鍵是確定A矩陣和B列陣，而這與列寫平衡方程有關(guān)。

2 Matlab 中帶參數(shù)代數(shù)方程的求解

Matlab 功能強(qiáng)大，除了可以進(jìn)行數(shù)值計(jì)算，還可以進(jìn)行符號推導(dǎo)。因此，如果某些靜力學(xué)問題沒有具體數(shù)值，也可以進(jìn)行求解。

案例 2：橫梁桁架結(jié)構(gòu)由橫梁AC和BC及五根細(xì)支撐桿組成，所受載荷及尺寸如圖 8 所示。求1，2，3 桿的內(nèi)力。

圖8 剛體系統(tǒng)

從理論力學(xué)教學(xué)的角度，希望學(xué)生適當(dāng)?shù)厝》蛛x體，但是有一定的技巧，解出的答案是(具體分析過程略)

如果采用Matlab 處理，則全部拆開，對節(jié)點(diǎn)和剛體分別列寫平衡方程，然后獲得A和B矩陣。為節(jié)省篇幅只畫出D節(jié)點(diǎn)和AC桿件的受力圖，見圖 9 和圖 10。

圖9 D 節(jié)點(diǎn)受力圖

圖10 AC 桿件受力圖

對D節(jié)點(diǎn)列寫水平和豎直方向力的平衡方程，有

對AC部件列寫水平和豎直方向力的平衡方程，再對A點(diǎn)取矩，有

其他桿件和節(jié)點(diǎn)的平衡方程類似，最后合在一起，把未知數(shù)放在方程一側(cè)，把已知載荷有關(guān)的量放在方程另一側(cè)，寫為AX=B的形式，有

利用X=inv(A)*B，可以求出解析解，整個程序的源代碼如圖 11。源代碼中用 syms 命令來定義符號，變量涉及到符號運(yùn)算時都需要先定義；simplify是化簡命令，可以自動化簡、合并表達(dá)式，例如simplify((cos(y))?2+(sin(y))?2)會自動化簡為 1；disp中的char 表示字符串。解的結(jié)果見圖12。

圖11 求解帶參數(shù)代數(shù)方程的源代碼

對比一下，可以看出Matlab 符號推導(dǎo)得到的前3 個解與傳統(tǒng)方法得到的式(6) 相同。

如果關(guān)心A矩陣的逆是怎樣的形式，可以單獨(dú)運(yùn)行inv(A)，圖13 為屏幕截圖。

圖12 解答表達(dá)式截圖

為了驗(yàn)證矩陣求逆是否正確，可以查看inv(A)?A 的結(jié)果，的確顯示為單位矩陣。

3 小結(jié)

理論力學(xué)教學(xué)應(yīng)注意基本概念、基本思路和基本方法，而具體繁瑣的計(jì)算工作可以交給數(shù)學(xué)軟件，這樣可以讓學(xué)生掌握最一般的靜力學(xué)分析、計(jì)算方法。

本篇介紹了Matlab 求解靜力學(xué)問題的方法，核心的函數(shù)是 inv(矩陣求逆)，其他相關(guān)的函數(shù)包括syms(定義符號)、simplify(符號化簡)和 disp(在屏幕上顯示)。利用這些函數(shù)，可以完成靜力學(xué)中代數(shù)方程的數(shù)值求解及帶參數(shù)的符號推導(dǎo)。

圖13 帶參數(shù)的A 矩陣的逆

數(shù)值計(jì)算看起來輸入的工作量較大，卻是一種通用的方法，關(guān)鍵的是可以獲得系統(tǒng)所有部件的受力 (包含數(shù)值解和符號解)，為后續(xù)進(jìn)一步分析打下了基礎(chǔ)(如強(qiáng)度分析、結(jié)構(gòu)優(yōu)化、失穩(wěn)等等)，也直接為后續(xù)的結(jié)構(gòu)力學(xué)打下了基礎(chǔ)。

傳統(tǒng)的建模方法，都需要針對具體問題，按一定的步驟推導(dǎo)才能得到靜力學(xué)或動力學(xué)方程。每遇到一個新的問題，由于系統(tǒng)的結(jié)構(gòu)不一樣，要按相同的步驟重復(fù)一遍。是否有一種方法可建立一個適合于任意系統(tǒng)的一般公式，只要把系統(tǒng)的最基本的一些參數(shù)，如剛體數(shù)目、連接類型、連接點(diǎn)位置、受力情況等帶入公式，就可以展開得到系統(tǒng)的動力學(xué)方程？本質(zhì)上就是系統(tǒng)的A和B矩陣如何生產(chǎn)，能否自動生成？從圖論的角度引入通路矩陣和聯(lián)通矩陣后，可以自動獲得系統(tǒng)的A和B矩陣[4]，而這正是一些商業(yè)力學(xué)軟件(例如Adams)的處理思路。也就是說，Matlab 處理力學(xué)問題，為學(xué)生打開了一個通向處理實(shí)際復(fù)雜問題的窗口。

可以想象，如果學(xué)生掌握了數(shù)值求解和符號推導(dǎo)，只要是靜力學(xué)能處理的問題，都可以很快獲得全部的解答，也許未來的靜力學(xué)題目可能就要換一種問法了，例如：已知兩岸的距離為100 m，如果要架一座桁架橋，要求最大承重是G，每根桿的最大受力為S，單位重量的桿件價格為p，市場有如下幾種尺寸的桿件可供選擇···。提出你的橋梁設(shè)計(jì)方案，如何在滿足約束條件下成本最低？這樣的題目既和實(shí)際接近，又把傳統(tǒng)的解題變?yōu)樵O(shè)計(jì)和優(yōu)化問題了，而這正是目前傳統(tǒng)力學(xué)教育所缺乏的。