Vibration energy harvester
A dual-mass piezoelectric harvester in a car suspension, driven by the roughness of the road.
- Name:
vibration_energy_harvester - Field:
engineering/circuits - Formalism: ordinary differential equations
- References: Trzaska 2023
Factory
@register_model(
field='engineering/circuits',
tags=('energy harvesting', 'piezoelectric', 'vibration', 'vehicles', 'ode'),
references=('Trzaska 2023',),
defaults={'until': 10.0, 'dt': 0.01},
observables=lambda system: [ObservableSpec('voltage', _piezo_voltage)],
)
def vibration_energy_harvester(
*,
m1: float = 40.0,
m2: float = 400.0,
k1: float = 2.0e5,
k2: float = 2.0e4,
b: float = 1500.0,
A: float = 0.01,
omega: float = 9.42,
n: float = 10.0,
d: float = 4.0e-10,
C: float = 1.0e-8,
) -> System:
"""Return a dual-mass piezoelectric harvester in a car suspension, driven by the roughness of the road.
The unsprung and sprung masses follow `m1 y1'' + b (y1' - y2') + k2 (y1 - y2) + k1 (y1 - z(t)) = 0` and
`m2 y2'' + b (y2' - y1') + k2 (y2 - y1) = 0`, with the road profile `z(t) = A sin(ω t)`. A lever with arm ratio
`n` presses a piezoelectric bar of constant `d` and capacitance `C`, whose charge is `n d k2 (y2 - y1)` and whose
voltage is the observable `voltage`. The damper `b` stands for the energy the bar absorbs.
Args:
m1:
The unsprung mass.
m2:
The sprung mass.
k1:
The stiffness of the tyre.
k2:
The stiffness of the suspension spring.
b:
The damping coefficient of the harvester.
A:
The amplitude of the road profile.
omega:
The angular frequency of the road profile.
n:
The ratio of the arms of the lever.
d:
The piezoelectric constant of the bar.
C:
The capacitance of the bar.
Returns:
system:
The model, with the displacements `harvester.y1` and `harvester.y2`, and the observable `voltage`.
Examples:
>>> from skcomplex.models.engineering.circuits import vibration_energy_harvester
>>> vibration_energy_harvester(b=1000.0).elements[0].b
1000.0
"""
harvester = PiezoHarvesterElement('harvester', m1=m1, m2=m2, k1=k1, k2=k2, b=b, A=A, omega=omega, n=n, d=d, C=C)
suspension = SpringDamperInteraction(
'suspension',
target='harvester',
states=('y1', 'y2'),
stiffness='k2',
damping='b',
divisor={'y1': 'm1', 'y2': 'm2'},
)
return System(elements=[harvester], interactions=[suspension])
Components
class PiezoHarvesterElement(Element):
"""Two-mass suspension of a car whose spring presses a piezoelectric bar as the road shakes it.
The element holds the tyre between the unsprung mass and the road, and the model joins the two masses by a
`SpringDamperInteraction` of the stiffness `k2` and the damping `b`.
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
y1: Annotated[object, StateSpec(order=2)] = 0.0,
y1_d1: float = 0.0,
y2: Annotated[object, StateSpec(order=2)] = 0.0,
y2_d1: float = 0.0,
m1: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 40.0,
m2: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 400.0,
k1: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 200000.0,
k2: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 20000.0,
b: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1500.0,
A: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
omega: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 9.42,
n: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
d: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 4e-10,
C: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 1e-08,
) -> None:
super().__init__(label, visual=visual)
self.y1 = y1
self.y1_d1 = y1_d1
self.y2 = y2
self.y2_d1 = y2_d1
self.m1 = m1
self.m2 = m2
self.k1 = k1
self.k2 = k2
self.b = b
self.A = A
self.omega = omega
self.n = n
self.d = d
self.C = C
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the acceleration of the unsprung mass by the tyre on the road, the suspension apart."""
return [RateLaw(s.y1, -s.k1 * (s.y1 - s.A * sin(s.omega * s.t)) / s.m1)]
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
harvester.y1: diff(y1_harvester(t), t, t) == (A_harvester*sin(omega_harvester*t) - y1_harvester(t))*k1_harvester/m1_harvester - k2_harvester*(y1_harvester(t) - y2_harvester(t))/m1_harvester - b_harvester*(diff(y1_harvester(t), t) - diff(y2_harvester(t), t))/m1_harvester
harvester.y2: diff(y2_harvester(t), t, t) == k2_harvester*(y1_harvester(t) - y2_harvester(t))/m2_harvester + b_harvester*(diff(y1_harvester(t), t) - diff(y2_harvester(t), t))/m2_harvester