@LiuYongJie 2016-10-23T10:21:59.000000Z 字数 11693 阅读 406

Ex_ 05 - The Trajectory of a Cannon Shell (2.9)

author: Liu Yongjie 2014301020094

Abstract

In this homework several problems involving the motion of objects through the atmosphere are considered. Still, the problems are all described by ordinary differential equations in which initial values are given and can be solved using Euler method. This report solved cannon trajectory problem, in which air drag and reduced air density at high altitudes are considered.

Key Words: Projectile Motion; Euler Method; Air Drag; Air Density; Cannon Shell

Background

In realistic projectile motion1, we will fail to to obtain a realistic description if we do not include Air resistance must be included if we are to obtain a realistic description of the problem.
In general, air resistance can be written in the fairly form:

²

$F_drag≈-B_1V-B_2V²$
and the ......... term dominates for most objects at any reasonable velocity. So we can use the fact that an object must push air in front it out of the way to move through the atmosphere.

²

$F_drag≈\frac{1}{2}CρAv²$
C is the drag coefficient, ρ is air density, A is the front area of the object.
Another important piece of physics we should consider is air density

$ρ=ρ_0exp(\frac{-y}{y_0})$ .
where

$α≈6.5×10^{-3}K/m , T_0$ is the sea level temperature, and the exponent α≈2.9 for air.
With density correction:

$ρ=ρ_0exp(\frac{-y}{y_0})$
where

$y_0=\frac {k_BT} {mg}≈1.0×10^4m$ and

³

$ρ_0=1.225kg/m³$ is the density at sea level.
This isothermal model of atmosphere is perhaps not realistic, since we know that the air temperature can vary quitea bit over altitude changes of a few kilometers. This leads to the second simplest model adiabatic model:

$ρ=ρ_0({1-\frac {ay} {T_0}})^α$
where a≈

$6.5×10^{-3}K/m \;and\; α≈2.5$ for air. Here

$T_0$ is the sea level temperature (in K)≈300K.
Since the drag force due to air resistance is proportional to the density, so for both models we have

$F_{drag}^{*}=\frac {ρ} {ρ_0}F_drag(y=0)$
Now we can get for isothermal ideal gas, the equations of motion become

$\frac {d^2x} {dt^2}=-\frac {B_2} {m}vv_xe^{\frac {-y} {y_0}}$

$\frac {d^2y} {dt^2}=-g-\frac {B_2} {m}vv_ye^{\frac {-y} {y_0}}$
Also, for adiabatic approximation, the equations of motion become

$\frac {d^2x} {dt^2}=-\frac {B_2} {m}vv_x({1-\frac {ay} {T_0}})^α$

$\frac {d^2y} {dt^2}=-g-\frac {B_2} {m}vv_y({1-\frac {ay} {T_0}})^α$
On top of that, we can use several numerical approaches including the simple Euler method and the Runge-Kutta method to yield the numerical solutions to these problems.

content

code:


import math
import matplotlib.pyplot as plt
# import modules
g=9.8
B2m=0.00004
y_zero=10000
a=0.0065
T0=300
alpha=2.5
# set constants
class cannon0:
    "the simplest model with no air drag, no air density variance, no probability distribution"
    # initialize variables    
    def __init__(self,v0,theta,yFinal=0):
        self.x0=0
        self.y0=0
        self.yFinal=yFinal
        self.v0=v0
        self.Theta=theta
        self.theta=theta*math.pi/180
        self.vx0=self.v0*math.cos(self.theta)
        self.vy0=self.v0*math.sin(self.theta)
        self.dt=0.01
        return None
    # external force other than gravity, no force in simplest case        
    def F(self,vx,vy,y=1):
        return 0,0
    # calculate the flying trajactory
    def fly(self):
        self.X=[self.x0]
        self.Y=[self.y0]
        self.Vx=[self.vx0]
        self.Vy=[self.vy0]
        self.T=[0]
        while not (self.Y[-1]<self.yFinal and self.Vy[-1]<0):
            newVx=self.Vx[-1]+self.F(vx=self.Vx[-1],vy=self.Vy[-1],y=self.Y[-1])[0]*self.dt
            newVy=self.Vy[-1]-g*self.dt+self.F(self.Vx[-1],self.Vy[-1])[1]*self.dt
            self.Vx.append(newVx)
            self.Vy.append(newVy)
            meanVx=0.5*(self.Vx[-1]+self.Vx[-2])
            meanVy=0.5*(self.Vy[-1]+self.Vy[-2])
#            meanV=math.sqrt(meanVx**2+meanVy**2) # not used in Cannon0 because there is no air drag
            newX=self.X[-1]+meanVx*self.dt
            newY=self.Y[-1]+meanVy*self.dt
            self.X.append(newX)
            self.Y.append(newY)
        # fix the final landing coordinate        
#        r=-self.Y[-2]/self.Y[-1]
        self.X[-1]=((self.Y[-2]-self.yFinal)*self.X[-1]+(self.yFinal-self.Y[-1])*self.X[-2])/(self.Y[-2]-self.Y[-1])
        self.Y[-1]=self.yFinal
        return 0
    # get the final distance shells can reach
    def distance(self):
        return self.X[-1]
    def height(self):
        return max(self.Y)
    # represent trajectory 
    def plot(self, color):
        plt.plot(self.X,self.Y,color,label="$%dm/s$,$%d\degree$, no air drag"%(self.v0,self.Theta))
        return 0
class cannon1(cannon0):
    "the second simplest model with no air drag under constant air density, no probability distribution"    
    # external force other than gravity        
    def F(self,vx,vy,y=1):
        vl=math.sqrt(vx**2+vy**2)
        self.Fx=-B2m*vx*vl
        self.Fy=-B2m*vy*vl
        return self.Fx,self.Fy
    def plot(self, color):
        plt.plot(self.X,self.Y,color,label="$%dm/s$,$%d\degree$, uniform air drag"%(self.v0,self.Theta))
        return 0
class cannon2(cannon0):
    "the model concerning ISOTHERMAL air density variance but no probability distribution"
    def F(self,vx,vy,x=1,y=1):
        vl=math.sqrt(vx**2+vy**2)
        self.Fx=-B2m*vx*vl*math.exp(-y/y_zero)
        self.Fy=-B2m*vy*vl*math.exp(-y/y_zero)
        return self.Fx,self.Fy
    def plot(self, color):
        plt.plot(self.X,self.Y,color,label="$%dm/s$,$%d\degree$, isothermal air drag"%(self.v0,self.Theta))
        return 0
class cannon3(cannon0):
    "the model concerning ADIABATIC air density variance but no probability distribution"    
    def F(self,vx,vy,y=1):
        vl=math.sqrt(vx**2+vy**2)
        self.Fx=-B2m*vx*vl*(1-a*y/T0)**alpha
        self.Fy=-B2m*vy*vl*(1-a*y/T0)**alpha
        return self.Fx,self.Fy
    def plot(self, color):
        plt.plot(self.X,self.Y,color,label="$%dm/s$,$%d\degree$, adiabatic air drag"%(self.v0,self.Theta))
        return 0
# select the angle casting the largest distance
Distance=[]
for i in range(90):
    A=cannon1(600,i)
    A.fly()
    newDistance=A.distance()
    Distance.append(newDistance)
maxDistance=max(Distance)
maxAngle=Distance.index(maxDistance)
print maxAngle
sub1=plt.subplot(221)
A=cannon0(600,45)
A.fly()
sub1.plot(A.X,A.Y,'r-',label='$45\degree$')
A1=cannon0(600,35)
A1.fly()
sub1.plot(A1.X,A1.Y,'r--',label='$35\degree$')
A2=cannon0(600,55)
A2.fly()
sub1.plot(A2.X,A2.Y,'r:',label='$55\degree$') 
sub1.set_title('No air drag')   
sub1.legend(loc="upper right",frameon=False)
plt.xlabel("Distance [m]")
plt.ylabel("Hieght [m]")  
sub2=plt.subplot(222)
B=cannon1(600,40)
B.fly()
sub2.plot(B.X,B.Y,'g-',label='$40\degree$')    
B1=cannon1(600,30)
B1.fly()
sub2.plot(B1.X,B1.Y,'g--',label='$30\degree$')
B2=cannon1(600,50)
B2.fly()
sub2.plot(B2.X,B2.Y,'g:',label='$50\degree$') 
sub2.set_title('Uniform Air Drag')
sub2.legend(loc="upper right",frameon=False)
plt.xlabel("Distance [m]")
plt.ylabel("Hieght [m]")  
sub3=plt.subplot(223)
C=cannon2(600,45)
C.fly()
sub3.plot(C.X,C.Y,'b-',label='$45\degree$')
C1=cannon2(600,35)
C1.fly()
sub3.plot(C1.X,C1.Y,'b--',label='$35\degree$')
C2=cannon2(600,55)
C2.fly()
sub3.plot(C2.X,C2.Y,'b:',label='$55\degree$') 
sub3.set_title('Isothermal Air Drag')
sub3.legend(loc="upper right",frameon=False)
plt.xlabel("Distance [m]")
plt.ylabel("Hieght [m]")  
sub4=plt.subplot(224)
D=cannon3(600,42)
D.fly()
sub4.plot(D.X,D.Y,'y-',label='$42\degree$')
D1=cannon3(600,32)
D1.fly()
sub4.plot(D1.X,D1.Y,'y--',label='$32\degree$')
D2=cannon3(600,52)
D2.fly()
sub4.plot(D2.X,D2.Y,'y:',label='$52\degree$') 
sub4.set_title('Adiabatic Air Drag')
sub4.legend(loc="upper right",frameon=False)
plt.xlabel("Distance [m]")
plt.ylabel("Hieght [m]")   
plt.show()
#print A.distance()
#print B.distance()
#print C.distance()
#print D.distance()

result

# -*- coding: utf-8 -*-
import math
import numpy as np
from matplotlib import pyplot as plt 
from matplotlib import animation  
g = 9.8
b2m = 4e-5
r = []
X = []
Y = []
class flight_state:
    def __init__(self, _x = 0, _y = 0, _vx = 0, _vy = 0, _t = 0):
        self.x = _x
        self.y = _y
        self.vx = _vx
        self.vy = _vy
        self.t = _t
class cannon:
    def __init__(self, _fs = flight_state(0, 0, 0, 0, 0), _dt = 0.01):
        self.cannon_flight_state = []
        self.cannon_flight_state.append(_fs)
        self.dt = _dt
    def next_state(self, current_state):
        global g
        next_x = current_state.x + current_state.vx * self.dt
        next_vx = current_state.vx
        next_y = current_state.y + current_state.vy * self.dt
        next_vy = current_state.vy - g * self.dt
        #print next_x, next_y
        return flight_state(next_x, next_y, next_vx, next_vy, current_state.t + self.dt)
    def shoot(self):
        while not(self.cannon_flight_state[-1].y < 0):
            self.cannon_flight_state.append(self.next_state(self.cannon_flight_state[-1]))
        r = - self.cannon_flight_state[-2].y / self.cannon_flight_state[-1].y
        self.cannon_flight_state[-1].x = (self.cannon_flight_state[-2].x + r * self.cannon_flight_state[-1].x) / (r + 1)
        self.cannon_flight_state[-1].y = 0
    def show_trajectory(self):
        x = []
        y = []
        for fs in self.cannon_flight_state:
            x.append(fs.x)
            y.append(fs.y)
        X.append(x)
        Y.append(y)
#        if n == 1:
#            plt.subplot(111)
#            line, = plot(x,y, label = labe)
#            xlabel(r'$x(m)$', fontsize=16)
#            ylabel(r'$y(m)$', fontsize=16)
#            text(40767, 14500, 'initial speed: 700m/s\n' + 'firing angle: 45' + r'$^{\circ}$', color='black')
#            title('Trajectory of cannon shell')#
#            ax = gca()
#            ax.spines['right'].set_color('none')
#            ax.spines['top'].set_color('none')
#            ax.xaxis.set_ticks_position('bottom')
#            ax.yaxis.set_ticks_position('left')
#            ax.set_xlim(0, 60000)
#            ax.set_ylim(0, 18000)
#        #show()
class drag_cannon(cannon):
    def next_state(self, current_state):
        global g, b2m
        v = math.sqrt(current_state.vx * current_state.vx + current_state.vy * current_state.vy)
        next_x = current_state.x + current_state.vx * self.dt
        next_vx = current_state.vx - b2m * v * current_state.vx * self.dt
        next_y = current_state.y + current_state.vy * self.dt
        next_vy = current_state.vy - g * self.dt - b2m * v * current_state.vy * self.dt
        #print next_x, next_y
        return flight_state(next_x, next_y, next_vx, next_vy, current_state.t + self.dt)
class adiabatic_drag_cannon(cannon):
    def next_state(self, current_state):
        global g, b2m
        factor = (1 - 6.5e-3 * current_state.y / 288.15) ** 2.5
        v = math.sqrt(current_state.vx * current_state.vx + current_state.vy * current_state.vy)
        next_x = current_state.x + current_state.vx * self.dt
        next_vx = current_state.vx - factor * b2m * v * current_state.vx * self.dt
        next_y = current_state.y + current_state.vy * self.dt
        next_vy = current_state.vy - g * self.dt - factor * b2m * v * current_state.vy * self.dt
        #print next_x, next_y
        return flight_state(next_x, next_y, next_vx, next_vy, current_state.t + self.dt)
class isothermal_drag_cannon(cannon):
    def next_state(self, current_state):
        global g, b2m
        factor = math.exp(-current_state.y / 1e4)
        v = math.sqrt(current_state.vx * current_state.vx + current_state.vy * current_state.vy)
        next_x = current_state.x + current_state.vx * self.dt
        next_vx = current_state.vx - factor * b2m * v * current_state.vx * self.dt
        next_y = current_state.y + current_state.vy * self.dt
        next_vy = current_state.vy - g * self.dt - factor * b2m * v * current_state.vy * self.dt
        #print next_x, next_y
        return flight_state(next_x, next_y, next_vx, next_vy, current_state.t + self.dt)
speed = 700
theta = np.arange(30., 56., 0.1)
v_x = [speed * math.cos(i * math.pi / 180) for i in theta]
v_y = [speed * math.sin(i * math.pi / 180) for i in theta]
def wahaha():
    b = []
    for i in range(len(theta)):
        b.append(adiabatic_drag_cannon(flight_state(0, 0, v_x[i], v_y[i], 0)))
        #labe2 = str(theta[i]) + r'$^{\circ}$'
        b[i].shoot()
#        xx = [b[i].cannon_flight_state[j].x for j in range(len(b[i].cannon_flight_state))]
#        yy = [b[i].cannon_flight_state[j].y for j in range(len(b[i].cannon_flight_state))]
#        X.append(xx)
#        Y.append(yy)
        r.append(b[i].cannon_flight_state[-1].x)
        b[i].show_trajectory()
        #legend(loc='upper left', frameon=False)
wahaha()
p = r.index(max(r))
print theta[p], max(r),type(max(r))
fig = plt.figure()  
ax = plt.axes(xlim=(0, 40000), ylim=(0, 18000))  
#labe = str(theta[i]) + r'$^{\circ}$'
line, = ax.plot([], [],lw = 2,label = 'adiabatic model' ,color = 'red') 
angle_text = ax.text(24167, 14400, '')
maxrange_text = ax.text(24167, 12400, '')
plt.xlabel(r'$x(m)$', fontsize=16)
plt.ylabel(r'$y(m)$', fontsize=16)
plt.title('Trajectory of cannon shell') 
plt.legend(loc='upper left', frameon=False)
# initialization function: plot the background of each frame  
def init():  
    line.set_data([], [])  
    angle_text.set_text('')
    maxrange_text.set_text('')
    return line, angle_text, maxrange_text
# animation function.  This is called sequentially  
# note: i is framenumber  
def animate(i):  
    x = X[i]  
    y = Y[i]
    line.set_data(x, y)  
    pct = float(X[i][-1])/float(max(r))*100
    angle_text.set_text('initial speed: 700m/s\n' + 'firing angle: %s' % theta[i] + r'$^{\circ}$' + '\nrange: %f %%' % pct)
    if i > p:
        maxrange_text.set_text('max range: %s' % max(r))
        maxrange_text.set_color(color='red')
    return line,  angle_text, maxrange_text
fram = len(theta)
# call the animator.  blit=True means only re-draw the parts that have changed.  
anim=animation.FuncAnimation(fig, animate, init_func=init, frames=fram, interval=30, blit=True) 
#anim.save('basic_animation.mp4', fps=30, extra_args=['-vcodec', 'libx264'])  
plt.show()
#show()

result

Conclusion

Through the first class progrem, We analysed the projectile motion of the cannon shell. We drawed the conclusion:
(1) Object orientation progreming can be more efficient compared to normal progreming.
(2) Big angle shell suffer a more serious influence not only from the airdrag but also from the air density change.
(3) With all the factors considered the angle 42.5 gained the farest diatance.

Thanks to