每日算法 — 使用java实现空竹：抛体运动学与角动量守恒模拟

引言

空竹是中国传统的民间玩具，通过两根抖杆和一根线绳使竹木制的空竹在空中高速旋转并做出各种抛接动作。本文用Java实现空竹运动的物理仿真系统，讲解抛体运动学建模与角动量守恒定律在旋转体动力学中的应用，模拟空竹的抛出、旋转、回落与重新接住的完整过程。

问题分析

空竹的结构特征

空竹可分为单轮空竹和双轮空竹，核心结构包括：
– 轮盘：提供转动惯量的主体，质量集中在轮缘
– 轴心：连接抖绳的细轴，也是旋转中心
– 风哨/风孔：高速旋转时产生哨声，同时受空气阻力影响

运动阶段划分

空竹的一次完整抛接动作包含四个阶段：
1. 加速旋转：双手交替拉动抖绳，通过摩擦力使空竹加速旋转
2. 抛射阶段：双手同时上提，空竹脱离绳约束做斜抛运动
3. 空中旋转：仅受重力与空气阻力矩作用，角速度逐渐衰减
4. 回落接住：空竹落回绳上，需匹配线速度与旋转速度

核心算法一：抛体运动学建模

空竹被抛出后的质心运动遵循斜抛运动规律，考虑空气阻力时的运动方程：

Fx = -k * vx * |v|   (空气阻力与速度平方成正比，方向相反)
Fy = -mg - k * vy * |v|

使用龙格-库塔四阶方法（RK4）数值积分求解微分方程，获得高精度的轨迹模拟。

核心算法二：角动量守恒与阻力矩

空竹在空中旋转时，主要受到空气阻力矩的作用：

τ = -c * ω * |ω|   (阻力矩与角速度方向相反)
I * dω/dt = τ      (转动定律)

其中 I 为空竹的转动惯量。对于轮盘状物体，I ≈ m * r²。

完整Java实现

package diabolo.physics;

/**
 * 空竹物理仿真系统
 * 
 * 核心功能：
 * 1. 抛体运动学：重力 + 空气阻力的RK4数值积分
 * 2. 旋转动力学：角动量守恒与空气阻力矩衰减
 * 3. 绳约束检测：判断空竹是否处于绳上
 * 4. 抛出-接住循环：模拟完整的抛接动作
 */
public class DiaboloSimulator {

    // 物理常量
    private static final double GRAVITY = 9.81;        // 重力加速度 m/s²
    private static final double AIR_DENSITY = 1.225;   // 空气密度 kg/m³
    private static final double TIME_STEP = 0.001;     // 仿真步长 1ms

    // 空竹参数
    private final double mass;           // 质量 kg
    private final double radius;         // 轮盘半径 m
    private final double axisLength;     // 轴长 m
    private final double dragCoeff;      // 空气阻力系数
    private final double momentDragCoeff; // 旋转阻力系数

    // 转动惯量 (轮盘模型: I = m * r²)
    private final double momentOfInertia;

    // 运动状态
    private double x, y;          // 质心位置 (m)
    private double vx, vy;        // 质心速度 (m/s)
    private double omega;         // 角速度 (rad/s)
    private double theta;         // 姿态角 (rad)
    private boolean onString;     // 是否在绳上

    // 绳参数
    private final double stringLength;  // 绳长 m
    private final double handSpacing;   // 两手间距 m
    private double handY;              // 手的高度 m

    public DiaboloSimulator(double mass, double radius, double axisLength,
                            double dragCoeff, double momentDragCoeff,
                            double stringLength, double handSpacing) {
        this.mass = mass;
        this.radius = radius;
        this.axisLength = axisLength;
        this.dragCoeff = dragCoeff;
        this.momentDragCoeff = momentDragCoeff;
        this.stringLength = stringLength;
        this.handSpacing = handSpacing;
        this.momentOfInertia = mass * radius * radius;
        reset();
    }

    public void reset() {
        x = 0;
        y = 0.5;  // 初始高度
        vx = 0;
        vy = 0;
        omega = 0;
        theta = 0;
        onString = true;
        handY = 1.0;
    }

    /**
     * 通过抖绳加速空竹旋转
     * @param pullForce 拉力 N
     * @param pullDuration 拉动持续时间 s
     */
    public void accelerateByPulling(double pullForce, double pullDuration) {
        if (!onString) {
            System.out.println("警告：空竹不在绳上，无法加速");
            return;
        }

        // 拉力通过绳与轴的摩擦产生力矩
        // 力矩 = 拉力 × 轴半径（简化模型）
        double torque = pullForce * (axisLength / 2);
        double angularAccel = torque / momentOfInertia;

        // 模拟拉动过程
        double steps = pullDuration / TIME_STEP;
        for (int i = 0; i < steps; i++) {
            omega += angularAccel * TIME_STEP;
            // 同时绳的轻微上提产生小幅弹跳
            y += 0.001 * Math.sin(i * 0.1);
        }

        System.out.printf("加速完成：角速度 = %.2f rad/s (%.2f RPM)%n",
            omega, omega * 60 / (2 * Math.PI));
    }

    /**
     * 执行抛出动作
     * @param throwVx 水平抛出速度 m/s
     * @param throwVy 垂直抛出速度 m/s
     */
    public void throwDiabolo(double throwVx, double throwVy) {
        if (!onString) {
            System.out.println("空竹已在空中");
            return;
        }
        onString = false;
        vx = throwVx;
        vy = throwVy;
        System.out.printf("抛出空竹：v=(%.2f, %.2f) m/s, ω=%.2f rad/s%n", vx, vy, omega);
    }

    /**
     * RK4四阶龙格-库塔积分步进
     * 求解二维运动方程（含空气阻力）
     */
    public void rk4Step() {
        // 状态向量: [x, y, vx, vy, omega, theta]
        double[] state = {x, y, vx, vy, omega, theta};
        double[] k1 = derivatives(state);
        double[] k2 = derivatives(add(state, multiply(k1, TIME_STEP / 2)));
        double[] k3 = derivatives(add(state, multiply(k2, TIME_STEP / 2)));
        double[] k4 = derivatives(add(state, multiply(k3, TIME_STEP)));

        // 加权平均: y_{n+1} = y_n + (k1 + 2*k2 + 2*k3 + k4) * h/6
        double[] delta = multiply(add(add(k1, multiply(k2, 2)), 
            add(multiply(k3, 2), k4)), TIME_STEP / 6);

        x += delta[0];
        y += delta[1];
        vx += delta[2];
        vy += delta[3];
        omega += delta[4];
        theta += delta[5];

        // 角速度衰减（空气阻力矩）
        applyRotationalDrag();
    }

    /**
     * 计算状态向量的导数
     * 返回 [dx/dt, dy/dt, dvx/dt, dvy/dt, domega/dt, dtheta/dt]
     */
    private double[] derivatives(double[] s) {
        double vx_ = s[2];
        double vy_ = s[3];
        double omega_ = s[4];

        double v = Math.sqrt(vx_ * vx_ + vy_ * vy_);
        double dragFx = 0, dragFy = 0;

        if (v > 1e-6) {
            // 空气阻力 Fd = 0.5 * ρ * Cd * A * v²
            double area = Math.PI * radius * radius;
            double dragMag = 0.5 * AIR_DENSITY * dragCoeff * area * v * v;
            dragFx = -dragMag * vx_ / v;
            dragFy = -dragMag * vy_ / v;
        }

        double ax = dragFx / mass;
        double ay = -GRAVITY + dragFy / mass;

        // 旋转阻力矩 (简化模型)
        double torqueDrag = -momentDragCoeff * omega_ * Math.abs(omega_);
        double alpha = torqueDrag / momentOfInertia;

        return new double[]{vx_, vy_, ax, ay, alpha, omega_};
    }

    /**
     * 应用旋转阻力导致的角速度衰减
     */
    private void applyRotationalDrag() {
        // 阻力矩导致的角速度指数衰减
        double decayRate = momentDragCoeff / momentOfInertia;
        omega *= Math.exp(-decayRate * TIME_STEP);
    }

    /**
     * 检测空竹是否回落到绳上
     * 绳的形状简化为两固定点间的悬链线或直线
     */
    public boolean checkCatch() {
        if (onString) return true;

        // 简化的绳检测：判断空竹y坐标是否低于绳高度
        // 绳在两手间呈V字形，中心最低点
        double stringCenterY = handY - Math.sqrt(
            Math.pow(stringLength / 2, 2) - Math.pow(handSpacing / 2, 2));

        // 空竹落回绳的条件：高度接近绳高度，水平位置在两手之间
        double halfSpan = handSpacing / 2;
        if (y <= stringCenterY + radius && Math.abs(x) < halfSpan && vy < 0) {
            onString = true;
            vy = -vy * 0.3; // 弹性碰撞，能量损失
            System.out.printf("接住空竹！位置=(%.3f, %.3f), 角速度=%.2f rad/s%n", x, y, omega);
            return true;
        }
        return false;
    }

    /**
     * 模拟一次完整的抛接过程
     */
    public void simulateTossAndCatch(double throwVx, double throwVy) {
        reset();

        // 阶段1：加速到足够转速
        System.out.println("===== 阶段1：加速旋转 =====");
        accelerateByPulling(5.0, 0.5);

        // 阶段2：抛出
        System.out.println("\n===== 阶段2：抛出 =====");
        throwDiabolo(throwVx, throwVy);

        // 阶段3：空中飞行
        System.out.println("\n===== 阶段3：空中飞行 =====");
        double maxHeight = y;
        double flightTime = 0;
        int printInterval = 100; // 每100步打印一次

        while (!onString && y > 0) {
            rk4Step();
            flightTime += TIME_STEP;
            if (y > maxHeight) maxHeight = y;

            checkCatch();

            if ((int)(flightTime / TIME_STEP) % printInterval == 0) {
                System.out.printf("t=%.3fs: pos=(%.3f, %.3f), v=(%.2f, %.2f), ω=%.1f rad/s%n",
                    flightTime, x, y, vx, vy, omega);
            }
        }

        // 阶段4：结果统计
        System.out.println("\n===== 抛接统计 =====");
        System.out.printf("飞行时间: %.3f s%n", flightTime);
        System.out.printf("最大高度: %.3f m%n", maxHeight);
        System.out.printf("水平位移: %.3f m%n", x);
        System.out.printf("落回时角速度: %.2f rad/s%n", omega);
        System.out.printf("是否成功接住: %s%n", onString ? "是" : "否（空竹落地）");
    }

    // 向量工具方法
    private double[] add(double[] a, double[] b) {
        double[] r = new double[a.length];
        for (int i = 0; i < a.length; i++) r[i] = a[i] + b[i];
        return r;
    }

    private double[] multiply(double[] a, double scalar) {
        double[] r = new double[a.length];
        for (int i = 0; i < a.length; i++) r[i] = a[i] * scalar;
        return r;
    }

    // Getters
    public double getX() { return x; }
    public double getY() { return y; }
    public double getOmega() { return omega; }
    public boolean isOnString() { return onString; }

    // ==================== 测试入口 ====================
    public static void main(String[] args) {
        // 典型空竹参数
        // 质量约200g，轮径约13cm，轴长约5cm
        DiaboloSimulator diabolo = new DiaboloSimulator(
            0.20,   // mass kg
            0.065,  // radius m
            0.05,   // axis length m
            0.47,   // drag coefficient (近似球体)
            0.001,  // rotational drag coefficient
            1.5,    // string length m
            0.5     // hand spacing m
        );

        System.out.println("===== 空竹物理仿真系统 =====\n");

        // 低抛：垂直上抛 3 m/s
        System.out.println("【测试1】低抛（vy=3.0 m/s）");
        diabolo.simulateTossAndCatch(0.5, 3.0);

        System.out.println("\n" + "=".repeat(50) + "\n");

        // 高抛：垂直上抛 6 m/s
        System.out.println("【测试2】高抛（vy=6.0 m/s）");
        diabolo.simulateTossAndCatch(0.3, 6.0);

        System.out.println("\n" + "=".repeat(50) + "\n");

        // 斜抛：带水平速度的抛接
        System.out.println("【测试3】斜抛（vx=2.0, vy=4.0 m/s）");
        diabolo.simulateTossAndCatch(2.0, 4.0);
    }
}

复杂度分析

扩展与优化

1. 绳的悬链线模型：用双曲余弦函数建模抖绳的真实形状，提升接住检测精度
2. 三维扩展：增加z轴运动，模拟绕身、绕腿等复杂空竹动作
3. 声音模拟：根据角速度和风孔面积计算哨声音频频率
4. 实时渲染：将物理状态同步到JavaFX或LibGDX渲染管线
5. 参数调优：通过实验数据拟合更精确的空气阻力系数

总结

本文用Java实现了空竹运动的物理仿真系统，将经典力学中的抛体运动与转动定律通过数值方法进行了精确模拟。RK4积分方法保证了轨迹计算的高精度，而角动量守恒的引入让旋转衰减更加真实。这个案例展示了如何将大学物理中的抽象公式转化为可运行的代码，是理解数值积分与刚体动力学的优秀入门项目。