Wheel and rail contact
A rail segment on a sleeper, loaded by the wheels of passing bogies.
- Name:
wheel_rail_contact - Field:
engineering/power_systems - Formalism: ordinary differential equations
- References: Trzaska 2023
Factory
@register_model(
field='engineering/power_systems',
tags=('energy', 'vibration', 'railways', 'pulses', 'ode'),
references=('Trzaska 2023',),
defaults={'until': 60.0, 'dt': 0.05},
)
def wheel_rail_contact(
*,
F: float = 1000.0,
T: float = 10.0,
width: float = 1.0,
spacing: float = 3.0,
k2: float = 0.0,
Mr: float = 60.0,
Ms: float = 156.0,
) -> System:
"""Return a rail segment on a sleeper, loaded by the wheels of passing bogies.
The rail and the sleeper are lumped masses joined by the pad, a spring and a damper, and the sleeper rests on the
ballast, `Mr z'' + br (z' - zs') + kr (z - zs) = F(t)` and
`Ms zs'' + br (zs' - z') + bs zs' + kr (zs - z) + ks zs + k2 (alpha zs^2 + β zs^3) = 0`. The load is two rectangular
pulses of force `F` per period `T`, one per wheel of a bogie, `spacing` apart. The nonlinear ballast term, with
`alpha = 1.75` and `β = -1`, is off for `k2 = 0`.
Args:
F:
The force of a wheel.
T:
The period of the passing bogies.
width:
The duration of the load of one wheel.
spacing:
The time between the two wheels of a bogie.
k2:
The stiffness of the nonlinear ballast term.
Mr:
The mass of the rail segment.
Ms:
The mass of the sleeper.
Returns:
system:
The model, with the deflections `rail.z` of the rail and `rail.zs` of the sleeper, and their velocities.
Examples:
>>> from skcomplex.models.engineering.power_systems import wheel_rail_contact
>>> wheel_rail_contact(k2=500.0).elements[0].k2
500.0
"""
rail = RailSegmentElement('rail', F=F, T=T, width=width, spacing=spacing, k2=k2, Mr=Mr, Ms=Ms)
pad = SpringDamperInteraction(
'pad',
target='rail',
states=('z', 'zs'),
stiffness='kr',
damping='br',
divisor={'z': 'Mr', 'zs': 'Ms'},
)
return System(elements=[rail], interactions=[pad])
Components
class RailSegmentElement(Element):
"""Rail segment on a sleeper, two masses joined by a spring and a damper and loaded by passing wheels.
The element holds the wheel load and the ballast under the sleeper, and the model joins the rail and the sleeper by
a `SpringDamperInteraction` of the stiffness `kr` and the damping `br`.
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
z: Annotated[object, StateSpec(order=2)] = 0.0,
z_d1: float = 0.0,
zs: Annotated[object, StateSpec(order=2)] = 0.0,
zs_d1: float = 0.0,
Mr: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 60.0,
Ms: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 156.0,
kr: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 6000.0,
br: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 300.0,
ks: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1500.0,
bs: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 520.0,
k2: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
alpha: Annotated[object, ParameterSpec()] = 1.75,
beta: Annotated[object, ParameterSpec()] = -1.0,
F: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1000.0,
T: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
width: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
spacing: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 3.0,
) -> None:
super().__init__(label, visual=visual)
self.z = z
self.z_d1 = z_d1
self.zs = zs
self.zs_d1 = zs_d1
self.Mr = Mr
self.Ms = Ms
self.kr = kr
self.br = br
self.ks = ks
self.bs = bs
self.k2 = k2
self.alpha = alpha
self.beta = beta
self.F = F
self.T = T
self.width = width
self.spacing = spacing
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the accelerations of the rail by the wheels and of the sleeper by the ballast, the pad apart."""
return [
RateLaw(s.z, RailSegmentElement._wheel_load(s) / s.Mr),
RateLaw(
s.zs,
(-s.bs * derivative(s.zs, 1) - s.ks * s.zs - s.k2 * (s.alpha * s.zs**2 + s.beta * s.zs**3)) / s.Ms,
),
]
@staticmethod
def _wheel_load(s: SimpleNamespace) -> Expression:
"""Return the force of the two wheels of a bogie: two pulses per period, `spacing` apart."""
return pulse_train(s.t, s.T, s.width, s.F) + pulse_train(s.t, s.T, s.width, s.F, start=s.spacing)
class SpringDamperInteraction(GradientInteraction):
"""Linear spring and damper between two coordinates `x1` and `x2` of order 2 of one element.
The spring is the gradient of the energy `k (x1 - x2)^2 / 2`, and the damper adds `-b (x1' - x2')` to the first
coordinate and its opposite to the second, both divided by the divisor of every coordinate, such as its own mass.
The stiffness and the damping are named parameters of the element.
Args:
label:
The label of the interaction.
target:
The path of the element with the two coordinates.
states:
The names of the two coordinates.
stiffness:
The name of the parameter of the stiffness `k`.
damping:
The name of the parameter of the damping `b`.
divisor:
The divisor of both coordinates, or a mapping from every coordinate to its own divisor.
Raises:
ValidationError:
The states are not two names, the stiffness or the damping is not a name, or a divisor is invalid.
Examples:
>>> from skcomplex.models.engineering.common.mechanics import SpringDamperInteraction
>>> SpringDamperInteraction('suspension', target='car', states=('y1', 'y2')).stiffness
'k'
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
states: Sequence[str] = ('x1', 'x2'),
stiffness: str = 'k',
damping: str = 'b',
divisor: str | Callable[[SimpleNamespace], object] | Mapping[str, str] | None = None,
) -> None:
super().__init__(label, target=target, states=states, divisor=divisor)
self.stiffness = stiffness
self.damping = damping
def _check(self) -> None:
"""Raise an error for states that are not two names, or a stiffness or damping that is not a name.
Raises:
ValidationError:
The states are not two names, the stiffness or the damping is not a name, or a divisor is invalid.
"""
super()._check()
if len(self._get_names()) != N_STATES:
error_msg = f'{self.label!r} joins the states {self.states!r}. Give two states.'
raise ValidationError(error_msg)
for name in (self.stiffness, self.damping):
if not (isinstance(name, str) and name.isidentifier()):
error_msg = f'{self.label!r} reads the parameter {name!r}. Give the name of a parameter.'
raise ValidationError(error_msg)
def _get_potential(self, namespace: SimpleNamespace) -> Expression:
"""Return the energy `k (x1 - x2)^2 / 2` of the spring.
Args:
namespace:
The namespace of the symbols the energy can use.
Returns:
energy:
The energy of the spring.
"""
first, second = (getattr(namespace.target, name) for name in self._get_names())
return getattr(namespace.target, self.stiffness) * (first - second) ** 2 / 2
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rates of the spring and of the damper of both coordinates.
Args:
namespace:
The namespace of the symbols the laws can use.
Returns:
laws:
The rate laws of the two coordinates.
"""
laws = super()._get_laws(namespace)
member = namespace.target
names = self._get_names()
first, second = (getattr(member, name) for name in names)
force = -getattr(member, self.damping) * (derivative(first, 1) - derivative(second, 1))
for name, state, sign in ((names[0], first, 1), (names[1], second, -1)):
divisor = self._get_divisor(member, name)
laws.append(RateLaw(state, sign * force if divisor is None else sign * force / divisor))
return laws
Equations
rail.z: diff(z_rail(t), t, t) == -kr_rail*(z_rail(t) - zs_rail(t))/Mr_rail - br_rail*(diff(z_rail(t), t) - diff(zs_rail(t), t))/Mr_rail + (where(both(t >= 0.0, mod(t, T_rail) < width_rail), F_rail, 0) + where(both(t >= spacing_rail, mod(-spacing_rail + t, T_rail) < width_rail), F_rail, 0))/Mr_rail
rail.zs: diff(zs_rail(t), t, t) == kr_rail*(z_rail(t) - zs_rail(t))/Ms_rail + br_rail*(diff(z_rail(t), t) - diff(zs_rail(t), t))/Ms_rail - ((beta_rail*zs_rail(t)^3 + alpha_rail*zs_rail(t)^2)*k2_rail + ks_rail*zs_rail(t) + bs_rail*diff(zs_rail(t), t))/Ms_rail