cf. David Darmofal. 16.901 Computational Methods in Aerospace Engineering, Spring 2005. (Massachusetts Institute of Technology: MIT OpenCourseWare), http://ocw.mit.edu (Accessed 22 May, 2016). License: Creative Commons BY-NC-SA

Lecture 1 Numerical Integration of Ordinary Differential Equations: An Introduction



In [2]:

    
import sympy
m_p, g = sympy.symbols('m_p g', real=True)



In [3]:

    
t = sympy.Symbol('t',real=True)



In [4]:

    
u = sympy.Function('u')(t)
D = sympy.Function('D')(u)



In [26]:

    
SphFreeFall = sympy.Eq(m_p*sympy.diff(u,t),m_p*g - D)
sympy.pprint( SphFreeFall )









    



    d                         
m_p⋅──(u(t)) = g⋅m_p - D(u(t))
    dt

For low speeds,



In [18]:

    
rho_g, mu_g = sympy.symbols("rho_g mu_g",real=True) 
# density and dynamic viscosity of the atmosphere,



In [19]:

    
a = sympy.symbols("a",real=True)



In [20]:

    
Re = 2*rho_g*u*a/mu_g



In [23]:

    
C_D = 24./Re + 6./(1. + sympy.sqrt(Re)) + 0.4



In [24]:

    
Dlow = 0.5*rho_g*sympy.pi*a**2*u**2*C_D



In [25]:

    
Dlow









    Out[25]:





0.5*pi*a**2*rho_g*(0.4 + 6.0/(sqrt(2)*sqrt(a*rho_g*u(t)/mu_g) + 1.0) + 12.0*mu_g/(a*rho_g*u(t)))*u(t)**2



In [28]:

    
sympy.pprint( SphFreeFall.subs( D, Dlow) )









    



    d                   2     ⎛                 6.0               12.0⋅μ_g ⎞  
m_p⋅──(u(t)) = - 0.5⋅π⋅a ⋅ρ_g⋅⎜0.4 + ───────────────────────── + ──────────⎟⋅u
    dt                        ⎜             ____________         a⋅ρ_g⋅u(t)⎟  
                              ⎜            ╱ a⋅ρ_g⋅u(t)                    ⎟  
                              ⎜      √2⋅  ╱  ──────────  + 1.0             ⎟  
                              ⎝         ╲╱      μ_g                        ⎠  

2           
 (t) + g⋅m_p



In [ ]: