Transient Numerical Simulations for Hydrodynamic Derivatives Predictions of an Axisymmetric Submersible Vehicle

In this study, a method for obtaining of hydrodynamic derivatives by numerically replicating the Planar Motion Mechanism (PMM) tests of an axisymmetric submersible model is demonstrated. The numerical simulations of PMM tests are regarded as transient due to the movement of the model in the discretized computational domain thus causing mesh deformation. To accommodate the sway and yaw oscillation motions of the model, the entire computational domain is divided into three zones namely rotating, inner and outer zone. Multi-block structured grid is generated with finer resolution in the proximity of the model to capture the boundary-layer flows. Non-conformal fluid interfaces are used to connect the three zones. Commercial CFD Solver FLUENT is used to simulate the flow characteristics while the dynamic mesh capability included in the software is applied to handle the mesh deformation during the movement of the model. In order to verify the CFD method, 6:1 prolate spheroid is used as it can be idealized as an axisymmetric submersible model. The CFD results of added mass derivatives of the model show very close agreement when compared with the theoretical values. The present study is an attempt towards developing an economical CFD method for evaluating the hydrodynamic derivatives of submersible platforms such as submarines, torpedoes and autonomous underwater vehicles during early design stages.


INTRODUCTION
As the technology is improving, the missions being expected for, both in defense and scientific research domain, submersibles platforms such as torpedoes, submarines, and autonomous undersea vehicles are becoming highly complicated and challenging (National Research Council, 2004).In order to meet the expectations of these missions, the performance of the newly designed vehicles in terms of speed, endurance, maneuverability, controls, etc., need to be superior to the existing ones.Understanding the hydrodynamic characteristics of a submersible is the prerequisite for the analysis of new designs which in turn involves a detailed investigation of differential equations, which govern its motion (Feldman, 1979).In these equations, hydrodynamic forces and moments appear as coefficients or derivatives, which are classified as static, rotary and acceleration derivatives.The static derivatives are due to linear velocity of the vehicle, rotary are due to angular velocity and acceleration derivatives are due to linear and angular accelerations.An accurate and reliable assessment of these derivatives is essential for cost effective yet advanced design configurations of the submersible vehicles.
The need to identify the hydrodynamic derivatives has led to the evolution of numerous captive model tests such as Straight-Line Tests (SLTs); Rotating Arm Tests (RATs) and using Planar Motion Mechanism (PMM) etc, Comstock (1967).However, performing the captive model tests has proven to be expensive and time-consuming especially during the preliminary stages of the design phase.With the rapid development of computational technology and increased maturity in numerical methods, it is now feasible to consider CFD (Computational Fluid Dynamics) as attractive and low cost tool suitable for deriving the hydrodynamic derivatives of the submersibles by replicating the experimental procedures (Phillips et al., 2010).The present study demonstrates the CFD predictions of hydrodynamic derivatives by replicating the forced oscillation motions of an axisymmetric submersible model as in PMM tests.As opposed to SLTs and RATs that directly supplies static and rotary coefficients only, PMM tests can be used to determine all the hydrodynamic coefficients in the equations of motion (Newman, 1977).The numerical approach consists of solving the flowfield using Unstaedy Reynolds Averaged Navier Stokes Equations, oscillating the model, controlling the mesh quality during the The organization of the stutdy is as follows; after this introduction, next section provides brief description of calculating hydrodynamic derivatives through PMM tests.After that the discussion moves onto the details of numerical approach covering geometric details of model, discretization of fluid domain, dynamic mesh algorithm, CFD solver settings, etc.The numerically computed results are then analyzed and compared.In the last section, conclusions regarding the CFD approach are drawn.

PMM TESTS PROCEDURE
The technique of conducting PMM tests was first described by Gertler (1959) who referred particularly to its use with submarine model.The essence of PMM tests is that the model is forced oscillated with predetermined amplitude and frequency whilst being towed with constant velocity down the tank.
Oscillatory motion in sway: Pure sway PMM tests consists of imposing an oscillatory motion in y direction with a constant forward speed such that model always remains parallel to the tank centerline as shown in the Fig. 1 In this way the angular velocity and acceleration are always zero in body coordinates during the course of the motion.The displacement ,velocity and acceleration can be expressed as: (3) where, and represents amplitude and frequency of the oscillation motion.While moved harmonically up and down in flow direction, an angle of attack is induced on the body.During one cycle of the motion, the velocity and acceleration vary between a positive and negative extreme value, and hence amplitude of force and moment are exerted on the model.These forces appears in sway and yaw equations of motion as: (4) (5) Since 0, therefore, motion equations will be simplified as: The force and moment can be decomposed into components in phase and out of phase with the displacement .From pure sway equation, the acceleration being sine function is in phase whereas the velocity being cosine function 90° out of phase with the displacement .The in phase component of the force and moment are directly related to acceleration and therefore can be used to compute acceleration deriveatives and .Similarly, the force and moment that are out of phase with displacement ( , will yield velocity derivatives and .Now when force and moment are in phase ( 0 and ),the associated acceleration derivatives will be computed as: For out of phase condition ( and 0), force and moment correpond to velocity derivative as: Oscillatory motion in yaw: For pure yawing motion, it is necessary that 0, so that the velocity of the model must always be tangential to the path of the model.This test enables us to calculate rotary hydrodynamic derivatives , , and .Figure 2 illuatrates the pure yaw motion,in which the resulatnt velocity of the body is directed along the tangent of the model.In this way the vertical component of the velocity and hence acceleration remain zero in body coordinates during the entire cycle.The motion parameter in case of pure yawing can be expressed as: where, is the amplitude of the angular oscillations.For pure yaw motion, the sway and yaw equations of motion will be simplified as: Decomposing and into components in phase ( 0) and out of phase ( 0) with the angular displacement , following hydrodynamic derivatives will obtained: Geometric model: Prolate spheroids are widely used to approximate the shape of an axisymmetric underwater vehicle.Although the geometrical shape of prolate  (2002).For the aforementioned reasons, the 6:1 prolate spheroid has been selected for the computation of the hydrodynamic derivatives in the present study.The coordinate system adopted in this study is such that the flow direction is in the positive x-axis, y points to the upward vertical direction, and the x y plane makes the vertical symmetry plane.The origin of the coordinate system is placed at the center of the spheroid as shown in the Fig. 3.The numerical simulations are conducted for 6:1 prolate spheroid model with Reynolds number approximately 4.5x10 ,which is based on the free stream velocity U ο and the model-length L.
Added mass derivatives: Prolate spheroid is surface of revolution generated by rotating the ellipse about its major axis as shown in the Fig. 3.The equation of ellipsoid is given as: For a prolate spheroid ( / 1, ), the eccentricity of the sections through the axis of symmetry is defined as: Other parameters are: The mass of the volume of fluid, with density , displaced by the spheriod and its moment of inertia ) are expressed as: Lamb's factors (Lamb, 1945) are represented as: With these definitions, the expressions of added mass and inertia derivatives for prolate spheroid will be Fossen (1994): COMPUTATIONAL DETAILS Domain decomposition: The numerical simulations of PMM tests are regarded as transient due to the movement of the body in the discretized computational domain thus causing mesh deformation.For this reason, the dynamic mesh capability of FLUENT software is utilized in this study.The dynamic mesh algorithm can maintain the domain grid qualities during the course of the movement of body using the three different grid updating schemes (Layering, spring smoothing and local re-meshing).For the rigid body moves in the domain; it is possible to make the grid near the body move with the body while using the spring smoothing and local remeshing techniques to maintain the grid qualities.However, fine unstructured grids and smalltime steps are often required for accuracy and stability, which increases the computational time (Zhang et al., 2010).An alternative method is to replicate the motion of the body with the fluid domain split into an inner and outer zone.The outer zone remains fixed in space while the inner zone containing the hull moves to create the motion induced by a PMM.In this way, the mesh in the inner zone remains locked in position relative to the motion of the vessel.This method is proven to very promising as high quality structured grids and larger time steps can be used.The only drawback of this Fig.4: Splitting up computational domain into different zones method is that it can be used for the Heave/Sway PMM tests where the motion of the body is purely linear (Xianzhoa and Yumin, 2010).In order to facilitate the simulations for pure yaw PMM test while maintaining the benefits of above method, the inner zone is further divided into two sub-zones as depicted in the Fig. 4. The rotating zone around the prolate spheroid has a semi spherical shape created to accommodate rotation motion and it is connected to inner subzone using nonconformal grid connection method.In this way, the rotation motion of the spheroid model can be handled dynamically without mesh distortion.
The oscillation motions on the spheroid model are imposed by compiling a User Defined Function (UDF) written in C++ programming language.The dynamic mesh UDF is implemented using the DEFINE-CG-MOTION macro.For pure sway motion, the rotating and inner zones move together with model in linear direction specified by UDF and mesh will remain unaffected.The mesh of the outer zone will be deformed in order to accommodate the motion of the spheroid model.In case of pure yaw motion, the spheroid model and rotating zone move together in rotation, at the same time, they also have translation motion with inner subzone.As in the pure sway motion, the mesh of the outer zone will be deformed.In this way, the pure yaw motion will be achieved without any mesh distortion in the inner zone.

Dynamic mesh method:
In all above-mentioned cases, the value of amplitude of the oscillation is taken very small ( = 0.04 m), because it represents the perturbation.For such small-scale body motions, the solver recommends localized smoothing method to update the mesh of the outer domain.In the smoothing method, the edges between any two mesh nodes are idealized as a network of interconnected springs.A displacement at a given boundary node will generate a force proportional to the displacement along all the springs connected to the node.Using Hooke's Law, the force on a mesh node can be written as: where, ∆ and ∆ are the displacements of node and its neighbor , is the number of neighboring nodes connected to node , and is the spring constant between node and its neighbor .
(27) At equilibrium, the net force on a node due to all the springs connected to the node must be zero.This condition results in an iterative equation such that: Since displacements are known at the boundaries (after boundary node positions have been updated), Eq. ( 28) is solved using a Jacobi sweep on all interior nodes.At convergence, the positions are updated such that: where, 1 and are used to denote the position at next and current time-steps respectively.

Grid generation and boundary conditions:
High quality multi-block grid is generated using FLUENT's pre-processor GAMBIT.Structured curvilinear grids (O-O topology) are employed with finer resolution in the proximity of the model.The first cell height for desired wall Δy is estimated from the following empirical correlations based on flat-plate boundary layer theory (Schlichting, 1966): where, Δ = First cell height (distance from the wall surface) / , = The wall shear stress Kinematic viscosity of sea water Now: From the above equation, first cell height is located such that wall Δ is kept in a range 30~100.Figure 5 shows that the mesh in the rotating zone which surrounds the spheroid model is much denser as compared to the other two regions.Figure 6 presents the enlarged view of the mesh in the vicinity of the model displaying the resolution of boundary layer.In total, the final mesh consists of 377,320 hexahedral cells with 395776 nodes.In Fig. 5 also shown are the boundary conditions imposed: at the inlet face and outer boundaries velocity inlet boundary condition is applied; static pressure is set at the outlet; no-slip condition is imposed at the surface of spheroid model and on the central plane symmetry boundary condition is applied.Outer boundaries of the computational domain are positioned far enough to have no effect on motion of the spheroid model.

Solver settings:
Series of unsteady RANS CFD simulations are carried out using the FLUENT's double precision segregated pressure based solver.For convection and diffusion terms, second order upwind scheme is adopted whereas SIMPLEC algorithm is chosen for velocity-pressure coupling.For convergence criteria, residuals of continuity, velocities and turbulence model is selected due to its superior performance track records for complex boundary layer flows under adverse pressure gradient and separation (Kim et al., 2003).The simulations of PMM tests are performed with three values of 1.256 rad/s, 1.885 rad/s and 2.512 rad/s.The number of time steps for each cycle is

RESULTS AND DISCUSSION
Planar sway motion simulations: Figure 7 illustrates typical velocity contours around the model during the planar sway motion.Due to the oscillations, wake generated behind the model represents oscillatory pattern .The time history of the sway force Y and moment N for three frequencies are plotted in Fig. 8 and 9, respectively.The steady state curves for three cycles are obtained after some instability at initial time steps.The computed forces and moments are decomposed into in phase and out of phase components with respect to displacement y.For numerical simulations, there will be no contribution of inertial forces and moment since only hydrodynamic forces are computed by integrating pressures on the body (Saout and Ananthakrishnan, 2011).Therefore,the equations 8 and 9 are reduced as: Using the above equations ,the hydrodynamic force and moment derivatives are calculated and converted into non-dimensional form (S. N.A.M.E, 1964).Table 1 shows the comparison between the added mass derivatives from the CFD simulations and theoretical   Planar yaw motion simulations: Figure 10 illustrates the velocity contours with the wake pattern behind the model during the planar yaw motion.The plots of the sway force Y and moment N for three values are shown in Fig. 11 and 12, respectively.Again instabilities in forces and moments are observed at the start of the simulation.However, after they reach a stable state the values are measured for three complete cycles.The forces and moments are then decomposed into in phase and out of phase components with respect to angular displacement .Filtering out contribution of inertial forces and moments, equations 15 and 16 are reduced as: , , Using the above equations ,the computed value of hydrodynamic force and moment derivatives in nondimensional form are given in Table 2. Comparing the inertia derivatives from the simulations and theoretical calculations, the maximum difference is about 6.5% which can be regarded as within the range of acceptable limits.

CONCLUSION
The current study investigates the use of CFD in the prediction of hydrodynamic derivatives for an axisymmetric submersible model.For this purpose, the captive model testing using planar motion mechanism is briefly described and numerically simulated using the capabilities of commercial CFD solver FLUENT.Forced sway and yaw oscillation motions are produced by implementing user defined functions.During the oscillation motions of the model, FLUENT solves the flow field by applying unsteady RANS equations while the deformation of the mesh during the model motion is controlled by dynamic spring smoothing technique.In order to maintain the quality of the grid during the motion the entire computational domain is divided into rotating, inner and outer zones.The three zones are meshed with high quality structured grids and connected via non conformal fluid interfaces.To get verification of the CFD method, prolate spheroid model with 6 to 1 aspect ratio is used due to its quality to have same shape and flow characteristics as an axisymmetric submersible vehicle.Comparison of the computed results of the added mass derivatives obtained from pure sway and yaw oscillation motion reveals very good agreement with the theoretical value.Based on this study, it is concluded that the CFD method is well capable of economically evaluating the hydrodynamic derivatives of submersible platforms such as submarines, torpedoes and autonomous underwater vehicles.

Fig. 1 :
Fig. 1: Oscillation motion in Sway movement of the model and recording the computed forces and moments.The organization of the stutdy is as follows; after this introduction, next section provides brief description of calculating hydrodynamic derivatives through PMM tests.After that the discussion moves onto the details of numerical approach covering geometric details of model, discretization of fluid domain, dynamic mesh algorithm, CFD solver settings, etc.The numerically computed results are then analyzed and compared.In the last section, conclusions regarding the CFD approach are drawn.

Fig. 3 :
Fig. 3: Prolate spheroid with major and minor axes spheroid is simple but in maneuvering its flow field carries over qualitatively to torpedoes, submarines, AUVs like stagnation flow, cross-flow separation, Boundary layer detachment, etc. Constantinescu et al.(2002).For the aforementioned reasons, the 6:1 prolate spheroid has been selected for the computation of the hydrodynamic derivatives in the present study.The coordinate system adopted in this study is such that the flow direction is in the positive x-axis, y points to the upward vertical direction, and the x y plane makes the vertical symmetry plane.The origin of the coordinate system is placed at the center of the spheroid as shown in the Fig.3.The numerical simulations are conducted for 6:1 prolate spheroid model with Reynolds number approximately 4.5x10 ,which is based on the free stream velocity U ο and the model-length L.

Fig. 8 :
Fig. 8: Time history of sway force for three cycles

Fig. 10 :
Fig. 10: Typical velocity contours around the model during yaw motion