2.1.1 Skin frication drag coefficient estimation

The skin friction drag coefficient (C _D)_sf is calculated using the following formula.

(1) \begin{align} \left( {{C_D}_{SF}} \right) = {C_f} \times \frac{{{S_{wetted}}}}{{{S_{referenc}}}} \end{align}

The reference area of the shell projectile is 0.0189m². For calculating the wetted surface area of the projectile, the model is divided into four sections. The Sections 1, 2 and 4 are considered as truncated cone, and the Section 3 is considered as cylinder respectively as shown in Fig. 3. The formula for a truncated cone lateral surface area is $F = \pi\!\left( {{r_{1{\rm{\;}}}} + {r_2}} \right)\sqrt {{{\left( {{r_{1{\rm{\;}}}} + {r_2}} \right)}^2} + {h^2}} $ and surface area is $F + \pi\!\left( {{r_{1{\rm{\;}}}} + {r_2}} \right)$ . Table 2 shows the surface area for all segmentation as per the design concerned.

Figure 3.

Segmentation of shell projectile.

Table 2.

Surface area for each segmentation

The total wetted surface area is calculated from addition of four surface areas ${\rm{\;}}{S_{wetted}} = 0.264{{\rm{m}}^2}$ .

Estimated wetted surface area by using CAD software SolidWorks is ${S_{wetted}} = 0.2763{\rm{\;}}{{\rm{m}}^2}\;$ (surface wetted area is very near to the theoretical value).

(2) \begin{align} {C_f} = {\rm{\;}}{\left( {{C_f}_{turbulent}} \right)_c} + \frac{{{x_1}}}{c}{\left( {{C_f}_{laminar}} \right)_{x1}} - \frac{{{x_1}}}{c}{\left( {{C_f}_{turbulent}} \right)_{x1}} \end{align}

Consider flow (M = 0.7) over shell projectile with a boundary layer extending from the leading edge over the distance ${x_1}$ to the transition point and a turbulent boundary extending over the distance ${x_1}$ from the transition point to the trailing edge. The critical Reynolds number is 5 × 10⁵, Reynolds number based on projectile length is 1.1 × 10⁷, (where c = 0.701m (Length of projectile). This locates the transition point relative to the chord length. So, transition length from leading edge is ${x_1}$ = 0.03m and from trailing edge is ${x_2}$ = 0.671m. For laminar flow, with Reynolds number based on ${x_1}$ ,

(3) \begin{align} {\left( {{C_f}_{laminar}} \right)_{x1}} = \frac{{1.328}}{{\sqrt {R{e_{x1}}} }}{\left( {1 + 0.12M_\infty ^2} \right)^{ - 0.12}} = 0.0018 \end{align}

For turbulent flow, with Reynolds number based on ${x_1}$ ,

(4) \begin{align} {\left( {{C_f}_{turbulent}} \right)_{x1}} = \frac{{0.455}}{{{{\left( {lo{g_{10}}R{e_{x1}}} \right)}^{2.58}}}}{\left( {1 + 0.21M_\infty ^2} \right)^{ - 0.32}} = 0.0049 \end{align}

For turbulent flow, with Reynolds number based on chord C,

(5) \begin{align} {\left( {{C_f}_{turbulent}} \right)_c} = \frac{{0.455}}{{{{\left( {lo{g_{10}}R{e_{x1}}} \right)}^{2.58}}}}{\left( {1 + 0.21M_\infty ^2} \right)^{ - 0.32}} = 0.0028 \end{align}

At, ${M_\infty }$ = 0.7 skin friction coefficient ${C_f}$ is 0.0026; skin friction drag coefficient ${\left( {{C_D}} \right)_{SF}}$ is 0.0394, as shown in Table 3.

2.1.2 Wave drag coefficient estimation

The bluntness ratio $'b'$ is the ratio between diameter of spherical blunting $'d'\;$ to maximum diameter of body $'D'$ as given in the following Equation (Reference Mirzaei, Karimi and Vaziri6).

(6) \begin{align} {\rm{\;}}b = \frac{d}{{D{\rm{\;}}}}{\rm{\;}} = {\rm{\;}}0.087 \end{align}

Table 3.

Coefficient of skin friction drag at various Mach number

Table 4.

Coefficient of wave drag at various Mach number

Forebody fineness ratio is the ratio between lengths of forebody ${l_f}$ to maximum diameter of body $'D'$ .

(7) \begin{align} f = \frac{{{l_f}}}{{D{\rm{\;}}}}{\rm{\;}} = {\rm{\;}}2.44 \end{align}

ESDU data Ref [20] gives wave drag coefficient for the various bluntness ratio and fineness ratio. But plots are available only for forebody fineness ratio of = 2 & 3. Therefore, to get the ${C_{DW}}$ values for interpolation technique was used. The calculated values of wave drag coefficient ( ${C_{DW}}$ ) are presented in the following Table 4.

Plots are available only for forebody fineness ratio 2 and 3, so getting ${C_D}$ values for $f = 2\; \& \;f = 3$ then interpolation both values for getting $\,f = 2.44$ values.

2.1.3. Base-drag coefficient estimation

The base-drag coefficients are estimated as shown in Table 5, using diameter ratio $\left(\dfrac{{{D_B}}}{{{D_M}}}\right)$ and boattail angle (β) using the ESDU data. Diameter ratio is ratio between base diameters ‘D_B’ to maximum diameter of body $'D'$ of the shell configuration. Boattail angle of M107 155mm shell projectile is β = 8^o. The ESDU data provides the total after body pressure drag coefficient ${C_D}_a$ . For subsonic and supersonic regions provides ${C_D}_a$ directly by using diameter ratio $\left(\dfrac{{{D_B}}}{{{D_M}}}\right)$ and boattail angle (β). For transonic region ${C_D}_a$ is sum of base pressure drag coefficient ${C_D}_b$ and boattail pressure drag coefficient ${C_D}_\beta $ .

(8) \begin{align} {C_D}_a{\rm{\;}} = {\rm{\;}}{C_D}_b + {C_D}_\beta \end{align}

Table 5.

Coefficient of base drag at various Mach numbers

Table 6 shows the cumulative result for shell projectile’s drag calculation, in terms of skin friction drag, base drag, wave drag and the total drag.

Table 6.

Coefficient of total drag at various freestream Mach numbers

Figure 4.

Computational domain.

2.2.1 Domain selection

Domain simplification is an extremely normal methodology and whenever utilised wisely, gives precise and quicker outcomes. The domain wrapping the projectile guarantees to embody the arrangement of shocks at the forebody and shedding the vortex at the rear. Figure 4(a) shows the projectile domain as used in the simulation software. Frontal section till the band is termed as forebody, followed by the rotating band, bleed intake with IWTB hole, boattail and base cavity. Figure 4(b) shows axisymmetric model domain (for revolution bodies) with upstream, top and bottom boundary of 3D. The downstream domain is extended up to 7D. The axisymmetric model is used to study for freestream Mach number ranging from subsonic to sonic region. From the numerical simulation, the wake region is observed up to 1.25D and its effect is seen up to 2.5D for all subsonic and sonic simulations. Therefore, the downstream domain selected at 7D is far behind the wake region. The top and bottom domain at 3D is sufficient to capture all details of flow field, because the angle-of-attack is zero degrees for all numerical simulations; moreover the disturbance is not seen beyond 2D in top and bottom surfaces.

3D model domain shown in Fig. 4(c) is used for freestream Mach number 1.5 and 2, to obtain the wake region more prominently. The distance from the inflow to the leading edge of the round, the round base to the outerflow and body wall to the far field are 0·15D, 6D and 11D, respectively, based on Ref. (Reference Jiajan, Chue, Nguyen and Yu5). This domain selection is confirmed and justified based on the shock standoff distance calculation made as given in Refs (Reference Lobb21 and Reference Shankara and Sreekanth22).

The following equation as given in Ref. (Reference Jiajan, Chue, Nguyen and Yu5) predicts the shock standoff distance for a sphere of diameter D. Where the ratio of the freestream density ${\rho _\infty }$ to the after-normal-shock density ${\rho _2}$ on the central stagnation streamline is used in the following relation from Ref [Reference Saiprakash, Senthil Kumar, Kadam Sunil, Singh Prakash, Shanmugam and Balu23]:

(9) \begin{align} \frac{\delta }{D} = 1.41\frac{{{\rho _\infty }{\rm{\;}}}}{{{\rho _2}{\rm{\;\;}}}} \end{align}

(10) \begin{align} \frac{\delta }{D} = 1 - {\left[ {1 - \frac{4}{{{{\left( {\gamma + 1} \right)}^2}}}\!\left( {\frac{1}{{{M^2}}} + \frac{{\gamma - 1}}{2}} \right)} \right]^{1/2}} \end{align}

From the above equations, it is observed that the shock standoff distance for Mach 1.5 is 0.07D and for Mach 2.0 is 0.05D. It is well known that in supersonic stream the disturbance will not propagate in upstream. Hence the selected domain size of 0.15D is more than sufficient and justified. For downstream domain selection, the wake region is seen up to 2.D and the effect of wake is seen up to the region of 5 to 6D. Therefore, the downstream domain is fixed at 11D, which is far behind the wake region.

Table 7 shows the CFD models used for the simulation. Initially 2D axisymmetric models were used to get the faster solution and to optimise the configuration. After optimisation to confirm the results and to compare with the experiment result 3D models were used in the CFD simulation.

Table 7.

Configuration for CFD simulation

2.2.2 Meshing condition

All models were created using SolidWorks, model meshing is carried out using ANSYS Workbench for axi-symmetric geometries and 3D model meshing is done using ICEM CFD. Figure 5 shows the unstructured extra-fine mesh created using ICEM CFD. Numerical simulation for all the above domain is carried out using ANSYS Fluent. Algebraic multi-grid (AMG) method for convergence of the implicit density-based solver, shear stress transport k-ω model was used. By keeping y⁺ value to 1, first layer thickness is found to be approx. 0.44nm for Mach number 2. Mesh quality for orthogonal quality is 0.269; and maximum aspect ratio is 36.1356. The governing equations are solved in the fluent software compressible fluid flow are as follows:

Figure 5.

3D model mesh in ICEM CFD.

Continuity Equation:

(11) \begin{align} \frac{{\delta \rho }}{{\delta t}} + \frac{{\delta \!\left( {\rho {u_i}} \right)}}{{\delta {x_i}}} = 0 \end{align}

Momentum Equation:

(12) \begin{align} \frac{{\delta \!\left( {\rho {u_i}} \right)}}{{\delta t}} + \frac{{\delta \!\left( {\rho {u_i}{u_j}} \right)}}{{\delta {x_j}}} = - \frac{{\delta P}}{{\delta {x_i}}} + \frac{\delta }{{\delta {x_i}}}{\tau _{ij}} \end{align}

Energy Equation:

(13) \begin{align} \frac{\delta }{{\delta t}}\left[ {\rho \!\left( {e + \frac{{{V^2}}}{2}} \right)} \right] + \frac{\delta }{{\delta {x_i}}}\left[ {\rho {u_{j{\rm{\!}}}}\left( {e + \frac{{{V^2}}}{2}} \right) + P + {q_j} - {u_i}{\tau _{ij}}} \right] = 0 \end{align}

Note that $\;\;u\;$ denotes instantaneous velocity, $V\;-$ the velocity modulus $,\;\rho $ – the gas density $,\;P$ – the gas pressure, ${q_j}\;$ – the heat flux, and ${\tau _{ij\;}}$ – the viscous stress tensor. In addition, the perfect gas equation of state was considered, the two equations mentioned in Menter’s k- $\omega$ SST model Ref. (24) is given below

(14) \begin{align} \frac{{\delta k}}{{\delta t}} + {u_i}\frac{{\delta k}}{{\delta {x_j}}} = \frac{1}{\rho }{p_k}\!\left( {\frac{{{M_\infty }}}{{Re}}} \right) - {\beta ^/}{\rm{k\omega }}\!\left( {\frac{{Re}}{{{M_\infty }}}} \right) + \frac{1}{\rho }\frac{\delta }{{\delta {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _T}}}{{{\sigma _k}}}} \right)\frac{{\delta k}}{{\delta {x_j}}}} \right]\left( {\frac{{{M_\infty }}}{{Re}}} \right) \end{align}

(15) \begin{align} \frac{{\delta \omega }}{{\delta t}} + {u_j}\frac{{\delta k}}{{\delta {x_j}}} & = \frac{1}{\rho }{p_\omega }\!\left( {\frac{{{M_\infty }}}{{Re}}} \right) - \beta {{\rm{\omega }}^2}\!\left( {\frac{{Re}}{{{M_\infty }}}} \right) + \frac{1}{\rho }\frac{\delta }{{\delta {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _T}}}{{{\sigma _\omega }}}} \right)\frac{{\delta \omega }}{{\delta {x_j}}}} \right]\left( {\frac{{{M_\infty }}}{{Re}}} \right) \nonumber \\ & \quad + 2\!\left( {1 - F} \right)\frac{1}{{{\sigma _\omega }_2}}\frac{1}{\omega }\frac{{\delta k}}{{\delta {x_j}}}\frac{{\delta \omega }}{{\delta {x_j}}}\!\left( {\frac{{{M_\infty }}}{{Re}}} \right) \end{align}

Where, ${\mu _T} = \min\!\left[ {\frac{{\rho k}}{\omega },\frac{{{a_1}\rho k}}{{\Omega {F_2}}}\left( {\frac{{{\mathop{\rm Re}\nolimits} }}{M}} \right)} \right]$ ; $p{}_k = {\mu _T}{\Omega ^2}$ ; $p{}_\omega = \gamma \rho {\Omega ^2}$ ; ${C_\mu } = {\beta ^/} = 0.09$ ; ${F_1} = \tanh \!\left( {{\Gamma ^4}} \right)$ ; $\Gamma = \min \!\left[ {\max\!\left( {{\Gamma _1},{\Gamma _3}} \right),{\Gamma _2}} \right]$ ; ${F_2} = \tanh\!\left( {{\Pi ^2}} \right)$ ; $\Pi = \max \!\left( {2{\Gamma _3},{\Gamma _1}} \right)$ ; ${\Gamma _1} = \frac{{500v}}{{{d^2}\omega }}{\left( {\frac{{{M_\infty }}}{{{\mathop{\rm Re}\nolimits} }}} \right)^2}$ ; ${\Gamma _2} = \frac{{4\rho k}}{{{d^2}{\sigma _{{\omega _2}}}\!\left( {C{D_{k - \omega }}} \right)}}$ ; ${\Gamma _3} = \frac{{\sqrt k }}{{{C_\mu }d\omega }}{\left( {\frac{{{M_\infty }}}{{{\mathop{\rm Re}\nolimits} }}} \right)^2}$ ; $C{D_{k - \omega }} = \max\! \left( {\rho \frac{2}{{{\sigma _{{\omega _2}}}\omega }}\frac{{\partial k}}{{\partial {x_j}}}\frac{{\partial \omega }}{{\partial {x_j}}},1 \times {{10}^{20}}} \right)$ .

Table 8.

GCI on the fine grid

Table 8 shows the grid convergence index $\left( {GCI} \right)\;$ study for the model Ref. (Reference Roache25). $GCI$ is a standardised way to report grid convergence quality. Following are the summary of the equations and steps that are followed during the grid convergence study:

1. 1. By keeping constant refinement ratio, i.e. $r\; = \;1.2$ , the coarse, medium and fine meshes were run in the ansys software. The coefficient of base drag is used to indicate grid convergence.

2. 2. The order of convergence, $p$ is explained in Equation (16), in terms of ${f_1},{\rm{\;}}{f_2}{\rm{\;and\;}}{f_3}.$

(16) \begin{align} p{\rm{\;}} = {\rm{\;}}\frac{{ln\!\left( {\frac{{{f_3} - {f_2}}}{{{f_2} - {f_1}}}} \right)}}{{ln{\rm{\;}}r}} \end{align}

1. 3. Richardson extrapolation method is used to predict the value at $H = 0$ , Where ${f_{H = 0{\rm{\;}}}}$ is the continuum value at zero grid spacing.

(17) \begin{align} {f_{H = 0{\rm{\;}}}} \cong {f_1} + {\rm{\;}}\frac{{{f_1} - {f_2}}}{{{r^p} - 1}} \end{align}

1. 4. The factor of safety is recommended to be $Fs = 3.0$ for comparisons of two grids and $Fs = 1.25$ for comparisons over three or more grids. Hence, $Fs$ considered here is 1.25. Then grid convergence index $\left( {GCI} \right)$ for the medium and fine refinement levels is,

(18) \begin{align} GC{I_{Fine}}{\rm{\;}} = {\rm{\;}}\frac{{{F_s}\left| {{\delta _{ij}}} \right|}}{{\left( {{r^p} - 1} \right)}} \end{align}

1. 5. It should be ensured that grids are in the asymptotic range of convergence by checking.

(19) \begin{align} GCI\; = \;\frac{{GC{I_{23}}}}{{{r^p}GC{I_{12}}}} \end{align}

Grid independence study was based on $GCI$ criteria. $GC{I_{23}}{\approx}{r^p}\;GC{I_{12}}$ . Based on the iteration values extracted from software, the solutions were checked to find whether they fall under the asymptotic range of convergence. To check the solution, Equation (19) was taken into account. The values are plotted in the Fig. 6. Overall GCI for k-ω SST model is closer to 1 as shown in Table 8, therefore for further analysis extra-fine mesh is considered.

Figure 6.

Detailed grid study.