Simulation of Bouncing Balls | #2: Bouncing in Boxes [gnuplot]

YouTube

Simulation [gnuplot]

Model of a Bouncing Ball

Equations of motion $$\begin{eqnarray} \label{eq:eom2} \left\{ \begin{array}{l} \frac{dX}{dt}&=&V_X\\ \frac{dY}{dt}&=&V_Y\\ \frac{d{V}_{X}}{dt}&=&-g\sin\alpha\\ \frac{d{V}_{Y}}{dt}&=&-g\cos\alpha \end{array} \right. \end{eqnarray}$$ Inelastic collision $$\begin{eqnarray} V_{X}^{\rm{out}}&=-\color{red}{e}V_{X}^{\rm{in}}\\ V_{Y}^{\rm{out}}&=-\color{red}{e}V_{Y}^{\rm{in}} \end{eqnarray}$$ Rotation matrix $$\begin{equation} R(\alpha) = \left[ \begin{array}{cc} \cos\alpha & -\sin\alpha \\ \sin\alpha & \cos\alpha \end{array} \right] \end{equation}$$ Rotation using the matrix $$\begin{eqnarray} \left[\begin{array}{c} x \\ y\end{array}\right] &=& R(\alpha) \left[\begin{array}{c}X \\ Y \end{array}\right] &=& \left[ \begin{array}{c} X\cos\alpha - Y\sin\alpha \\ X\sin\alpha + Y\cos\alpha \end{array} \right]\\ \left[\begin{array}{c} v_x \\ v_y\end{array}\right] &=& R(\alpha) \left[\begin{array}{c}V_X \\ V_Y \end{array}\right] &=& \left[ \begin{array}{c} V_X\cos\alpha - V_Y\sin\alpha \\ V_X\sin\alpha + V_Y\cos\alpha \end{array} \right]\\ \end{eqnarray}$$

Case 1: $\alpha=20^\circ$

PLT file ①

# Setting ----------------------------------------
reset
set term gif animate delay 5 size 1280, 720
set output "BallBox.gif"
set nokey
unset grid
 
A = -250        # A<=x<=B
B = 250
C = -60            # C<=y<=D
D = 180
set xr[A:B]
set yr[C:D]
set xl "{/TimesNewRoman:Italic=24 x}"
set yl "{/TimesNewRoman:Italic=24 y}"
set tics front font 'Times New Roman,20'
set xtics 50
set ytics 50
set size ratio -1
 
# Parameter ----------------------------------------
m = 1.5             # mass [kg]
g = 9.8             # gravitational acceleration [m/s2]
r = 4.0             # Radius of the ball
V = 40              # velocity [m/s]
a = 20*pi/180       # Slope angle [rad] (-90<a<90)
b = 70*pi/180       # Projection angle [rad]
e = 0.8             # The coefficient of restitution [-]
ep = 2              # Epsilon
 
dt = 0.001          # Step [s]
dh = dt/6.0
cut = 200           # Decimation
dis = 800           # Start to disappear
 
N = 5               # The number of balls
cnt = 0             # The number of balls being framed out
array X[N]          # Position of balls
array Y[N]
array VX[N]         # Velocity of balls
array VY[N]
array x[N]          # Position of balls after a rotation
array y[N]
array vx[N]         # Velocity of balls after a rotation
array vy[N]
 
# Color of balls
array color[N] = ["royalblue", "red", "orange", "green", "black"]
                  
# Border line
Lu = 130
Ld = 0
Lr = 150
Ll = -150
 
# Goal
if(a>0){
    GOALx = Ll + r + ep
} else { if(a<0) {
    GOALx = Lr - r - ep
} else { if(a==0) {
    # Not define GOALx
}}}
 
GOALy = Ld + r + ep
 
#Rotation matrix R(a)
# R[1]=R11, R[2]=R12, R[3]=R21, R[4]=R22
array R[4] = [cos(a), -sin(a), sin(a), cos(a)]
 
 
# Functions ----------------------------------------
# Equations of Motion
f1(x, y, vx, vy) = vx          # dx/dt
f2(x, y, vx, vy) = vy          # dy/dt
f3(x, y, vx, vy) = -g*sin(a)   # dvx/dt
f4(x, y, vx, vy) = -g*cos(a)   # dvy/dt
 
# 4th order Runge-Kutta (Define rk_i(x, y, vx, vy))
do for[i=1:4]{
    rki = "rk"
    fi  = "f".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, fi, fi, fi, fi)
    eval rki
}
 
# X->x, Y->y (x=R(a)X)
x(a, b) = R[1]*a + R[2]*b
y(a, b) = R[3]*a + R[4]*b
 
# x->X, y->Y (X=R(-a)x)
X(a, b) = R[4]*a - R[2]*b
Y(a, b) = -R[3]*a + R[1]*b
 
# Time
Time(t) = sprintf("{/TimesNewRoman:Italic t} = %.1f", t)
 
 
# Plot ----------------------------------------
# Initial Value
t  = 0
b_inc = 36               # Increment for b
do for[j=1:N]{
    X[j]  = X(0, 60)
    Y[j]  = Y(0, 60)
    VX[j] = V*cos(b + (b_inc*(j-int(N/2))) * pi/180)
    VY[j] = V*sin(b + (b_inc*(j-int(N/2))) * pi/180)
    
    # Vector rotation (Rotation matrix R(a))
    x[j]  = x(X[j] , Y[j] )
    y[j]  = y(X[j] , Y[j] )
    vx[j] = x(VX[j], VY[j])
    vy[j] = y(VX[j], VY[j])
}
 
# Draw initiate state for 70 steps
do for [i=1:70] {
    # Time
    set label 1 Time(t) left at graph 0.05, 0.93 font 'TimesNewRoman, 25' front
    
    # Balls
    do for[j=1:N]{
        set obj j circ at x[j], y[j] size r fc rgb color[j] fs solid noborder
    }
    
    # Draw ground and balls
    if(a != 0){
        plot tan(a)*(x-x(0, Ld)) + y(0, Ld) with filledcurve x1 lc rgb 'gray50', \
             tan(a)*(x-x(0, Lu)) + y(0, Lu) with filledcurve x2 lc rgb 'gray50', \
             -1/tan(a)*(x-x(Ll, 0)) + y(Ll, 0) with filledcurve y1 lc rgb 'gray50', \
             -1/tan(a)*(x-x(Lr, 0)) + y(Lr, 0) with filledcurve y2 lc rgb 'gray50'
    } else {
        plot Ld with filledcurve x1 lc rgb 'gray50', \
             Lu with filledcurve x2 lc rgb 'gray50', \
             ((x<Ll)||(x>Lr) ? Lu : 1/0) with filledcurve y1=Ld lc rgb 'gray50'
    }
}
 
# Update until all of balls are framed out
i = 0               # The number of loops
while(1){
    i = i +1
    t = t + dt
    set label 1 Time(t)
 
    do for[j=1:N]{
        # 4th order Runge-Kutta
        temp_X  = X[j]  + rk1(X[j], Y[j], VX[j], VY[j])
        temp_Y = Y[j]  + rk2(X[j], Y[j], VX[j], VY[j])
        temp_VX = VX[j] + rk3(X[j], Y[j], VX[j], VY[j])
        temp_VY = VY[j] + rk4(X[j], Y[j], VX[j], VY[j])
        
        X[j]=temp_X; Y[j]=temp_Y; VX[j]=temp_VX; VY[j]=temp_VY
        
        # Judge whether balls bounce or not
        if(Y[j] > Lu-r){
            Y[j]  = Lu-r
            VY[j] = -e*VY[j]
        }
        if(Y[j] < Ld+r){
            Y[j]  = Ld+r
            VY[j] = -e*VY[j]        
        }
        if(X[j] > Lr-r){
            X[j]  = Lr-r
            VX[j] = -e*VX[j]
        }
        if(X[j] < Ll+r){
            X[j]  = Ll+r
            VX[j] = -e*VX[j]
        }
        
        # Vector rotation (Rotation matrix R(a))
        x[j]  = x(X[j] , Y[j] )
        y[j]  = y(X[j] , Y[j] )
        vx[j] = x(VX[j], VY[j])
        vy[j] = y(VX[j], VY[j])
 
        # Update objects
        set obj N*i+j circ at x[j], y[j] size r fc rgb color[j] fs solid noborder
        # Make old objects trajectory of the ball 
        set obj N*(i-1)+j at x[j], y[j] size 0.1
        
        # Start to disappear
        if(i>=dis){
            unset object 1+N*(i-dis)+j
        }
    }
    
    # Decimate and plot
    if(i%cut==0){
        replot
    }
    
    # Count balls meeting a condition
    do for[j=1:N]{
        if(a > 0){
            if(X[j]<GOALx && Y[j]<GOALy){
                cnt = cnt + 1   # whether a ball stands still or not
            }
        } else { if(a < 0){
            if(X[j]>GOALx && Y[j]<GOALy){
                cnt = cnt + 1   # whether a ball stands still or not
            }
        } else {
            if(Y[j]<GOALy){
                cnt = cnt + 1   # whether a ball is grounded or not
            }
        }}
    }
    
    # Exit from the loop when all of balls meet the condition
    if(cnt == N){
        break
    } else {
        cnt =0
    }
}
 
set out

GIF file ①

$V=40\ \mathrm{m/s},\ a=20^\circ,\ b=70^\circ,\ e=0.8,\ \varepsilon=0.1,\\L_u=1320,\ L_d=0,\ L_r=150,\ L_l=-150$

Case 2: $\alpha=45^\circ$

PLT file ②

A = -200
B = 200
C = -200
D = 200

V = 40
a = 45*pi/180
b = 50*pi/180
e = 0.8
ep = 2

Lu = 120
Ld = -120
Lr = 120
Ll = -120

GIF file ②

$V=40\ \mathrm{m/s},\ a=45^\circ,\ b=50^\circ,\ e=0.8,\ \varepsilon=2,\\L_u=120,\ L_d=-120,\ L_r=120,\ L_l=-120$

Case 3: $\alpha=0^\circ$

GIF file ④

$V=40\ \mathrm{m/s},\ a=0^\circ,\ b=85^\circ,\ e=0.8,\ \varepsilon=2,\\L_u=130,\ L_d=0,\ L_r=150,\ L_l=-150$

Extra

Case 4: $\alpha=5^\circ$

GIF file ③

$V=60\ \mathrm{m/s},\ a=5^\circ,\ b=85^\circ,\ e=0.5,\ \varepsilon=2,\\L_u=130,\ L_d=0,\ L_r=150,\ L_l=-150$