PAN Yu-cun （潘雨村）， ZHOU Qi-dou （周其斗）， ZHANG Huai-xin （張懷新）

1. State Key Laboratory of Ocean Engineering， Shanghai Jiao Tong University， Shanghai 200240， China

2. Department of Naval Architecture， Naval University of Engineering， Wuhan 430033， China，E-mail： pyc_navy@163.com

Numerical*simulation of rotating arm test for prediction of submarine rotary derivatives

PAN Yu-cun （潘雨村）1，2， ZHOU Qi-dou （周其斗）2， ZHANG Huai-xin （張懷新）1

1. State Key Laboratory of Ocean Engineering， Shanghai Jiao Tong University， Shanghai 200240， China

2. Department of Naval Architecture， Naval University of Engineering， Wuhan 430033， China，E-mail： pyc_navy@163.com

（Received June 18， 2013， Revised November 6， 2014）

The numerical method is used for predicting the rotary-based hydrodynamic coefficients of a submarine. Unsteady RANS simulations are carried out to numerically simulate the rotating arm test performed on the SUBOFF submarine model. The dynamic mesh method is adopted to simulate the rotary motions. From the hydrodynamic forces and moments acting on the submarine at different angular velocities， the rotary derivatives of the submarine can be derived. The computational results agree well with the experimental data. The interaction between the sail tip vortex and the cross flow in the hull boundary layer is discussed， and it is shown that the interaction leads to the “out-of-plane” loads acting on the submarine.

submarine maneuverability， hydrodynamic coefficients， rotating arm test， dynamic mesh

Introduction

The rotary derivative of a submarine is one of the most important hydrodynamic parameters， and it significantly affects the maneuverability and the dynamic stability of the vehicle. The rotating arm experiment is an effective method to determine the rotating-related hydrodynamic coefficients. The submarine model is fixed to a rotating arm， while the radius of rotation and the angular velocity can be adjusted systematically. The transverse force and the yawing moment acting on the hull at different angular velocities are measured. Consequently， by analyzing these results， the rotary derivatives of the model can be derived. Moreover， it is believed that the rotary derivatives derived from the rotating arm test are generally more accurate than those from the planar motion mechanism （PMM）experiment.

With the advances of the computational fluid dynamics （CFD）， much effort is put in the numerical simulations for manoeuvering purposes［1-6］. The CFD method can be considered as a “numerical rotating arm basin”， which can be used to predict the forces and moments directly from the flow field around a submarine model. Gregory［7］put a deformed body in a rectilinear flow to investigate the flow separation over a body of revolution in a steady turning state. The total force on the curved body can account for the case of rotation within 5% of deviation. The results for the moments see a difference of 20% for /=5R L and a difference of 100% for /=3R L. Sung et al.［8］simulated the flow around a turning submarine named ONR Body 1. The deviations between the computed forces and moments and those of experiments were within 20%. Zhang et al.［9］performed a computational study of the Series 58， SUBOFF and DRDC STR bare hulls undergoing steady turning maneuvering. They found that the rotation increases the lateral force and reduces the yawing moment relative to a hull in a pure translation at equivalent drift angles. Hu and Lin［10］computed the hydrodynamic coefficients of an autonomous underwater vehicle SMAL01 based on an added momentum source method.

In some investigations［8-10］， the submarine rota-tion was considered in two ways： first， at the inlet of the computational domain， the incident flow varies linearly with the constant turning rate from the center of rotation， second， the Navier-Stokes （NS） equations were adjusted to a body-fixed frame of reference. In doing so， the unsteady flow can be treated as a steady problem. The vehicle could remain stationary in the control volume. It is convenient for meshing and numerical computations. But there are some drawbacks：both the NS equations and the turbulence models should be modified in a rotating coordinate system. Further， in the presence of the background rotation，the solid wall and far field boundary conditions should be treated carefully， which can significantly affect both accuracy and convergence. The gradients of the Cartesian velocity are set to zero， and the pressure is obtained by a non-reflecting condition as suggested by Sung［8］. Since in the physical experiment， the hull rotates through undisturbed fluid， it is intuitive to adopt the inertial frame of reference in solving the NS equations for the submarine motion. The artificial pressure gradient［11］is avoided in simulating the flow field in such a way， and results might be closer to the physical reality. The difficulty of this method lies in the moving boundary. To overcome this difficulty， the dynamic mesh method or the overset mesh method was purposed. Carrica et al.［12］simulated the flow field of a surface combatant in the steady turning state and the PMM test using an overset grid method. Pan et al.［13］applied the dynamic mesh method to simulate the PMM experiment performed on the SUBOFF submarine model.

The aim of the present study is to explore the possibility of developing a numerical method to predict the rotary derivatives of a submarine. The virtual rotating arm basin experiments are conducted using the unsteady RANS solver in an inertial reference system，and the moving boundary of the vehicle is taken into account by the dynamic mesh method. The resultant forces and moments are then post–processed to compute the rotating-related coefficients. The flow field around the rotating submarine is also discussed.

1. Numerical model

1.1 Description of the model test

The target studied in this paper is the SUBOFF model. The entity model is a body of revolution， without bow planes， with a sail， two horizontal planes and two vertical rudders， and a ring wing supported by four struts in an “X” configuration. The overall length of the SUBOFF model is 4.356 m， while the length between the perpendicular edges is 4.261 m， and the maximum diameter is 0.508 m.

Generally， the 6 DOF motion of the submarine is described using two coordinate systems. The first is a right-handed， body-fixed coordinate system， with its origin O at a point 2.013 m aft of the forward perpendicular edge on the hull centerline， which is prescribed as the center of the buoyancy （CB）. The -x axis points upstream， the -yaxis points starboard and the -zaxis points downward.

Fig.1 Principal earth-fixed and body-fixed coordinate systems

The second coordinate system， an inertial reference frame， is used to define the motions of the first coordinate system， as shown in Fig.1. In this earthfixed coordinate system， the position of the vehicle’s CB is then expressed in ξ， η， ζ coordinates. The orientation of the body-fixed coordinate system is described by Euler angles ψ（yaw））， θ（pitch）， φ（roll）.

The origin E of the earth-fixed coordinate system is located at the center of the virtual basin， and the model is mounted to a virtual rotating arm which revolves about the axis Eζ. In Fig.1， viewing from above， the model rotates clockwise with a steady angular velocity， r， that implies that the model turns to the starboard. The model?s -yaxis coincides with the arm， while -xaxis and -zaxis are kept normal to the arm. Thus， the transverse velocity component of the model’s CB is always zero and the longitudinal velocity component is equal to its linear speed.

1.2 Governing equations

Numerical simulations are performed with the CFD software Ansys Fluent. The flow around the vehicle is modeled using the incompressible， Reynolds averaged Navier-Stokes （RANS） equations：

where t is the time，iu are the time averaged velo-city components in Cartesian coordinates xi（ i=1，2， 3）， ρ is the fluid density，iF are the body forces，P is the time averaged pressure， μ is the viscous coefficient，iu′ are the fluctuating velocity components.

The finite volume method is employed to discretize the governing equations with the second-order upwind scheme. The semi-implicit method for the pressure-linked equations （SIMPLE） is used for the pressure-velocity coupling. In order to allow the closure of the time averaged Navier-Stokes equations， various turbulence models were introduced to provide an estimation of the Reynolds stress tensor

In Sung?s work［8］， the predictions about a submarine model ONR Body-1 in a steady turning state with different turbulence models were compared. Moreover，the simulations of the SUBOFF model at the incidence was also studied. In most cases， the realizable k-ε model provides better agreements with the experimental results. Therefore， this model is also chosen in the present paper.

2. Computational method

2.1 Computational domain and boundary conditions

The computational domain is a cylindrical volume without an inlet or outlet. The outer boundaries are two body lengths， 2L， away from the hull. These boundaries are either of constant radius or of constant depth so that the flow is nominally tangential to the flow direction. A static pressure opening condition is applied on these boundaries， which allow for the entrainment without a predetermined flow direction， so that the flow can either enter or leave the boundary as dictated by the solution［9］. The no-slip boundary condition is applied to the hull surfaces.

Fig.2 Surface grid of the model

2.2 Mesh definition

The fluid domain is divided into three regions： an inner region， an outer region， and an intermediate layer between them. In the inner region， the multiblock hexahedral grids are used and refined in critical regions such as the boundary layer， the trailing edge of the sail and the control surface， as shown in Fig.2. The coarse hexahedral elements are used in the outer region， which are adequate for the far field. The intermediate layer is covered by unstructured tetrahedral grids， which can be conveniently re-meshed in the case of the boundary deformation， as shown in Fig.3. Such a hybrid mesh strategy allows a fine grid resolution to resolve the viscous flow features near the hull，while saves the computation cost in the outer region where the flow has small variations. The geometry modeling and the grid generation are done by using the Gambit software.

Fig.3（a） Meshes in three sub-regions

Fig.3（b） The near-body meshes in three sub-regions

To numerically simulate the moving boundary of the vehicle， the dynamic mesh method is used. In the transient simulation， the outer domain remains stationary in space， while the inner region containing the SUBOFF model rotates around the -zaxis. It should be noted that the mesh in the inner region remains locked in a position relative to the motion of the vessel. Hence， the mesh of the intermediate layer is deformed to accommodate the motion of the inner region. The new node locations are updated at each time step according to the calculation of the user defined function （UDF）， while the overall mesh topology is maintained.

2.3 Numerical experimental parameters

In the rotating arm experiment， the model rotatesat a constant linear speed， while the radius of rotation and the angular velocity are adjusted accordingly. The non-dimensional angular velocity， r′， is related to the length of the model， L， and the linear velocity U as follows：

the linear velocity is， in turn， related to the angular velocity， r， and the towing radius， R， by

therefore， the non-dimensional angular velocity is typically expressed as

As is known， one of the difficulties in the rotating arm tests is the determination of the forces at r=0. Since the linear velocity is constant， if the angular velocity approaches zero， the radius would be infinity. Due to the restriction on the diameter of a physical basin， the non-dimensional angular velocity r′ generally ranges from 0.2 to 0.5.

Fig.4（a） Lateral forces coefficient Y′ at different angular velocities

Fig.4（b） Yawing moment coefficient N′ at different angular velocities

As for the numerical rotating arm basin， the restriction on the arm length is not a primary obstacle. In this paper， the linear velocity of the model keeps as 4 m/s （with the Reynolds number of 1.693×107based on the vehicle length）， while the radius of rotation varies from 60 m to 20 m， the corresponding non-dimensional angular velocity is about 0.07-0.2.

The hydrodynamic force and moment acting on the turning model are non-dimensionalized as follows：

Figure 4 shows the lateral force and the yawing moment varying with the angular velocity. The values of the rotary derivatives are determined from the slopes of the curves about the force and moment coefficients versus the angular velocity at =0r′， as follows：

3. Results and discussions

3.1 Verification and validation

In the present study， the verification and validation （V&V） method proposed by Stern et al.［14，15］is performed for the rotary derivativesrY′ andrN′. The benchmark data are the results of Roddy’s PMM experiment［16］.

Table1 Results forrY′ and rN′ with different grids

3.1.1 Grid convergence

To inspect the grid convergence， a series of systematically refined grids are generated. In view of the limited computational resource， different refinement ratios are adopted in different regions. Since the nearwall grid has a dominant influence on the precision of the CFD simulation， in the inner region， the grid refinement is conducted with a ratio of2 in three dire-ctions， while in the intermediate layer and the outer region， the refinement ratio is less than2. The fine， medium and coarse grids have 1.341×107， 7.15×106and 4.31×106grid cells， respectively. The+y value，i.e.， the non-dimensional normal distance to the wall from the first grid， is approximately 30 for the fine grid， 43 for the medium grid and 60 for the coarse grid.

Table 1 summarizes the rotary derivativesrY′andrN′ computed with different grids. The solutions corresponding to the three grids are presented as1S， S2and S3. The differences of the solutions ε for medium-fine and coarse-medium solutions are defined as

The convergence ratioGR is defined as

Table 2 shows the grid convergence ratioGR，the order of accuracyGP， the correction factorGC，and the grid uncertaintyGU.

Table2 Grid convergence verification ofrY′ and rN′

Since 0＜RG＜1， the monotonic convergence is obtained forrN′. Thus， the generalized Richardson extrapolation［14］can be used to estimate the grid uncertainty GU. The estimates for the errorGδ*and the order of accuracyGP are defined as

The correction factor is defined as

And，gestp is an estimate for the limiting order of accuracy as the spacing size tends to zero. In the present paper， the expected order of the flow solver is the second. Here， CG=0.28indicates that the leading-order term over-predicts the error. The grid uncertaintyGU can be estimated as

ForrY′， the convergence ratioGR is negative，which implies thatrY′ is oscillatorily converged. Under this condition， the grid uncertainty is estimated simply by bounding the error based on the oscillation maximumsUS and minimumsLS， i.e.，

The values of the grid uncertaintyGU （1.7%1S，4.4%1S） for rY′ and rN′ are reasonable in view of the overall number of grid points used.

3.1.2 Time step convergence

For the time step convergence study， three different time steps are used to obtain solutions with Dt= 0.003318， 0.004694， and 0.006637 （non-dimensionalized with /L V）.

Table3 Results forrY′ andrN′ with different time steps

Table4 Time step convergence verification ofrY′ andrN′

Tables 3 shows the solutions of the three time steps， Table 4 shows the verification parameters for the time step study. The definitions of the time step convergence parameters are similar to the grid convergence parameters. The subscript “G” indicates “Grid”，and the subscript “T” indicates “Time step”.

For bothrY′ andrN′， theTR values satisfy 0＜RT＜1. That is， the solution converges monotonically with the time step refining. Hence， the numerical verification is demonstrated for this case. Reasonably small levels for the time step uncertaintyTU （0.97%，2.04%） are predicted for the coefficientsrY′ andrN′，respectively. The uncertainty levels are smaller than those in the grid study， which suggests that the system is less sensitive to the time step size.

Table5 Validation ofrY′ and rN′

3.1.3 Validation

From the above results， the rotary derivatives derived from the CFD agree well with those from the experiments［16］. It can be preliminarily concluded that the CFD calculation is reliable for estimating the rotary derivatives of a fully appended submarine by simulating the rotating arm test.

3.2 Flow field

The longitudinal distribution of the lateral force Y′ on the model in a typical turning state（r′=0.15） is shown in Fig.5. It is necessary to point out that the lateral velocity varies along the hull. In the present study， the angle of drift at the CB is set to zero. However， the local drift angle （）β λ varies with its longitudinal location over the entire length of the submarine， as shown in Fig.6. Here， λ is the distance between the local point and the forward perpendicular edge. At the forward perpendicular edge， where =0λ， the local drift angle approaches its negative peak， and this angle increases as its location moves toward the stern. Thus， the CB is the critical location where the drift angle changes from negative to positive.

Fig.5 Longitudinal distribution of lateral force Y′ on the SUBOFF model

Fig.6 Drift angle variation in longitudinal locations

As is known， the sail is located between /=Lλ 0.2 and 0.3. In Fig.5， for the region in front of the sail，that is， the location moves ahead from /=0.2Lλ to λ/L =0， the absolute value of the cross flow velocity increases， thus the absolute value of the lateral force increases accordingly. However， when /0.03Lλ＜，the windward projected area sharply decreases， thus the absolute value of Y′ drops.

At the point of /=0.2Lλ， the sail creates a high lift， which brings the lateral force Y to its minimum value. For the points about /=0.89Lλ and λ/L =0.95， they correspond to the control surfaces and the duct. These appendages have relatively large windward projected areas， and they have relatively larger lateral velocities. Thus， it is not surprising that the maximum lateral force Y is created here.

To help understand the mechanism that generates the forces and moments on the body， the flow fieldaround the submarine in =0.15r′ is studied.

Fig.7 The flow field around the SUBOFF submarine

Figure 7 shows the velocity contours around the submarine as it travels in the rotating arm test simulation. The Q criterion， i.e.， based on the positive second invariant of the velocity gradient tensor， is used to define the type and the location of the vortex core［17］. Here， the vortical structures are shown as isosurfaces of =50Q.

During the submarine’s maneuvering， the sail encounters a lateral flow and works as a finite span wing. The lift develops on the sail， as a consequence， a strong vortex is shed from the sail tip， and it is convected downstream along the hull. The sail tip vortex rolls the low velocity wake flow from the sail， and interacts with the boundary layer in the cross flow around the hull.

Fig.8 Flow field at /=0.23Lλ for SUBOFF model at =r′0.15

Another typical feature of the flow field is the horseshoe vortex， which is generated around the sailhull junction， and it is rolled and elongated as it travels along the side face of the sail.

Figure 8 shows an example of the velocity and pressure fields in a cross plane， looking toward the model’s bow. The plane of λ/L =0.232is just behind the leading edge of the sail. Here， a pair of horseshoe vortices are shown at both sides of the sail-hull junction. The local drift angle is negative， that is， the starboard of the body is expected to be on the windward side. Hence， the lateral flow creates a pressure difference between the starboard and the port， what is more， the presence of the sail mightily amplifies this difference. Therefore， the negative lateral force （）Y-is produced.

However， for the aft body， the direction of the cross flow changes its direction. The plane of /=Lλ 0.58 is a cross section behind the CB. The portside of the body becomes the windward side， and the lateral force induced by the pressure difference points to the starboard. （i.e.， Y）. The horseshoe vortices are also affected by the cross flow velocity and move towards the starboard.

Fig.9 Flow field at /=0.58Lλ for SUBOFF model at =r′0.15

Another point of interest is that the hull in the horizontal motion is acted by a vertical force， called the“out-of-plane” loads. In Fig.9（a）， the direction of thevortical flow shedding from the sail tip is consistent with that of the cross flow in the hull boundary layer. And the merging of the two flows results in a higher velocity on the upper hull in comparison to that on the lower hull. According to the Bernoulli?s theorem， the relatively lower pressure on the upper hull surface is produced， as shown in Fig.9（b）， which consequently results in an upward normal force （）Z- and a bow pitch down moment.

4. Conclusions

A numerical procedure for the prediction of the rotary derivative of a submarine is proposed， based on solving the unsteady incompressible Reynolds-averaged Navier-Stokes equations in an inertial frame of reference system. The force and the moment on the SUBOFF model during the rotating arm test are successfully calculated. The prediction of the rotary coefficients of the submarine model enjoys an acceptable level of accuracy.

The longitudinal distribution of the lateral force on the model in a turning state is investigated. The local drift angle varies with its longitudinal location，which significantly affects the lateral force. From the CFD solution， the plots of the velocity and pressure distributions at several longitudinal locations are created and analyzed.

Another noticeable phenomenon is the interaction between the vortex shed from the sail and the hull boundary layer. For a submarine in the horizontal rotating arm test， this interaction results in a upward force on the aft hull and consequently a bow pitch down moment. This flow physics is illustrated and explained with the help of the computed flow field.

［1］ PAN Zi-ying， WU Bao-shan and SHEN Hong-cui. Research of CFD application in engineering estimation of submarine maneuverability hydrodynamic forces［J］. Journal of Ship Mechanics, 2004， 8（5）： 42-51（in Chinese）.

［2］ SIMONSEN C. D.， STERN F. RANS maneuvering simulation of Esso Osaka with rudder and a body-force propeller［J］. Journal of Ship Research, 2005， 49（2）：98-120.

［3］ LI Ying-hua， WU Bao-shan and ZHANG Hua. Research on application of prediction for unsteady maneuvering motion of underwater vehicle by dynamic mesh technique in CFD［J］. Journal of Ship Mechanics, 2010，14（10）： 1100-1108（in Chinese）.

［4］ PANG Yong-jie， YANG Lu-cun and LI Hong-wei et al. Approaches for predicting hydrodynamic characteristics of submarine objects［J］. Journal of Harbin Engineering University, 2009， 30（8）： 903-908（in Chinese）.

［5］ PHILLIPS A. B.， TURNOCK S. R. and FURLONG M. Evaluation of manoeuvring coefficients of a self-propelled ship using a blade element momentum propeller model coupled to a Reynolds averaged Navier Stokes flow solver［J］. Ocean Engineering, 2009， 36： 1217-1225.

［6］ RACINE B. J.， PATERSON E. G. CFD-based method for simulation of marine-vehicle maneuvering［C］. 35th AIAA Fluid Dynamics Conference and Exhibit. Toronto， Ontario， Canada， 2005.

［7］ GREGORY P. Flow over a body of revolution in a steady turn［D］. Doctoral Thesis， Melbourne， Australian：University of Melbourne， 2006.

［8］ SUNG C. H.， JIANG M. Y. and RHEE B. et al. Validation of the flow around a turning submarine［C］. The 24th Symposium on Naval Hydrodynamics. Fukuoka，Japan， July， 2002.

［9］ ZHANG J. T.， MAXWELL J. A. and GERBER A. G. et al. Simulation of the flow over axisymmetric submarine hulls in steady turning［J］， Ocean Engineering, 2013，57： 180-196.

［10］ HU Z. Q.， LIN Y. Computing the hydrodynamic coefficients of underwater vehicles based on added momentum sources［C］. The 18th International Offshore and Polar Engineering Conference. Vancouver，Canada， 2008， 451-456.

［11］ RAY A.， SINGH S. and SESHADRI V. CFD as a concept design tool for manoeuvring hydrodynamics of Underwater vehicles［C］. Warship 2010: Advanced Technologies in Naval Design and Construction. London，UK， 2010.

［12］ CARRICA P. M.， WILSON R. V. and NOACK R. W. et al. A dynamic overset， single-phase level set approach for viscous ship flows and large amplitude motions and maneuvering［C］. The 26th Symposium on Naval Hydrodynamics. Rome， Italy， 2006.

［13］ PAN Yu-cun， ZHANG Huai-xin and ZHOU Qi-dou. Numerical prediction of submarine hydrodynamic coefficients using CFD simulation［J］. Journal of Hydrodynamics, 2013， 25（1）： 840-847.

［14］ STERN F.， WILSON R. V. and COLEMAN H. W. et al. Comprehensive approach to verification and validation of CFD simulations-Part 1： Methodology and procedures［J］. Journal of Fluids Engineering, 2001， 123（4）：793-802.

［15］ WILSON R. V.， STERN F. and COLEMAN H. W. et al. Comprehensive approach to verification and validation of CFD simulations-Part 2： Application for RANS simulation of a cargo/container ship［J］. Journal of Fluids Engineering, 2001， 123（4）： 803-810.

［16］ RODDY R. F. Investigation of the stability and control characteristics of several configurations of the DARPA SUBOFF model （DTRC model 5470） from captivemodel experiments［R］. DTRC/SHD 1298-08， 1990， 1-108.

［17］ LIU Zhi-hua.， XIONG Ying and WANG Zhan-zhi et al. Numerical simulation and experimental study of the new method of horseshoe vortex control［J］. Journal of Hydrodynamics， 2010， 22（4）： 572-581.

10.1016/S1001-6058（15）60457-7

* Project supported by the National Natural Science Foundation of China （Grant No. 11272213）.

Biography: PAN Yu-cun （1980- ）， Male， Ph. D. Candidate，Lecturer

Corresponding ahthor: ZHANG Huai-xin，E-mail： hxzhang@sjtu.edu.cn