Springs and spring-like motion are essential components in various aspects of physics simulations. They can be used for anything from creating cool animations to cloth physics and collision.

Overview

1. What you'll need
2. The theory
3. Implementation
4. Linear vs nonlinear springs
5. References

What you'll need

To explore springs in both rigidbody and particle physics, you should have the following components:

For your particle system:

• A particle class containing attributes such as mass and position.
• The ability to apply forces, set velocities, and perform integration.

For rigidbody physics:

• Rigidbodies equipped with collision shapes (which can be simple or complex) mass, and a transform.
• The capability to apply forces at specific locations and perform integration.

Note that for the latter you’ll probably need a more advanced integration algorithm, such as semi-implicit Euler. This is due to things like damping not always being applied correctly, causing the object to e.g. spin forever. Without going into too much detail over the difference between explicit and semi-implicit Euler, here’s a short code snippet of the integration & applying force.

void Cinkes::CRigidBody::AddForceAtPoint(const CVector3& a_ForceToAdd, const CVector3& a_Point)
{
	m_Force += a_ForceToAdd;
	CVector3 rel = (a_Point - GetTransform().getOrigin());
	m_Torque += a_ForceToAdd.Cross(rel);
}
void Cinkes::CRigidBody::Integrate(CScalar a_T)
{
	if (m_InverseMass < 0.001f)
		return;

	// Calculate linear acceleration from force inputs.
	m_Acceleration = m_Force * m_InverseMass;

	// Calculate angular acceleration from torque inputs.
	CVector3 angularAcceleration =
		m_InverseIntertiaTensorWorld * (m_Torque);

	// Adjust velocities
	// Update linear velocity from both acceleration and impulse.
	m_Velocity += m_Acceleration * a_T;

	// Update angular velocity from both acceleration and impulse.
	m_AngularVelocity += angularAcceleration * a_T;

	// Impose drag.
	m_Velocity *= CUtils::Pow(m_LinearDamping, a_T);
	m_AngularVelocity *= CUtils::Pow(m_AngularDamping, a_T);

	// Adjust positions
	// Update linear position.
	m_Transform.setOrigin(m_Transform.getOrigin() + m_Velocity * a_T);

	// Update angular position.
	CQuaternion q = m_Transform.getQuaternion();
	CVector3 temp = m_AngularVelocity * a_T;
	CQuaternion o = CQuaternion(temp[0], temp[1], temp[2], 0);
	o *= 0.5;
	q += o;
	m_Transform.setBasis(CMat3x3(q)); 

	// Normalise the orientation, and update the matrices with the new
	// position and orientation
	SetInverseInertiaTensorWorld();

	// Clear accumulators.
	ClearForces();

}

The theory

Newton's law

You’ve probably heard of Newton’s third law by now: every action has an equal and opposite reaction. Newton's third law plays a crucial role in understanding the behavior of springs. When a spring is compressed or stretched, it exerts a force in one direction. However, according to Newton's third law, there is an equal and opposite force exerted by the spring in the opposite direction. This force is what allows the spring to bounce back or elongate.
For instance, when you compress a spring by pushing on it, the spring pushes back with the same amount of force. This balanced force creates a restoring force that tries to return the spring to its original shape. Similarly, when you stretch a spring, it pulls back with an equal and opposite force, attempting to bring the spring back to its relaxed state.

Hooke’s law

The working of springs, at least in terms of physics programming is described by Hooke’s law. According to it, the force exerted by a spring is directly proportional to the amount it is stretched or compressed.
This relationship can be expressed as
F = kx
where F represents the force, k is the spring constant (a measure of the spring's stiffness), and x is the displacement from the spring's equilibrium position.
Hooke's law allows us to quantify the relationship between the applied force and the resulting deformation of a spring. As we stretch or compress a spring, the force it exerts increases proportionally. This linearity allows for predictable behavior, making Hooke's law a valuable tool for calculating and understanding spring forces in various situations (we’ll talk about nonlinear springs at the end of the article).
When we introduce damping to the equation of a spring-mass system, we consider the effects of a resisting force that opposes the motion of the system. This resistance is known as damping and is typically caused by factors such as air resistance or friction.

Damping

If you implement the above equation, you might notice that the spring will oscillate forever, never reaching equilibrium (the state where the force applied to the object, e.g. gravity is equal to the force applied by the spring and therefore the object stops moving).
Damping can be incorporated into the equation of motion through the addition of a damping coefficient. This coefficient accounts for the velocity of the object and introduces a proportional damping force. The resulting equation becomes:
F = kx - bv
where F represents the total force, k is the spring constant, x is the displacement from equilibrium, b is the damping coefficient, and v is the velocity of the object. The damping coefficient, b, determines the strength of the damping force. Higher values of b correspond to stronger damping, which leads to the object reaching equilibrium faster.

And that’s it for the theory! Now let’s look at some code.

Implementation

You can connect pretty much anything with springs: a point and an object, two objects, two points, whatever your heart desires. For the purposes of this article, we will connect one static point to one rigidbody (or particle).

The spring class

Other than the basic data that we talked about in the previous section, we also need to add a connection point. This is important because the spring is connected to the object at a point on the rigidbody. This will create torque, which will introduce angular movement and therefore rotate the object.

class CSpring
{
public:
	CRigidBody* m_Body[2]{};
	CVector3 m_Point[2];
	CScalar m_SpringConstant;
	CScalar m_RestLength;
	CScalar m_DampeningConstant;
	CScalar m_CubicStiffness;
	ESPRING_TYPE m_Type;
};

The force generator

The Cinkes library works based on force generators. These provide a modular way to add and remove forces such as gravity easily. First, let’s start with processing the data received from the spring class.

	void UpdateForce(void* a_Body) override
	{
		m_Current = static_cast(a_Body);
		
		//one end is some type of object
		if (m_Current->GetBody1() != nullptr)
		{
			m_Positions[0] = m_Current->GetBody1()->GetTransform().getOrigin() + m_Current->GetPoint1();
			if (m_Current->GetBody1()->GetType() == EOBJECT_TYPE::TYPE_RIGID) { m_Velocities[0] = m_Current->GetBody1()->GetLinearVelocity(); }
		}
		else
		{
		//end is only a point in space
			m_Positions[0] = m_Current->GetPoint1();
		}
		//other is an object
		if (m_Current->GetBody2() != nullptr)
		{
			m_Positions[1] = m_Current->GetBody2()->GetTransform().getOrigin() + m_Current->GetPoint2();
			if (m_Current->GetBody2()->GetType() == EOBJECT_TYPE::TYPE_RIGID) { m_Velocities[1] = m_Current->GetBody2()->GetLinearVelocity(); }
		}
		else
		{
		//other end is a point
			m_Positions[1] = m_Current->GetPoint2();

		}

		//calculate velocity if they are rigidbodies
		if (m_Current->GetBody1() != nullptr) CalculateFirst();
		if (m_Current->GetBody2() != nullptr) CalculateSecond();
	}

This is fairly self-explanatory. We calculate the spring forces acting on the objects separately since we can’t always be sure that they are rigidbodies (remember, they can just be points in space). This might open up some room for optimization, but for simplicity’s sake we’ll keep it separate.

Now let's take a look at the actual implementations:

void CalculateFirst() const
{
	CVector3 displacement = m_Positions[1] - m_Positions[0];
	CScalar distance = displacement.Length();
	displacement.Normalize();
	
	// Calculate the spring force based on Hooke's Law
	//We multiply by the displacement to give the force a direction (all other values are scalar)
	CVector3 springForce =  displacement * m_Current->GetSpringConstant() * (distance - m_Current->GetRestLength());

	// Calculate the damping force based on damping coefficient
	//We don't need to worry about the rest length, so there is no reason to multiply by the direction
	CVector3 dampingForce =  (m_Velocities[1] - m_Velocities[0]) * m_Current->GetDampeningConstant();

	// Return the total force applied at the attachment point
	//m_Positions[0] is a point in world space
	m_Current->GetBody1()->AddForceAtPoint(springForce + dampingForce,m_Positions[0]);
}

Note that depending on your implementation of AddForceAtPoint() function either needs coordinates in local or world space. As we’ve seen in section 1, when actually applying the torque (using the cross product), local space coordinates need to be used.

Linear vs. nonlinear springs

While in 90% of cases Hooke’s law will be enough for you, there are cases where you might want to go as realistic as possible. An example for a linear spring is, for example a metal oscillator. A nonlinear spring is, for example a rubber band.
There are various equations/algorithms/mathematical models to simulate nonlinear springs. These range from adding another coefficient to completely new, complex algorithms. A commonly used method for rubber-like springs is the Mooney-Rivlin model, which goes like so: F = c₁(I₁ - 3) + c₂(I₂ - 3)
where F represents the force exerted by the spring, c₁ and c₂ are material constants, and I₁ and I₂ are the first and second invariants of the deformation tensor, respectively. The first invariant (I₁) is a measure of the volume change of the material, while the second invariant (I₂) characterizes the deviatoric (shear) deformation.

By incorporating these invariants into the equation, the Mooney-Rivlin model captures complex nonlinear behaviors that arise in real-world materials. In code it would look something like this:

CScalar calculateNonlinearSpringForce(CScalar a_C1, CScalar a_C2) {
    // Calculate the displacement between the two attachment points
    CScalar displacementX = m_Positions[1][0] - m_Positions[0][0];
    CScalar displacementY = m_Positions[1][1] - m_Positions[0][1];
    CScalar displacementZ = m_Positions[1][2] - m_Positions[0][2];

    // Calculate the Euclidean distance between the attachment points
    CScalar deformation = std::sqrt(displacementX * displacementX + displacementY * displacementY + displacementZ * displacementZ);

    // Calculate the force based on the Mooney-Rivlin model equation
    CScalar I1 = deformation - 3.0;
    CScalar I2 = deformation * deformation - 3.0;
    CScalar force = a_C1 * I1 + a_C2 * I2;
	m_Current->GetBody1()->AddForceAtPoint(force, m_Positions[0]);
    return force;
}

It’s worth noting again, that unless you’re planning to make a very advanced library (such as Bullet or PhysX), you’ll be fine without nonlinear springs. This is mostly here for the sake of completion.

And that’s it! Now you can make catapults and bouncy balls to your heart’s content!

References:

• Book: Game Physics Engine Development by Ian Millington
• My own implementation on Github