运动算例基础

摘要

运动算例是计算机辅助工程（CAE）和三维设计软件中用于模拟机械系统运动行为的关键工具。本文深入剖析了运动算例的三种核心类型——动画、基本运动与Motion分析，从原理、实现机制到应用场景进行了全面对比。通过详细的代码示例和工程实践，帮助工程师在机构设计、干涉检查、载荷分析等环节做出合理选择。无论你是机械设计新手还是资深工程师，本文都能为你提供从理论到实践的完整知识体系。

---

1. 引言

在机械产品设计流程中，运动仿真占据着至关重要的地位。传统的物理样机测试成本高昂、周期漫长，而运动算例技术允许工程师在数字环境中验证机构的运动逻辑、动力特性以及系统响应。然而，许多工程师在面对"动画"、"基本运动"和"Motion分析"这三个术语时往往感到困惑——它们看起来相似，实则服务于完全不同的工程需求。

想象一下这样一个场景：你需要设计一个四连杆机构。如果只是为了让客户看到机构的大致运动姿态，使用动画功能即可；但若需要精确计算连杆上的受力情况，则必须借助Motion分析。本文将从技术底层出发，系统性地解析这三种算例类型的本质差异，并提供可复用的代码示例，帮助你建立清晰的选型决策框架。

---

2. 运动算例技术全景

2.1 运动仿真的数学基础

在深入具体类型前，我们先建立必要的理论背景。运动算例本质上是求解多体系统动力学问题，其核心方程可表示为：

对于刚体系统，牛顿-欧拉方程：

M(q) * q̈ + C(q, q̇) + G(q) = τ + J^T * λ

其中：

• M(q)：质量矩阵
• q̈：广义加速度
• C(q, q̇)：科里奥利与离心力项
• G(q)：重力项
• τ：驱动力/力矩
• J：约束雅可比矩阵
• λ：拉格朗日乘子（约束反力）

三种算例类型本质上是对上述方程不同层级的简化求解：

算例类型	求解策略	是否考虑质量/惯性	是否计算力/力矩
动画	纯运动学插值	否	否
基本运动	简化的动力学	是（仅部分）	是（近似值）
Motion分析	完整动力学	是（完整）	是（精确值）

2.2 时间积分与数值稳定性

Motion分析通常采用隐式积分器（如HHT算法），而基本运动使用显式积分器。这导致了两者在处理刚性系统时的稳定性差异：

# 显式与隐式积分对比（简化示例）
import numpy as np

def explicit_euler(M, C, K, F, dt, q, v, steps):
    """显式欧拉积分 - 用于基本运动"""
    for i in range(steps):
        a = np.linalg.solve(M, F - C @ v - K @ q)
        q = q + v * dt
        v = v + a * dt
    return q, v

def implicit_newmark(M, C, K, F, dt, q, v, beta=0.25, gamma=0.5):
    """Newmark隐式积分 - 用于Motion分析"""
    # 预测步
    q_pred = q + v * dt + (0.5 - beta) * a * dt**2
    v_pred = v + (1 - gamma) * a * dt
    
    # 修正步（迭代求解）
    for _ in range(10):  # 通常需要Newton-Raphson迭代
        residual = M @ a + C @ v_pred + K @ q_pred - F
        J = M + gamma * dt * C + beta * dt**2 * K
        delta_a = np.linalg.solve(J, -residual)
        a += delta_a
        q = q_pred + beta * dt**2 * delta_a
        v = v_pred + gamma * dt * delta_a
    return q, v

---

3. 动画：视觉驱动的运动模拟

3.1 技术原理

动画是运动算例中最基础的层级，它不求解任何动力学方程，而是通过关键帧插值来驱动部件运动。其核心算法是：

class KeyframeAnimation:
    """基于关键帧的动画引擎"""
    def __init__(self):
        self.keyframes = {}  # {time: {component: transform_matrix}}
        
    def add_keyframe(self, time, component_states):
        """添加关键帧"""
        self.keyframes[time] = component_states
        
    def interpolate(self, t):
        """在时间t进行插值"""
        # 找到相邻关键帧
        times = sorted(self.keyframes.keys())
        t_prev = max([ti for ti in times if ti <= t])
        t_next = min([ti for ti in times if ti >= t])
        
        if t_prev == t_next:
            return self.keyframes[t_prev]
        
        # 线性插值因子
        alpha = (t - t_prev) / (t_next - t_prev)
        
        # 对每个组件进行插值
        result = {}
        for component in self.keyframes[t_prev]:
            # 使用SLERP进行旋转插值，线性插值进行平移
            prev_mat = self.keyframes[t_prev][component]
            next_mat = self.keyframes[t_next][component]
            result[component] = self.slerp_transform(prev_mat, next_mat, alpha)
        
        return result
    
    def slerp_transform(self, mat1, mat2, alpha):
        """球形线性插值（用于旋转矩阵）"""
        # 提取旋转部分
        R1 = mat1[:3, :3]
        R2 = mat2[:3, :3]
        
        # 转换为四元数
        q1 = self.rotation_matrix_to_quaternion(R1)
        q2 = self.rotation_matrix_to_quaternion(R2)
        
        # SLERP
        dot = np.dot(q1, q2)
        if dot < 0:
            q2 = -q2
            dot = -dot
        
        theta = np.arccos(np.clip(dot, -1, 1))
        if theta < 1e-6:
            q_interp = q1
        else:
            q_interp = (np.sin((1-alpha)*theta) * q1 + 
                       np.sin(alpha*theta) * q2) / np.sin(theta)
        
        # 组合平移
        t1 = mat1[:3, 3]
        t2 = mat2[:3, 3]
        t_interp = (1-alpha) * t1 + alpha * t2
        
        # 重新组装矩阵
        R_interp = self.quaternion_to_rotation_matrix(q_interp)
        result = np.eye(4)
        result[:3, :3] = R_interp
        result[:3, 3] = t_interp
        return result

3.2 应用场景与局限

适用场景：

• 产品展示动画（无需物理真实性）
• 机构运动姿态预览
• 装配体爆炸视图
• 早期概念设计沟通

典型局限：

• 无法检测干涉（除非手动设置）
• 不考虑惯性效应
• 速度/加速度无物理意义
• 不能计算受力

3.3 工程实践示例

在SolidWorks中创建动画的典型流程：

# 通过API创建动画（伪代码）
def create_simple_animation(sw_app):
    """在SolidWorks中创建旋转动画"""
    model = sw_app.ActiveDoc
    anim_mgr = model.GetAnimationManager()
    
    # 创建动画算例
    anim_study = anim_mgr.CreateStudy("旋转展示", swAnimationType_e.swAnimation)
    
    # 添加第一个关键帧（0秒）
    key1 = anim_study.AddKeyframe(0.0)
    key1.SetComponentTransform("Crank", np.eye(4))  # 初始位置
    
    # 添加第二个关键帧（5秒，旋转360度）
    key2 = anim_study.AddKeyframe(5.0)
    rot_matrix = create_rotation_matrix_z(2*np.pi)
    key2.SetComponentTransform("Crank", rot_matrix)
    
    # 生成动画
    anim_study.Generate()

---

4. 基本运动：轻量级动力学引擎

4.1 技术原理与简化假设

基本运动在动画基础上增加了简化的动力学计算，其核心假设包括：

1. 刚体假设：所有部件视为刚性（无变形）
2. 恒定质量分布：不考虑质心变化
3. 线性摩擦模型：使用库仑摩擦简化
4. 无约束反力计算：不求解拉格朗日乘子

其求解器实现：

class BasicMotionSolver:
    """基本运动求解器 - 简化动力学"""
    def __init__(self, gravity=np.array([0, -9.81, 0])):
        self.gravity = gravity
        self.bodies = {}  # {id: {mass, inertia, position, velocity}}
        self.joints = []  # 运动副列表
        self.forces = []  # 施加的力/力矩
        
    def add_body(self, body_id, mass, inertia, initial_pose):
        """添加刚体"""
        self.bodies[body_id] = {
            'mass': mass,
            'inertia': inertia,
            'pose': initial_pose,
            'velocity': np.zeros(6),  # [vx, vy, vz, wx, wy, wz]
            'force': np.zeros(6)
        }
    
    def add_spring(self, body1, body2, stiffness, damping, rest_length):
        """添加弹簧力"""
        self.forces.append({
            'type': 'spring',
            'bodies': (body1, body2),
            'k': stiffness,
            'c': damping,
            'L0': rest_length
        })
        
    def compute_forces(self):
        """计算所有作用力（简化版）"""
        for body_id, body in self.bodies.items():
            # 重力
            body['force'][:3] = body['mass'] * self.gravity
            
            # 弹簧力
            for spring in self.forces:
                b1, b2 = spring['bodies']
                if body_id in (b1, b2):
                    p1 = self.bodies[b1]['pose'][:3, 3]
                    p2 = self.bodies[b2]['pose'][:3, 3]
                    delta = p2 - p1
                    dist = np.linalg.norm(delta)
                    direction = delta / dist
                    
                    # 弹簧力：F = k*(L-L0) + c*v
                    v_rel = (self.bodies[b2]['velocity'][:3] - 
                            self.bodies[b1]['velocity'][:3])
                    v_proj = np.dot(v_rel, direction)
                    
                    force_mag = (spring['k'] * (dist - spring['L0']) + 
                                spring['c'] * v_proj)
                    
                    if body_id == b1:
                        body['force'][:3] += force_mag * direction
                    else:
                        body['force'][:3] -= force_mag * direction
    
    def step(self, dt):
        """单步积分（显式欧拉）"""
        self.compute_forces()
        
        for body_id, body in self.bodies.items():
            # 线加速度
            accel = body['force'][:3] / body['mass']
            
            # 角加速度（简化：假设惯性矩阵为球张量）
            angular_accel = body['force'][3:6] / body['inertia']
            
            # 更新速度
            body['velocity'][:3] += accel * dt
            body['velocity'][3:6] += angular_accel * dt
            
            # 更新位置
            body['pose'][:3, 3] += body['velocity'][:3] * dt
            
            # 更新旋转（使用小角度近似）
            omega = body['velocity'][3:6]
            delta_rot = self.skew_symmetric(omega * dt)
            body['pose'][:3, :3] = body['pose'][:3, :3] @ (np.eye(3) + delta_rot)
    
    def skew_symmetric(self, v):
        """构造反对称矩阵"""
        return np.array([
            [0, -v[2], v[1]],
            [v[2], 0, -v[0]],
            [-v[1], v[0], 0]
        ])

4.2 应用场景与局限性

适用场景：

• 机构干涉检测（考虑基本物理）
• 电机选型初步估算
• 机械手运动规划验证
• 简单弹簧-阻尼系统模拟

局限性：

• 无法精确计算约束反力
• 不考虑柔性体变形
• 摩擦模型过于简化
• 高速运动时数值不稳定

4.3 工程实践示例

在SolidWorks中使用基本运动进行凸轮机构分析：

def setup_cam_follower_basic_motion(sw_app):
    """设置凸轮从动件的基本运动分析"""
    model = sw_app.ActiveDoc
    study_mgr = model.GetStudyManager()
    
    # 创建基本运动算例
    study = study_mgr.CreateStudy("凸轮分析", swStudyType_e.swBasicMotion)
    
    # 设置重力
    study.SetGravity(np.array([0, -9.81, 0]))
    
    # 添加旋转马达（驱动凸轮）
    motor = study.AddMotor("Cam", swMotorType_e.swRotaryMotor)
    motor.SetConstantSpeed(60)  # 60 RPM
    
    # 添加弹簧（从动件回程）
    spring = study.AddSpring("Follower", "Ground", 
                             stiffness=1000,  # N/m
                             damping=10,      # N·s/m
                             free_length=0.05) # 初始长度
    
    # 设置结果输出
    result = study.GetResult()
    result.AddPlot("Follower_Displacement", swPlotType_e.swDisplacement)
    result.AddPlot("Spring_Force", swPlotType_e.swForce)
    
    # 运行仿真
    study.Run(5.0, 100)  # 5秒，100帧

---

5. Motion分析：精确的多体动力学

5.1 完整的动力学建模

Motion分析是运动算例的巅峰，它求解完整的约束多体动力学方程：

class MotionAnalysisSolver:
    """完整Motion分析求解器"""
    def __init__(self):
        self.bodies = []  # 刚体列表
        self.constraints = []  # 约束条件
        self.forces = []  # 力/力矩
        self.motors = []  # 驱动
        self.t = 0.0
        
    def add_revolute_joint(self, body1, body2, point, axis):
        """添加旋转副约束"""
        self.constraints.append({
            'type': 'revolute',
            'bodies': (body1, body2),
            'point': point,
            'axis': axis
        })
        
    def assemble_system(self):
        """组装系统方程"""
        n_bodies = len(self.bodies)
        n_constraints = len(self.constraints)
        
        # 构建质量矩阵（6n x 6n）
        M = np.zeros((6*n_bodies, 6*n_bodies))
        for i, body in enumerate(self.bodies):
            M[6*i:6*i+3, 6*i:6*i+3] = body['mass'] * np.eye(3)
            M[6*i+3:6*i+6, 6*i+3:6*i+6] = body['inertia']
        
        return M
    
    def compute_constraint_jacobian(self, q):
        """计算约束雅可比矩阵"""
        n_c = len(self.constraints)
        n_q = 6 * len(self.bodies)
        J = np.zeros((n_c, n_q))
        
        for i, constr in enumerate(self.constraints):
            b1, b2 = constr['bodies']
            p = constr['point']
            
            # 旋转副的约束方程导数
            # 位置约束：r1 + R1 * s1 - r2 - R2 * s2 = 0
            # 方向约束：R1 * axis1 - R2 * axis2 = 0
            
            # 简化示例：只计算位置约束部分
            J[i, 6*b1:6*b1+3] = np.eye(3)
            J[i, 6*b2:6*b2+3] = -np.eye(3)
            
        return J
    
    def solve_dynamics(self, dt, total_time):
        """求解完整动力学（使用隐式积分）"""
        n_q = 6 * len(self.bodies)
        q = np.zeros(n_q)  # 广义坐标
        qd = np.zeros(n_q)  # 广义速度
        M = self.assemble_system()
        
        results = []
        while self.t < total_time:
            # 计算当前状态
            J = self.compute_constraint_jacobian(q)
            
            # 计算外力
            F_ext = self.compute_external_forces(q, qd)
            
            # 求解带约束的动力学方程
            # [M, J^T; J, 0] * [qd_dot; lambda] = [F_ext; -J_dot*qd]
            
            # 构建系统矩阵
            n_c = len(self.constraints)
            A = np.block([
                [M, J.T],
                [J, np.zeros((n_c, n_c))]
            ])
            
            # 右端项
            rhs = np.concatenate([F_ext, -J @ qd])
            
            # 求解
            solution = np.linalg.solve(A, rhs)
            qdd = solution[:n_q]
            lambdas = solution[n_q:]