第五次作业

---

作业内容

考虑风阻，以及空气密度随高度增加而减小这两个因素，绘制出导弹的轨迹（Exercise2.6）。

思路

欧勒法：

坐标

$$x_{i+1}=x_i+v_{xi}Δt$$
$$y_{i+1}=y_i+v_{yi}Δt$$
阻力
$$F_{drag,x}=-(B_2/m)vv_{xi}$$
$$F_{drag,y}=-(B_2/m)vv_{yi}$$
其中$v=\sqrt{v_{xi}^2+v_{yi}^2}$
速度
$$v_{x,i+1}=v_{x,i}+F_{drag,x}Δt$$
$$v_{y,i+1}=v_{y,i}+(-g+F_{drag,x})Δt$$
其中，$B_2$正比于$\rho$，取$g=9.8m/s^2$,炮弹发射处$B_{2,0}/m=4*10^{-5}m^{-1}.$

---

等温过程

$$\rho=\rho_0e^{-y/y_0}$$ ,故 $$B_2=B_{2,0}e^{-y/y_0}$$
其中$y$为竖直方向位移。取$y_0=k_BT/mg=1.0*10^4m.$

代码

import math
import pylab as pl
class cannon_shell():
    def __init__(self, velocity=700, fixing_angle=9*math.pi/36, position_x=0,
            position_y=0,drag_constant=0.00004,       #drag_constant is B2/m
            acceleration_of_gravity=9.8, decay_index=10000, time_step=0.1):
                self.theta=fixing_angle*180/math.pi
                self.x = [position_x]
                self.y = [position_y]
                self.v = [velocity]
                self.v_x = [velocity*math.cos(fixing_angle)]
                self.v_y = [velocity*math.sin(fixing_angle)]
                self.B2_m = [drag_constant]
                self.g = acceleration_of_gravity
                self.y_0 = decay_index
                self.dt = time_step
                self.t = [0]
    def run(self):
        while(self.y[-1]>=0):
            self.x.append(self.x[-1]+self.v_x[-1]*self.dt/1000)
            self.y.append(self.y[-1]+self.v_y[-1]*self.dt/1000)
            self.v_x.append(self.v_x[-1]-self.B2_m[-1]*self.v[-1]*self.v_x[-1]*self.dt)
            self.v_y.append(self.v_y[-1]-self.g*self.dt-self.B2_m[-1]*self.v_y[-1]*self.v[-1]*self.dt)
            self.B2_m.append(self.B2_m[-1]*math.exp(-self.y[-1]/self.y_0))
            self.v.append(math.sqrt(pow(self.v_x[-1],2)+pow(self.v_y[-1],2)))
            self.t.append(self.t[-1]+self.dt)
    def show_results(self):
        plot1=pl.plot(self.x,self.y)
        pl.title('cannon shell trajectory')
        pl.xlabel('x($km$)')
        pl.ylabel('y($km$)')
        pl.xlim(0,25)
        pl.ylim(0,10)
        pl.show
a = cannon_shell()
a.run()
a.show_results()

若不考虑空气密度随高度的变化，则程序将简化：

import math
import pylab as pl
class cannon_shell():
    def __init__(self, velocity=700, fixing_angle=9*math.pi/36, position_x=0,
            position_y=0,drag_constant=0.00004,       #drag_constant is B2/m
            acceleration_of_gravity=9.8, time_step=0.1):
                self.theta=fixing_angle*180/math.pi
                self.x = [position_x]
                self.y = [position_y]
                self.v = [velocity]
                self.v_x = [velocity*math.cos(fixing_angle)]
                self.v_y = [velocity*math.sin(fixing_angle)]
                self.B2_m = drag_constant
                self.g = acceleration_of_gravity
                self.dt = time_step
                self.t = [0]
    def run(self):
        while(self.y[-1]>=0):
            self.x.append(self.x[-1]+self.v_x[-1]*self.dt/1000)
            self.y.append(self.y[-1]+self.v_y[-1]*self.dt/1000)
            self.v_x.append(self.v_x[-1]-self.B2_m*self.v[-1]*self.v_x[-1]*self.dt)
            self.v_y.append(self.v_y[-1]-self.g*self.dt-self.B2_m*self.v_y[-1]*self.v[-1]*self.dt)
            self.v.append(math.sqrt(pow(self.v_x[-1],2)+pow(self.v_y[-1],2)))
            self.t.append(self.t[-1]+self.dt)
    def show_results(self):
        plot1=pl.plot(self.x,self.y,'--')
        pl.title('cannon shell trajectory')
        pl.xlabel('x($km$)')
        pl.ylabel('y($km$)')
        pl.xlim(0,25)
        pl.ylim(0,10)
        pl.show
a = cannon_shell()
a.run()
a.show_results()

运行结果

1. 考虑密度随高度变化与不考虑的定性比较：（虚线为未考虑，实线为考虑；初速度为700m/s，角度为45°）
   结果图1
2. 考虑空气密度随高度变化的情况下，改变炮弹发射角度的值，在同一坐标中作图：（速度不变，角度从30°到65°，每隔5°取一个次）
   结果图2

结论

1. 考虑空气密度随高度的变化，则炮弹将比不考虑时飞得更高更远。
2. 所选的角度值中，40°飞行距离最远。