Processing math: 100%
Home / YouTube /

Simulation of Bouncing Balls | #1: Bouncing on Slopes [gnuplot]

Saturday, February 16, 2019

gnuplot Mechanics Simulation of Bouncing Balls YouTube

t f B! P L

Contents[hide]

YouTube

Simulation [gnuplot]

Model of a Bouncing Ball


Equations of motion \begin{eqnarray} \label{eq:eom1} \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{equation} V_{Y}^{\rm{out}}=-\color{red}{e}V_{Y}^{\rm{in}} \end{equation}
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}

BallBouncing.plt

  1. # Setting ----------------------------------------
  2. reset
  3. set term gif animate delay 5 size 1280, 720
  4. set output "BallBouncing.gif"
  5. set nokey
  6. set grid
  7.  
  8. L = 250            # -L<=x<=L
  9. A = -40            #  A<=y<=B
  10. B = 160
  11. set xr[-L:L]
  12. set yr[A:B]
  13. set xl "{/TimesNewRoman:Italic=24 x}"
  14. set yl "{/TimesNewRoman:Italic=24 y}"
  15. set tics font 'Times New Roman,20'
  16. set xtics 50
  17. set ytics 40
  18. set size ratio -1
  19.  
  20. # Parameter ----------------------------------------
  21. g = 9.8                # gravitational acceleration [m/s2]
  22. V = 30              # velocity [m/s]
  23. r = 4.0             # Radius of the ball
  24. a = 5*pi/180        # Slope angle [rad]
  25. b = -85*pi/180      # Projection angle [rad]
  26. dt = 0.001          # Step [s]
  27. dh = dt/6.0
  28. cut = 150           # Decimation
  29. dis = 800           # Start to disappear
  30.  
  31. N = 5               # The number of balls
  32. array e[N]          # The coefficient of restitution
  33. array X[N]          # Position of balls
  34. array Y[N]
  35. array VX[N]         # Velocity of balls
  36. array VY[N]
  37. array x[N]          # Position of balls after a rotation
  38. array y[N]
  39. array vx[N]         # Velocity of balls after a rotation
  40. array vy[N]
  41. array color[N] = ["royalblue", "red", "orange", "green", "black"]
  42. cnt = 0             # The number of balls being framed out
  43.  
  44. # Rotation matrix R(a)
  45. # R[1]=R11, R[2]=R12, R[3]=R21, R[4]=R22
  46. array R[4] = [cos(a), -sin(a), sin(a), cos(a)]
  47.  
  48. # Functions ----------------------------------------
  49. # Equations of Motion
  50. f1(x, y, vx, vy) = vx          # dx/dt
  51. f2(x, y, vx, vy) = vy          # dy/dt
  52. f3(x, y, vx, vy) = -g*sin(a)   # dvx/dt
  53. f4(x, y, vx, vy) = -g*cos(a)   # dvy/dt
  54.  
  55. # 4th order Runge-Kutta (Define rk_i(x, y, vx, vy))
  56. do for[i=1:4]{
  57.     rki = "rk"
  58.     fi  = "f".sprintf("%d", i)
  59.     rki = rki.sprintf("%d(x, y, vx, vy) = (\
  60.         k1 = %s(x, y, vx, vy),\
  61.         k2 = %s(x + dt*k1/2., y + dt*k1/2., vx + dt*k1/2., vy + dt*k1/2.),\
  62.         k3 = %s(x + dt*k2/2., y + dt*k2/2., vx + dt*k2/2., vy + dt*k2/2.),\
  63.         k4 = %s(x + dt*k3, y + dt*k3, vx + dt*k3, vy + dt*k3),\
  64.         dh * (k1 + 2*k2 + 2*k3 + k4))", i, fi, fi, fi, fi)
  65.     eval rki
  66. }
  67.  
  68. # (X,Y) ->(x,y) (x=R(a)X)
  69. x(X, Y) = R[1]*X + R[2]*Y
  70. y(X, Y) = R[3]*X + R[4]*Y
  71.  
  72. # (x,y) -> (X,Y) (X=R(-a)x)
  73. X(x, y) = R[4]x - R[2]*y
  74. Y(x, y) = -R[3]*x + R[1]*y
  75.  
  76. # Time
  77. Time(t) = sprintf("{/TimesNewRoman:Italic t} = %.1f s", t)
  78.  
  79. # Plot ----------------------------------------
  80. # Initial Value
  81. t  = 0
  82. b_inc = 3
  83.  
  84. do for[i=1:N]{
  85.     X[i]  = X(150, 60)
  86.     Y[i]  = Y(150, 60)
  87.     VX[i] = V*cos(b + (b_inc*(i-int(N/2))) * pi/180)
  88.     VY[i] = V*sin(b + (b_inc*(i-int(N/2))) * pi/180)
  89.     e[i] = 0.3+(1.0-0.3)/N*i                        # float division
  90.     
  91.     # Vector rotation (Rotation matrix R(a))
  92.     x[j]  = x(X[j] , Y[j] )
  93.     y[j]  = y(X[j] , Y[j] )
  94.     vx[j] = x(VX[j], VY[j])
  95.     vy[j] = y(VX[j], VY[j])  
  96. }
  97.  
  98. # Draw initiate state for 70 steps
  99. do for [i=1:70] {
  100.     # Time
  101.     set label 1 Time(t) left at graph 0.01, 0.93 font 'TimesNewRoman, 25'
  102.     
  103.     # Balls
  104.     do for[j=1:N]{
  105.         set obj j circ at x[j], y[j] size r fc rgb color[j] fs solid noborder
  106.     }
  107.     
  108.     # Draw ground and balls
  109.     plot tan(a)*x with filledcurve x1 lc rgb 'gray50'
  110. }
  111.  
  112. # Update until all of balls are framed out
  113. i = 0                     # The number of loops
  114. while(1){
  115.     i = i +1
  116.     t = t + dt
  117.     set label 1 Time(t)
  118.  
  119.     do for[j=1:N]{
  120.         # 4th order Runge-Kutta
  121.         temp_X  = X[j]  + rk1(X[j], Y[j], VX[j], VY[j])
  122.         temp_Y  = Y[j]  + rk2(X[j], Y[j], VX[j], VY[j])
  123.         temp_VX = VX[j] + rk3(X[j], Y[j], VX[j], VY[j])
  124.         temp_VY = VY[j] + rk4(X[j], Y[j], VX[j], VY[j])
  125.         
  126.         X[j]=temp_X; Y[j]=temp_Y; VX[j]=temp_VX; VY[j]=temp_VY
  127.         
  128.         # Judge whether balls bounce or not
  129.         if(Y[j]<r){
  130.             Y[j]  = r
  131.             VY[j] = -e[j]*VY[j]      # inelastic collision
  132.         }
  133.         
  134.         # Vector rotation (Rotation matrix R(a))
  135.         x[j]  = x(X[j] , Y[j] )
  136.         y[j]  = y(X[j] , Y[j] )
  137.         vx[j] = x(VX[j], VY[j])
  138.         vy[j] = y(VX[j], VY[j])
  139.  
  140.         # Update objects
  141.         set obj N*i+j circ at x[j], y[j] size r fc rgb color[j] fs solid noborder
  142.         # Make old objects trajectory of the ball 
  143.         set obj N*(i-1)+j at x[j], y[j] size 0.1
  144.  
  145.         # Start to disappear
  146.         if(i>=dis){
  147.             unset object N*(i-dis)+j
  148.         }
  149.     }
  150.     
  151.     # Decimate and plot
  152.     if(i%cut==0){
  153.         replot
  154.     }
  155.     
  156.     # Count the number of balls being framed out
  157.     do for[j=1:N]{
  158.         if(x[j]<-L-r){
  159.             cnt = cnt + 1
  160.         }
  161.     }
  162.     
  163.     # Exit from the loop when all balls being framed out
  164.     if(cnt == N){
  165.         break
  166.     } else {
  167.         cnt =0
  168.     }
  169. }
  170.  
  171. set out

BallBouncing.gif

V=30\mathrm{m/s},\ a=5^\circ,\ b=-85^\circ,\ e_i=0.3+0.14i,\ b_{\mathrm{inc}}=3
t=0: X_j=200\cos\left(a\right)+60\sin\left(a\right),\ Y_{j}=-200\sin\left(a\right)+60\cos\left(a\right)\\ \left(t=0: x_{j}=200,\ y_{j}=60\right)

BallBounding2.gif

V=30\mathrm{m/s},\ a=15^\circ,\ b=-75^\circ,\ e=0.7, b_{\mathrm{inc}}=3
t=0: X_j=200\cos\left(a\right)+120\sin\left(a\right),\ Y_{j}=-200\sin\left(a\right)+120\cos\left(a\right)\\ \left(t=0: x_{j}=200,\ y_{j}=120\right)

BallBounding3.gif

V=30\mathrm{m/s},\ a=5^\circ,\ b=15^\circ,\ e=0.8, b_{\mathrm{inc}}=9
t=0: X_j=-0.8L,\ Y_{j}=0.25B

Extra

BallBounding4.gif

V=30\mathrm{m/s},\ a=5^\circ,\ b=-75^\circ,\ e=0.8, b_{\mathrm{inc}}=5
t=0: X_j=60\sin\left(a\right),\ Y_{j}=60\cos\left(a\right)\\ \left(t=0: x_{j}=0,\ y_{j}=60\right)

Search This Blog

Translate

QooQ