Computation of Axisymmetric Turbulent Viscous Flow Around Sphere

Axisymmetric turbulent viscous flow around sphere is computed using finite volume method based on Reynolds-averaged Navier-Stokes (RANS) equations. Two-dimensional axisymmetric flow solver has been used to analyze flow at Reynolds number of 5×10. Spalart-Allmaras (S-A) and shear stress transport (SST) k-ω turbulence models are used to capture turbulent viscous flow. The numerical results in terms of the skin friction coefficient, pressure coefficient and drag coefficient at that Reynolds number have been shown either graphically or in the tabular form and velocity vectors have been displayed graphically. The computed results show good agreement with published experimental measurements.


Introduction
Applications of computational fluid dynamics (CFD) to the maritime industry continue to grow as this advanced technology takes advantage of the increasing speed of computers.In the last two decades, different areas of incompressible flow modeling including grid generation techniques, solution algorithms and turbulence modeling, and computer hardware capabilities have witnessed tremendous development.In view of these developments, computational fluid dynamics (CFD) can offer a cost-effective solution to many problems in underwater bodies.However, effective utilization of CFD for marine hydrodynamics depends on proper selection of turbulence model, grid generation and boundary resolution.
Turbulence modeling is still a necessity as even with the emergence of high performance computing since analysis of complex flows by direct numerical simulations (DNS) is untenable.The peer approach, the large-Eddy simulation (LES), still remains expensive.Hence, simulation of underwater hydrodynamics continues to be based on the solution of the Reynolds-averaged Navier-Stokes (RANS) equations.Various researchers used turbulence modeling to simulate flow around axisymmetric bodies since late seventies.Patel and Chen [1] made an extensive review of the simulation of flow past axisymmeric bodies.Choi and Chen [2] gave calculation method for the solution of RANS equation, together with k-ε turbulence model.Sarkar et al. [3] used a low-Re k-ε model of Lam and Bremhorst [4] for simulation of flow past underwater axisymmetric bodies.
A considerable amount of research work has been published on flow around the axisymmetric sphere.The basic structure of the flow past a sphere has been experimentally investigated using a variety of approaches, including flow visualization by Achenbach [5], Taneda [6], Bakic [7] etc.Recent time-accurate computations of laminar and turbulent flow around spheres using different methods are reported by many researchers, among them the work of Gregory [8], Kalro [9] and Sun and Chwang [10] are remarkable.In this present study, Spalart-Allmaras (S-A) and shear stress transport (SST) k-ω turbulence models are used to simulate fully turbulent flow past underwater axisymmetric sphere.

Theoretical formulation
For the incompressible flow past an axisymmetric underwater body, the continuity equation in cylindrical co-ordinate is given by: where x is the axial coordinate, r is the radial coordinate, u is the axial velocity and v is the radial velocity.The source term S m is the mass added to the continuous phase from the dispersed second phase and any user-defined sources.S m is taken as zero here since single phase is considered.Also, the axial and radial momentum equations are given by: where, p = static pressure, μ = molecular viscosity, ρ = density, F x and F r are external body forces and . F x and F r are taken as zero since no external force is considered in this study.

The Spalart-Allmaras (S-A)
The Spalart-Allmaras turbulence model used in this study is a simple one-equation model [11] that solves a modeled transport equation for the turbulent viscosity.This model is designed for wall-bounded flows and gives good results for boundary layers subjected to adverse pressure gradients, much like the flow fields encountered in this study.Although the original Spalart-Allmaras model requires that the viscous-affected region of the boundary layer be properly resolved through the use of a fine mesh inside the boundary layer, the model has been modified for its implementation so that wall functions are used when the mesh resolution is not sufficiently fine near object surfaces.The fact that the S-A model is a one-equation model with relatively lax grid density requirements further enhances its suitability for this particular study since, for the computer platform used, maximum computational efficiency is critical.The transported variable in the Spalart-Allmaras model, − ν , is identical to the turbulent kinem tic viscosity except in the near-wall (viscous-affected) region.The transport equation for where, G ν is the production of turbulent viscosity and Y ν is the destruction of turbulent viscosity that occurs in the near-wall region due to wall blocking and viscous damping.σ ν and C b2 are constants and v is the molecular kinematic viscosity.S ν is a user-defined source term.Note that the turbulence kinetic energy k is not calculated in the Spalart-Allmaras model.
To obtain the modified turbulent viscosity, ν, for the Spalart-Allmaras model from the turbulence intensity, I and length scale, l, the following equation can be used: where, (here, L = sphere diameter) In this model the constants are considered as:

The shear-stress transport (SST) k-ω model
The SST k-ω turbulence model is a two-equation eddy-viscosity model developed by Menter [12] to effectively blend the robust and accurate formulation of the k-ω model in the near-wall region with the free-stream independence of the k-e model in the far field.To achieve this, the k-e model is converted into a k-ω formulation.The SST k-ω model is similar to the standard k-ω model, but includes the following refinements: These features make the SST k-ω model more accurate and reliable for a wider class of flows (e.g., adverse pressure gradient flows, airfoils, transonic shock waves) than the standard k-ω model.
The shear-stress transport (SST) k-ω model is so named because the definition of the turbulent viscosity is modified to account for the transport of the principal turbulent shear stress.The use of a k-ω formulation in the inner parts of the boundary layer [13] makes the model directly usable all the way down to the wall through the visous sub-layer, hence the SST k-ω model can be used as a Low-Re turbulence model without any extra damping functions.The SST formulation also switches to a k-ε behaviour in the free-stream and thereby avoids the common k-ω problem that the model is too sensitive to the inlet freestream turbulence properties.It is this feature that gives the SST k-ω model an advantage in terms of performance over both the standard k-ω model and the standard k-ε model.Other modifications include the addition of a cross-diffusion term in the ω equation and a blending function to ensure that the model equations behave appropriately in both the near-wall and far-field zones.
Transport equations for the SST k-ω model are given by: In these equations,

Boundary conditions
Since the geometry of an axisymmetric body is, in effect, a half body section rotated about an axis parallel to the free stream velocity, the bottom boundary of the domain is modeled as an axis boundary.Additionally, the left and top boundaries of the domain are modeled as 'velocity inlet', the right boundary is modeled as an 'outflow boundary', and the surface of the body itself is modeled as a 'wall'.

Viscous drag
The viscous drag of a body is generally derivable from the boundary-layer flow either on the basis of the local forces acting on the surface of the body or on the basis of the velocity profile of the wake far downstream.The local hydrodynamic force on a unit of surface area is resolvable into a surface shearing stress or local skin friction tangent to the body surface and a pressure p normal to the surface.The summation over the whole body surface of the axial components of the local skin friction and of the pressure gives, respectively, the skin-friction drag D f and the pressure drag where, r w is the radius from the axis to the body surface, α is the arc length along the meridian profile, and x e is the total arc length of the body from nose to tail.The sum of the two drags then constitutes the total viscous drag, D or D=D f +D p The drag coefficient, C D and the pressure coefficient, C p based on some appropriate reference area A are given by: where, p ∞ is pressure of free stream and U ∞ is free stream velocity.

Methodology
The computational domain is extended ten times the sphere diameter in fore and aft of the sphere.The region also extended to ten times the sphere diameter in the vertical direction from the edge.It is ensured that in the selected model, the numerical results would be accurate and that the problem would be solvable in a reasonable amount of time.
For the purposes of grid construction, the computational domain for sphere model is divided into two regions: the boundary layer region and the free stream region.Dividing the domain in this fashion is a common practice in problems where the effects of the viscous boundary layer that forms on the body are expected to significantly affect the flow field and where enhanced grid resolution in the vicinity of the boundary layer is important.The boundary layers are attached to the spheres and the direction of the boundary layer grids is defined such that the grids extended into the interior of the domains.Based on prior experience with numerical simulations involving boundary layers and the expected growth of the boundary layer meridionally along the sphere, both boundary layer meshes are approximately 3 cm in height.Increasing the number of rows in the boundary layer meshes only served to vary cell density, and did not change the total height of the mesh.Finally, the growth factors are chosen to increase the resolution of the meshes at the base of the boundary layers (where flow parameter gradients are largest) while still maintaining high grid resolution, low cell skewness at the top of the boundary layers, and a total boundary layer mesh thickness of approximately 3 cm.Low skewness is important to ensure similar cell proportions between outer boundary layer cells and neighboring free stream region cells.The boundary layer grid parameters for the axisymmetric sphere models are shown in Table 1.If the growth factor is not listed in the following tables, it would be considered as unity.Meshing of the free stream regions took place in two steps.First, the edges of the regions are meshed, and then, using the edge meshes, the interiors of the regions (or faces) are meshed.Since boundary layer meshing has already been performed, only the axis boundary, inlet, outlet, and top edges have to be meshed.Comparatively coarse meshes are specified on the exterior (inlet, outlet, and top) boundaries due to the lack of large flow property fluctuations (and thus low grid densities) in those regions.For better control of edge node spacing, the bottom boundary is constructed in multiple sections.Grading is necessary to ensure a smooth transition between the relatively small cell sizes near the boundary layer grids and the relatively large cell sizes on the outer edges of the domains.Table 2 shows the node spacing on the edges of the domains for each edge node distribution.For the purpose of grid construction, the computational domain is divided into three faces: Middle face, Front face and Rear face.At first, the edges of the faces are meshed, and then, using the edge meshes, the interiors of the faces are meshed.The node spacing on the edges of the domain for each node distribution is given in the Table 3.
Once the edges are meshed, the interior of the domains need to be meshed using automatic face mesh generation scheme.The meshing scheme that is chosen is pave meshing scheme.The pave scheme creates an unstructured grid of mesh elements, which is particularly desirable for its applicability to a wide range of face geometries, its ability to deal with irregularly shaped interiors, and its ease of use.There is no restriction on mesh node spacing imposed by the pave scheme since only triangular face elements are used.More cells are constructed near the surface of the sphere to tackle the high velocity gradient in the boundary layer region of the viscous flow.Fig. 1(a) shows the grid for the axisymmetric sphere, which is symmetric about the axis of rotation.Also, Fig. 1

Sphere (Wall) Axis
A finite volume method [14,15] is employed using a commercial software FLUENT to obtain a solution of the Reynolds averaged Navier-Stokes equations.The coupling between the pressure and velocity fields is achieved using PISO algorithm [14].A second order upwind scheme is used for the convection and the central-differencing scheme for diffusion terms.

Results and Discussion
Computed results for axisymmetric turbulent flow around sphere are compared to Achenbach's experimental data [5] for transcitical flow at Re = 5×10 6 .According to the experimental observation of Achenbach, the flow around sphere can be classified into four regions depending on the Reynolds number.In the subcritical region (Re < 3×10 5 ) the drag coefficient is namely independent of Reynolds number.The critical region (3×10 5 < Re < 4×10 5 ) is characterized by a rapid drop of the drag coefficient.The minimum being reached at the critical Re = 3.7×10 5 , with further increase in Reynolds number, C D slightly increases again which is known as supercritical flow (4×10 5 < Re< 2×10 6 ) and it seems that the curve is going to reach another maximum.The transition from supercritical to transcritical (Re > 2×10 6 ) is rather floating.
Achenbach's experiment is performed at a transcritical Reynolds number where the flow is considered fully turbulent and thus can be directly compared to the turbulent computational models.In the present study, the turbulent flow is simulated using Spalart-Allmaras (S-A) and shear stress transport k-ω turbulence model at Re = 5×10 The predicted pressure coefficient over the surface of the sphere is shown in Fig. 2. The computed results are very close to the Achenbach's experimental values [5].Fig. 3 shows computed values of skin friction coefficient over the sphere.In this case, the computed skin friction coefficient curves does not track well with Achenbach's data forward of the separation point.However, the general trends of the curves are the same.
The accuracy of skin friction coefficient prediction in numerical simulations is highly dependent on the accurate resolution of the turbulent boundary layer near the surface of the body.Accurate calculation of near-wall effects requires an extremely fine mesh in that region.Since boundary layer separation arises due to pressure variations, accurate separation point predictions are dependent on accurate pressure calculations, which require a less fine mesh than skin friction calculations.
However, the discrepancies between actual and computed C f curves are not expected to greatly affect the reliability of the total drag prediction since skin friction drag accounts only for 10 percent of the total drag in this case.Here, the computed C p = 0.  4 shows the angular position of separation points as well as the percentage of difference from experimental values.The numerical predictions of separation point matched Achenbach's experimental data well.Table 5 shows values of drag coefficient predicted by two turbulence models and also experimental values measured by Achenbach.Achenbach measured drag coefficients from 0.09 to 0.18 as Reynolds number varied from 4x10 5 to 5x10 6 .This implies that drag is Reynolds number-dependent in the supercritical to transcritical range, and that the small discrepancy between Achenbach's drag coefficient and the numerically computed drag coefficient may be due to Reynolds number mismatch more than inaccuracies in the numerical method.Achenbach also a.The standard k-ω model and the transformed k-e model are both multiplied by a blending function and both models are added together.b.The blending function is designed to be one in the near-wall region, which activates the standard k-ω model, and zero away from the surface, which activates the transformed k-e model.c.The SST model incorporates a damped cross-diffusion derivative term in the ω equation.d.The definition of the turbulent viscosity is modified to account for the transport of the turbulent shear stress.e.The modeling constants are different.

kG
represents the generation of turbulence kinetic energy due to mean velocity gradients, G ω represents the generation of ω, Г k and Г ω represent the effective diffusivity of k and ω, respectively, Y k and Y ω represent the dissipation of k and ω due to turbulence, D ω represents the cross-diffusion term, S k and S ω are user-defined source terms.

~
D p which for a body of revolution in axisymmetric flow become ; Fig. 1(a).Axisymmetric sphere unstructured grid with boundary conditions.

Fig. 2 .
Plot of pressure coefficient on the surface of sphere at Re = 5×10 6 .
Fig. 4(a) shows the velocity vectors around sphere.The separated region and vortex shedding are clearly visible in the close up view near wall as shown in Fig. 4(b).Table4shows the angular position of separation points as well as the percentage of difference from experimental values.The numerical predictions of separation point matched Achenbach's experimental data well.Table5shows values of drag coefficient predicted by two turbulence models and also experimental values measured by Achenbach.

Fig. 4 .
Fig. 4. (a) Velocity vectors around sphere, (b) close up view that shows the separation point.

Table 1 .
Boundary layer parameters of axisymmetric sphere grids.

Table 2 .
The node spacing of sphere on the edges of the domains for each edge node distribution.

Table 4 .
Angle of separation for axisymmetric turbulent flow around sphere.

Table 5 .
Drag coefficient for axisymmetric turbulent flow around sphere.