如果在理论力学教学中引入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中代数方程的求解
静力学问题的求解一般可以化为代数方程的求解,代数方程一般可以写为
$\begin{eqnarray} AX=B \end{eqnarray}$
其中$A$是$n\times n$阶的矩阵,由系统的位置、尺寸等参数构成,$X$是$n\times1$阶的列阵,由系统中待求解的未知数构成,$B$是$n\times 1$阶的列阵,由系统中已知载荷、尺寸等参数构成。具体内容和形式见下面案例中的具体表达式。
列写力和力矩的平衡方程是理论力学教学中的重点。一旦有了平衡方程就可以获得式(1)中的$A$矩阵和$B$列阵,而Matlab处理矩阵运算特别方便,其求解格式为
$\begin{eqnarray} {\rm X=inv(A)^\ast B} \end{eqnarray}$
其中inv是Matlab中矩阵求逆的函数[2],运行后就能直接解出系统中所有的未知力。
案例1：图示桁架系统中(图1),$ABC$是正三角形,边长为1 m,$DEF$也是正三角形,且$\angle ACD=\angle BAE=\angle CBF=15^{\circ}$,水平力$P=10$ N,垂直力$Q=20$ N,求1,2,3杆的内力[3]。
图1
从理论力学教学的角度,希望学生采用特殊截面法,把1,2,3杆截断,把三角形$DEF$"挖出来"(图2),把$DEF$看作刚体,三个未知数正好可以求解(求解过程略)。但是这个特殊截面对很多同学而言有一定的难度,不容易想到。而且每个桁架问题都可能有特殊性,求解时需要经验和技巧。
图2
而采用节点法就不需要什么技巧：将所有的杆件都编号(图3),全部拆开,设杆件受拉为正,对各节点列写平衡方程(为节省篇幅只以外部$A$节点和内部$E$节点为例,见图4和图5)。
图3
图4
图5
对每个节点,根据水平和竖直方向力的平衡方程,分别有
$\begin{eqnarray} \label{eq1} \left.\begin{array}{l} S_{1} \cdot \cos 15^{\circ}+S_{7} +S_{9} \cdot \cos 60^{\circ}+X_{A} =0 \\ S_{1} \cdot \sin 15^{\circ}+S_{9} \cdot \sin 60^{\circ}+Y_{A} =0 \\ \end{array}\right\} \end{eqnarray}$
$\begin{eqnarray} \left.\begin{array}{l} -S_{1} \!\cdot\! \cos 15^{\circ}-S_{4} \!\cdot\! \cos 60^{\circ}+ S_{5} \!\cdot\! \cos 60^{\circ}=0 \\ -S_{1} \!\cdot\! \sin 15^{\circ}+S_{4} \!\cdot\! \sin 60^{\circ}+ S_{5} \!\cdot\! \sin 60^{\circ}-Q=0 \\ \end{array}\right\} \end{eqnarray}$
其他节点的平衡方程类似,最后合在一起,写为$AX=B$的形式,有
$\begin{eqnarray} \label{eq3} &&\left[\begin{array}{*{20}c} {\cos 15^{\circ}} & 0 & 0 & 0 & 0 & 0 & 1 & 0 & {\cos 60^{\circ}} & 1 & 0 & 0 \\ {\sin 15^{\circ}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {\sin 60^{\circ}} & 0 & 1 & 0 \\ 0 & {-\cos 45^{\circ}} & 0 & 0 & 0 & 0 & {-1} & {-\cos 60^{\circ}} & 0 & 0 & 0 & 0 \\ 0 & {\sin 45^{\circ}} & 0 & 0 & 0 & 0 & 0 & {\sin 60^{\circ}} & 0 & 0 & 0 & 1 \\ 0 & 0 & {-\cos 75^{\circ}} & 0 & 0 & 0 & 0 & {\cos 60^{\circ}} & {-\cos 60^{\circ}} & 0 & 0 & 0 \\ 0 & 0 & {-\sin 75^{\circ}} & 0 & 0 & 0 & 0 & {-\sin 60^{\circ}} & {-\sin 60^{\circ}} & 0 & 0 & 0 \\ 0 & 0 & {\cos 75^{\circ}} & {\cos 60^{\circ}} & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & {\sin 75^{\circ}} & {-\sin 60^{\circ}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ {-\cos 15^{\circ}} & 0 & 0 & {-\cos 60^{\circ}} & {\cos 60^{\circ}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ {-\sin 15^{\circ}} & 0 & 0 & {\sin 60^{\circ}} & {\sin 60^{\circ}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & {\cos 45^{\circ}} & 0 & 0 & {-\cos 60^{\circ}} & {-1} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & {-\sin 45^{\circ}} & 0 & 0 & {-\sin 60^{\circ}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array} \right] \\&& \left[\begin{array}{*{20}c} {S_{1}} \\ {S_{2}} \\ {S_{3}} \\ {S_{4}} \\ {S_{5}} \\ {S_{6}} \\ {S_{7}} \\ {S_{8}} \\ {S_{9}} \\ {X_{A}} \\ {Y_{A}} \\ {Y_{B}} \\ \end{array} \right]\!=\!\left[\begin{array}{*{20}c} 0 \\ 0 \\ 0 \\ 0 \\ {-P} \\ 0 \\ 0 \\ 0 \\ 0 \\ Q \\ 0 \\ 0 \\ \end{array} \right] \end{eqnarray}$
在Matlab中运行X$=$inv(A)*B,就得到所有杆件以及$A$和$B$铰链处的力,具体桁架问题求解程序源代码见图6。源代码中clc是清除屏幕;Matlab中表示矩阵很方便,例如[1, 2, 3]是$3\times1$的行阵,而[1; 2; 3]是$1\times 3$的列阵;zeros(12, 12)是生成$12\times 12$的零矩阵,里面元素全是0;$A(i,j)$表示$A$矩阵中第$i$行第$j$列的元素;在屏幕上显示的格式为disp(可以显示特定的文字),所以用num2str命令把具体数值转换为符号。如果想 更简单些,X=inv(A)*B后面不写分号";"就能直接显示结果(如显示为3.4509,而不是$S_1=-3.4509$ N)。
图6
图7
因此使用Matlab求解静力学问题,关键是确定$A$矩阵和$B$列阵,而这与列写平衡方程有关。
2 Matlab中带参数代数方程的求解
Matlab功能强大,除了可以进行数值计算,还可以进行符号推导。因此,如果某些静力学问题没有具体数值,也可以进行求解。
案例2：横梁桁架结构由横梁$AC$和$BC$及五根细支撑杆组成,所受载荷及尺寸如图8所示。求1,2,3杆的内力。
图8
从理论力学教学的角度,希望学生适当地取分离体,但是有一定的技巧,解出的答案是(具体分析过程略)
$\begin{eqnarray} \left.\begin{array}{l} S_{1} =\dfrac{\sqrt 2}2(P+3qa)\\ S_{2} =-\dfrac12(P+3qa)\\ S_{3} =\dfrac12(P+3qa)\\ \end{array}\right\} \end{eqnarray}$
如果采用Matlab处理,则全部拆开,对节点和刚体分别列写平衡方程,然后获得$A$和$B$矩阵。为节省篇幅只画出$D$节点和$AC$杆件的受力图,见图9和图10。
图9
图10
$\begin{eqnarray} \left.\begin{array}{l} -\dfrac{\sqrt 2}2S_{1} +S_{3} =0 \\ \dfrac{\sqrt 2}2S_{1} +S_{2} =0 \\ \end{array}\right\} \end{eqnarray}$
对$AC$部件列写水平和竖直方向力的平衡方程,再对$A$点取矩,有
$\begin{eqnarray} \left.\begin{array}{l} \dfrac{\sqrt 2}2S_{1} +X_{A} +X_{C} =0 \\ -\dfrac{\sqrt 2}2S_{1} -S_{2} +Y_{A} +Y_{C} -P-qa=0 \\ -S_{2} a+2Y_{C} a-Pa-qa\cdot \dfrac32a=0 \\ \end{array}\right\} \end{eqnarray}$
其他杆件和节点的平衡方程类似,最后合在一起,把未知数放在方程一侧,把已知载荷有关的量放在方程另一侧,写为$AX=B$的形式,有
$\begin{eqnarray} \label{eq4} \left[ {{\begin{array}{c@{\quad\ }c@{\quad\ }c@{\quad\ }c@{\quad\ }c@{\quad\ }c@{\quad\ }c@{\quad\ }c@{\quad\ }c@{\quad\ }c} -\dfrac{\sqrt2}2 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \dfrac{\sqrt 2}2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & {-1} & 0 & \dfrac{\sqrt 2}2 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & \dfrac{\sqrt 2}2 & 0 & 0 & 0 & 0 & 0 \\ \dfrac{\sqrt 2}2 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 \\ -\dfrac{\sqrt2}2 & {-1} & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 \\ 0 & {-a} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {2a} \\ 0 & 0 & 0 & 0 & -\dfrac{\sqrt2}2 & 0 & 0 & 0 & {-1} & 0 \\ 0 & 0 & 0 & {-1} & -\dfrac{\sqrt2}2 & 0 & 0 & 1 & 0 & {-1} \\ 0 & 0 & 0 & a & 0 & 0 & 0 & 0 & 0 & {2a} \\ \end{array} }} \right]\left[ {{\begin{array}{*{20}c} {S_{1}} \\ {S_{2}} \\ {S_{3}} \\ {S_{4}} \\ {S_{5}} \\ {X_{A}} \\ {Y_{A}} \\ {Y_{B}} \\ {X_{C}} \\ {Y_{C}} \\ \end{array} }} \right]\mbox{=}\left[ {{\begin{array}{*{20}c} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ {P+qa} \\ {Pa+\dfrac32qa^{2}} \\ 0 \\ {qa} \\ {-\dfrac32qa^{2}} \\ \end{array} }} \right] \end{eqnarray}$
利用X=inv(A)*B,可以求出解析解,整个程序的源代码如图11。源代码中用syms命令来定义符号,变量涉及到符号运算时都需要先定义;simplify是化简命令,可以自动化简、合并表达式,例如simplify((cos(y))^2+(sin(y))^2)会自动化简为1;disp中的char表示字符串。解的结果见图12。
图11
图12
对比一下,可以看出Matlab符号推导得到的前3个解与传统方法得到的式(6)相同。
如果关心$A$矩阵的逆是怎样的形式,可以单独运行inv(A),图13为屏幕截图。
图13
为了验证矩阵求逆是否正确,可以查看inv(A)$^*$ A的结果,的确显示为单位矩阵。
3 小结
理论力学教学应注意基本概念、基本思路和基本方法,而具体繁琐的计算工作可以交给数学软件,这样可以让学生掌握最一般的静力学分析、计算方法。
本篇介绍了Matlab求解静力学问题的方法,核心的函数是inv(矩阵求逆),其他相关的函数包括syms(定义符号)、simplify(符号化简)和disp(在屏幕上显示)。利用这些函数,可以完成静力学中代数方程的数值求解及带参数的符号推导。
数值计算看起来输入的工作量较大,却是一种通用的方法,关键的是可以获得系统所有部件的受力(包含数值解和符号解),为后续进一步分析打下了基础(如强度分析、结构优化、失稳等等),也直接为后续的结构力学打下了基础。
传统的建模方法,都需要针对具体问题,按一定的步骤推导才能得到静力学或动力学方程。每遇到一个新的问题,由于系统的结构不一样,要按相同的步骤重复一遍。是否有一种方法可建立一个适合于任意系统的一般公式,只要把系统的最基本的一些参数,如刚体数目、连接类型、连接点位置、受力情况等带入公式,就可以展开得到系统的动力学方程?本质上就是系统的$A$和$B$矩阵如何生产,能否自动生成?从图论的角度引入通路矩阵 和联通矩阵后,可以自动获得系统的$A$和$B$矩阵[4],而这正是一些商业力学软件(例如Adams)的处理思路。也就是说,Matlab处理力学问题,为学生打 开了一个通向处理实际复杂问题的窗口。
可以想象,如果学生掌握了数值求解和符号推导,只要是静力学能处理的问题,都可以很快获得全部的解答,也许未来的静力学题目可能就要换一种问法了,例如：已知两岸的距离为100 m,如果要架一座桁架桥,要求最大承重是$G$,每根杆的最大受力为$S$,单位重量的杆件价格为$p$,市场有如下几种尺寸的杆件可供选择$\cdots$。提出你的桥梁设计方案,如何在满足约束条件下成本最低?这样的题目既和实际接近,又把传统的解题变为设计和优化问题了,而这正是目前传统力学教育所缺乏的。
参考文献
高云峰. 理论力学辅导与习题集. 北京: 清华大学出版社, 2003
[本文引用: 3]
甘勤涛, 程政田, 胡仁喜. Matlab数学计算与工程分析——从入门到精通. 北京: 机械工业出版社, 2019
[本文引用: 3]
李俊峰, 张雄. 理论力学. 北京: 清华大学出版社, 2010
[本文引用: 3]
贾书惠. 刚体动力学. 北京: 高等教育出版社, 1987
[本文引用: 1]
3
2003
... (1)在静力学中,以往对于复杂系统的受力分析通常要适当取分离体,有时需要高度的技巧[1];同时由于传统计算能力的限制,往往只要求解出某些部件的受力;如果采用Matlab处理,可以采用统一的处理方式,把系统全部拆开,快速求出所有部件的受力,对系统的整体和各部件受力有更全面的了解. ...
... 在Matlab中运行X$=$inv(A)*B,就得到所有杆件以及$A$和$B$铰链处的力,具体桁架问题求解程序源代码见图6.源代码中clc是清除屏幕;Matlab中表示矩阵很方便,例如[1, 2, 3]是$3\times1$的行阵,而[1; 2; 3]是$1\times 3$的列阵;zeros(12, 12)是生成$12\times 12$的零矩阵,里面元素全是0;$A(i,j)$表示$A$矩阵中第$i$行第$j$列的元素;在屏幕上显示的格式为disp(可以显示特定的文字),所以用num2str命令把具体数值转换为符号.如果想 更简单些,X=inv(A)*B后面不写分号";"就能直接显示结果(如显示为3.4509,而不是$S_1=-3.4509$ N). ...
... [1; 2; 3]是$1\times 3$的列阵;zeros(12, 12)是生成$12\times 12$的零矩阵,里面元素全是0;$A(i,j)$表示$A$矩阵中第$i$行第$j$列的元素;在屏幕上显示的格式为disp(可以显示特定的文字),所以用num2str命令把具体数值转换为符号.如果想 更简单些,X=inv(A)*B后面不写分号";"就能直接显示结果(如显示为3.4509,而不是$S_1=-3.4509$ N). ...
3
2019
... 其中inv是Matlab中矩阵求逆的函数[2],运行后就能直接解出系统中所有的未知力. ...
... 在Matlab中运行X$=$inv(A)*B,就得到所有杆件以及$A$和$B$铰链处的力,具体桁架问题求解程序源代码见图6.源代码中clc是清除屏幕;Matlab中表示矩阵很方便,例如[1, 2, 3]是$3\times1$的行阵,而[1; 2; 3]是$1\times 3$的列阵;zeros(12, 12)是生成$12\times 12$的零矩阵,里面元素全是0;$A(i,j)$表示$A$矩阵中第$i$行第$j$列的元素;在屏幕上显示的格式为disp(可以显示特定的文字),所以用num2str命令把具体数值转换为符号.如果想 更简单些,X=inv(A)*B后面不写分号";"就能直接显示结果(如显示为3.4509,而不是$S_1=-3.4509$ N). ...
... ; 2; 3]是$1\times 3$的列阵;zeros(12, 12)是生成$12\times 12$的零矩阵,里面元素全是0;$A(i,j)$表示$A$矩阵中第$i$行第$j$列的元素;在屏幕上显示的格式为disp(可以显示特定的文字),所以用num2str命令把具体数值转换为符号.如果想 更简单些,X=inv(A)*B后面不写分号";"就能直接显示结果(如显示为3.4509,而不是$S_1=-3.4509$ N). ...
3
2010
... 案例1：图示桁架系统中(图1),$ABC$是正三角形,边长为1 m,$DEF$也是正三角形,且$\angle ACD=\angle BAE=\angle CBF=15^{\circ}$,水平力$P=10$ N,垂直力$Q=20$ N,求1,2,3杆的内力[3]. ...
... 在Matlab中运行X$=$inv(A)*B,就得到所有杆件以及$A$和$B$铰链处的力,具体桁架问题求解程序源代码见图6.源代码中clc是清除屏幕;Matlab中表示矩阵很方便,例如[1, 2, 3]是$3\times1$的行阵,而[1; 2; 3]是$1\times 3$的列阵;zeros(12, 12)是生成$12\times 12$的零矩阵,里面元素全是0;$A(i,j)$表示$A$矩阵中第$i$行第$j$列的元素;在屏幕上显示的格式为disp(可以显示特定的文字),所以用num2str命令把具体数值转换为符号.如果想 更简单些,X=inv(A)*B后面不写分号";"就能直接显示结果(如显示为3.4509,而不是$S_1=-3.4509$ N). ...
... ; 3]是$1\times 3$的列阵;zeros(12, 12)是生成$12\times 12$的零矩阵,里面元素全是0;$A(i,j)$表示$A$矩阵中第$i$行第$j$列的元素;在屏幕上显示的格式为disp(可以显示特定的文字),所以用num2str命令把具体数值转换为符号.如果想 更简单些,X=inv(A)*B后面不写分号";"就能直接显示结果(如显示为3.4509,而不是$S_1=-3.4509$ N). ...
1
1987
... 传统的建模方法,都需要针对具体问题,按一定的步骤推导才能得到静力学或动力学方程.每遇到一个新的问题,由于系统的结构不一样,要按相同的步骤重复一遍.是否有一种方法可建立一个适合于任意系统的一般公式,只要把系统的最基本的一些参数,如刚体数目、连接类型、连接点位置、受力情况等带入公式,就可以展开得到系统的动力学方程?本质上就是系统的$A$和$B$矩阵如何生产,能否自动生成?从图论的角度引入通路矩阵 和联通矩阵后,可以自动获得系统的$A$和$B$矩阵[4],而这正是一些商业力学软件(例如Adams)的处理思路.也就是说,Matlab处理力学问题,为学生打 开了一个通向处理实际复杂问题的窗口. ...
