Vibro-impact Model as a Dynamical System
Overview
The Vibro-impact Model represents the motion of an object inside a vibrating cylindrical capsule in a nonsmooth dynamical system, which can broadly be classified into two types: Hard and Soft.
Model 1
Let’s refer to the Cylindrical Capsule for convenience as a capsule. In the illustrations below, $s$ is the length of the capsule, $M$ is the mass of the capsule, $m_{b}$ is the mass of the object, and it is assumed that the object’s motion, represented by $M \gg m_{b}$, does not affect the capsule’s movement.
The capsule is assumed to perform harmonic oscillations in the direction of $F_{\text{ext}}$ with a frequency of $\omega$ per time unit $\tau$. Therefore, in the absence of collision, the capsule’s displacement $Z$, given the amplitude of vibration $A_{0}$ and the derivative with respect to $\tau$, $'$, can be represented by the following differential equation: $$ M Z '' = - A_{0} \cos \omega \tau $$ $g$ is the gravitational acceleration, and $\beta$ is the incline of the capsule. Consequently, in the absence of collision, the displacement of the object $Y$ can be expressed as: $$ Y '' = - g \sin \beta $$ Looking in the direction of $F_{\text{ext}}$, when the coordinate of the object is $Y$ and the coordinate of the capsule is $Z$, their relative displacement $X := Y - Z$ with respect to $A := A_{0} / M$ can be expressed as: $$ X '' = - g \sin \beta + A \cos \omega \tau $$
Hard Impact
The Hard Impact Model reflects the impact by reversing the direction when the object touches the bounds of the capsule. $$ \begin{align*} X '' =& - g \sin \beta + A \cos \omega \tau & \text{if } \left| x \right| < {{ s } \over { 2 }} \\ X ' \left( \tau \right) =& - r X ' \left( \tau - d\tau \right) & \text{if } \left| x \right| = {{ s } \over { 2 }} \end{align*} $$ Here, $d\tau \approx 0$ represents a very short time interval. When abstracted to a non-physical system, it can be expressed in terms of the derivative of $t$ with respect to $\dot{}$ as: $$ \begin{align*} \ddot{x} =& - \gamma + \cos \pi t & \text{if } \left| x \right| < {{ d } \over { 2 }} \\ \dot{x} \left( t \right) =& - r \dot{x} \left( t - dt \right) & \text{if } \left| x \right| = {{ d } \over { 2 }} \end{align*} $$
Soft Impact
The Soft Impact Model simulates the impact through a forced harmonic oscillation, where the direction changes as the object hits a spring-loaded plate. The relative displacement when not in contact with the spring $U := Y - Z$ and when in contact $U_{c}$ are considered separately. $$ \begin{align*} U '' =& - g \sin \beta + A \cos \omega \tau & , \text{for } \left| U \right| < {{ s } \over { 2 }} \\ U_{c} '' =& - g \sin \beta + A \cos \omega \tau - {{ k } \over { m_{b} }} \sign \left( U_{c} \right) \left( \left| U_{c} \right| - {{ s } \over { 2 }} \right) - {{ c } \over { m_{b} }} U_{c} ' & , \text{for } \left| U_{c} \right| \ge {{ s } \over { 2 }} \end{align*} $$ Abstracting to a non-physical system without distinguishing between $u$ and $u_{c}$, it can be expressed in terms of the derivative of $t$ with respect to $\dot{}$ as: $$ \ddot{u} = - \gamma + \cos \left( \pi t \right) + \begin{cases} 0 & , \text{for } \left| u \right| < {{ d } \over { 2 }} \\ - \kappa^{2} \sign (u) \left( \left| u \right| - {{ d } \over { 2 }} \right) - \mu \dot{u} & , \text{for } \left| u \right| \ge {{ d } \over { 2 }} \end{cases} $$ Here, $\sign$ represents the sign function.
Variables and Parameters
The model’s parameters are generally given as follows:
- $m_{b} = 0.0035$: The mass of the object.
- $A = 1.0$: A value related to the amplitude of vibration.
- $\gamma = g \sin \beta = 0$: A value acting vertically due to gravity.
- $\beta = 0.0$: The incline of the capsule.
- $g = 9.8$: The gravitational acceleration.
- $s = 0.3$: The length of the capsule.
- $d = 0.3$: A transformed value of the capsule’s length.
- $\omega = \pi$: The frequency of vibration.
- $\omega_{0}$: The system’s natural frequency.
- $r = 0.5$: In the Hard Impact Model, a value related to the velocity after collision, where the velocity changes to $r$ times its original right after the collision. At $r = 0.5$, the velocity is halved.
- $k = 560.0$: The elastic coefficient in the Soft Impact Model.
- $\kappa = 400.0$: A value related to the elastic coefficient.
- $\mu = 172.363$: A value related to the damping force in the Soft Impact Model.
When the model is transformed to non-dimensional, the transformation of variables is as follows: $$ \begin{align*} t :=& {{ \tau \omega } \over { \pi }} \\ x(t) :=& {{ X (\tau) \omega^{2} } \over { A \pi^{2} }} \\ u(t) :=& {{ U (\tau) \omega^{2} } \over { A \pi^{2} }} \end{align*} $$ The transformation of parameters is as follows: $$ \begin{align*} d :=& {{ s \omega^{2} } \over { A \pi^{2} }} \\ \gamma :=& {{ g \sin \beta } \over { A }} \\ \kappa :=& {{ \omega_{0} \pi } \over { \omega }} \\ \mu :=& {{ c \pi } \over { m \omega }} \end{align*} $$
Explanation
Trajectories
Animation for the Soft Impact Model when $d = 0.15$. The $x$-axis represents the relative displacement $u$, and the $y$-axis represents the velocity $v = \dot{u}$. When reaching $\pm d /2 = \pm 0.075$, the velocity changes abruptly along with a strong impact:
Typically, the parameter of interest is $d$ related to the capsule’s length. Let’s see how the system changes as $d$ changes in the Soft Impact Model. The red solid line in the right figure corresponds to the capsule’s boundary, and it’s noticeable that the object’s displacement slightly exceeds the boundary.
- At $d = 0.30$:
- At $d = 0.20$:
- At $d = 0.16$:
- At $d = 0.10$:
Bifurcation
Drawing a bifurcation diagram while decreasing the parameter $d$ from $0.3$ to $0.1$ reveals the bifurcation shown above. Characteristic windows become prominent below $d \approx 0.12$.
Code
The following is Julia code used to reproduce the images in this post.
Bifurcation Diagram
const κ = 400.0
const μ = 172.363
using ProgressBars
packages = [:DataFrames, :CSV, :Plots, :LaTeXStrings]
@time for package in ProgressBar(packages)
@eval using $(package)
end
println(join(packages, ", "), " loaded!")
function RK4(f::Function,v::AbstractVector, h=10^(-2))
V1 = f(v)
V2 = f(v + (h/2)*V1)
V3 = f(v + (h/2)*V2)
V4 = f(v + h*V3)
return v + (h/6)*(V1 + 2V2 + 2V3 + V4), V1
end
function factory_soft(idx::Int64, d::Number)
d2 = d/2
function soft(tuv)
t, u, v = tuv
impact = ifelse(abs(u) < d2, 0, -(κ^2)*sign(u)*(abs(u)-d2) - μ*v)
ṫ = 1
u̇ = v
v̇ = cospi(t) + impact
return [ṫ, u̇, v̇]
end
dt = 10^(-5); tend = 50
t_ = 0:dt:tend
ndatapoints = round(Int64, tend/(10dt))
len_t_ = length(t_)
traj = zeros(6, len_t_+1)
x = [0.0, 0.05853, 0.47898]
dx = soft(x)
traj[1:3, 1] = x
for t in 1:length(t_)
x, dx = RK4(soft, x, dt)
if t ≥ ndatapoints
traj[4:6, t] = dx
traj[1:3, t+1] = x
end
end
traj = traj[:, (end-ndatapoints):(end-1)]'
data = DataFrame(traj,
["t", "u", "v", "dt", "du", "dv"])
return data
end
d_range = 0.1:0.0001:0.3
schedule = DataFrame(idx = eachindex(d_range), d = d_range)
xdots = Float64[]; ydots = Float64[]
for dr in ProgressBar(eachrow(schedule))
data = factory_soft(dr.idx, dr.d)
idx = [false; diff(abs.(data.u) .> (dr.d/2)) .> 0]
sampled = data.v[idx]
append!(xdots, fill(dr.d, length(sampled)))
append!(ydots, sampled)
end
@time a1 = scatter(xdots, ydots,
xlabel = L"d", ylabel = "Impact velocity",
label = :none, msw = 0, color = :black, ms = 0.5, alpha = 0.5, size = (700, 300))
png(a1, "soft_bifurcation")
Animation
Plots.gr()
anim = @animate for tk in ProgressBar(1:100:nrow(_data))
plot(_data.u[1:10:tk], _data.v[1:10:tk],
xlabel = "position", ylabel = "velocity", legend = :none,
xlims = (-0.2, 0.2), ylims = (-0.6, 0.6), alpha = (1:10:tk)/tk)
end
mp4(anim, "soft.mp4")
Dimitri Costa, Rachel Kuske, Daniil Yurchenko; Qualitative changes in bifurcation structure for soft vs hard impact models of a vibro-impact energy harvester. Chaos 1 October 2022; 32 (10): 103120. https://doi.org/10.1063/5.0101050 ↩︎