Simulation of Bouncing Balls | #4: Bouncing on Explicit Functions (2) [gnuplot]

Sunday, April 7, 2019

gnuplot Mechanics Simulation of Bouncing Balls YouTube

t f B! P L

YouTube

Simulation

Model of a Bouncing Ball


Impulse-momentum theorem \begin{equation} \label{eq:im4} m\boldsymbol{v}'=m\boldsymbol{v}+\boldsymbol{I}=m\boldsymbol{v}+I\boldsymbol{n} \end{equation} Collision \begin{eqnarray} \left|\boldsymbol{v}'\cdot\boldsymbol{n}\right|=e\left|\boldsymbol{v}\cdot\boldsymbol{n}\right| \end{eqnarray} Calculating I \begin{eqnarray} I &=& \left| \left(m\boldsymbol{v}'\right) \cdot \boldsymbol{n} \right| + \left| \left(m\boldsymbol{v}\right) \cdot \boldsymbol{n} \right| \nonumber\\ &=& m\left(\left|\boldsymbol{v}' \cdot \boldsymbol{n} \right| + \left| \boldsymbol{v} \cdot \boldsymbol{n} \right|\right) \nonumber\\ &=& m\left(e\left|\boldsymbol{v}\cdot\boldsymbol{n}\right|+\left|\boldsymbol{v}\cdot\boldsymbol{n}\right|\right) \nonumber\\ &=& -m\left(e+1\right) \left( \boldsymbol{v} \cdot \boldsymbol{n} \right) \label{eq:I} \end{eqnarray} Velocity after the collision (inserting (\ref{eq:I}) into (\ref{eq:im4})) \begin{eqnarray} \boldsymbol{v}'=\boldsymbol{v}-(e+1)\left(\boldsymbol{v}\cdot\boldsymbol{n}\right)\boldsymbol{n} \end{eqnarray}

Wall function \begin{equation} W(x, y)=y-f(x) \end{equation} Inequality constraint \begin{equation} W(x, y)\geq 0 \end{equation} Unit normal vector \begin{equation} \boldsymbol{n}=\frac{\nabla W}{\left|\nabla W\right|} \end{equation} Collision detection \begin{equation} \boldsymbol{v}'=\boldsymbol{v}-\left(e+1\right)\left(\boldsymbol{v}\cdot \frac{\nabla W}{\left|\nabla W\right|}\right)\frac{\nabla W}{\left|\nabla W\right|} \end{equation}

Case 1

\begin{eqnarray*} g(x) &=& x\left(x^3-28x^2-420\right)\left(10+\sin 8x\right)\\ f(x) &=& \frac{1}{10^5}\times g\left(\frac{5x}{7}+8\right)-5 \end{eqnarray*}

PLT file

# Setting --------------------
reset
set term gif animate delay 5 size 1280, 720
set output"fxy03.gif"

set nokey
set sample 10000
set xr[-35:35]
set yr[-15:15]
set xl "{/TimesNewRoman:Italic=24 x}"
set yl "{/TimesNewRoman:Italic=24 y}"
set tics font 'Times New Roman,20'
set xtics 5
set ytics 5

set size ratio -1       # ratio = yrange / xrange 
set grid

N = 5
array x[N]
array y[N]
array vx[N]
array vy[N]
array c[N] = ["royalblue", "red", "orange", "green", "black"]

# Parameter --------------------
r     = 0.3             # Radius of the ball
g     = 9.8             # Gravitational acceleration
dis   = 500             # Start to disappear
e     = 0.95            # Elasticity coefficient
dt    = 0.001           # Time step
dh    = dt/6
cut   = 80              # Decimation
ep = 1e-4               # for collision detection
time1 = 2               # Stop time4
lim1  = time1/dt/cut

time2 = 40              # Time limit40
lim2  = time2/dt

# Functions --------------------
# Wall
g(x) = x**2
f(x) = g(x/7.)-13.
W(x, y) = y-f(x)
 
# Partial derivative of f(x)
Wx(x) = -(f(x-2*dt)-8*f(x-dt)+8*f(x+dt)-f(x+2*dt)) / (12*dt) # 4th order central difference formula
Wy(y) = 1
                    
# n = (fx, fy) / sqrt(fx^2+fy^2)
d(x, y)  = sqrt(x**2+y**2)
nx(x, y) = Wx(x)/d(Wx(x), Wy(y))
ny(x, y) = Wy(y)/d(Wx(x), Wy(y))

# Inner product of v and n
vn(x,y,vx,vy) = vx*nx(x,y)+vy*ny(x,y)
I(x,y,vx,vy) = -(1+e)*vn(x,y,vx,vy)

# Equations of Motion
rk1(x,y,vx,vy) = vx         # dx/dt
rk2(x,y,vx,vy) = vy         # dy/dt
rk3(x,y,vx,vy) = 0          # dvx/dt
rk4(x,y,vx,vy) = -g         # dvy/dt

# 4th order Runge-Kutta (Define RK_i(x, y, vx, vy))
do for[i=1:4]{
    RKi = "RK"
    rki  = "rk".sprintf("%d", i)
    RKi = RKi.sprintf("%d(x, y, vx, vy) = (\
        k1 = %s(x, y, vx, vy),\
        k2 = %s(x + dt*k1/2., y + dt*k1/2., vx + dt*k1/2., vy + dt*k1/2.),\
        k3 = %s(x + dt*k2/2., y + dt*k2/2., vx + dt*k2/2., vy + dt*k2/2.),\
        k4 = %s(x + dt*k3, y + dt*k3, vx + dt*k3, vy + dt*k3),\
        dh * (k1 + 2*k2 + 2*k3 + k4))", i, rki, rki, rki, rki)
    eval RKi
}

# Label of displaying parameters
# Elastic
e(e) = sprintf("{/TimesNewRoman:Italic e} = %.2f", e)
# Time
Time(t) = sprintf("{/TimesNewRoman:Italic t} = %.1f s", t)


# Plot --------------------
# Initiate value
t = 0.0                             # Time

do for[j=1:N]{
    x[j]  = 13.*(j-3)-2             # Position
    y[j]  = 10.
    vx[j] = 0.                      # Velocity
    vy[j] = 0. 
}
    
# Draw initial state for lim1 steps
do for[i=1:lim1]{
    # Time
    set label 1 Time(t) font 'Times New Roman, 24' at graph 0.03, 0.03
    set label 2 e(e) font 'Times New Roman, 24' at graph 0.023, 0.08
    
    # Ball
    do for[j=1:N]{
        set object j+1 circle at x[j], y[j] size r fc rgb c[j] fs solid front
    }
    
    plot f(x) lc rgb 'black'
}

# Update
do for[i=1:lim2]{
    # Time
    t = t + dt
    set label 1 Time(t)
    
    # Balls
    do for[j=1:N]{
        # 4th order Runge-Kutta
        tmp_x  = x[j]  + RK1(x[j], y[j], vx[j], vy[j])
        tmp_y  = y[j]  + RK2(x[j], y[j], vx[j], vy[j])
        tmp_vx = vx[j] + RK3(x[j], y[j], vx[j], vy[j])
        tmp_vy = vy[j] + RK4(x[j], y[j], vx[j], vy[j])
        
        # Collision detection
        if(W(x[j], y[j]) >= 0 && W(tmp_x, tmp_y) < 0){
            # vec = x_[i+1] - x_[i]
            vec_x = tmp_x - x[j]
            vec_y = tmp_y - y[j]
            nor_v = d(vec_x, vec_y)
            
            # Find collision point
            while(W(x[j], y[j]) >= 0){
                x[j] = x[j] + vec_x / nor_v * ep
                y[j] = y[j] + vec_y / nor_v * ep
            }
            
            x[j] = x[j] - vec_x / nor_v * ep
            y[j] = y[j] - vec_y / nor_v * ep
            
            # Calculate bounce direction
            size = d(tmp_x-x[j], tmp_y-y[j])
            tmp_x = x[j] - (1+e) * size * vn(x[j],y[j],vec_x,vec_y) / nor_v * nx(x[j], y[j])
            tmp_y = y[j] - (1+e) * size * vn(x[j],y[j],vec_x,vec_y) / nor_v * ny(x[j], y[j])
            
            # Collision
            tmp_vx = vx[j] + I(x[j], y[j], vx[j], vy[j]) * nx(x[j], y[j])
            tmp_vy = vy[j] + I(x[j], y[j], vx[j], vy[j]) * ny(x[j], y[j])
        }
        
        # Update position and velocity of a ball
        x[j]  = tmp_x
        y[j]  = tmp_y
        vx[j] = tmp_vx
        vy[j] = tmp_vy
        
        # Draw a ball
        set object N*i+j circ at x[j], y[j] size r fc rgb c[j] fs solid front
        set object N*(i-1)+j circ size r*0.01    # Old ball turns smaller

        # Start to disappear
        if(i>=dis){
            unset object N*(i-dis)+j
        }
    }
    
    # Decimate and draw
    if(i%cut==0){
        replot
    }
}

set out 

GIF file ①

$e=0.95$

Case 2

\[ f(x)=5\left[1+\sin\left(\frac{2x}{3}\right)\right]\cos\left(\frac{2x}{3}\right) \]

GIF file ②

$e=1.00$

Case 3

\begin{eqnarray*} g(x) &=& \left(\sqrt{\frac{\left|1-x\right|+1-x}{2}}+\frac{1}{4}\right)e^{-\left(1-x\right)^2}\\ & &+\frac{1}{40}\left[\left|1-2^4\left(5x-3\right)^4\right|+\left|1-2^9\left(5x-3\right)^4\right|+2-528\left(5x-3\right)^4\right]\\ f(x) &=& -13\times g\left(\frac{x}{13}\right) \end{eqnarray*}

GIF file ③

"Oppai Function" by Hanao (はなお) [1]
$e=0.95$

References

  1. 東大参戦!大学対抗おっぱい関数選手権ついにトップ10大学発表! (はなお) [YouTube in Japanese].

External links

▶︎ Hanao's YouTube channel https://www.youtube.com/channel/UCPyNsNSTUtywkekbDdCA_8Q
▶︎ Hanao's Twitter https://twitter.com/hanao87_0

Search This Blog

Translate

QooQ