@rrtcc 2016-10-16T14:20:40.000000Z 字数 9047 阅读 40

第五次作业

炮弹轨迹 极值 空气阻力

摘要

通过解决Exercise2.9，复现Finger2.5的结果，并进一步考虑,引入空气阻力，进行相关计算。

Exercise 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 Figure 2.5. Perform your calculation for different firing angles and determine the value of the angle that gives the maximum range.

背景

不考虑空气阻力

以下斜体部分摘自课本page25
xm-1111型中程炮弹
The Euler method we used to treat the problem in chapter 1 can easily be generalized to deal with motion in two spatial dimensions. To be specific, we can consider a projectile such as a shell shot by a cannon. If we ignore the air resistance, the equations of motion, which are again obtained from Newton's second law, can be written as
牛顿方程
We can write each of these second-order equations as two first-order differential equations

利用在上一次作业中所使用的一级欧拉近似法有以下方程
近似方程
于是可以利用python进行相关计算

考虑空气阻力

The Euler method can easily be generalized to deal with motion in two spatial dimensions.
If we ignore the air resistance, the equations of motion, which are again obtained from Newton's second law, can be written as

$\frac{d^2x}{dt^2}=0$
$\frac{d^2y}{dt^2}=-g$ where x and y are the horizontal and vertical coordinates of the projectile, and g is the acceleration due to the gravity. These are second-order differential equations.We can put them into four first-order ordinary differential equations as

$\frac{dx}{dt}=v_x$
$\frac{dv_x}{dt}=0$
$\frac{dy}{dt}=v_y$
$\frac{dv_y}{dt}=-g$ where $v_x$ and $v_y$ are the x and y components of the velocity.
To use the Euler method, we write each derivative in finite difference form

$x_{i+1}=x_i+v_{x,i}\Delta t$
$v_{x,i+1}=v_{x,i}$
$y_{i+1}=y_i+v_{y,i}\Delta t$
$v_{y,i+1}=v_{y,i}-g\Delta t$ Then we consider air resistance, we assume that the magnitude of the drag force on cannon shell is given by

$F_{drag}=-B_2v^2$ The components of the drag force are thus

$F_{drag,x}=-B_2vv_x$
$F_{drag,y}=-B_2vv_y$ Adding this force to the equations of motion leads to

$x_{i+1}=x_i+v_{x,i}\Delta t$
$v_{x,i+1}=v_{x,i}-\frac{B_2vv_{x,i}}{m}\Delta t$
$y_{i+1}=y_i+v_{y,i}\Delta t$
$v_{y,i+1}=v_{y,i}-g\Delta t-\frac{B_2vv_{y,i}}{m}\Delta t$ If we consider the reduced air density at high altitudes, there is a formula of density which depends on altitude

$\rho=\rho_0(1-\frac{ay}{T_0})^\alpha$ For air： $a≈6.5×10^{-3}K/m$ , $α≈2.5$
So the drag force at altitude can be written as

$F^*_{drag}=\frac{\rho}{\rho_0}F_{drag}(y=0)$ Then the differential equation of motion can be modified into

$x_{i+1}=x_i+v_{x,i}\Delta t$
$v_{x,i+1}=v_{x,i}-(1-\frac{ay}{T_0})^\alpha\frac{B_2vv_{x,i}}{m}\Delta t$
$y_{i+1}=y_i+v_{y,i}\Delta t$
$v_{y,i+1}=v_{y,i}-g\Delta t-(1-\frac{ay}{T_0})^\alpha\frac{B_2vv_{y,i}}{m}\Delta t$ The air resistance was included with $B_2/m=4×10^{-5}m^{-1}$ ,and the initial speed was 700 $m/s$ .

正文

在不考虑空气阻力的情况下，设计代码如下

import math  
import pylab as pl
class bicycle:
    def __init__(self, power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/6),\
       initial_Vy=50*math.sin(math.pi/6)):
        self.v = [initial_velocity]
        self.x=[initial_x]
        self.y=[initial_y]
        self.Vx=[initial_Vx]
        self.Vy=[initial_Vy]
        self.t = [0]
        self.m = mass
        self.p = force
        self.dt = time_step
    def run(self):
        while(self.y[-1]>=0):
            self.x.append(self.x[-1]+self.Vx[0]*self.dt)
            self.y.append(self.y[-1]+self.Vy[-1]*self.dt)
            self.Vy.append(self.Vy[-1]-9.8*self.dt)
        print("yn=",  self.y[-1])
        print("y(n-1)=",self.y[-2])
        print("xn=",self.x[-1])
        print("x(n-1)=",self.x[-2])
class diff_step_check(bicycle):    
    def show_result1(self,style='b'):
        pl.plot(self.x, self.y,style,label='30')
        pl.title('projectile motion without air resistance')
        pl.xlabel('x')
        pl.ylabel('y')
        pl.legend(loc='best')
a=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/6),\
       initial_Vy=50*math.sin(math.pi/6))
a.run()
a.show_result1('b')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,'k--',label='35')
        pl.title('projectile motion without air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/6*(7/6)),\
       initial_Vy=50*math.sin(math.pi/6*(7/6)))
b.run()
b.show_result2('g')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,'m-.',label='40')
        pl.title('projectile motion without air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
c=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/6*(4/3)),\
       initial_Vy=50*math.sin(math.pi/6*(4/3)))
c.run()
c.show_result2('r')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,'y-',label='45')
        pl.title('projectile motion without air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
d=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/4),\
       initial_Vy=50*math.sin(math.pi/4))
d.run()
d.show_result2('c')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,'b:',label='50')
        pl.title('projectile motion without air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/4*(10/9)),\
       initial_Vy=50*math.sin(math.pi/4*(10/9)))
b.run()
b.show_result2('m')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self,'r--',style,label='55')
        pl.title('projectile motion without air resistance')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=50*math.cos(math.pi/4*(11/9)),\
       initial_Vy=50*math.sin(math.pi/4*(11/9)))
b.run()
b.show_result2('y')
show()

结果如下图所示无空气阻力
课本结果如下
可以看出复现了课本的结果
考虑空气阻力，代码如下

import math  
import pylab as pl
class bicycle:
    def __init__(self, power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/6),\
       initial_Vy=700*math.sin(math.pi/6)):
        self.v = [initial_velocity]
        self.x=[initial_x]
        self.y=[initial_y]
        self.Vx=[initial_Vx]
        self.Vy=[initial_Vy]
        self.t = [0]
        self.m = mass
        self.p = power
        self.dt = time_step
    def run(self):
        while(self.y[-1]>=0):
            self.x.append(self.x[-1]+self.Vx[-1]*self.dt)
            self.Vx.append(self.Vx[-1]-math.exp(-self.y[-1]/(10**4))*4/(10**5)*self.v[-1]*self.Vx[-1]/self.m*self.dt)
            self.y.append(self.y[-1]+self.Vy[-1]*self.dt)
            self.Vy.append(self.Vy[-1]-9.8*self.dt-math.exp(-self.y[-1]/(10**4))*4/(10**5)*self.v[-1]*self.Vy[-1]/self.m*self.dt)
            self.v.append(math.sqrt(self.Vx[-1]**2+self.Vx[-1]**2))
        print("yn=",  self.y[-1])
        print("y(n-1)=",self.y[-2])
        print("xn=",self.x[-1])
        print("x(n-1)=",self.x[-2])
        r=-self.y[-2]/self.y[-1]
        xl=(self.x[-2]+r*self.x[-1])/(r+1)
        print("零点为",xl)
class diff_step_check(bicycle):    
    def show_result1(self,style='b'):
        pl.plot(self.x, self.y,style,label='30')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x')
        pl.ylabel('y')
        pl.legend(loc='best')
from pylab import *
a=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/6),\
       initial_Vy=700*math.sin(math.pi/6))
a.run()
a.show_result1('b')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,'k--',label='35')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/6*(7/6)),\
       initial_Vy=700*math.sin(math.pi/6*(7/6)))
b.run()
b.show_result2('g')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,'m-.',label='40')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
c=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/6*(4/3)),\
       initial_Vy=700*math.sin(math.pi/6*(4/3)))
c.run()
c.show_result2('r')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='45')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
d=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/4),\
       initial_Vy=700*math.sin(math.pi/4))
d.run()
d.show_result2('c')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,'b:',label='50')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/4*(10/9)),\
       initial_Vy=700*math.sin(math.pi/4*(10/9)))
b.run()
b.show_result2('m')
class diff_step_check(bicycle):
     def show_result2(self,style='b'):
        pl.plot(self.x, self.y,style,label='55')
        pl.title('projectile motion with air resistance and isothermal atmosphere')
        pl.xlabel('x(m)')
        pl.ylabel('y(m)')
        pl.legend(loc='best')
b=diff_step_check( power=10, mass=1, time_step=0.1,\
        initial_velocity=1 ,initial_x=0,initial_y=0,initial_Vx=700*math.cos(math.pi/4*(11/9)),\
       initial_Vy=700*math.sin(math.pi/4*(11/9)))
b.run()
b.show_result2('y')
show()

结果如下图所示有空气阻力
无法得出角度为何值时有最大值，故增加60度的曲线如下

故45到60度时取得最大值，可考虑缩小范围
最后结果如下
故最大值在52到54度之间

结论

复现了书上的结果
得到最大值在52到54度之间

致谢

感想吴帆帆同学的指导