Example - Explosion of a Fireworks Shell

A fireworks shell is launched from the origin (launch speed $v_0$, initial angle $\theta$) and travels along the usual parabolic path above flat ground. Had it failed to explode, it would have had the usual range

$$ R = \frac{2v_0^2}{g}\sin(\theta)\cos(\theta) $$

Right at the very peak of its path, it explodes into two pieces of equal mass. Just after the explosion, one of the pieces shoots straight upward with a speed $v_1$ after the explosion. What is the velocity of the second piece after the explosion and where is the center of mass of the system after both pieces land?


In [1]:
%matplotlib inline

In [2]:
from __future__ import division, print_function
from ivisual import *
from numpy import *
import matplotlib.pyplot as plt



In [27]:
m=1 #mass of shell
m1=m/2 #mass of piece 1
m2=m1 #mass of piece 2
g=9.8

v0=20 #initial speed of shell
theta=45*pi/180 #initial angle of shell
v=v0*vector(cos(theta), sin(theta), 0)
v1=v0*vector(cos(theta), sin(theta), 0)
v2=v0*vector(cos(theta), sin(theta), 0)

phi1=90*pi/180 #final angle of piece 1
v1mag=5 #final speed of piece 1

r=vector(0,0,0) #position of shell
r1=vector(0,0,0) #position of piece 1
r2=vector(0,0,0) #position of piece 2

#gravitational forces on the shell and pieces
gfield=vector(0,-g,0)
Fgrav=m*gfield
Fgrav1=m1*gfield
Fgrav2=m2*gfield

tpeak=v.y/g

t=0
dt=0.005

#initial net force on the shell and pieces
Fnet=Fgrav
Fnet1=Fgrav1
Fnet2=Fgrav2

#we need a boolean to change the conditions when the explosion occurs
explosion=False

#store the data
tlist=[]
x1list=[]
y1list=[]
vx1list=[]
x2list=[]
y2list=[]
vx2list=[]
xshell=[]
yshell=[]

#while particles are above the ground
while r1.y>-0.01 or r2.y>-0.01:

    # we only want to run this code once, just as soon as the particle crosses the peak
    if(t>tpeak and explosion==False):
        Fnet=vector(0,0,0)
        Fnet1=Fgrav1
        Fnet2=Fgrav2
        v1=v1mag*vector(cos(phi1),sin(phi1),0)
        v2=1/m2*vector(m*v.x-m1*v1.x, -m1*v1mag,0)

        #set velocity of shell to zero
        v=vector(0,0,0)
        
        #set explosion to True so that this code will only run once
        explosion=True

    #if the particle hits the ground, set the force and velocity to zero
    if(r1.y<-0.01):
        Fnet1=vector(0,0,0)
        v1=vector(0,0,0)

    #if the particle hits the ground, set the force and velocity to zero
    if(r2.y<-0.01):
        Fnet2=vector(0,0,0)
        v2=vector(0,0,0)

    #Model the motion of the shell and pieces
    v=v+Fnet/m*dt
    r=r+v*dt
    
    v1=v1+Fnet1/m1*dt
    r1=r1+v1*dt
    
    v2=v2+Fnet2/m2*dt
    r2=r2+v2*dt
    
    t=t+dt

    #store data
    tlist.append(t)
    x1list.append(r1.x)
    y1list.append(r1.y)
    vx1list.append(v1.x)
    x2list.append(r2.x)
    y2list.append(r2.y)
    vx2list.append(v2.x)
    xshell.append(r.x)
    yshell.append(r.y)

xcmlist=1/2*(array(x1list)+array(x2list))
ycmlist=1/2*(array(y1list)+array(y2list))

In [28]:
plt.title('path')
plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.plot(x1list,y1list,'m.',x2list,y2list,'y.')
plt.show()



In [29]:
plt.title('CM path')
plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.plot(xcmlist,ycmlist,'c.')
plt.show()



In [30]:
plt.title('x vs t for piece 1')
plt.xlabel('t (s)')
plt.ylabel('x (m)')
plt.plot(tlist,x1list,'b.')
plt.show()



In [31]:
plt.title('y vs t for piece 1')
plt.xlabel('t (s)')
plt.ylabel('y (m)')
plt.plot(tlist,y1list,'r.')
plt.show()



In [32]:
plt.title('x vs t for piece 2')
plt.xlabel('t (s)')
plt.ylabel('x (m)')
plt.plot(tlist,x2list,'b.')
plt.show()



In [33]:
plt.title('y vs t for piece 2')
plt.xlabel('t (s)')
plt.ylabel('y (m)')
plt.plot(tlist,y2list,'r.')
plt.show()



In [43]:
#calculate cm of pieces on ground
print("numerical results")
print("============================")
xcm=0.5*(x1list[size(x1list)-1]+x2list[size(x2list)-1])
print("x of CM (numerical) =",xcm)
h=max(yshell)
print("h (numerical)",h)


#theoretical calculations
print()
print("theoretical results")
print("============================")

R=2*v0*cos(theta)*v0*sin(theta)/g
print("R = ",R)

htheor=R/4*tan(theta)
print("h = ",htheor)
xcmtheor=R/2-v0/g*v1mag*cos(theta)+v0/g*cos(theta)*sqrt(v1mag*v1mag+2*g*htheor)
print("x of CM (theoretical) = ", xcmtheor)


numerical results
============================
x of CM (numerical) = 34.8603643125
h (numerical) 10.1687552982

theoretical results
============================
R =  40.8163265306
h =  10.2040816327
x of CM (theoretical) =  40.8163265306

In [ ]: