📂동역학

동역학계로써의 진동 충격 모델

개요

진동 충격 모델^{vibro-impact model}은 진동하는 원통형 관 속에 있는 물체의 움직임을 나타내는 넌스무스한 다이내믹 시스템으로, 크게 하드와 소프트^{hard, Soft} 두 가지 타입으로 나눌 수 있다.

모델 ¹

원통형 관^{cylindric0al 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}하게 바뀔 때 변수의 변환은 다음과 같다. $$ \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 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")

애니메이션

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