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


### Modified basic finner¶

The modified basic finner or Air Force modified finner missile (AFF) considers a generic missile configuration for which longitudinal stability coefficients are available for Mach numbers in a range of 0.6-2.5. The modified finner is a 2.5 caliber tangent-ogive cylinder forebody with 4 trapezoidal fins.

References

1. Bhagwandin and Sahu (2013) Numerical prediction of pitch damping stability derivatives for finned projectiles.
2. Samardzic et al. (2007) Some experimental results of subsonic derivative obtained in the T-38 wind tunnel by forced oscillation.
3. Murman (2005) A reduced-frequency approach for calculating dynamic derivatives.
4. Dunn (1989) Aeropredictive methods for missile analysis

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/affSchem.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.6
1 101325 293.15 0.8
2 101325 293.15 0.9
3 101325 293.15 1.1
4 101325 293.15 1.2
5 101325 293.15 1.5
6 101325 293.15 2.5

Free flight data 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 3.5 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 [12]:
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 mixing length, $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*u
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.6 -> k = 636.2 and omega = 2589584.3
For Mach = 0.8 -> k = 1131.0 and omega = 3452779.0
For Mach = 0.9 -> k = 1431.4 and omega = 3884376.4
For Mach = 1.1 -> k = 2138.3 and omega = 4747571.1
For Mach = 1.2 -> k = 2544.7 and omega = 5179168.5
For Mach = 1.5 -> k = 3976.1 and omega = 6473960.7
For Mach = 2.5 -> k = 11044.8 and omega = 10789934.4


For the non-slip wall the following wall functions were selected: kLowReWallFunction, nutkWallFunction, omegaWallFunction 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.144

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, and to resolve the viscous sublayer the $y^+$ should be less than 1.

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.6 -> Re = 4.1e+06 and y(y+=1) = 1.9e-06 m
For Mach = 0.8 -> Re = 5.5e+06 and y(y+=1) = 1.5e-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 = 1.2 -> Re = 8.2e+06 and y(y+=1) = 1.0e-06 m
For Mach = 1.5 -> Re = 1.0e+07 and y(y+=1) = 8.2e-07 m
For Mach = 2.5 -> Re = 1.7e+07 and y(y+=1) = 5.1e-07 m


### Results¶

Below, the normal force and pitching moment slopes are compared to numerical results from the Army Research Lab (ARL) as well as free flight data from Defence Research and Development Canada (DRDC). The slopes for two test flights SF and MF are available.

In [13]:
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():