作业八：质点弹簧系统

我突然感觉之前的文章内容太乱了，详略不得当！这篇文章尝试大道至简。

本文内容：构造质点弹簧系统。两种模拟方法实现：隐式/显式欧拉法 与 Verlet积分数值方法。

HW8:物理仿真之质点弹簧系统

构造Rope对象

在 rope.cpp 中，实现 Rope 绳索类的构造。该对象从起点start开始，到end结束，中间包含了 num_nodes个节点。每一个节点 mess 都有质量，所以称为质点。相邻两个质点之间的线段视为弹簧。在作业框架的构造函数中，传入了几个参数：

• start 和 end: 绳索的起点和终点

• num_nodes: 质点数量

• node_mass: 质点质量

• k: 弹簧弹性系数

• pinned_nodes: 指示了哪个质点是固定的

线性插值的公式：

$$\text{pos} = \text{start} + (\text{end} - \text{start}) \times \frac{i}{\text{num\_nodes} - 1.0}$$

其中，$i$表示第几个节点。

Rope::Rope(Vector2D start, Vector2D end, int num_nodes, float node_mass, float k, vector<int> pinned_nodes) {
    // Create a rope starting at `start`, ending at `end`, and containing `num_nodes` nodes.
    for(int i = 0; i < num_nodes; i++) {
        Vector2D currentPos = start + (end - start) * static_cast<double>(i) / (num_nodes - 1.0);
        Mass* newMass = new Mass(currentPos, node_mass, false);
        masses.push_back(newMass);
    }

    for(int i = 0; i < num_nodes - 1; i++) {
        Spring* newSpring = new Spring(masses[i], masses[i + 1], k);
        springs.push_back(newSpring);
    }

    for (auto &i : pinned_nodes) {
        masses[i]->pinned = true;
    }
}

完成后的效果是，可以看到屏幕上画出绳子，但它不发生运动。

胡克定律 与 隐式/显式欧拉法

首先使用胡克定律处理弹簧力（中间部分把两个位置向量的差单位化了，最右变括号内是弹簧的伸缩量）：

$$\mathbf{f}_{\mathbf{b} \rightarrow \mathbf{a}}=-k_s \frac{\mathbf{b}-\mathbf{a}}{\|\mathbf{b}-\mathbf{a}\|}(\|\mathbf{b}-\mathbf{a}\|-l)$$

遍历所有的弹簧，对弹簧两端的质点施加正确的弹簧力，并且将计算出的力加到 m1 质点上，并从 m2 质点上减去，保证力的作用和反作用平衡。

接下来是实现隐式欧拉法，即：

$$\begin{aligned} & \mathrm{F}=\mathrm{ma}\\ & v(t+1)=v(t)+a(t) \cdot dt\\ & x(t+1)=x(t)+v(t) \cdot dt &&\text{ // For explicit method}\\ & x(t+1)=x(t)+v(t+1) \cdot dt &&\text{ // For semi-implicit method} \end{aligned}$$

遍历所有的质点，对于每一个未固定的质点，计算其受到的总力并更新其速度和位置。包括重力、全局阻尼力（我定义为：dampingConstant）。这里的阻尼力我在Assignment中找不到更详细的说明了，这里直接做一个与速度反方向的小阻力，且这个阻力与速度成正比。

void Rope::simulateEuler(float delta_t, Vector2D gravity) {
    // Calculate forces due to springs using Hooke's Law
    for (auto &s : springs) {
        Vector2D displacement = s->m2->position - s->m1->position;
        float stretchAmount = displacement.norm() - s->rest_length;
        Vector2D springForce = s->k * displacement.unit() * stretchAmount;

        s->m1->forces += springForce;
        s->m2->forces -= springForce;
    }

    // Damping constant for global damping
    const float dampingConstant = 0.01;

    // Update forces, velocities, and positions for each mass
    for (auto &m : masses) {
        if (!m->pinned) {
            // Add gravitational force
            m->forces += gravity * m->mass;

            // Add damping force
            Vector2D dampingForce = -dampingConstant * m->velocity;
            m->forces += dampingForce;

            // Compute acceleration
            Vector2D acceleration = m->forces / m->mass;

            // Update velocity and position using implicit Euler
            m->velocity += acceleration * delta_t;
            m->position += m->velocity * delta_t;
            
            // //explicit  Euler  
            // m->position += m->velocity * delta_t;
            // m->velocity += acceleration * delta_t;
            // */
        }

        // Reset forces for next iteration
        m->forces = Vector2D(0, 0);
    }
}

• explicit欧拉方法需要很高的SPF才能够稳定

显式 Verlet

Verlet方法的优势在于其简单性和能量保守特性。其基本思想是用当前和前一时刻的位置，而不是速度，来估算下一时刻的位置。基本公式如下：

$$x(t+Δt)=2x(t)−x(t−Δt)+a(t)Δt^2$$

其中，

• $x(t+Δt) $是下一个时刻的位置。

• $x(t)$ 是当前时刻的位置。

• $x(t−Δt)$ 是前一个时刻的位置。

• $a(t)$ 是当前时刻的加速度。

• $Δt$ 是时间步长。

速度公式（本文不需要用到）：

$$v(t)=\frac{x(t+\Delta t)-x(t-\Delta t)}{2 \Delta t}$$

void Rope::simulateVerlet(float delta_t, Vector2D gravity) {
    for (auto &s: springs) {
        // TODO (Part 3): Simulate one timestep of the rope using explicit Verlet （solving constraints)
        Vector2D ab = s->m2->position - s->m1->position;
        float length = ab.norm();
        Vector2D forceDirection = ab.unit();
        Vector2D force = s->k * (length - s->rest_length) * forceDirection;

        s->m1->forces += force;
        s->m2->forces -= force;
    }
    for (auto &m: masses) {
        Vector2D accelerations = (gravity + (m->forces / m->mass));
        if (!m->pinned) {
            // TODO (Part 3.1): Set the new position of the rope mass
            Vector2D temp_position = m->position;
            // TODO (Part 4): Add global Verlet damping
            double damping_factor = 0.000001;// 阻尼
            m->position += (1 - damping_factor) * (m->position - m->last_position);
            m->position += accelerations * delta_t * delta_t;
            m->last_position = temp_position;
        }
        m->forces = Vector2D(0, 0);
    }
}

• 第一个循环：applySpringForces

• 第二个循环：updateMassPositions

蓝色是欧拉方法，绿色是Verlet方法。

Reference

1. GAMES101 Lingqi Yan

2. https://www.cnblogs.com/coolwx/p/15059551.html