📂動力学

力学系としての振動衝撃モデル

概要

振動衝撃モデル^{vibro-impact model}は、振動する円筒形カプセル内の物体の動きを表す非スムーズなダイナミックシステムで、主にハード^hardとソフト^softの2つのタイプに分けられます。

モデル ¹

円筒形カプセル^{cylindrical Capsule}は便宜上カプセルと呼びます。以下に登場するイラストでは、$s$はカプセルの長さ、$M$はカプセルの質量、$m_{b}$は物体の質量で、$M \gg m_{b}$であるため、物体の動きがカプセルの動きに影響を与えないと仮定します。

カプセルは$F_{\text{ext}}$の方向に時間単位$\tau$に基づいて振動数$\omega$の調和振動を行うと仮定します。これにより、衝突がない場合のカプセルの変位^displacement $Z$は振動の振幅^amplitude $A_{0}$と$\tau$に対する微分$'$について、以下の微分方程式で表すことができます。 $$M Z '' = - A_{0} \cos \omega \tau$$ $g$は重力加速度、$\beta$はカプセルの傾きです。これにより、衝突がない場合の物体の変位$Y$は $$Y '' = - g \sin \beta$$ です。$F_{\text{ext}}$方向に見たときの物体の座標を$Y$、カプセルの座標を$Z$とすると、それらの相対変位$X := Y - Z$は$A := A_{0} / M$について、以下のように表すことができます。 $$X '' = - g \sin \beta + A \cos \omega \tau$$

ハードインパクト

ハードインパクトモデル^{hard Impact model}は、物体がカプセルの端に触れたときに方向が反転する方法で衝撃を反映します。 $$\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*}$$ ここで$d\tau \approx 0$は非常に短い時間間隔を表します。物理的な意味を持たないように抽象的なシステムに変換すると、$t$に対する微分$\dot{}$について、以下のようになります。 $$\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 model}は、物体がスプリングを付けた板にぶつかって方向が変わることを強制調和振動を通じて反映します。スプリングに触れていないときの相対変位$U := Y - Z$とスプリングに触れているときの相対変位$U_{c}$が別々に分けられます。 $$\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*}$$ 物理的な意味を持たないように抽象的なシステムに変換すると、$u$と$u_{c}$を別々に区別せずに、$t$に対する微分$\dot{}$について、以下のようになります。 $$\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}$$ ここで$\sign$は符号関数です。

変数とパラメータ

基本的にモデルのパラメータは以下のように与えられます：

• $m_{b} = 0.0035$: 物体の質量です。
• $A = 1.0$: 振動の強さに関連する値です。
• $\gamma = g \sin \beta = 0$: 重力によって垂直方向に作用する値です。
  • $\beta = 0.0$: カプセルの傾きです。
  • $g = 9.8$: 重力加速度です。
• $s = 0.3$: カプセルの長さです。
  • $d = 0.3$: カプセルの長さが変換された値です。
• $\omega = \pi$: 振動数です。
  • $\omega_{0}$: システムの固有振動数です。
• $r = 0.5$: ハードインパクトモデルで、衝突後の速度に関連する値で、衝突直後に速度は元の$r$倍に変わります。$r = 0.5$の場合、速度が半分になります。
• $k = 560.0$: ソフトインパクトモデルで、弾性係数です。
  • $\kappa = 400.0$: 弾性係数に関連する値です。
  • $\mu = 172.363$: ソフトインパクトモデルで、減衰力に関連する値です。

モデルが非次元化^{non-dimensional</sup}

されるときの変数の変換は以下のようになります。 $$\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*}$$ パラメータの変換は以下のようになります。 $$\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*}$$

説明

トラジェクトリ

ソフトインパクトモデルで$d = 0.15$の場合のアニメーションです。$x$軸は相対変位$u$、$y$軸は速度$v = \dot{u}$です。$\pm d /2 = \pm 0.075$に達すると、強い衝撃とともに速度が変わります：

通常、私たちが関心を持つパラメータはカプセルの長さに関連する$d$です。ソフトインパクトモデルで$d$が変化することによってシステムがどのように変化するか見てみましょう。右側の図の赤い実線はカプセルの境界に対応し、よく見ると物体の変位がわずかに境界の外に出ていることがわかります。

• $d = 0.30$の場合：
• $d = 0.20$の場合：
• $d = 0.16$の場合：
• $d = 0.10$の場合：

バイフケーション

パラメータ$d$を$0.3$から$0.1$まで減らしながらバイフケーションダイアグラムを描くと、上のようなバイフケーションが見られます。$d \approx 0.12$以下で特徴的なウィンドウが顕著です。

コード

以下は、このポストで使用された図を再現するジュリアのコードです。

バイフケーションダイアグラム

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")

アニメーション

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")