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，市场有如下几种尺寸的杆件可供选择···。提出你的桥梁设计方案，如何在满足约束条件下成本最低？这样的题目既和实际接近，又把传统的解题变为设计和优化问题了，而这正是目前传统力学教育所缺乏的。