《理论力学》课程教学资源（文献资料）Matlab求解理论力学问...一）刚体系统及桁架受力问题_高云峰

力学与实践第43卷第2期2021年4月教育研地Matlab求解理论力学问题系列（一）刚体系统及桁架受力问题高云峰1)（清华大学航天航空学院，北京100084）专题：如果在理论力学教学中引入Matlab，根据经验，只需要三次课，就可以让学生掌握代数方程和微分(1)静力学专题1篇：刚体系统及桁架的受力问方程的数值求解、符号推导、动画演示等，让学生对题（看重介绍Matlab中代数方程的数值求解和符号理论力学问题的理解有飞跃式的提升；而教学中某求解)：些解题技巧性的内容则可以压缩，保持总学时不变。(2）运动学专题1篇：典型机构的运动分析（着具体来说：重介绍Matlab中非线性方程组的求解、动画显示，(1)在静力学中，以往对于复杂系统的受力分析如何对运动方程求导数)；通常要适当取分离体，有时需要高度的技巧[1]；同时(3)动力学专题2篇：单摆和椭圆摆的运动和周由于传统计算能力的限制，往往只要求解出某些部期（着重介绍Matlab中微分方程的数值求解、计算件的受力：如果采用Matlab处理，可以采用统一的可靠性、根据数据的快速傅里叶变换求周期）、乒乓处理方式，把系统全部拆开，快速求出所有部件的受球滚动问题（着重介绍Matlab中分段积分的处理方法，以及与分段对应的积分中断点问题)：力，对系统的整体和各部件受力有更全面的了解。(②)在运动学中，以往分析系统运动时，强调求(4)综合运用专题1篇：数据转换问题（着重介特定时刻或特定位置某点或刚体的速度和加速度，而绍在不同坐标系中看到结果，包括运动和动力学问系统的整体运动特点、某些点的运动轨迹有时难以想题)。象；而采用Matlab处理，可以求出系统任意点或刚通过这几篇文章，可以让学生们了解、熟悉Mat-体在任意时刻的速度和加速度等运动量，特别是其lab，大大提高解决问题的能力。画图和动画演示功能，可以快速直观地显示系统的下面首先从静力学开始。根据教学经验，学生整个运动过程、给出任意点的运动轨迹。在静力学中对于刚体系统和桁架的受力分析感到相(3)在动力学中，以往绝大部分问题都只能列写对比较困难，通常要适当拆开，否则解不出来。而动力学方程，通常没有解析解，传统数学分析的方法Matlab可以采用统一的方法求解，降低了解题的技也用不上，系统丰富复杂的动力学现象很难从方程巧，但是得到的解答更全面。中看出；而采用Matlab处理，可以获得系统整个运1Matlab中代数方程的求解动过程中的受力、速度和加速度等量，还可以快速直静力学问题的求解一般可以化为代数方程的求观地演示系统的运动过程。考虑到目前理论力学教学中对于数值计算、符解，代数方程一般可以写为号推导很少介绍，为此专门准备系列理论力学教学AX=B(1)文章，每篇介绍1~2个典型的理论力学问题及如何其中A是n×n阶的矩阵，由系统的位置、尺寸等利用Matlab进行处理。系列文章具体计划分为如下本文于2020-06-01收到。1) E-mail: gaoyunfeng@tsinghua.edu.cn引用格式：高云峰.Matlab求解理论力学问题系列（一）刚体系统及桁架受力问题.力学与实践，2021，43(2)：256-26）Gao Yunfeng.Force analysis of rigid body system and trusses.Mechanics in Engineering,2021, 43(2): 256-26(C)1994-2023ChinaAcademicJournalElectronicPublishingHouse.Allrightsreserved.http://www.cnki.net

第 43 卷 第 2 期 力 学 与 实 践 2021 年 4 月 Matlab 求解理论力学问题系列 (一) 刚体系统及桁架受力问题 高云峰 1) (清华大学航天航空学院, 北京 100084) 如果在理论力学教学中引入 Matlab，根据经验， 只需要三次课，就可以让学生掌握代数方程和微分 方程的数值求解、符号推导、动画演示等，让学生对 理论力学问题的理解有飞跃式的提升；而教学中某 些解题技巧性的内容则可以压缩，保持总学时不变。 具体来说： (1) 在静力学中，以往对于复杂系统的受力分析 通常要适当取分离体，有时需要高度的技巧[1]；同时 由于传统计算能力的限制，往往只要求解出某些部 件的受力；如果采用 Matlab 处理，可以采用统一的 处理方式，把系统全部拆开，快速求出所有部件的受 力，对系统的整体和各部件受力有更全面的了解。 (2) 在运动学中，以往分析系统运动时，强调求 特定时刻或特定位置某点或刚体的速度和加速度，而 系统的整体运动特点、某些点的运动轨迹有时难以想 象；而采用 Matlab 处理，可以求出系统任意点或刚 体在任意时刻的速度和加速度等运动量，特别是其 画图和动画演示功能，可以快速直观地显示系统的 整个运动过程、给出任意点的运动轨迹。 (3) 在动力学中，以往绝大部分问题都只能列写 动力学方程，通常没有解析解，传统数学分析的方法 也用不上，系统丰富复杂的动力学现象很难从方程 中看出；而采用 Matlab 处理，可以获得系统整个运 动过程中的受力、速度和加速度等量，还可以快速直 观地演示系统的运动过程。 考虑到目前理论力学教学中对于数值计算、符 号推导很少介绍，为此专门准备系列理论力学教学 文章，每篇介绍 1 ∼ 2 个典型的理论力学问题及如何 利用 Matlab 进行处理。系列文章具体计划分为如下 专题： (1) 静力学专题 1 篇：刚体系统及桁架的受力问 题 (着重介绍 Matlab 中代数方程的数值求解和符号 求解)； (2) 运动学专题 1 篇：典型机构的运动分析 (着 重介绍 Matlab 中非线性方程组的求解、动画显示， 如何对运动方程求导数)； (3) 动力学专题 2 篇：单摆和椭圆摆的运动和周 期 (着重介绍 Matlab 中微分方程的数值求解、计算 可靠性、根据数据的快速傅里叶变换求周期)、乒乓 球滚动问题 (着重介绍 Matlab 中分段积分的处理方 法，以及与分段对应的积分中断点问题)； (4) 综合运用专题 1 篇：数据转换问题 (着重介 绍在不同坐标系中看到结果，包括运动和动力学问 题)。 通过这几篇文章，可以让学生们了解、熟悉 Mat￾lab，大大提高解决问题的能力。 下面首先从静力学开始。根据教学经验，学生 在静力学中对于刚体系统和桁架的受力分析感到相 对比较困难，通常要适当拆开，否则解不出来。而 Matlab 可以采用统一的方法求解，降低了解题的技 巧，但是得到的解答更全面。 1 Matlab 中代数方程的求解 静力学问题的求解一般可以化为代数方程的求 解，代数方程一般可以写为 AX = B (1) 其中 A 是 n × n 阶的矩阵，由系统的位置、尺寸等 本文于 2020–06–01 收到。 1) E-mail: gaoyunfeng@tsinghua.edu.cn 引用格式: 高云峰. Matlab 求解理论力学问题系列 (一) 刚体系统及桁架受力问题. 力学与实践, 2021, 43(2): 256-261 Gao Yunfeng. Force analysis of rigid body system and trusses. Mechanics in Engineering, 2021, 43(2): 256-261

第2期257高云峰：Matlab求解理论力学问题系列（一）刚体系统及架受力问题1Ss参数构成，X是n×1阶的列阵，由系统中待求解的未知数构成，B是n×1阶的列阵，由系统中已知载荷、尺寸等参数构成。具体内容和形式见下面案例中的具体表达式。列写力和力矩的平衡方程是理论力学教学中的-YA重点。一旦有了平衡方程就可以获得式（1)中的ASC矩阵和B列阵，而Matlab处理矩阵运算特别方便，其求解格式为+Q图2特殊截面法X = inv(A)*B(2)其中inv是Matlab中矩阵求逆的函数[2]，运行后就能直接解出系统中所有的未知力。案例1：图示桁架系统中（图1)，ABC是正三角形，边长为1m，DEF也是正三角形，且LACD=/BAE=ZCBF=15°,水平力P=10N,E垂直力Q=20N，求1，2，3杆的内力[3]。7+QCF图3所有杆件编号Y.1S,D1>S-+STX.AEAB图4A节点受力图+QS:图1桁架系统S,从理论力学教学的角度，希望学生采用特殊截1面法，把1，2，3杆截断，把三角形DEF“挖出来”（图2)，把DEF看作刚体，三个未知数正好可以求+Q解（求解过程略）。但是这个特殊截面对很多同学而图5E节点受力图言有一定的难度，不容易想到。而且每个桁架问题都可能有特殊性，求解时需要经验和技巧。对每个节点，根据水平和竖直方向力的平衡方而采用节点法就不需要什么技巧：将所有的杆程，分别有件都编号（图3)，全部拆开，设杆件受拉为正，对各S1cos15°+S?+Sg:cos60°+XA=0(3)节点列写平衡方程（为节省篇幅只以外部A节点和S1·sin15°+Sg·sin60°+YA=0内部E节点为例，见图4和图5)。(C)1994-2023ChinaAcademicJournalElectronicPublishingHouse.All rightsreserved.http://www.cnki.net

第 2 期 高云峰：Matlab 求解理论力学问题系列 (一) 刚体系统及桁架受力问题 257 参数构成，X 是 n × 1 阶的列阵，由系统中待求解 的未知数构成，B 是 n × 1 阶的列阵，由系统中已知 载荷、尺寸等参数构成。具体内容和形式见下面案例 中的具体表达式。 列写力和力矩的平衡方程是理论力学教学中的 重点。一旦有了平衡方程就可以获得式 (1) 中的 A 矩阵和 B 列阵，而 Matlab 处理矩阵运算特别方便， 其求解格式为 X = inv(A)∗B (2) 其中 inv 是 Matlab 中矩阵求逆的函数[2]，运行后就 能直接解出系统中所有的未知力。 案例 1：图示桁架系统中 (图 1)，ABC 是正 三角形，边长为 1 m，DEF 也是正三角形，且 ∠ACD = ∠BAE = ∠CBF = 15◦，水平力 P = 10 N， 垂直力 Q = 20 N，求 1，2，3 杆的内力[3]。 C A B E F D    Q P 图 1 桁架系统 从理论力学教学的角度，希望学生采用特殊截 面法，把 1，2，3 杆截断，把三角形 DEF“挖出来” (图 2)，把 DEF 看作刚体，三个未知数正好可以求 解 (求解过程略)。但是这个特殊截面对很多同学而 言有一定的难度，不容易想到。而且每个桁架问题都 可能有特殊性，求解时需要经验和技巧。 而采用节点法就不需要什么技巧：将所有的杆 件都编号 (图 3)，全部拆开，设杆件受拉为正，对各 节点列写平衡方程 (为节省篇幅只以外部 A 节点和 内部 E 节点为例，见图 4 和图 5)。 E F D    S3 S1 S2 Q 图 2 特殊截面法 E F D A B C    Q       图 3 所有杆件编号 A  XA YA S S S   图 4 A 节点受力图 E  S S S   Q 图 5 E 节点受力图 对每个节点，根据水平和竖直方向力的平衡方 程，分别有 S1 · cos 15◦ + S7 + S9 · cos 60◦ + XA = 0 S1 · sin 15◦ + S9 · sin 60◦ + YA = 0 } (3)

258力学与实践2021年第43卷其他节点的平衡方程类似，最后合在一起，写为-S1·cos15°-S4·cos60°+S5·cos60°=0(4)-Si·sin15°+S4sin60°+Ss-sin60°-Q=0JAX=B的形式，有00cos15°0010cos6000100S100sin15000000sin 6000S200010COS45°0000COs6000000S30-1000000sin 45°00sin60°001S400000cos 60°-P-COs75°0-cos60°000S5000 sin 75°000sin60°sin60°S60000(5)=00COs75°cos60°010000000S700sin 75° sin 60°000000000S80-cos15°0-cos60°COs6000000000Sg0sin15000sin60°00Qsin60°00000XA000cOS 60°00000YA0COS 45°-10000-sin60°000000LYB]0 sin 45°0L在Matlab中运行X=inv(A)*B，就得到所有杆i行第j列的元素；在屏幕上显示的格式为disp(可件以及A和B铰链处的力，具体桁架问题求解程序以显示特定的文字），所以用num2str命令把具体数源代码见图6。源代码中clc是清除屏幕；Matlab中值转换为符号。如果想更简单些，X=inv(A）*B后面表示矩阵很方便，例如[1，2.3]是3×1的行阵，而不写分号“：”就能直接显示结果（如显示为3.4509，[1;2;3]是1×3的列阵：zeros(12,12是生成12×12而不是Si=-3.4509N)。的零矩阵，里面元素全是0；A(ii)表示A矩阵中第其解答的截图见图7。%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%桁架系统受力分析+符号推导%%969696969696969696969696969696969696969696969696969696969696969696clc;%清除屏幕%开始工作P=10:Q=20:%载荷赋值%后面用到的角度theta15=15*pi/180;theta45=45*pi/180;theta60=60*pi/180:theta75=75*pi/180:%角度转换为弧度形成B列阵B=[0;0;0;0;-P:0;0;0:0;Q:0;0];%B列阵形成A矩阵A=zeros(12,12);%A矩阵先清零A(1,1)=cos(theta15);A(1,7)=1;A(1,9)=cos(theta60);A(1,10)=1;%A中元素%省略A(12,2)=-sin(theta45);A(12,5)=-sin(theta60):%A中元素一解代数方程%X=inv(A)*B:%一在屏幕上显示disp（[s1=，num2str（x（1)），N"）%显示s1的结果%省略不符结果disp（YB=，num2str（X（12)），NJ）%显示YB的结果图6桁架问题源代码(C)1994-2023ChinaAcademic JournalElectronicPublishingHouse.All rightsreserved.http://www.cnki.net

258 力 学 与 实 践 2021 年 第 43 卷 −S1 ·cos 15◦ − S4 ·cos 60◦ + S5 ·cos 60◦ = 0 −S1 ·sin 15◦ + S4 ·sin 60◦ + S5 ·sin 60◦ − Q = 0 } (4) 其他节点的平衡方程类似， 最后合在一起， 写为 AX = B 的形式，有                            cos 15◦ 0 0 0 0 0 1 0 cos 60◦ 1 0 0 sin 15◦ 0 0 0 0 0 0 0 sin 60◦ 0 1 0 0 − cos 45◦ 0 0 0 0 −1 − cos 60◦ 0 0 0 0 0 sin 45◦ 0 0 0 0 0 sin 60◦ 0 0 0 1 0 0 − cos 75◦ 0 0 0 0 cos 60◦ − cos 60◦ 0 0 0 0 0 − sin 75◦ 0 0 0 0 − sin 60◦ − sin 60◦ 0 0 0 0 0 cos 75◦ cos 60◦ 0 1 0 0 0 0 0 0 0 0 sin 75◦ − sin 60◦ 0 0 0 0 0 0 0 0 − cos 15◦ 0 0 − cos 60◦ cos 60◦ 0 0 0 0 0 0 0 − sin 15◦ 0 0 sin 60◦ sin 60◦ 0 0 0 0 0 0 0 0 cos 45◦ 0 0 − cos 60◦ −1 0 0 0 0 0 0 0 − sin 45◦ 0 0 − sin 60◦ 0 0 0 0 0 0 0                                                       S1 S2 S3 S4 S5 S6 S7 S8 S9 XA YA YB                            =                            0 0 0 0 −P 0 0 0 0 Q 0 0                            (5) 在 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 桁架问题源代码

第2期259高云峰：Matlab求解理论力学问题系列（一）刚体系统及桁架受力问题S1=-3.4509N如果采用Matlab处理，则全部拆开，对节点和S2=-9.4281N刚体分别列写平衡方程，然后获得A和B矩阵。为S3=12.879N节省篇幅只画出D节点和AC杆件的受力图，见S4=14.3647N图9和图10。S5=7.698NS6=-10.5157NStSaS7=13.5912NS8=13.849NS9=-0.51567NSaXA=-10N?YA=1.3397NDYB=18.6603N图9D节点受力图图7全部解答的截图因此使用Matlab求解静力学问题，关键是确定YAYA矩阵和B列阵，而这与列写平衡方程有关。2Matlab中带参数代数方程的求解Matlab功能强大，除了可以进行数值计算，还XeXAA可以进行符号推导。因此，如果某些静力学问题没有C具体数值，也可以进行求解。案例2：横梁桁架结构由横梁AC和BC及五c根细支撑杆组成，所受载荷及尺寸如图8所示。求S.1，2，3杆的内力。图10AC杆件受力图P0对D节点列写水平和竖直方向力的平衡方程，有A?V2S1 + S3 = 02(7)VSi + S2= 02D对AC部件列写水平和竖直方向力的平衡方程，再对A点取矩，有图8刚体系统V2Si+XA+Xc=0从理论力学教学的角度，希望学生适当地取分2V2离体，但是有一定的技巧，解出的答案是（具体分析(8)Si-S2+YA+Yc-P-qa=02过程略)3-S2a+2Yca-Pa-qa·=a= 022S1:(P +3qa)2其他杆件和节点的平衡方程类似，最后合在一(6)S2:(P + 3qa)起，把未知数放在方程一侧，把已知载荷有关的量放S3=(P + 3qa)在方程另一侧，写为AX=B的形式，有(C)1994-2023ChinaAcademicJournalElectronicPublishingHouse.All rightsreserved.http://www.cnki.net

第 2 期 高云峰：Matlab 求解理论力学问题系列 (一) 刚体系统及桁架受力问题 259 图 7 全部解答的截图 因此使用 Matlab 求解静力学问题，关键是确定 A 矩阵和 B 列阵，而这与列写平衡方程有关。 2 Matlab 中带参数代数方程的求解 Matlab 功能强大，除了可以进行数值计算，还 可以进行符号推导。因此，如果某些静力学问题没有 具体数值，也可以进行求解。 案例 2：横梁桁架结构由横梁 AC 和 BC 及五 根细支撑杆组成，所受载荷及尺寸如图 8 所示。求 1，2，3 杆的内力。 A P q D C B E a a a a a      图 8 刚体系统 从理论力学教学的角度，希望学生适当地取分 离体，但是有一定的技巧，解出的答案是 (具体分析 过程略) S1 = √ 2 2 (P + 3qa) S2 = − 1 2 (P + 3qa) S3 = 1 2 (P + 3qa)    (6) 如果采用 Matlab 处理，则全部拆开，对节点和 刚体分别列写平衡方程，然后获得 A 和 B 矩阵。为 节省篇幅只画出 D 节点和 AC 杆件的受力图，见 图 9 和图 10。 D S1 S2 S3 图 9 D 节点受力图 A C   P q YA YC XA XC S1 S2 图 10 AC 杆件受力图 对 D 节点列写水平和竖直方向力的平衡方程，有 − √ 2 2 S1 + S3 = 0 √ 2 2 S1 + S2 = 0    (7) 对 AC 部件列写水平和竖直方向力的平衡方程，再 对 A 点取矩，有 √ 2 2 S1 + XA + XC = 0 − √ 2 2 S1 − S2 + YA + YC − P − qa = 0 −S2a + 2YC a − P a − qa · 3 2 a = 0    (8) 其他杆件和节点的平衡方程类似，最后合在一 起，把未知数放在方程一侧，把已知载荷有关的量放 在方程另一侧，写为 AX = B 的形式，有

260力学与实践2021年第43卷V201000000020V200100000020V2S10000-1000002S2V2S300001000002S4V20000010010Ss2(9)P + qaXAV2-100001001YA32Pa + qa200000002aYB-a00XcV200000000-1Yc2qaV230000-1001-1290220000002a00a利用X=inv(A)*B，可以求出解析解，整个程序S1=(2°（1/2)*(P+3*a*q))/2S2=-P/2-（3*a*q）/2的源代码如图11。源代码中用syms命令来定义符S3=P/2+（3*a*q)/2号，变量涉及到符号运算时都需要先定义；simplifyS4=-P/2（3*a*q)/2是化简命令，可以自动化简、合并表达式，例如S5=(2°（1/2)*(P+3*a*q))/2simplify(cos(y)~2+(sin(y））2)会自动化简为1;dispXA=0中的char表示字符串。解的结果见图12。YA=(3*P)/4 + a*qYB=P/4+a*q%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%XC=-P/2-（3*a*g)/2%%刚体系统受力分析+符号推导：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%YC=P/4clc:%清除屏幕：开始工作%图12解答表达式截图symsPqaABX%定义符号一求B列阵5B=[0:0;0;0:0;P+q*a;P*a+3*q*a*a/2:0q*a:-3*q*a*a/2];为了验证矩阵求逆是否正确，可以查看inv(A)*求A矩阵%A的结果，的确显示为单位矩阵。%A矩阵中的元素A(1,1)=-sqrt(2)/2;A(1,3)=1:%省略部分A中的元素%A矩阵中的元素A(10,4)=a;A(10,10)=2*a;3小结一求解方程X=simplify(iny(A)*B) :%X=inv(A)B，然后化筒理论力学教学应注意基本概念、基本思路和基一在屏幕上显示结果本方法，而具体繁琐的计算工作可以交给数学软件，disp（'S1=，char(x(1)）J）%显示S1的结果%省略部分结果这样可以让学生掌握最一般的静力学分析、计算方disp（["YC=，char（x（10)）J）%显示YC的结果法。图11求解带参数代数方程的源代码本篇介绍了Matlab求解静力学问题的方法，核心的函数是inv(矩阵求逆），其他相关的函数包括对比一下，可以看出Matlab符号推导得到的前syms(定义符号）、simplify(符号化简）和disp(在屏幕3个解与传统方法得到的式（6)相同。上显示)。利用这些函数，可以完成静力学中代数方如果关心A矩阵的逆是怎样的形式，可以单独运行inv(A)，图13为屏幕截图。程的数值求解及带参数的符号推导。(C)1994-2023 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net

260 力 学 与 实 践 2021 年 第 43 卷                                       − √ 2 2 0 1 0 0 0 0 0 0 0 √ 2 2 1 0 0 0 0 0 0 0 0 0 0 −1 0 √ 2 2 0 0 0 0 0 0 0 0 1 √ 2 2 0 0 0 0 0 √ 2 2 0 0 0 0 1 0 0 1 0 − √ 2 2 −1 0 0 0 0 1 0 0 1 0 −a 0 0 0 0 0 0 0 2a 0 0 0 0 − √ 2 2 0 0 0 −1 0 0 0 0 −1 − √ 2 2 0 0 1 0 −1 0 0 0 a 0 0 0 0 0 2a                                                             S1 S2 S3 S4 S5 XA YA YB XC YC                       =                                0 0 0 0 0 P + qa P a + 3 2 qa2 0 qa − 3 2 qa2                                (9) 利用 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(在屏幕 上显示)。利用这些函数，可以完成静力学中代数方 程的数值求解及带参数的符号推导

第2期高云峰：Matlab求解理论力学问题系列（一）刚体系统及桁架受力问题261[-2°(1/2)/2,2°(1/2)/2,-2(1/2)/2,2°(1/2)/2,0,0,2°(1/2)/(2*a),0,0,-2(1/2)/(2*a))[1/2,1/2,1/2,-1/2,0, 0,-1/(2*a),0, 0,1/(2*a))[1/2,1/2,-1/2,1/2,0,0,1/(2*a),0,0,-1/(2*a))[-1/2,-1/2,-1/2,1/2,0,0,-1/(2*a),0, 0,1/(2*a))[2°(1/2)/2,2°(1/2)/2,2(1/2)/2,2°(1/2)/2,0,0,2(1/2)/(2*a),0,0,-2(1/2)/(2*a)）[1,0,1,010, 1, 0,0,1, 0,[3/4,-1/4,-1/4,1/4,0,1,-1/(4*a),0,0,-1/(4*a))[1/4,1/4,1/4,3/4,0,0,1/(4*a),0, 1,1/(4*a))T-1/2,-1/2,-1/2,1/2,0,0,-1/(2*a),-1,0,1/(2*a))[1/4,1/4,1/4,-1/4,0,0,1/(4*a)，0,0,1/(4*a))图13带参数的A矩阵的逆数值计算看起来输入的工作量较大，却是一种可以想象，如果学生掌握了数值求解和符号推导，只要是静力学能处理的问题，都可以很快获得全通用的方法，关键的是可以获得系统所有部件的受力（包含数值解和符号解），为后续进一步分析打下部的解答，也许未来的静力学题目可能就要换一种了基础（如强度分析、结构优化、失稳等等），也直接问法了，例如：已知两岸的距离为100m，如果要架为后续的结构力学打下了基础。一座桁架桥，要求最大承重是G，每根杆的最大受力传统的建模方法，都需要针对具体问题，按一定为S，单位重量的杆件价格为p，市场有如下几种尺的步骤推导才能得到静力学或动力学方程。每遇到寸的杆件可供选择。提出你的桥梁设计方案，如一个新的问题，由于系统的结构不一样，要按相同的何在满足约束条件下成本最低？这样的题目既和实步骤重复一遍。是否有一种方法可建立一个适合于际接近，又把传统的解题变为设计和优化问题了，而任意系统的一般公式，只要把系统的最基本的一些这正是目前传统力学教育所缺乏的。参数，如刚体数目、连接类型、连接点位置、受力情况等带入公式，就可以展开得到系统的动力学方程？参考文献本质上就是系统的A和B矩阵如何生产，能否自动生成？从图论的角度引入通路矩阵和联通矩阵后，1高云峰.理论力学辅导与习题集。北京：清华大学出版社，20032甘勤涛，程政田，胡仁喜Matlab数学计算与工程分析——从可以自动获得系统的A和B矩阵[4]，而这正是一入门到精通.北京：机械工业出版社，2019些商业力学软件（例如Adams）的处理思路。也就是3李俊峰，张雄。理论力学。北京：清华大学出版社，2010说，Matlab处理力学问题，为学生打开了一个通向4贾书惠.刚体动力学.北京：高等教育出版社，1987处理实际复杂问题的窗口。(贵任编辑：胡漫)(C)1994-2023ChinaAcademic JournalElectronicPublishingHouse.Allrightsreserved.http://www.cnki.net

第 2 期 高云峰：Matlab 求解理论力学问题系列 (一) 刚体系统及桁架受力问题 261 图 13 带参数的 A 矩阵的逆 数值计算看起来输入的工作量较大，却是一种 通用的方法，关键的是可以获得系统所有部件的受 力 (包含数值解和符号解)，为后续进一步分析打下 了基础 (如强度分析、结构优化、失稳等等)，也直接 为后续的结构力学打下了基础。 传统的建模方法，都需要针对具体问题，按一定 的步骤推导才能得到静力学或动力学方程。每遇到 一个新的问题，由于系统的结构不一样，要按相同的 步骤重复一遍。是否有一种方法可建立一个适合于 任意系统的一般公式，只要把系统的最基本的一些 参数，如刚体数目、连接类型、连接点位置、受力情 况等带入公式，就可以展开得到系统的动力学方程？ 本质上就是系统的 A 和 B 矩阵如何生产，能否自 动生成？从图论的角度引入通路矩阵和联通矩阵后， 可以自动获得系统的 A 和 B 矩阵 [4]，而这正是一 些商业力学软件 (例如 Adams) 的处理思路。也就是 说，Matlab 处理力学问题，为学生打开了一个通向 处理实际复杂问题的窗口。 可以想象，如果学生掌握了数值求解和符号推 导，只要是静力学能处理的问题，都可以很快获得全 部的解答，也许未来的静力学题目可能就要换一种 问法了，例如：已知两岸的距离为 100 m，如果要架 一座桁架桥，要求最大承重是 G，每根杆的最大受力 为 S，单位重量的杆件价格为 p，市场有如下几种尺 寸的杆件可供选择 · · ·。提出你的桥梁设计方案，如 何在满足约束条件下成本最低？这样的题目既和实 际接近，又把传统的解题变为设计和优化问题了，而 这正是目前传统力学教育所缺乏的。 参 考 文 献 1 高云峰. 理论力学辅导与习题集. 北京: 清华大学出版社, 2003 2 甘勤涛, 程政田, 胡仁喜. Matlab 数学计算与工程分析 —— 从 入门到精通. 北京: 机械工业出版社, 2019 3 李俊峰, 张雄. 理论力学. 北京: 清华大学出版社, 2010 4 贾书惠. 刚体动力学. 北京: 高等教育出版社, 1987 (责任编辑: 胡 漫)