作业八:质点弹簧系统

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

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

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

构造Rope对象

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

  • startend: 绳索的起点和终点

  • num_nodes: 质点数量

  • node_mass: 质点质量

  • k: 弹簧弹性系数

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

线性插值的公式:

其中,$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;
    }
}

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

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

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

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

接下来是实现隐式欧拉法,即:

遍历所有的质点,对于每一个未固定的质点,计算其受到的总力并更新其速度和位置。包括重力、全局阻尼力(我定义为: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) $是下一个时刻的位置。

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

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

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

  • $Δ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