简介面向机械原理课程学习者与课程设计者的平面六杆机构运动分析大作业围绕角位移、角速度、角加速度以及E点位移、速度、加速度等运动变量展开适合机械类本科生完成课程作业、复盘数值求解思路时参考。压缩包共1个文件为单一doc文档约670KB内含题目说明与原始参数表L1至L6、α、xG、yG等尺寸、封闭图形建模与矩阵推导过程、程序流程以及用C语言编写的源程序清单涉及math.h、agaus.c、dnetn.c等数学与数值解算库。作业按原动件0到360度、37个计算点求解通过结构体存储各运动变量并将结果输出到num-output.txt读者可据此还原计算流程、对照公式与代码并绘制运动曲线图与轨迹曲线。目前已有268人学习下载。1. 平面六杆机构的运动分析为什么最后都得落到数值解法上37 个采样点、4 个从动件、位移速度加速度三类量一份机械原理大作业把解析建模、数值方法和 C 语言三件事捆成了一条链。麻烦在于给定第 1-B 组的杆长 L124.0、L2105.6、L265.0、L367.5、L487.5、L534.4、L625.0以及 60° 的装配夹角 α 之后从动件的角位移 θ2、θ3、θ5、θ6 并没有可以直接代入的显式表达式——两个封闭环把四个未知角耦合成一组非线性方程。想画出 θ2~θ6 随 θ1 变化的曲线绕不开牛顿迭代求位置、高斯消元求速度与加速度这条路径。这份文档适合正在做平面六杆机构运动分析的人对照复现也适合想复习解析建模加数值求解完整链路的人。2. 两份闭环方程的建立与牛顿迭代求角位移2.1 从封闭图形到四个非线性方程机构的运动学约束来自两个封闭环。第一个环绕 O2 出发L1 L2 L3 L4其中 O4 固定L4 是机架矢量。第二个环从 O2 经 B、C、D 一直连到 G 点L1 L2 L5 L6 AG这里的 AG 是把机架参考点 G 的绝对坐标XG153.5YG41.7也纳入闭合条件。把两个矢量环分别向 X、Y 轴投影就得到dnetnf里的四个残差方程自变量 x[0]~x[3] 依次是 θ2、θ3、θ5、θ6/* 第一个闭环O2-B-C-O4 */ y[0] L[1]*cos(p-theta[0]) L[2]*cos(x[0]) - L[3]*cos(x[1]) - L[4]; y[1] L[1]*sin(p-theta[0]) L[2]*sin(x[0]) - L[3]*sin(x[1]); /* 第二个闭环C-D-E-GG 点坐标折算进右端 */ y[2] L[3]*cos(x[1]) L[0]*cos(x[0]-Alpha) L[5]*cos(x[2]) - L[6]*cos(x[3]) - XG L[4]; y[3] L[3]*sin(x[1]) L[0]*sin(x[0]-Alpha) L[5]*sin(x[2]) - L[6]*sin(x[3]) - YG;注意数组下标的对应关系L[0]65.0 是 L2L[2]105.6 才是 L2这一处最容易写错。L2 与 L2 之间存在 60° 的固定夹角所以它的方向角写成x[0]-Alpha而不是 x[0]。G 点坐标不参与求导作为常数项放在右端也因此出现-XGL[4]这种看起来别扭但正确的写法——这是把 G 点相对某个中间参考系的坐标换算回 O 点坐标系的结果。2.2 牛顿迭代的骨架与差商雅可比四个残差方程没有解析解标准做法是牛顿法。dnetn(4, eps, t, h, x, k)的调用约定是方程个数 4、收敛精度 eps1e-7、数值微分步长 t0.1 与 h0.1、解向量 x既是初值也是出口、最大迭代次数 k100。因为dnetnf只给了残差雅可比矩阵靠中心差商近似即可/* 用中心差商构造 4x4 雅可比J[i][j] dF_i / dx_j */ for (j 0; j n; j) { double xj x[j]; x[j] xj h; dnetnf(y1, x, n); /* F(x h*e_j) */ x[j] xj - h; dnetnf(y2, x, n); /* F(x - h*e_j) */ x[j] xj; /* 立即还原避免污染下一列 */ for (i 0; i n; i) J[i][j] (y1[i] - y2[i]) / (2.0 * h); }得到 J 之后解线性方程 J·Δx -F再令 x x Δx反复直到max|Δx| eps。参数 h 取 0.1 是个折中太大则差商误差 O(h²) 明显太小则 y1 与 y2 相减出现有效位抵消。对角度的弧度值来说0.1 意味着约 5.7°量级上够用。2.3 初值给多粗才收敛源程序里的初值是x[4] {26.23, 49.75, 87.16, 37.25}乘Angle转弧度。对照 θ10 时收敛结果 θ20.656023 rad≈37.6°、θ31.267191 rad≈72.6°初值偏差超过 10°照样一次收敛说明机构在这一位形附近的雅可比条件数还不错。变量初值度θ10 收敛值度偏差θ226.2337.5911.36θ349.7572.6022.85θ587.16132.3245.16θ637.25110.8373.58但够宽不等于随便给。牛顿法的收敛域依赖具体位形θ1 越接近死点同一个初值的迭代次数越多甚至可能跳到另一个装配分支。源码里每算完一个点都getchar()停一下正是为了盯住循环变量 i迭代次数有没有异常放大。3. 速度与加速度方程的 4×4 矩阵装配3.1 一阶求导得到速度方程把四个位置方程对时间求一次导得到的线性方程组写成 A·ω bA 就是位置方程对四个角的偏导矩阵ω 是待求的 ω2、ω3、ω5、ω6b 收集已知的 ω1 驱动项。装配顺序与dnetnf的方程顺序一一对应a[0][0] -L[2]*sin(p-theta[1]); /* dF0/dθ2 */ a[0][1] L[3]*sin(p-theta[2]); /* dF0/dθ3 */ a[0][2] 0.; /* θ5 不在第一个闭环里 */ a[0][3] 0.; /* θ6 同理 */ a[1][0] L[2]*cos(p-theta[1]); a[1][1] -L[3]*cos(p-theta[2]); a[2][0] -L[0]*sin(p-theta[1]-Alpha); /* L2别写成 L2 */ a[2][1] -L[3]*sin(p-theta[2]); a[2][2] -L[5]*sin(p-theta[3]); /* θ5 的方向角 */ a[2][3] L[6]*sin(p-theta[4]); /* θ6 的符号相反 */ a[3][0] L[0]*cos(p-theta[1]-Alpha); a[3][1] L[3]*cos(p-theta[2]); a[3][2] L[5]*cos(p-theta[3]); a[3][3] -L[6]*cos(p-theta[4]);右端向量只有前两行有值因为 ω1 是输入b[0] L[1]*sin(p-theta[0])*w1; b[1] -L[1]*cos(p-theta[0])*w1; b[2] 0.; b[3] 0.;符号规律很简单位置方程里某项是Lk*cos(θk)对 θk 求导就是-Lk*sin(θk)*ωk是-Lk*cos(θk)则求导后-Lk*sin的系数翻正。把 θ1 的导数移到右端时L1*sin(θ1)*ω1原样留在 b[0]。装配完直接agaus(a,b,4)即可。3.2 二阶求导得到加速度方程再对时间求一次导A 矩阵与速度方程完全一致位置方程的二阶偏导仍落在同一组偏导上变的只有右端 b。右端多出离心项 ω²以及驱动角加速度项——这里 ω1 恒定所以角加速度驱动项为 0b[0] L[2]*cos(p-theta[1])*w2*w2 - L[3]*cos(p-theta[2])*w3*w3 w1*w1*L[1]*cos(p-theta[0]); /* α10 时的驱动余项 */ b[1] L[2]*sin(p-theta[1])*w2*w2 - L[3]*sin(p-theta[2])*w3*w3 w1*w1*L[1]*sin(p-theta[0]); b[2] L[0]*cos(p-theta[1]-Alpha)*w2*w2 L[3]*cos(p-theta[2])*w3*w3 L[5]*cos(p-theta[3])*w5*w5 - L[6]*cos(p-theta[4])*w6*w6;每一项的来源都是把含 α 的项留在左边把含 ω² 的项挪到右边符号移动时跟着变号。写代码时最容易漏的是 w1*w1*L[1]*cos(theta[0])这一项——它来自 θ1 的二阶导数因为 ω1 是常数、α10只有 ω1² 这一半留下来。3.3 agaus 的调用约定与必须重装矩阵agaus(a, b, n)用高斯消元带主元选取解 n 元线性方程组返回值非 0 表示求解成功解写在 b 里。三个关键点参数含义注意an×n 系数矩阵会被就地LU分解破坏b右端向量出口变成解向量n阶数这里是 4因为 a 被就地破坏加速度那一步必须把 a 的 16 个元素重新赋值一遍源码里确实重装了一次完整的 a。还有一点if(agaus(a,b,4)!0)是判断成功的条件返回 0 说明主元过小、矩阵接近奇异此时 b 里的解不可信应该跳过不写入结构体而不是继续用。4. E 点轨迹与 C 程序结构4.1 motion 结构体与索引映射整个程序的数据载体是一个结构体数组struct motion mot[37]每个元素装一个采样点的全部结果struct motion { int theta1; /* 原动件角度单位度 */ double theta[5]; /* θ1,θ2,θ3,θ5,θ6 —— 注意跳过 θ4 */ double w[4]; /* ω2,ω3,ω5,ω6 */ double alpha[4]; /* α2,α3,α5,α6 */ double XYe[2]; /* E 点位置 */ double Ve[3]; /* E 点速度分量 合速度 */ double ae[3]; /* E 点加速度分量 合加速度 */ };数组下标与物理量的映射是这份代码里最容易读错的地方。θ4 是机架不参与求解所以theta[]里跳过了它w 和 alpha 都只存从动件下标 0~3 对应 θ2、θ3、θ5、θ6。写后处理脚本时如果按下标即编号理解theta[3] 会被误当成 θ4整条曲线就错位了。4.2 37 点主循环与结果落盘主循环从 θ10° 走到 360°步长 10°用 n36 保证端点闭合。每轮先调用牛顿迭代更新 x再把 x 写回结构体然后装配速度、加速度矩阵for (n 0, p mot; n 36; n, p) { p-theta1 n * 10; /* 度写文件用 */ p-theta[0] n * 10 * Angle; /* 弧度计算用 */ i dnetn(4, eps, t, h, x, k); /* x 被就地更新为本次解 */ for (m 0; m 4; m) p-theta[m1] x[m]; /* θ2,θ3,θ5,θ6 */ /* …装配 a、b 求 w再装配 a、b 求 alpha… */ fprintf(fp, %d\t, p-theta1); for (m 0; m 4; m) fprintf(fp, %lf\t, p-theta[m]); for (m 0; m 3; m) fprintf(fp, %lf\t, p-w[m]); for (m 0; m 3; m) fprintf(fp, %lf\t, p-alpha[m]); for (m 0; m 1; m) fprintf(fp, %lf\t, p-XYe[m]); for (m 0; m 2; m) fprintf(fp, %lf\t, p-Ve[m]); for (m 0; m 2; m) fprintf(fp, %lf\t, p-ae[m]); fprintf(fp, \n); }输出文件num-output.txt每行 22 列制表符分隔可以直接丢进 Excel 或 Python 画图。x 被dnetn就地更新这一点很重要下一轮循环的初值自动继承上一轮的解这比每点都用同一组粗初值稳健得多尤其在 θ1 接近 340°~350° 那段。4.3 E 点位移、速度、加速度的合成E 点挂在 CD 杆的延长端位置由 D 点经 L5 和 G 点经 L6 共同决定p-XYe[0] XG L[6]*cos(p-theta[4]) - L[5]*cos(p-theta[3]); p-XYe[1] YG L[6]*sin(p-theta[4]) - L[5]*sin(p-theta[3]); p-Ve[0] -L[6]*sin(p-theta[4])*p-w[3] L[5]*sin(p-theta[3])*p-w[2]; p-Ve[1] L[6]*cos(p-theta[4])*p-w[3] - L[5]*cos(p-theta[3])*p-w[2]; p-Ve[2] sqrt(p-Ve[0]*p-Ve[0] p-Ve[1]*p-Ve[1]);注意速度里的下标p-w[2]是 ω5对应 θ5p-w[3]是 ω6对应 θ6正是 4.1 节强调过的映射。加速度表达式多出 ω² 项p-ae[0] -L[6]*cos(p-theta[4])*w6*w6 - L[6]*sin(p-theta[4])*alpha6 L[5]*cos(p-theta[3])*w5*w5 L[5]*sin(p-theta[3])*alpha5; p-ae[1] -L[6]*sin(p-theta[4])*w6*w6 L[6]*cos(p-theta[4])*alpha6 L[5]*sin(p-theta[3])*w5*w5 - L[5]*cos(p-theta[3])*alpha5; p-ae[2] sqrt(p-ae[0]*p-ae[0] p-ae[1]*p-ae[1]);这一组公式用对位置表达式连续求导就能得到cos 一次导变 -sin·ω再导一次变 -cos·ω² - sin·α注意两项符号不要混。5. 收敛失败与曲线突变初值、步长与数值验证5.1 340°~350° 的加速度尖峰从哪来翻结果表会发现两处异常。θ1340° 时 α5-17.5054、α6-25.2675θ1350° 时 α5 突然翻成 17.5277、α6364.8278 量级的跳变速度表里 ω5 从 -3.03617 跳到 -4.53764 再到 -2.67285。这不是程序写错了而是机构进入了接近死点的位形雅可比 A 在主元上趋近于零agaus求出的解对输入误差极度敏感几何上的微小变化被放大成加速度尖峰。看到这种形态先检查dnetn返回的迭代次数 i 有没有骤增再确认agaus是否返回 0。如果 i 超过 30 或者agaus返回 0就不要把该点的结果当作有效解写进文件。5.2 用上一点解当前初值源码里的做法已经是最省事的稳健策略把 x 声明在循环外每轮dnetn就地更新下一轮直接当初值。如果换成每点都重新给{26.23,49.75,87.16,37.25}在 340° 附近几乎必然迭代失败或跳分支。改法只有一行double x[4] {26.23*Angle, 49.75*Angle, 87.16*Angle, 37.25*Angle}; /* 循环内不再重新赋值让 dnetn 覆写 x 即可 */代价是一旦某个点解错误差会顺着 θ1 的正方向一路传染。折中方案是每算 10 个点强制回一次粗初值或者在i 30时回退到上一轮的解重新迭代。5.3 三个可复现的验证动作第一闭环残差校验。把每点的 θ2、θ3、θ5、θ6 代回dnetnf看四个 y 是否都在 1e-6 以内超过就说明该点没收敛。第二数值微分对拍。用相邻两点的角位移差分除以角度步长再乘 ω1和解析求得的 ω2 比较# data.txt 为 num-output.txt第 1 列是 θ1度第 3 列是 θ2弧度 import numpy as np d np.loadtxt(num-output.txt) theta1, theta2 d[:, 0], d[:, 2] w1 1.0 d_theta1 np.deg2rad(10.0) w2_num np.gradient(theta2, d_theta1) * w1 # 中心差分 print(np.max(np.abs(w2_num[1:-1] - d[1:-1, 7]))) # 与解析 ω2 对比中心差分精度 O(h²)10° 步长下差值通常落在 1e-3 量级如果某点差到 1e-1 以上就说明该点迭代质量差。第三周期闭合检查。θ10° 与 360° 是同一位置把两行的四个从动角做差如果差值在 1e-9 量级上闭合说明整轮 37 个点的迭代没有跳装配分支如果差了 0.1 rad 以上回去看 5.2 节的初值继承链断在哪一步。本文还有配套的精品资源点击获取
