Ex_06 2.10 Consider cases in which the target is higher and lower than the cannon.

author: Liu Yongjie 2014301020094

---

The title

2.10 Generalize the program developed for the previous problem so that it can deal with situations in which the target is at a different altitude than the cannon. Consider cases in which the target is higher and lower than the cannon. Also investigate how the minimum firing velocity required to hit the target varies as the altitude of the target is varied.

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.

code

import time
import math
import matplotlib.pyplot as plt
def correct(string,y_t):           ##### delete all the data of (x,y) where y<0
    x_record = string[0][-1] 
    y_record = y_t
    while True:
        if (string[1][-1] < y_t):
            if (string[1][-1] < y_t and string[1][-2] > y_t):
                x_record = string[0][-1] 
                y_record = string[1][-1]               
            del string[0][-1]
            del string[1][-1]
        else:
            string[0].append(string[0][-1]+(y_t - string[1][-1])*(string[0][-1] - x_record)/(string[1][-2] - y_record))
            string[1].append(y_t)
            break
    return string[0],string[1]
def calculate(v,theta):
    vx = []
    vy = []
    vx.append( v * math.cos(theta*math.pi/180))
    vy.append( v * math.sin(theta*math.pi/180))
    dt = 0.1
    t = []
    x = []
    y = []
    t.append(0)
    x.append(0)
    y.append(0)
    i = 1
    while (y[-1] >= 0):
        x.append(x[i - 1] + vx[i - 1]*dt)
        vx.append(vx[i - 1] - dt*0.00004*math.sqrt(vx[i - 1]**2+vy[i - 1]**2)*vx[i - 1]*(1-(0.0065*y[i - 1])/280)**2.5)
        y.append(y[i - 1] + vy[i - 1]*dt)
        vy.append(vy[i - 1]-9.8*dt -dt*0.00004*math.sqrt(vx[i - 1]**2+vy[i - 1]**2)*vy[i - 1]*(1-(0.0065*y[i - 1])/280)**2.5)
        t.append(t[i - 1] + dt)
        i+= 1
    r = y[-2] / y[-1]
    x[-1] = (x[-2] + r * x[-1]) / (r + 1)
    y[-1] = 0
    return x, y
def find_maxheight(string):
    for i in range(len(string[1])):
        if ( string[1][0] < string[1][i] ):
            string[1][0] = string[1][i]
    return string[1][0]
def scan_angle(v,theta, x_t, y_t, degree):
    ran_a = [20, 5, 2, 0.8, 0.2]
    delta = [1, 0.5, 0.1, 0.05, 0.01]
    theta = theta - ran_a[degree]
    theta_record = []
    x = []
    data = [[],[]]
    for i in range(int(ran_a[degree]*2/delta[degree]+1)):
        data = calculate(v,theta)[:]
        if ( find_maxheight(data) < y_t):
            theta = theta + delta[degree]
            x.append(100000000)
            theta_record.append(theta)
        else:
            data = correct(data,y_t)[:]
            x.append(data[0][-1])
            theta_record.append(theta)
            theta = theta + delta[degree]
    for j in range(len(x)):
        if ( abs(x[0] - x_t) > abs(x[j] - x_t)):
            x[0] = x[j]
            theta_record[0] = theta_record[j]
    return theta_record[0],x[0]   ##### debug
def scan_v(v,theta, x_target, y_target, degree):
    ran_v = [100, 10, 2, 0.8, 0.2]
    delta = [10, 2, 0.5, 0.1, 0.02]
    v = v - ran_v[degree]
    x = []
    theta_record = []
    v_record = []
    for i in range(int(ran_v[degree]*2/delta[degree]+1)):
        record_v = scan_angle(v,theta, x_target, y_target, degree)[:]
        x.append(record_v[1])
        theta_record.append(record_v[0])
        v_record.append(v)
        v = v + delta[degree]
    for j in range(len(x)):
        if ( abs(x[0] - x_target) > abs(x[j] - x_target)):
            x[0] = x[j]
            v_record[0] = v_record[j]
            theta_record[0] = theta_record[j]
    print v_record
    print theta_record
    print x
    return v_record[0],theta_record[0]
def deep_scan(v, theta, x_target, y_target, precision):
    degree = 0
    record = [[],[]]
    for i in range(precision):
        record = scan_v(v,theta,x_target,y_target,degree)[:]
        v = record[0]
        theta = record[1]
        degree = degree + 1
    return record
def judge_hitting(string,x_t,y_t):
    alpha = 0
    for i in range(len(string[0])):
        if ( x_t - 8 < string[0][-1] < x_t + 8 ):
            alpha = 1
    if (alpha == 1):
        print 'Successfully hitted'
        return 1
    else:
        print 'Miss'
        return 0
def main():
    start = time.clock()
    x_target = 10000
    y_target = 2000
#    theta0 = 180*math.atan(y_target/x_target)/math.pi
    data_record = deep_scan(700,60,x_target,y_target,5)   ### 3 for grade 1, 4 for grade 2 and 5 for grade 3
    v = data_record[0]
    theta = data_record[1]
    cannon_record = correct(calculate(v, theta),y_target)
    x_max = cannon_record[0][-1]
    judge_hitting(cannon_record,x_target,y_target)
    print '%f ' % x_max
    print '%f ' % theta
    print '%f ' %v
    end = time.clock()
    print "read: %f s" % (end - start)
    plt.figure(figsize = (8,6))
    plt.title('Trajectory of cannon shell')
    plt.xlabel('x(m)')
    plt.ylabel('y(m)')
    plt.plot(x_target,y_target,'k*',linewidth = 10,label='Target')
    plt.plot(cannon_record[0],cannon_record[1],label= 'Trajectroy')
    plt.legend(loc = 'upper left')
    plt.savefig('chapter2.png',dpi = 144)
    plt.show()
main()

Click to view the results

Hitting point (10000, 2000), Angle: $75.50^。$, Speed: 610.48m/s

Conclusion

This assignment allows us to learn how to solve high-order differential equation via transferring it into first-order diferential equations. Also, air density at differential altitude plays an important role in the trejectory of cannon shell.

thanks to