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

高云峰

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

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

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

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

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

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

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

(2) 運動學專題1 篇：典型機構的運動分析(著重介紹Matlab 中非線性方程組的求解、動畫顯示，如何對運動方程求導數)；

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

(4) 綜合運用專題1 篇：數據轉換問題(著重介紹在不同坐標系中看到結果，包括運動和動力學問題)。

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

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

1 Matlab 中代數方程的求解

靜力學問題的求解一般可以化為代數方程的求解，代數方程一般可以寫為

其中A是n×n階的矩陣，由系統的位置、尺寸等參數構成，X是n×1 階的列陣，由系統中待求解的未知數構成，B是n×1 階的列陣，由系統中已知載荷、尺寸等參數構成。具體內容和形式見下面案例中的具體表達式。

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

其中inv 是Matlab 中矩陣求逆的函數[2]，運行后就能直接解出系統中所有的未知力。

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

圖1 桁架系統

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

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

圖2 特殊截面法

圖3 所有桿件編號

圖4 A 節點受力圖

圖5 E 節點受力圖

對每個節點，根據水平和豎直方向力的平衡方程，分別有

其他節點的平衡方程類似，最后合在一起，寫為AX=B的形式，有

在Matlab 中運行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 命令把具體數值轉換為符號。如果想更簡單些，X=inv(A)*B 后面不寫分號“；”就能直接顯示結果(如顯示為3.4509，而不是S1=?3.4509 N)。

其解答的截圖見圖7。

圖6 桁架問題源代碼

圖7 全部解答的截圖

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

2 Matlab 中帶參數代數方程的求解

Matlab 功能強大，除了可以進行數值計算，還可以進行符號推導。因此，如果某些靜力學問題沒有具體數值，也可以進行求解。

案例 2：橫梁桁架結構由橫梁AC和BC及五根細支撐桿組成，所受載荷及尺寸如圖 8 所示。求1，2，3 桿的內力。

圖8 剛體系統

從理論力學教學的角度，希望學生適當地取分離體，但是有一定的技巧，解出的答案是(具體分析過程略)

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

圖9 D 節點受力圖

圖10 AC 桿件受力圖

對D節點列寫水平和豎直方向力的平衡方程，有

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

其他桿件和節點的平衡方程類似，最后合在一起，把未知數放在方程一側，把已知載荷有關的量放在方程另一側，寫為AX=B的形式，有

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

圖11 求解帶參數代數方程的源代碼

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

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

圖12 解答表達式截圖

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

3 小結

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

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

圖13 帶參數的A 矩陣的逆

數值計算看起來輸入的工作量較大，卻是一種通用的方法，關鍵的是可以獲得系統所有部件的受力 (包含數值解和符號解)，為后續進一步分析打下了基礎(如強度分析、結構優化、失穩等等)，也直接為后續的結構力學打下了基礎。

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

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