In [1]:
import pandas as pd
import numpy as np
import matplotlib as plt
import math
from IPython.display import Image, HTML, display
from glob import glob

Basic finner

The basic finner or Army-Navy generic finner missile (ANF) considers a generic missile configuration for which longitudinal stability coefficients are available for Mach numbers in a range of 0.5-4.5

References

  1. Allen et al. (2017) Forced motions design for aerodynamic identification and modelling of a generic missile configuration.
  2. Bhagwandin and Sahu (2013) Numerical prediction of pitch damping stability derivatives for finned projectiles.

Below, a schematic representation of the generic configuration is shown:

In [2]:
imagesList=''.join( ["<img style='width: 410px; margin: 10px; float: left; border: 0px solid black;' src='%s' />" % str(s) 
                 for s in sorted(glob('screenShot/geo/anfSchem.png')) ])
display(HTML(imagesList))

Initial conditions

Assuming the ideal gas approximation, the material properties are

In [3]:
mP = pd.read_csv('./input/matProp.csv')
mP
Out[3]:
Gas constant Specific heat ratio
0 287.058 1.4

and the specified freestream inlet conditions are

In [4]:
iCI = pd.read_csv('./input/intialCond.csv')
iCI
Out[4]:
Static pressure Static temperature Mach number
0 101325 293.15 0.5
1 101325 293.15 0.9
2 101325 293.15 1.1
3 101325 293.15 2.5
4 101325 293.15 4.5

Experimental results for the static aerodynamic coefficients, namely, the axial force, the normal force slope and pitching moment slope at 0$^\circ$, are available.

Grid

For the purpose of the analysis, an unstructured, cut-cell Cartesian mesh with approximately 2 million cells is created using the hexPress utility. The viscous sublayer is resolved and an expansion ratio of 1.2 is used in the boundary layer. To resolve the sublayer a non-dimensional first cell height or $y^+$ of less than 1 is required. The farfield boundaries are placed 50 characteristic lengths away. The surface mesh is shown below

In [5]:
imagesList=''.join( ["<img style='width: 650px; margin: 10px; float: left; border: 0px solid black;' src='%s' />" % str(s) 
                 for s in sorted(glob('screenShot/mesh/mesh*.png')) ])
display(HTML(imagesList))

Turbulence model

For the analysis the "Standard"-Menter SST implementation is prescibed. If a turbulent intensity, $I$, of 1% and mixinglength, $L$, of 1x10$^{-5}$ m are prescibed, the corresponding farfield $k$ and $\omega$ values can be calculated from the relations $$ k = \frac{3}{2} (I |U|)^2$$ $$ \omega = \frac{\sqrt{k}}{C_{\mu}^{0.25} L}$$ where $C_{\mu} = 0.09$

In [6]:
# Ideal gas density
# rho = iCI.iloc[0][0]/(mP.iloc[0][0]*iCI.iloc[0][1])
c = math.sqrt(mP.iloc[0][1]*mP.iloc[0][0]*iCI.iloc[0][1])
for index, row in iCI.iterrows():
    M = iCI.iloc[index][2]
    u = M*c
    k = 3/2*(0.01*u)**2
    omega = math.sqrt(k)/(math.pow(0.9,0.25)*1e-5)
    print("For Mach = {} -> k = {:.1f} and omega = {:.1f}".format(M,k,omega))
For Mach = 0.5 -> k = 4.4 and omega = 215798.7
For Mach = 0.9 -> k = 14.3 and omega = 388437.6
For Mach = 1.1 -> k = 21.4 and omega = 474757.1
For Mach = 2.5 -> k = 110.4 and omega = 1078993.4
For Mach = 4.5 -> k = 357.9 and omega = 1942188.2

For the non-slip wall the following wall functions are selected: kqRWallFunction, omegaWallFunction, nutkWallFunction and alphatWallFunction

To non-dimensionalise the integrated aerodynamic forces the following relations are employed to compute the body normal force $$ c_n = \frac{F_N}{\frac{1}{2} \rho u_i u_i A}$$ and pitching moment $$ c_{m} = \frac{M_m}{\frac{1}{2} \rho u_i u_i A c}$$ where $A$ and $c$ are the reference area and the centre of rotation respectively.

In [7]:
rV = pd.read_csv('./input/refValues.csv')
rV
Out[7]:
Length Area Centre of rotation
0 0.03 0.000707 0.165

The freestream viscosity is computed using Sutherland's Law $$\mu = \mu_{0} \left( \frac{T}{T_{0}} \right)^{3/2}\frac{T_{0} + S}{T + S}$$

In [8]:
# Sutherland
mu0 = 1.716e-5
T0  = 273.15
S   = 110.4
C1  = 1.458e-6
mu  = mu0*np.power(iCI.iloc[0][1]/T0,1.5)*(T0 + S)/(iCI.iloc[0][1] + S)
# Ideal gas density
rho = iCI.iloc[0][0]/(mP.iloc[0][0]*iCI.iloc[0][1])
# print(rho)
# Kinematic viscosity
nu  = mu/rho

The corresponding Reynolds number is computed using the relation $$\mathrm{Re} = \frac{\rho u L}{\mu}$$ from which the first cell height can be computed using the $y^+$ approximation. When wall functions are employed, a $y^+$ between 30 and 150 is required.

In [9]:
c = math.sqrt(mP.iloc[0][1]*mP.iloc[0][0]*iCI.iloc[0][1])
for index, row in iCI.iterrows():
    M = iCI.iloc[index][2]
    u = M*c
#     Re = rho*u*rV.iloc[0][0]/mu
    len = 0.3
    Re = rho*u*len/mu
    if (Re > 1e9):
        print ("******************************************************************************")
        print ("WARNING: The Schlichting skin-friction correlation is only valid for Re < 10e9")
        print ("******************************************************************************")
    # Schlichting skin-friction correlation
    Cf    = np.power(2*np.log10(Re)-0.65,-2.3)
    # Wall shear stress
    tauW  = 0.5*Cf*rho*u*u
    # Friction velocity
    uStar = np.sqrt(tauW/rho)
    # Wall distance
    yPlus = 1
    y     = yPlus*mu/(rho*uStar) 
    print("For Mach = {} -> Re = {:.1e} and y(y+={}) = {:.1e} m".format(M,Re,yPlus,y))
For Mach = 0.5 -> Re = 3.4e+06 and y(y+=1) = 2.2e-06 m
For Mach = 0.9 -> Re = 6.2e+06 and y(y+=1) = 1.3e-06 m
For Mach = 1.1 -> Re = 7.5e+06 and y(y+=1) = 1.1e-06 m
For Mach = 2.5 -> Re = 1.7e+07 and y(y+=1) = 5.1e-07 m
For Mach = 4.5 -> Re = 3.1e+07 and y(y+=1) = 2.9e-07 m
In [ ]:
 

Results

Below, the normal force and pitching moment slopes are compared to numerical results from the Army Research Lab as well as experimental data from Defence Research and Development Canada (DRDC). The slopes for both static wind tunnel tests (WT) and free flight (FF) are available.

In [10]:
from IPython.display import Image, HTML, display
from glob import glob
imagesList=''.join( ["<img style='width: 650px; margin: 5px; float: left; border: 0px solid black;' src='%s' />" % str(s) 
                 for s in sorted(glob('screenShot/results/aero*.png')) ])
display(HTML(imagesList))
In [11]:
from IPython.core.display import HTML

def css_styling():
    styles = open("./styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
Out[11]: