📂Dynamics

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 ¹

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 ProgressMeter
packages = [:DataFrames, :CSV, :Plots, :LaTeXStrings]
@time for package in 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[]
@showprogress for dr in 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")