@computationalphysics-2014301020090 2016-10-16T16:48:44.000000Z 字数 12640 阅读 243

The 5th homework

Chapter 2.2 problem 2.9:The trajectory of a cannon shell

The 5th homework
Chapter 2.2 problem 2.9:The trajectory of a cannon shell

在此输入正文

●Abstract

problem 2.9

Calculate the trajectory of our cannon shell including both air drag and the reduced air density at high altitudes so that you can reproduce the results in Figure2.5 . Perform your calculation for different firing angles and determine the value of angle the gives the maximum range.

●background

We investigate the effects of both air drag and the reduced air density at high altitudes in order to reproduce the results in textbook for L1. We generalize it to cover targets of different altitude and use Python to design the programs, which can realize the purpose of this assignment. This article deals with the ideas for these problems, the programs ,and the results.
Without air drag, the trajectory of a cannon shell would be a perfect parabola, and the maximum range occurs for a firing angle of 45°. However, when the air drag is considered, things are distinctly different: the maximum range would be a smaller one.
Since the cannon shell could reach considerable high altitudes, we must be aware of the reduced air density. In this problem, we solved both the isothermal and adiabatic cases, and generalize the problem to cover targets at different altitudes.
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

●Idea

●we assumed that the drag force for cannon shell is given by

$F_{drag}=-B_2v^2$
●according to Newton's second law,

$\frac {d^2x} {dt^2}=\frac {F_{drag}cosθ} {m}=\frac {F_{drag}} {m}\frac {v_x} {v}=-\frac {B_2} {m}vv_x$

$\frac {d^2y} {dt^2}=\frac {F_{drag}sinθ} {m}-g=\frac {F_{drag}} {m}\frac {v_y} {v}-g=-\frac {B_2} {m}vv_y-g$
●We also take the reduced air density at high altitudes into account. There are two mainstream models regard of this problem in statistical mechanics: isothermal ideal gas, and adiabatic approximation.
For isothermal ideal gas, the pressure p depends on altitudes according to

$p(y)=p(0)e^{\frac {-mgy} {k_BT}}$
●For an ideal gas the pressure is proportional to the destiny,so this leads to

$ρ=ρ_0e^{\frac {-y} {y_0}}$
where

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

$ρ_0$ is the density at sea level≈1.225

$Kg/m^3$ .
●For adiabatic approximation, the density depends on altitudes according to

$ρ=ρ_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.

●Code1

# -*- coding: utf-8 -*-
"""
Created on Sun Oct 16 16:57:48 2016
@author: 胡胜勇
"""
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()

●fifure

●For more intuitive observation about the relation between the and the angle,We introduce another list of code

●Code2

# -*- coding: utf-8 -*-
"""
Created on Sun Oct 16 20:11:38 2016
@author: 胡胜勇
"""
# -*- 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

Approximately，relations between the and the angle are as follows:

Condition	Max Angle	Max Range
no drag	$45.0 ^o$	50004.950 m
uniform air drag	$38.8^o$	22070.050 m
isothermal air drag	$45.9^o$	26621.572 m
adiabatic air drag	$43.8^o$	24641.195m

Pytion is really interesting as well as powerful,I need to learn more about it.

●Acknowledgments

Thanks to Yue Shaosheng,Wang Shibing