Two-mass spring-damper between walls
Two masses between two walls, tied to the walls and to each other by springs, with dampers and forces.
- Name:
two_mass_spring_damper - Field:
physics/mechanics - Formalism: ordinary differential equations
- References: Ryzhov et al. 2021
Factory
@register_model(
field='physics/mechanics',
tags=('mechanics', 'ode', 'oscillation', 'normal modes', 'forcing'),
references=('Ryzhov et al. 2021',),
defaults={'until': 30.0, 'dt': 0.05},
)
def two_mass_spring_damper(
x1: float = 1.2,
x2: float = 2.0,
*,
m1: float = 1.0,
m2: float = 1.0,
k1: float = 1.0,
k2: float = 1.0,
k: float = 1.0,
b1: float = 0.0,
b2: float = 0.0,
l1: float = 1.0,
l2: float = 1.0,
F1: float = 0.0,
F2: float = 0.0,
omega1: float = 1.0,
omega2: float = 1.0,
) -> System:
"""Return two masses between two walls, tied to the walls and to each other by springs, with dampers and forces.
The first mass rests at `x01 = l1` and the second at `x02 = l1 + l2`. They follow
`m1 x1'' = -k1 (x1 - x01) + k ((x2 - x1) - (x02 - x01)) - b1 x1' + F1 cos(ω1 t)` and
`m2 x2'' = -k ((x2 - x1) - (x02 - x01)) - k2 (x2 - x02) - b2 x2' + F2 cos(ω2 t)`.
Args:
x1:
The initial position of the first mass.
x2:
The initial position of the second mass.
m1:
The first mass.
m2:
The second mass.
k1:
The stiffness of the spring between the left wall and the first mass.
k2:
The stiffness of the spring between the second mass and the right wall.
k:
The stiffness of the spring between the masses.
b1:
The damping coefficient of the first mass.
b2:
The damping coefficient of the second mass.
l1:
The rest length of the first spring.
l2:
The rest length of the middle spring.
F1:
The amplitude of the periodic force on the first mass.
F2:
The amplitude of the periodic force on the second mass.
omega1:
The angular frequency of the force on the first mass.
omega2:
The angular frequency of the force on the second mass.
Returns:
system:
The model, with the masses `first` and `second` and their positions `first.x` and `second.x`.
Examples:
>>> from skcomplex.models.physics.mechanics import two_mass_spring_damper
>>> two_mass_spring_damper(l1=2.0).interactions[0].center
2.0
"""
rest1, rest2 = l1, l1 + l2
return System(
elements=[MassPointElement('first', x=x1, x_d1=0.0, m=m1), MassPointElement('second', x=x2, x_d1=0.0, m=m2)],
interactions=[
HarmonicPotentialInteraction('left_spring', target='first', k=k1, center=rest1),
HarmonicPotentialInteraction('right_spring', target='second', k=k2, center=rest2),
SpringInteraction('middle_spring', first='first', second='second', k=k, length=l2),
DampingInteraction('first_damper', target='first', gamma=b1),
DampingInteraction('second_damper', target='second', gamma=b2),
HarmonicForceInteraction('first_force', target='first', F=F1, omega=omega1),
HarmonicForceInteraction('second_force', target='second', F=F2, omega=omega2),
],
space=ContinuousSpace(n_dim=1),
)
Components
class MassPointElement(PointElement):
"""Point mass, a point whose position `x` has order 2, with its velocity `x_d1` and a mass `m`.
The forces on the point divide by the mass through their divisor `'m'`.
Args:
label:
The label of the point.
x:
The initial position. A number applies to every component.
x_d1:
The initial velocity. A number applies to every component.
visual:
How the point is drawn, or `None` for `default_visual_`.
m:
The mass.
Attributes:
default_visual_:
Points at the positions.
Examples:
>>> from skcomplex.models.physics.common.particles import MassPointElement
>>> MassPointElement('ball', x=[0.0, 1.0], m=2.0).m
2.0
"""
def __init__(
self,
label: str,
x: npt.ArrayLike = 0.0,
x_d1: npt.ArrayLike = 0.0,
visual: VisualSpec | None = None,
*,
m: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
super().__init__(label, x=x, x_d1=x_d1, visual=visual)
self.m = m
class HarmonicPotentialInteraction(PotentialEnergyInteraction):
"""Harmonic potential that pulls a particle towards a center.
The energy is `k * |x - center|**2 / 2`.
Args:
label:
The label of the interaction.
target:
The path of the particle.
k:
The stiffness.
center:
The center. A number applies to every component.
states:
The name of the vector state that holds the position of the particle.
divisor:
The name of the state or parameter that holds the mass of the particle.
Examples:
>>> from skcomplex.models.physics.common.particles import HarmonicPotentialInteraction
>>> HarmonicPotentialInteraction('trap', target='ball', k=2.0).k
2.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
k: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
center: Annotated[npt.ArrayLike, ParameterSpec(vector=True)] = 0.0,
states: str = 'x',
divisor: str = 'm',
) -> None:
super().__init__(label, target=target, states=states, divisor=divisor)
self.k = k
self.center = center
def _get_center(self, s: SimpleNamespace) -> object:
"""Return the center, the hook a subclass overrides to move it."""
return s.center
def _get_potential(self, s: SimpleNamespace) -> Expression:
"""Return `k * |x - center|**2 / 2`.
Args:
s:
The namespace of the symbols the energy can use.
Returns:
energy:
The potential energy as a symbolic expression.
"""
displacement = getattr(s.target, self._get_names()[0]) - self._get_center(s)
return s.k * displacement.dot_product(displacement) / 2
class SpringInteraction(PotentialEnergyInteraction):
"""Spring between two particles.
The energy is `k * (|x_first - x_second| - length)**2 / 2`, and its gradient moves both particles.
Args:
label:
The label of the interaction.
first:
The path of the first particle.
second:
The path of the second particle.
k:
The stiffness.
length:
The rest length.
states:
The name of the vector state that holds the positions of the particles.
divisor:
The name of the state or parameter that holds the masses of the particles.
Examples:
>>> from skcomplex.models.physics.common.particles import SpringInteraction
>>> SpringInteraction('spring', first='a', second='b', k=3.0, length=1.0).length
1.0
"""
def __init__(
self,
label: str,
*,
first: Annotated[str, RoleSpec()],
second: Annotated[str, RoleSpec()],
k: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
length: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
states: str = 'x',
divisor: str = 'm',
) -> None:
super().__init__(label, target=first, states=states, divisor=divisor)
self.first = first
self.second = second
self.k = k
self.length = length
def _get_potential(self, s: SimpleNamespace) -> Expression:
"""Return `k * (|x_first - x_second| - length)**2 / 2`.
Args:
s:
The namespace of the symbols the energy can use.
Returns:
energy:
The potential energy as a symbolic expression.
"""
gap = getattr(s.first, self._get_names()[0]) - getattr(s.second, self._get_names()[0])
return s.k * (sqrt(gap.dot_product(gap)) - s.length) ** 2 / 2
class DampingInteraction(RateInteraction):
"""Linear damping, a force `-gamma * v` proportional and opposite to the velocity `v` of a target.
The velocity is the first derivative of the position, a state of order 2, whose second derivative changes by
`-gamma * v / m` with the mass as divisor. The LAMMPS backend runs it as `fix viscous`.
Args:
label:
The label of the interaction.
target:
The path of the element or group whose velocity is damped.
gamma:
The damping coefficient.
state:
The name of the state of order 2 that holds the position of the target.
divisor:
The name of the state or parameter that holds the mass of the target.
Examples:
>>> from types import SimpleNamespace
>>> from sage.all__sagemath_symbolics import var
>>> from skcomplex.models.physics.common.particles import DampingInteraction
>>> target = SimpleNamespace(x=var('x'), m=2)
>>> damping = DampingInteraction('air', target='ball', gamma=0.1).check().prepare()
>>> [law.value for law in damping.run(SimpleNamespace(target=target, gamma=var('gamma'))).laws_]
[-1/2*gamma*derivative(x, 1)]
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
gamma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
state: str = 'x',
divisor: str = 'm',
) -> None:
super().__init__(label, target=target, state=state, divisor=divisor)
self.gamma = gamma
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the drag `-gamma * v`.
Args:
namespace:
The namespace of the symbols the rate can use.
Returns:
rate:
The drag, before the division by the mass.
"""
return -namespace.gamma * velocity_of(getattr(namespace.target, self.state))
class HarmonicForceInteraction(RateInteraction):
"""Periodic force `F cos(omega t)` on a target, the second derivative `F cos(omega t) / m` of its position.
It is the rate `F cos(omega t)` with the mass as divisor, applied to every component of the position.
Args:
label:
The label of the interaction.
target:
The path of the element or group the force acts on.
F:
The amplitude of the force.
omega:
The angular frequency of the force.
state:
The name of the state of order 2 that holds the position of the target.
divisor:
The name of the state or parameter that holds the mass of the target.
Examples:
>>> from skcomplex.models.physics.mechanics import HarmonicForceInteraction
>>> HarmonicForceInteraction('drive', target='mass', F=0.5, omega=2.0).omega
2.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
F: Annotated[float, ParameterSpec()] = 0.0,
omega: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
state: str = 'x',
divisor: str = 'm',
) -> None:
super().__init__(label, target=target, state=state, divisor=divisor)
self.F = F
self.omega = omega
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the force `F cos(omega t)`.
Args:
namespace:
The namespace of the symbols the rate can use.
Returns:
rate:
The force, before the division by the mass.
"""
return namespace.F * cos(namespace.omega * namespace.t)
Equations
first.x[0]: diff(x_first_0(t), t, t) == (center_left_spring_0 - x_first_0(t))*k_left_spring/m_first + k_middle_spring*(length_middle_spring - abs(x_first_0(t) - x_second_0(t)))*(x_first_0(t) - x_second_0(t))/(m_first*abs(x_first_0(t) - x_second_0(t))) + F_first_force*cos(omega_first_force*t)/m_first - gamma_first_damper*diff(x_first_0(t), t)/m_first
second.x[0]: diff(x_second_0(t), t, t) == (center_right_spring_0 - x_second_0(t))*k_right_spring/m_second - k_middle_spring*(length_middle_spring - abs(x_first_0(t) - x_second_0(t)))*(x_first_0(t) - x_second_0(t))/(m_second*abs(x_first_0(t) - x_second_0(t))) + F_second_force*cos(omega_second_force*t)/m_second - gamma_second_damper*diff(x_second_0(t), t)/m_second