Using reduced hydrodynamic models to accelerate the predictor–corrector convergence of implicit 6-DOF URANS submarine manoeuvring simulations Academic research paper on "Earth and related environmental sciences"

Accepted Manuscript

Using Reduced Hydrodynamic Models to Accelerate the Predictor-Corrector Convergence of Implicit 6-DOF URANS Submarine Manoeuvring Simulations

Mark C. Bettle, Andrew G. Gerber, George D. Watt

PII: DOI:

Reference: To appear in:

S0045-7930(14)00086-3

http://dx.doLorg/10.1016/j.compfluid.2014.02.023 CAF 2474

Computers & Fluids

An. International Joiwnal

Received Date: Revised Date: Accepted Date:

30 July 2013 6 February 2014 25 February 2014

Please cite this article as: Bettle, M.C., Gerber, A.G., Watt, G.D., Using Reduced Hydrodynamic Models to Accelerate the Predictor-Corrector Convergence of Implicit 6-DOF URANS Submarine Manoeuvring Simulations, Computers & Fluids (2014), doi: http://dx.doi.org/10.10167i.compfluid.2014.02.023

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Using Reduced Hydrodynamic Models to Accelerate the Predictor-Corrector Convergence of Implicit 6-DOF URANS Submarine Manoeuvring Simulations

Mark C. Bettlea*, Andrew G. Gerberb, George D. Watta

aDefence Research and Development Canada - Atlantic, P.O. Box 1012, Dartmouth, NS,

Canada B2Y 3Z7

bDepartment of Mechanical Engineering, University of New Brunswick, P.O. Box 4400, Fredericton, NB, Canada E3B 5A3

Abstract

An implicit predictor-corrector method is presented for the simultaneous integration of the six degrees-of-freedom (DOF) equations of motion for a manoeuvring submarine and the unsteady Reynolds-Averaged Navier Stokes (URANS) equations describing the vehicle hydrodynamics. The novel method uses coefficient-based hydrodynamic models for estimating the Jacobian matrix for Newton iteration. The method is applied to emergency rising and horizontal plane zig-zag manoeuvres. It is shown to converge faster at each timestep than under-relaxed fixed-point iteration with an optimum relaxation parameter. A simple model containing only primary linear hydrodynamic coefficients that are relatively easy to estimate or measure was found to be adequate for modelling the Jacobian matrix in these simulations.

Keywords: six DOF, manoeuvring, simulation, submarine, predictor-corrector method, URANS, coefficient-based hydrodynamic models

* Corresponding author Email address: Mark.Bettle@drdc-rddc.gc.ca (Mark C. Bettle )

Preprint submitted to Elsevier

3rd March 2014

Nomenclature

r,N<Z_

B = pVg

CB, CG d

Ed ,Ei EOM f

Tempoi Equatic

F¿ ,1 y , F

Fxy, Fzx: Fyz

Jbs Jf

Jo, JL, JLp, JAMo

K, M,N

12 x 12 coefficient matrix for the 12 vehicle EOM (Eq. (4)).

Inertial coupling vector in 6 DOF EOM (Eq. (2)). Buoyancy.

Backwards differentiation formula.

Second-order two-step backwards differentiation formula. Centres of buoyancy and gravity. Maximum hull diameter. Degrees of freedom.

discretization error; iterative error. ions of motion. External load vector in 6 DOF EOM (Eq. (2)). Non-linear function being solved for y, y.

Coefficient-based model for steady translational hydrodynamics. Hydrodynamic load vector, [XH, YH, ZH, KH, MH, NH]. Gravitational constant, generic function in Eq. (1). Right hand side of the 12 vehicle EOM (Eq. (4)). Moments of inertia in body axes. Products of inertia in body axes.

Propeller self-propulsion behind-the-boat advance ratio. Jacobian matrix for Newton iteration. Coefficient models for the Jacobian matrix (Table 1). Body axis moments.

Nc ,Nf

p, q,r R = UL/v s = [u, v, w,p, q, r]T t, At

u, v, w uf ,vf ,wf

velocity velocity

luid sub-it<

U = %/ u2 + v2 + w2

W = mg Wn

Wt = Wo - B

x,y,z xo, yo, zo xb,yB ,ZB

xf ,yf ,zf

XG,yG, ZG

X,Y,Z XP

Greek letters

Effect of propeller torque on rolling moment K. Overall length of the submarine hull. Mass within V.

6 x 6 vehicle mass and added mass matrices.

Number of correctors per timestep; number of fluid sub-iterations per corrector.

Body axis angular velocities. Reynolds number.

Vector of vehicle linear and angular velocity components. Time, timestep size. Body axis velocities.

Components of fluid velocities in fluid coordinate directions.

Overall speed of vehicle.

Unsteady Reynolds-Averaged Navier Stokes.

Volume of the external hydrodynamic envelope, including ballast tanks. Weight within V. Initial weight within V. eight trim. dy-fixed axes. Inertial (earth-fixed) axes.

Coordinates of CB (centroid of V ) in body axes. Fluid axes.

Coordinates of CG (centre of mass m) in body axes. Body axis forces.

Effect of propeller thrust on axial force X. Submarine state vector, [u, v, w,p, q, r, xo, yo, zo, 9,

¡3 = tan 1(—v/u) Angle of drift.

Sr ,0s

Rudder (for yaw control) and sternplane (for pitch control) deflections; direction is found from the right hand rule using body axes. Relaxation parameter for fixed-point iteration. Flow incidence, always positive. Dynamic viscosity of sea water. Kinematic viscosity of sea water. Density of sea water.

Yaw, pitch, and roll Euler angles giving body axes relative to inertial axes.

Yaw angle at which the rudder is reversed in zig-za

Subscripts AM

f o R T

dot ove

/er a sy

Added mass. A command. Control surface. Propeller. Fluid.

Initial condition.

ational component. anslational component.

manoeuvres.

Rota Tran

A dot over a symbol indicates differentiation with respect to time. >

Introduction

Realistic submarine manoeuvring predictions are needed to understand operational limitations and set safe operating limits. Manoeuvring simulations have traditionally used coefficient-based models for the hydrodynamic loads, such as the model of Gertler and Hagen [1], which was revised by Feldman [2]. These models have been successful at predicting the motions of conventional submarines during moderate, deeply-submerged manoeuvres. However, they have been found inadequate for extreme manoeuvres or non-conventional hull

forms [3]. Also, these standard models do not incorporate bottom clearance or free surface effects, or disturbance loads from passing ships. An understanding of these effects is becoming increasingly important due to a shift in focus towards more littoral operations [4], in which submarines operate in possibly congested seaways.

Computational fluid dynamics (CFD) based on the Reynolds-Averaged Navier Stokes (RANS) equations have been found to give relatively good predictions of hydrodynamic loads on manoeuvring submarine hull forms (e.g. [5, 6, 7, 8, 9, 10]) and submarine hulls with appendages (e.g. [4, 11, 12]). Early studies focused on the computation of hydrodynamic loads for various prescribed submarine motions in a manner analogous to captive-model experiments. These simulations can be used in place of (or in conjunction with) experiments to build up coefficient-based hydrodynamic models to be used in 6 DOF manoeuvring simulations. With the advent of faster parallel computers, it has become feasible to use a different approach in which the 6 DOF motion of the vehicle is evaluated as part of an unsteady RANS (URANS) simulation. With this method, the flow computation provides the hydrodynamic loads needed for the 6 DOF vehicle equations of motion and the solution for vehicle motion is used to update boundary conditions in the flow solution. This motion-coupled approach provides a more comprehensive accounting of the unsteady hydrodynamic loads because information is not lost in the process of making CFD results fit an empirical model. It can also be used to predict unsteady multi-body manoeuvring scenarios for which empirical models have not been developed and experiments are more expensive.

Work on motion-coupled, 6 DOF URANS simulations of submarine manoeuvres was pioneered by the Engineering Research Center at the Mississippi State University and the Applied Research Laboratory at the Pennsylvania State University in the development of the code UNCLE [3]. Pankajakshan et al. [13] used UNCLE to simulate control-surface induced overshoot manoeuvres of the ONR Body1 Radio Controlled model. Their simulations incorporate deflecting control surfaces and a rotating propeller within the RANS solution, which is solved on a

structured multi-block dynamic grid. The simulated submarine motion agreed very well with a free-swimming model experiment. A similar study conducted by Venkatesan and Clark [14] using the flow solver COMET also showed reasonable agreement with experiment. Six DOF unsteady RANS simulations by Bettle et al. [15] have also predicted the development of sudden rolling motions in a diesel electric submarine rising buoyantly to the surface in an emergency.

The code OVER-REL-TCURS was developed from a variant of UNCLE to handle the 6 DOF motions of multiple bodies using overset grids. Dreyer and Boger [16] used OVER-REL-TCURS to perform model-scale simulations of a tanker overtaking a submarine with different initial clearances and relative velocities. Simulation results for submarine pitch and displacement in the case where control surfaces were fixed agreed well with free-running model scale experiments.

Other researchers have applied similar methods to simulating the motion of other marine vehicles. One of the most impressive examples is the study by Carrica et al. [17] of a surface combatant, ONR Tumblelhome, broaching in waves. Simulations were conducted using the CFD code CFDShip-Iowa v4.0. Moving rudders and rotating propellers are explicitly incorporated in the CFD solution using dynamic overset grids and the free surface is modelled using a single-phase level set approach. A mesh with 21 million grid points and detached eddy turbulence modelling were used for the simulations, which required approximately 4 weeks of wall clock time on 195 processors running in parallel. Results for the ship's motion compared well with experiment.

Results for CFD-based 6 DOF manoeuvring simulations for marine applications are promising. However, the computational expense is several orders-of-magnitude greater than coefficient-based simulations and it is worthwhile to investigate methods for improving the efficiency of the simulations. The present work is focused on the efficiency of the simultaneous integration of the 6 DOF rigid body equations and the URANS equations in time. As described further in Section 2 below, the 6 DOF rigid body equations are a set of 6 non-linear, coupled, first order ordinary differential equations (ODE) relating the rate of

the c Rung

change in vehicle momentum to the external loads on it, having the form:

Ms = f (t, s, s) (1)

where t is time, s is a vector containing the 6 vehicle velocities (3 translations (u,v,w), 3 rotations (p,q,r)), M(t, s) is the mass matrix, and f (t, s, S) is the summation of the external load vector and an inertial coupling vector (defined in Eq. (2) below). The hydrodynamic loads in f are a function of s because of added mass effects, making Eq. (1) implicit. If (1) is to be solved in this form, an implicit numerical integration method is, in our experience, most reliable. This is the case when f is provided by a URANS solver. However, with coefficient based methods, such as that described by Watt [18], added mass effects in f are readily modelled as a linear function of ¡s:

f (t, s, s)= g(t, s)+MAs (1a)

where M^ is a 6 x 6 matrix of constant, geometry-dependent added mass coefficients. The M^ term can be moved to the LHS of (1) thereby allowing an explicit integration scheme to be used.

Implicit multistep methods based on the backward differentiation formula (BDF) and implicit Runge-Kutta methods have been identified as methods of choice for (1) [19, 20]. However, it is still possible to integrate (1) using explicit methods. Explicit Runge-Kutta methods were used for the 6 DOF equations in the code UNCLE [21] and its successor OVER-REL-TCURS [16]. On the other hand, implicit multi-step methods were used in CFDShip-Iowa [22, 23]. For urrent work, we briefly tried the explicit Heun's method (a second order unge-Kutta method) but found it unstable. Henceforth, our effort has been on a implicit method based on the second-order two-step backward differentiation formula, BDF2, which we have found to be robust for a wide variety of submarine manoeuvres. In addition to stability considerations, the BDF2 was selected to be consistent with the time integration scheme of the flow solver in order to minimize numerical errors in the fluid solution near the far-field boundaries.

In this paper, the efficiency of the BDF2 predictor-corrector method is in-

vestigated for two iterative schemes: 1) fixed-point iteration with an under-relaxation parameter and 2) Newton iteration using a coefficient-based hydro-dynamic model to estimate Jacobian matrix terms. Two unsteady submarine manoeuvres are simulated in order to test these methods: the emergency rising manoeuvre [24, 25, 18, 15] and a horizontal zig-zag manoeuvre.

2. Modelling

ic standa

2.1. Submarine Geometry

The hull-sail-tail (HST) configuration of DRDC's generic standard submarine geometry, shown in Figure 1, is used for all CFD simulations in this work. This submarine has an axisymmetric hull with a length-to-maximum diameter ratio of 8.75 and consists of three profile sections: a Riegels type D2 nose, a constant diameter midsection, and a parabolic tail. A full-scale length of L = 70 meters is considered for this work. The sail has a NACA0020 profile with a rectangular profile and a flat tip. The tail is a cruciform ('+') configuration of 4 identical NACA015 tailplanes with flat tips. The geometry is fully defined in Reference [26].

2.2. Experimental Data and the Coefficient-Based Model

Wind tunnel experiments were conducted with an almost identical geometry referred to as the DRDC-STR model (see note in Figure 1), with the Static test rig (STR) model test facility in the Institute for Aerospace Research (IAR) 9 meter wind tunnel in Ottawa [27, 26] using a 6 meter long model. STR data includes yaw sweeps of ±30 degrees with the HST model rolled in 30 degree increments at a Reynolds number based on length, R, up to 23 million; this data is presented in [25]. Propulsion experiments were also performed at incidence angles from —30 to 30 degrees at Reynolds numbers over 20 million [28, 29]. These data are the basis of a high-incidence, coefficient based model presented in [18]. Unsteady hydrodynamics are modelled using added mass coefficients estimated using replacement ellipsoids for each submarine component [30] and

Figure 1: DRDC standard submarine geometry (hull-sail-tail configuration). Reproduced from [26]. * 2.8836d for the DRDC-STR model.

rotational and control derivatives are estimated using semi-empirical methods. These models are used extensively in the present work and will be referred to throughout as simply 'the coefficient-based model'. Six DOF simulations that use the coefficient-based model for evaluating hydrodynamic loads in place of URANS will be referred to as 'coefficient-based simulations'.

2.3. Rigid Body Equations of Motion

CIn all simulations considered in this work, the submarine is modelled as a deeply-submerged rigid body far removed from other bodies and the free surface. The motion of the submarine is evaluated using the 6 DOF rigid-body equations of motion, expressed in a local body-axis coordinate system (x,y, z) that moves with the submarine, as shown in the Nomenclature. These equations have been described elsewhere [1]; however, they will be repeated to assist in describing the details of a new integration scheme used in this work that

incorporates coefficient-based modelling for improved convergence. The 6 DOF rigid body equations relate the rate of change in momentum of the submarine to the external forces acting on it as follows:

(f - b) (2)

where s = [u, v, w,p, q, r]T is the velocity state vector, f is the vector of external loads applied to the rigid body, and b is an inertial coupling vector that arises because the equations are formulated in the non-inertial body-axes coordinate system. The components of each term in Eq. (2) are

Mass matrix:

mz, —my,

as follow:

mzG -my

—mzG 0 mxG

myg -mxG 0

Ix I j-xy I J zx

I j-xy Iy I 1yz

I - zx I Iyz Iz

Inertial coupling vector:

m[wq — vr — x,(q2 + r2) + y,pq + z,pr] m[ur — wp — y,(r2 + p2) + z,qr + x, qp] m[vp — uq — z,(p2 + q2) + x,rp + y,rq] (Iz — Iy )qr + Izxpq + Iyz(r2 — q2) + Ixypr + m[y,(vp — uq) — z,(ur — wp)] (Ix — Iz )rp — Ixy qr + Izxp — r2) + Iyzqp + m[z,(wq — vr) — x,(vp — uq)] (Iy — Ix)pq — Iyzrp + Ixy (q2 — p2) + Izxrq + m[x,(ur — wp) — y,(wq — vr)]

External load vector:

-(W - B)sin0 + XH(t, s, s) (W - B) cos 0 sin $ + YH(t, s, S) (W - B) cos 0 cos $ + ZH(t, s, S) (yG W - yBB) cos 0 cos $ - (zGW - zBB) cos 0 sin $ + KH (t, s, s) -(xGW - xBB) cos0 cos $ - (zGW - zBB) sin0 + MH(t, s (xGW - xB B) cos 0 sin $ + (yGW - yBB) sin0 + NH (t, s, s

where m is submarine mass, Ix,Iy, Iz are the moments of inertia, and Ixy, Izx, Iyz are the products of inertia. The centre of mass^(CG) and centre of buoyancy (CB) are located in the body axis coordinate system by (xG,yG,zG) and (xB,yB,zB), respectively. The six rows of the load vector f correspond to the following components: the axial force X in the x-direction, the lateral force Y in the y-direction, the normal force Z in the z-direction, the rolling moment K about the x-axis, the pitching moment M about the y-axis, and the yawing moment N about the z-axis. There are three types of forces and moments acting on the submarine: weight W, hydrostatic buoyancy B, and hydrodynamic forces and moments FH = (XH ,YH ,ZH ,KH ,MH ,NH). These forces and moments are transformed to the body axes directions based on the roll, $, pitch, 0, and yaw, ^ angles of the submarine (defined below). The weight of the submarine includes the flood water in the ballast tanks, which is assumed to move with the same rigid body motion as the submarine. The buoyancy is equal to the volume e submarine contained in the external hydrodynamic envelope multiplied the density of seawater. Its point of application is at the centre of buoyancy, which is located at the centroid of the appended submarine. The equations used to evaluate FH are presented in Sections 2.4 and 2.5.

With the exception of the modelling assumptions that go into the evaluation of the hydrodynamic loads, the above equations are exact for a rigid body with constant mass properties. In the simulation of the emergency rising manoeuvre in this work, the mass properties of the submarine change with time because

water is blown out of the ballast tanks using compressed air. The mathematical model used to simulate ballast tank blowing is presented in [18]. The submarine mass properties are updated in a quasi-steady manner each time step in the simulations presented herein. That is, additional terms described by Strumpf [31] modelling mass exit velocity (requiring ballast tank opening sizes and locations) and rates of change of mass and moments of inertia are neglected in the momentum change being modelled in Eq. (2). Mackay [32] found that these terms had only a small effect in a ballast blowing rising manoeuvre, probably because ballast tank flood water mass is a small percentage of combined vehicle and added mass and the mass change occurs relatively slowly.

A full description of the submarine state includes its position and orientation in inertial space. The 6 DOF momentum equations above do not, by themselves, provide this information. An inertial (Earth-fixed) coordinate system needs to be defined for this purpose. In this work, the inertial coordinate system is oriented with the x0 axis pointing North, the y0 axis to the East, and the z0 axis towards the centre of the Earth. The position of the submarine is given by (xo,yo,zo) and the orientation is defined using roll ($), pitch (0), and yaw (—) Euler angles. These angles are defined as follows: if the initial orientation of the body axes is co-incident with the inertial coordinate system, the final orientation is found by 1) yawing by the angle - about the z-axis, 2) pitching about the intermediate y-axis by the angle 0, and 3) rolling about the x-axis by the angle $, where the order is important. A set of six auxiliary kinematic relations transform the submarine's body-axes velocities to rates of change of orientation and inertial position [1]:

xo =u cos 0 cos — + v (sin $ sin 0 cos — - cos $ sin —)

+ w (sin $ sin — + cos$ sin 0 cos —) (3a)

y0 =u cos 0 sin — + v (cos $ cos — + sin $ sin 0 sin —)

+ w (cos $ sin 0 sin — - sin $ cos —)

¿0 = — u sin 0 + v cos 0 sin $ + w cos 0 cos $

<p =p + (r cos $ + q sin $) tan 0

0 =q cos $ — r sin $ r cos $ + q sin $

(3d) (3e) (3f)

Eqs. (2) and (3) can be combined together to form a set of 12 first order, ordinary differential equations that have the following implicit form:

Ay = g(t, y, y )

I implicit

where A is a 12 by 12 coefficient matrix that is composed of the 6 by 6 mass matrix in the upper left corner and a 6 by 6 identity matrix in the lower right corner and zeros elsewhere, g is the RHS of Eqs. (2) and (3), and y is the submarine state vector:

P, q,r,xo,yo,zo,$, 0,

2.4- Evaluating Hydrodynamic Loads

In this work, the hydrodynamic loads FH on the right hand sides of the first six equations of Eq. (4) are evaluated primarily using the URANS equations. However, some coefficient-based modelling is retained to account for the effect of the propeller and control surface deflections in order to reduce the complexity and computational requirements for the CFD solution (which is orders of magnitude more computationally intensive than solving a coefficient-based model). These loads are superimposed on hydrodynamic loads obtained from a CFD solution for a submarine with fixed tail fins and no propeller as discussed in [15].

To summarize, the hydrodynamic loads are broken down as follows:

XC + XP + XCFD

Yc + Ycfd Zc + ZcFD Kp + KcFD Mc + McFD Nc + Ncfd

where FC = [XC, YC, ZC, KC, MC, NC ]T are the modelled control forces due to tail fin deflection, XP and KP are the modelled propeller thrust and torque, and Fcfd = [Xcfd,Ycfd,Zcfd,Kcfd,Mcfd,Ncfd]t are the hydrodynamic loads obtained from the URANS solution. There is no roll-control (KC = 0) because the tail fins do not deflect differentially. The tail fin and propeller algebraic models are a straightforward function of y and computationally fast to evaluate.

2.5. Fluid Equations of Motion

A large computational fluid domain is created around the DRDC standard hull-sail-tail submarine geometry as shown in Figure 2 to discretize the URANS equations for the 6 DOF simulations. The far-field external boundaries were placed at least 2.5 submarine lengths away from the submarine based on a preliminary study that showed that an increase in the minimum clearance from 2.5L to 7L results in less than 0.1% change in hydrodynamic load predictions. domain is filled with a fully-structured hexahedral mesh with a total of x 106 cells; a near body cross-section through the symmetry plane is shown in Figure 3. The first node spacing above the submarine boundary is such that the average y+ values are around 50 for a full-scale Reynolds number of 175 x 106; wall functions are thus employed in the boundary layer. It was shown in a previous work that good agreement with experiment is obtained with a close variant of this mesh (StrCrs in [12]) for the same y+ values over a range of low to moderate angles of drift and angles of attack.

4.46 x

ogy of the 4.46 X 106 hexahedral

Figure 2: Fluid domain with 3 cross-sections showing thi cell mesh.

In a previous study [15], we solved the fluid equations in the body frame of reference and used a combination of spatially varying boundary conditions and momentum source terms to account for apparent linear, angular, centripetal, and Coriolis fluid accelerations due to frame motion. It was not used in the present work because the custom executable of the commercial code that applied the rotational momentum source terms became unavailable. An alternative approach is taken whereby the URANS equations are solved in a frame-of-reference that translates with the submarine but does not rotate with respect to the in-al frame. The fluid xf,yf,zf coordinate directions are always aligned with inertial xo,yo, z0 coordinate directions and the origins of the fluid and body coordinate systems are coincident. Submarine rotation is handled by rotating the submarine boundary relative to the fluid frame via mesh motion. This requires the mesh to deform throughout the domain, as illustrated in Figure 3, because the external boundaries are made stationary relative to the fluid frame,

ertial the in

Figure 3: Surface mesh on the submarine and a cross-section through the symmetry plane of the mesh in the vicinity of the submarine a) at t = 0 and b) after the submarine has pitched up in the emergency rising manoeuvre (S9, t = 30 s).

a choice imposed by a solver restriction on mesh motion at some boundaries1. Although mesh deformation adds complexity, it simplifies the treatment of ex-

ternal boundary conditions. A separate study by Bettle [12] showed that the moving mesh approach used here produces the same result as the body frame approach in [15]. The proposed predictor-corrector integration scheme is independent of the method used to handle 6 DOF vehicle motion in the U solution.

An arbitrary Lagrangian-Eulerian (ALE) formulation of the continuity and URANS equations is used to account for the motion of the computational mesh. In tensor notation we have:

dp + d (p (Uj - umj)) = 0

dt dxj

d(pUi) + d (pUi (Uj - Umj)) = _ dp + _d_ dt dxj dxi dxj

dUi dUj dxj dxi

—pu'uj +Fb

The mean and fluctuating components of the instantaneous fluid velocities (in the fluid frame of reference) are denoted as Ui and ui, respectively, p is the density, p is viscosity, and P is the mean component of the instantaneous pressure. The velocity of the mesh is denoted by Umj. The fluid velocity is referred to by its three components in subsequent sections using the subscript f to distinguish them from submarine velocity, i.e.: Ui = uf ,vf ,wf for i = 1, 2, 3.

The Reynolds stress term, puiuj, arises due to Reynolds averaging of the Navier-Stokes equation to account for turbulence. As in any RANS simulation, a turbulence model is needed for this term in order to close the equations. In this work, Menter's shear stress transport (SST) model [34] is used as it is a well-established model that has been found to model flow separation in adverse pressure gradients well [34, 33]. The fc-w formulation allows the model to be used stably down through the viscous sublayer and the switch to the k-e model

1ANSYS CFX v12.1 documentation indicates that "solution errors will be introduced if any of the specified [mesh] motion is normal to the boundary faces" for "inlets", "outlets", and "opening" boundaries (Chapter 2 of [33]).

eliminates the problem of the k-u model being overly sensitive to free-stream conditions. The transport equations for the SST turbulence model are written

in ALE form as follows:

d(pk) d (pk (Uj - Umj)) = _d_

dt dxj dxi

P akJ dxi

+ Pk - ¡3'pku (8)

d (pu) + d (pu (Uj - Umj)) = _d_ dt dxj dxj

+ pA d—

a,,, dxj

(1 - Wp—dd^dr + kPk - 3ipu2 (9)

ak — dxj dxj k

where k is the turbulent kinetic energy per unit mass, Pk is the shear production of turbulence, u is the turbulence frequency, Fi is a blending function, and a1, 3i, 3', a,, and ak are closure coefficients. The unknown Reynolds stress tensor,

puiuj, is modelled in the SST turbulence model with the following equation for eddy viscosity:

Vt aik nr.,

vt = — = --^ttt (10)

p max (aiw, SF2)

where a1 is a constant, S is an invariant measure of the strain rate and F2 is a

blending function.

The body force on the right hand side of Eq. (7), FB, is a momentum source term that accounts for the linear acceleration of the fluid domain relative to the inertial frame:

FB = -P^) (11)

where R0 is the position vector locating the submarine with respect to the inertial coordinate system. When the fluid coordinate system is always aligned w,hthe„1ertla1co°[d,„ate8y8temas,8theca8eforthep[ese„t.8tudy,t,

acceleration in Eq. (11) can be written as follows:

dt2 i= xoi + yoj + zok (11a)

where (i, j, k) are unit vectors in the (x0,y0,z0) (and thus also the (xf ,yf ,zf )) directions, respectively. Note that FB contains additional terms for domain rotation when the body frame of reference is used, as described in [15].

The mesh deformation is computed in this work by solving the following partial differential equation for mesh displacement:

dxj ( dxj ^ ( )

where xi = Xi — x0 is the mesh displacement and r is a diffusion coefficient, analogous to a spring stiffness, that can be changed to affect the way the deformation is diffused into the domain. In regions of high stiffness (large r), the mesh moves as a rigid body and there is very little distortion. For this reason, the stiffness is generally set high near the surface of solid boundaries to preserve the mesh quality in these regions.

The hydrodynamic loads, Fqfd, are obtained by integrating the URANS solution for pressure and shear stress distributions on the submarine boundary.

Csolut: equa

2.5.1. Boundary Conditions

The submarine boundary is considered to be a hydrodynamically smooth surface and is given a no-slip boundary condition that enforces the fluid velocity at the boundary to be equal to the boundary velocity.

The boundary conditions for the outer boundaries are set based on the assumption that the ocean is undisturbed at that distance from the submarine. Under the further assumptions that the submarine is deeply submerged in a quiescent ocean (no waves or current), this implies that the inertial-frame fluid velocity and the pressure are zero on these boundaries. In reality, the far-field pressure varies with depth, but this hydrostatic field is not included in the CFD solution because it is accounted for by the buoyancy terms in the rigid body ation of motion for the submarine. A combination of pressure-specified and velocity-specified boundary conditions are used for the 6 external boundaries as follows:

• Boundary —xf (behind submarine at t = 0): Outlet with zero average static pressure (pave = 0).

• Boundaries +xf, —zf (in front and above submarine at t = 0): Velocity

= 0 in the inertial frame as implemented using Eq. (13) below, turbulent intensity and eddy viscosity ratio set to 5%.

Boundaries +yf, -yf, +zf: Static pressure = 0. Velocity gradient perpen-

dicular to the boundary = 0.

The velocity on the +xf and -zf boundaries is set relative to the movi frame. It is set to the negative of the translational velocity of the s impose zero velocity relative to the inertial frame:

¡ubmarme to

U6 = - (Xoi + yoj + Zok) (13)

where Ub is the fluid velocity set on the boundary.

All transient simulations are initialized using results for a steady simulation with the submarine moving straight ahead in the xo direction with its initial velocity uo = X0o, where the subscript "o" denotes an initial condition. Specific conditions for the emergency rise and zig-zag are presented in Section 5.

3. Predictor Corrector Integration Method

A six DOF manoeuvring simulation can be considered an initial value problem in which the 12 non-linear coupled ordinary differential equations describing the vehicle motion state (Eq. (4)) are integrated in time, given a set of initial conditions. As discussed in the Introduction, an implicit predictor-corrector method based on the two-step, second order backward difference formula (BDF2) is selected. The same implicit method is used for coefficient-based simulations. In fact, the only difference between the coefficient and URANS simulations presented below is how the hydrodynamic loads FH are evaluated; the coefficient based simulations are extremely useful as a means of rapidly developing and verifying the algorithms eventually used for the URANS simulations, as well as for comparison purposes.

A schematic of the predictor-corrector method is shown in Figure 4. Note that subscripts used throughout this section indicate the time level while superscripts in square brackets indicate the iteration, j. Three main steps are

Initial Conditions

[ Start

Predictor

Evaluation

Corrector

-( J > mc?1

<r1 Yes

-{t > tend1 )

,----------Evaluation---------1

URANS Computation (Nf iterations)

Control Model

Propeller Model

Xp, Kf

Tailfin Model

Compute remaining terms on RHS of Eq. 4

Ballast Blowing Model

To Corrector

discrete

1: Predictc

3r-Correct(

Figure 4: Predictor-Corrector solution procedure.

performed at each discrete time level in the solution: a predictor step (P), an evaluation step (E), and a corrector step (C). The solution for the "new" time level, tj+i, begins with the predictor step, in which a prediction is made for the submarine state, y[+1, and state derivative, y[+1, using the known state and state derivative of the submarine at previous equally-spaced time levels ti,ti-i, ■■■,tj-k. The predictor acts as an initial condition for subsequent evaluation and corrector steps, which are performed iteratively NC times each time

O level. A useful mnemonic for describing the procedure is P(EC)Nc , which indicates that one predictor step is performed followed by NC Evaluation and Corrector steps at each time level.

The evaluation step consists of evaluating the coefficient matrix A and g(t, y, y) in Eq. (4). A summary of the evaluation of the individual components of g(t, y, y) using a combination of CFD and coefficient-based models is

shown in Figure 4. The evaluation step is by far the most computationally expensive step in the procedure and will be discussed in further detail in Section 3.5 below.

In the corrector step, the state of the submarine is "corrected" in a process of solving and integrating Eq. 4 with an implicit method described below.

3.1. Predictor

A predictor based on the assumption of constant submarine acceleration was found to provide a good initial guess at each time level. Since the first six components of y are accelerations whereas the remaining six components are velocities, they are treated separately in order to approximate constant acceleration. Acceleration components are held constant while velocity components are linearly extrapolated:

stant whil

y[+i[1 - 6]= yi (14a)

yi+i[7 - 12] = 2.i - yi-i (14b)

All 12 components of the predicted submarine state yi+]1 are obtained by using the above prediction for y i+1 with the BDF2 formula:

e predicte on for y[0]

▼L°] 41 , 2 . , . [0] n ^

yi+1 = 3yi - 3yi-1 + gAt yyi+1 (15)

where At is the time step size. Eqs. (14) and (15) constitute a predictor step. The BDF2 is normally associated with implicit time integration; however its use in the predictor is explicit because y[+1 is calculated directly using known

values of the submarine state in Eqs. (14). 3.2. Corrector

The objective of the corrector iterations is to converge on a solution for the integrated submarine state yi+1 at each fixed time level. Note that the implicit nature of Eq. (4) requires that a converged solution for yi+1 be obtained at each time level as well. The relationship between and yi+1 is fixed at each time

level by the integration method. As with the predictor step, the second order BDF2 is used, except this time it is applied in an implicit manner:

yj+i] = 4 y 1 y + 2Ai y j+1l (16) yi+i = 3yi - 3y— + 3Ai yi+i (16)

This relationship is used to eliminate the direct dependence of Eq. (4) on y so that it can be written as a direct function of only y at the fixed time level ii+1 so that at convergence:

[A (yi+i)] yi+i = g(yi+i) (17)

if if the sub

Note that the coefficient matrix A depends on the state of the submarine in scenarios with mass change due to ballast blowing.

The procedure at each corrector iteration is to first solve Eq. (17) for y j+Li], then substitute yi++i] into Eq. (16) to calculate yi++i]. Two commonly used methods for solving Eq. (17) are fixed-point iteration and Newton iteration. Both methods are evaluated in terms of their iterative efficiency for 6 DOF submarine manoeuvring simulations.

3.2.1. Fixed-Point Iteration

To use fixed-point iteration, which is also referred to as the method of successive substitutions, Eq. (17) must be placed in the following form:

y = F (y ) (18)

This can be done by setting

his can be do &

F (y ) = [A (y r1 g (y ) (18a)

xed-point iteration is defined by the following recursion relation:

y[J'+1] = F (yt^) (19)

The solution of Eq. (19) is initiated with an initial guess for y provided by the predictor. Fixed-point iteration is guaranteed to converge to a unique value for an arbitrary starting point if the function F (y) satisfies a Lipschitz condition

\\F (y) - F (y*)|| < L yy - y*y (20)

for all y, y*, where the Lipschitz constant L satisfies 0 < L < 1 ([35], pg. 12). Here, y* is the value of y that satisfies Eq. (18). If F (y) is a differentiable function, the above Lipschitz condition is equivalent to requiring the absolute value of the derivative dF/dy to be greater or equal to zero but less than 1 in the region of interest. For the rigid body equations of motion applied to the fully submerged DRDC standard submarine considered in this work, this Lipschitz condition is not met and the solution diverges with standard fixed-point iteration. The large unsteady hydrodynamic loads due to submarine acceleration, the so-called "added-inertia" forces, were found to be the primary cause. It is shown by Bettle [12] that if the ratio of the added mass to the mass of the submarine has a magnitude greater than 1 (in any of the six DOF), fixed-point iteration will diverge. Several components of added inertia are on the order of 1 for marine vessels immersed in water because the density of the surrounding fluid is approximately the same as the body itself.

Although standard fixed-point iteration was found to diverge, stable solutions can be obtained by using the following modification of standard fixed-point iteration:

(1 - Y) y[j] + (Y) F (y [j] ) (21)

ion parame

where 7 is a relaxation parameter that blends the old value of y with the function evaluation. When 7 is 1, Eq. (21) reduces to standard fixed-point iteration. An under-relaxation parameter in the range 0 < 7 < 1 dampens the oscillations that occur from iteration to iteration, and if 7 is set to a value smaller than ritical value, the oscillations are completely eliminated. However, making 7 mall will slow convergence due to excessive damping.

a critical v too small i

3.2.2. Newton Iteration

Newton iteration, also referred to as the Newton-Raphson method, is a well-known procedure for solving scalar non-linear equations of the form F (y) = 0. It uses the first derivative of the non-linear function, F', to guide convergence

according to the following recursive relation:

V[j] -

F(yj])

Eq. (22) is applied to the system of equations to be solved (Eq. (17)) as follows [35]:

y [j+1] = y j] -

J-1 yj]

y ¿+1

g (y j])

:obian matri:

and JF (y[j+1]) = (dF/dy)|y[i+1] is the 12 x 12 Jacobian matrix of F (y [j]) with respect to y. Rather than invert the Jacobian matrix, Eq. (23) is rearranged as follows for solution using LU decomposition

Ä y[j]

where Ay,

[j+1] _ y [j] ^¿+i y ¿+1

is the increment to be added to the old iterate to

correct it. Eq. (24) is evaluated at every iteration j for yj , which is used in Eq. (16) to update y[j+1] according to second order BDF2 integration. The solution of Eqs. (24) and (16) constitute a corrector step using Newton iteration. Note that the RHS of Eq. (24) is the imbalance between the RHS and LHS of the submarine equations of motion. At convergence, this imbalance goes to zero as does the correction Ay[j].

An advantage of Newton iteration over fixed-point iteration (or relaxed fixed-iteration) is that a second order convergence rate is achieved in the former only first order convergence is obtained in the latter. This benefit, however, is achieved at the expense of evaluating the Jacobian matrix each iteration. In practice a modified Newton iteration is often employed, in which the Jaco-bian matrix is not updated at every iteration, or even every time step, in order to save on computational effort. This can be done because the Jacobian matrix need not be exact, since at convergence Ay^ goes to zero; slight errors in JF affect the convergence rate but not the converged solution. Updates for JF can

point while

be delayed until convergence becomes slow or cannot be achieved in a time step. In the present work, JF is updated each corrector step because its evaluation is orders of magnitudes less expensive than the CFD solution.

3.3. Estimating the Jacobian Matrix Using Coefficient-Based Models

An approximation from theory can provide an affordable estimate for the Jacobian matrix needed for Newton iteration. In this work, the coefficient-based hydrodynamic model [18] for the DRDC-STR submarine provides a good estimate for JF. It should be noted that some of the most important terms in JF are added mass coefficients, which can be estimated from theory reasonably accurately. Also, the derivatives of most hydrodynamic loads in the coefficient-based model have simple analytical expressions.

In order to show the importance of different components making up JF, it is expanded as follows:

Jf = dF = + dAy[J] - dy (25)

dy dy dy

where F is given in Eq. (23). The term dAy[j] is non-zero only for simulations in which the submarine mass changes as a function of submarine state, as is the case for simulations with ballast blowing (the submarine mass properties are a function of depth, pitch angle, and time in the ballast blowing model). This term is neglected in this work because it adds significant complexity, which does not appear to be justified based on the excellent iterative convergence achieved without it. The last term in Eq. (25) can be further expanded as follows:

dg = F + OFrb + 0Fam + OFn + m + dFp + dFc (26)

dy dy dy dy dy dy dy dy

where Fg are hydrostatic loads (terms multiplied by weight W or buoyancy B in the external load vector f of Eq. (2)) and FRB accounts for the remainder of terms in Eqs. (2) and (3) that are explicit in y. The derivatives of Fg and FRB can be evaluated analytically directly from Eqs. (2) and (3) and are the same for any 6 DOF manoeuvring simulation that use these governing equations. The remaining terms in Eq. (26) represent different components of hydrodynamic

loads in the coefficient-based model: Fam are unsteady "added mass" loads accounting for linear and angular acceleration, while Fr, Ft, FP, and Fc are loads due to steady rotation, steady translation, propulsion, and control surface deflection, respectively.

Note that, with the exception of dFAM , the chain rule of differentiation to be applied when e« a!, co.m'ponents of f, because tlley are e: functions of y rather than y:

dg _ dg dyi+i

dyi+i dy.

where, for the second order BDF2 selected for integrati

explicit

integration (i

. (16)):

dy i+i 3"

Eq. (28) shows that most terms in the Jacobian matrix are timestep size dependent, with the notable exception of added mass coefficients.

Since the coefficient-based model is only an approximation to the true Jaco-bian matrix of hydrodynamic loads in the unsteady RANS simulations and there is uncertainty in the coefficients, the added complexity of incorporating all nonlinear and smaller secondary coefficients may not be worthwhile. As summarized in Table 1, three levels of simplifications are applied to the coefficient based model in an attempt to find a model that achieves a similar performance as the full model, J0, but is easier to implement. Model JL is obtained by linearising

the hydrodynamic model about zero incidence (© = tan-1 (%/v2 + w2fu) = 0) while retaining all added mass coefficients and all derivatives associated with control surface deflection and propulsion. The derivatives for steady translation are found for a given forward velocity u by taking the derivatives of Fuvw [18] with respect to u,v,w at zero incidence:

d(uk )

where F = X,Y, Z, K, M, N for and uk = u,v,w for k = 1, 2, 3. There are a total of 6 x 3 = 18 steady translation derivatives at zero incidence, but only

Hydrodynamic coefficients included in each model for JF

Fam Fr Ft Fc

All All F J- uvw All

Xuq Xuu ■ Xuw

Y Y up ur Y 1 uv

All Zuq Z / Z u ^uw All

K K up ur K uv

Muq Muu■ Muw

Nup ■ Nur Nuv i

Y- Y Y- 1 V > 1 p> 1 r Y Y up ur Yuv

Zw ■ Zq Zuq Zuw

K K K K K up ur Kuv 0

Mw, Mq Muq Muw

Nv , Np, Nr Nup ■ Nur Nuv

Same as JLp 0 0 0

Table 1: Hydrodynamic coefficients included for each model of Jp. All models contain the hydrostatic/gravitational component (Fg) and the component due to terms explicit in y in Eq. (2), Frb . The baseline case, Jo, contains all hydrodynamic coefficients in [18] and all numerical values are presented there.

the 9 shown for JL in Table 1 are non-zero for a submarine with one plane of symmetry.

Model JLp is the same as JL except derivatives associated with control surface deflection and propulsion are omitted and only the primary coefficients that are easiest to estimate or measure experimentally are retained. The simplest model, Jamo, retains only the primary added mass matrix for the hydrodynamic component of the Jacobian matrix.

The effect of uncertainty in hydrodynamic coefficients on iterative convergence is investigated in Section 6 by perturbing the hydrodynamic coefficients in JF by ±30%, a value selected based on a pessimistic estimate for the accuracy

with which primary coefficients can be estimated.

It should be noted that if hydrodynamic coefficients (or methods to estimate them) are not available, they can also be computed using steady and unsteady RANS simulations. Model JLp, for example, could be obtained by performing 11 preliminary RANS simulations before starting the 6 DOF simulation. These would include: one steady simulation at zero incidence and with the nominal forward velocity for the manoeuvre, 5 steady simulations with small perturbations (Av, Aw, Ap, Aq, Ar) superimposed on the nominal forward velocity , and 5 unsteady simulations with constant accelerations (V,W,W,p,r) to get added mass coefficients. This approach, left for future study, may have the advantage that the coefficients should approximate the Jacobian matrix in the 6 DOF RANS simulation well, especially if the same mesh, turbulence model, etc., are used for generating the coefficients as for the 6DOF simulation. The advantage of using an existing model or estimating coefficients with semi-empirical methods is that it does not require any significant computational effort.

3.4. Convergence

ng predict

One way of performing predictor-corrector integration is to allow iterations to continue until the change in submarine state from one iteration to the next falls below a specified tolerance. In this mode of "correcting to convergence", the total number of corrector iterations NC for a time step is not known in advance because it depends on the solution. The problem with this mode is that the solution can hang at a time step if the convergence criteria cannot be met. For this reason, it is common practice in predictor-corrector integration of ODE's to specify NC in advance, where NC is typically set to a small number

[35]. When using a small number of corrector iterations, the accuracy of the solution may be influenced by the predictor. An analysis by Lambert [35] shows that if the order of accuracy of the predictor is one less than the corrector, as is the case in this work, then if only one corrector step is performed (NC = 1), the principal local truncation error (PLTE) at each time level is affected by the PLTE for the predictor but the overall order of accuracy will be that of the

;ormed

oiution. 3rred to

varioi

corrector. If NC > 1 then both the PLTE and the overall order of accuracy are that of the corrector alone. At least two corrector steps are used in all URANS simulations performed in this work and thus the time integration is second order accurate, and the PLTE is that of the BDF2.

Note that the numerical solution of the fluid equations of motion perfo: dunng the evaluation step at each lM.io„ , a,so requ.res an .terative so, The fluid iterations are nested within the corrector iterations and are refe: as "coefficient-loop iterations"2, k. For each corrector loop, Nf coefficient-loop iterations are performed, making the total number of coefficient-loop iterations performed at a time level equal to NC x Nf.

3.5. Implementation in ANSYS CFX v12.1

The commercial software ANSYS-CFX v12.1 is used for solving the fluid equations of motion (Eqs. (6) - (12)). ANSYS CFX is a finite-volume CFD solver using element assembly of its equations and an implicit temporal discretization method. The hydrodynamic (uf, vf, wf, p) equations are solved as a single system using a coupled solver [36, 37] and the Additive Correction [38] Algebraic Multi-grid [39] procedure is used to accelerate the solution. Fully second order-accurate discretization is used for advection terms in the URANS equations for our simulations and the second order accurate backward difference formula is used for temporal discretization.

The predictor-corrector integration scheme for the 6 DOF submarine motion is implemented within the framework of the fluid flow solver as summarized gure 5. User Fortran subroutines are connected to the CFX solver at ious stages in the solution of the fluid equations for evaluating the auxiliary coefficient-based models, evaluating the predictor and corrector equations, and setting boundary conditions based on the state of the submarine.

2 Here "coefficient" refers to the coefficients in the discretized URANS equations and is unrelated to the "coefficient-based model" for hydrodynamic loads.

Figure 5: Computational solution procedure for URANS simulations.

4. Estimating Numerical Errors

CFD predictions contain both numerical errors and modelling error. Numerical error here is defined as the difference between the numerical solution and the exact solution of the governing equations, while modelling error is a measure of how well the governing equations represent reality. There are three main sources of numerical error: round-off error, iterative error, Ej, and discretization error, Ed ; see for example Eca and Hoekstra [40]. The main objective of the current work is to analyse the iterative convergence of the predictor-corrector integration scheme described in Section 3 when applied to simulating 6 DOF submarine manoeuvres. As a result, the primary verification work done in this study involves the estimation of iterative error and temporal discretization error for simulations conducted with a single 4.46 x 106 cell mesh. It was found in a previous work [12] that the hydrodynamic loads are well predicted with this level of refinement (the mesh used is a slightly modified version of StrCrs in [12]). To give an indication of the sensitivity to grid refinement, predictions for the six hydrodynamic load components at a drift angle of ¡3 = 20° change by less than 2.5% when the number of cells are increased by a factor of 2.8 (%/2 refinement in each direction) and even less for an additional %/2 refinement, such that the maximum change was 4.2% for an eight-fold increase in cells. These changes are small relative to the estimated 20% modelling error in SST RANS predictions of out-of-plane hydrodynamic loads at 3 = 20° [12]. Load predictions are even less sensitive to mesh refinement at smaller flow incidence angles, where results for all grids match experiment well (3 < 15°). The 4.46 x 106 cell mesh provides adequately low discretization errors for the purposes of the present work, without incurring excessive computational effort.

The iterative and temporal discretization errors are estimated for each unsteady 6 DOF simulation. A large set of coefficient-based simulations are conducted to analyse several aspects of the predictor-corrector method before performing computationally intensive URANS simulation. An advantage of the coefficient-based simulations is that they can be run to machine accuracy in

a reasonable amount of time. In this case, the iterative errors and temporal discretization errors can be evaluated directly as follows:

• The exact solution for submarine state, ye(t), is approximated by a solution with a timestep size, Ate, that is small enough for discretization errors to be negligible while also large enough for round-off errors to be insignificant.

• A solution with negligible iterative error, yc(t), is obtained by converging to round-off error each timestep.

• The temporal discretization error at a time level ti in a solution, EDi, is the difference between the solution, yi, and the exact solution ye for that time t.

i ETi, is

The L,

• The iterative error at a time level ti, E/i, is the difference between the solution, yi, and the fully-converged solution, yc, for that time t.

A similar approach is taken for the URANS simulations, except the simulations cannot be run to machine accuracy. In this case the discretization error for a solution is estimated by comparing against solution with At about one order of magnitude smaller, and iterative error is estimated by comparing with a solution that is converged at least an order of magnitude further at each timestep.

orm is used as an overall measure of the discretization and iterative errors for a given simulation. It is defined as the maximum error for all timesteps:

\\E=max(|Ei|), 1 < i < imax (29)

where Ei is the error at time step i and imax is the total number of time steps for the solution.

Double precision is used for all simulations in order to reduce round-off error as much as possible.

5. Manoeuvring Scenarios

5.1. Emergency Rising Manoeuvre Two types of manoeuvres are simulated to test the predictor-corrector method:

an emergency rising manoeuvre and a horizontal plane zig-zag. Emergency rising manoeuvres have been described by Watt [18], who simulated 9 limited and controlled rising scenarios using coefficient-based hydrodynamic models. In this manoeuvre, the submarine rises buoyantly to the surface by emptying the ballast tanks using compressed air. Small to medium sized submarines are known to develop large roll angles just before surfacing. The development of this roll instability has been previously modelled using unsteady motion-coupled RANS simulations with under-relaxed fixed-point iteration (see Section 3.2.1) [15]. In the present study, a rising scenario with a more extreme rolling motion is simulated with a focus on comparing the iterative performance of fixed-point iteration with Newton iteration employing coefficient-based models to estimate the Jacobian matrix. Larger time step sizes than used in [15] are considered to test the robustness of the method.

The simulation modelled is S9 from [18] using the same initial conditions and control parameters, as summarized in Table 2. At t = 0, the submarine is moving with straight-and-level flight at a speed of 3 m/s and depth of 100 meters. The predicted hydrodynamic loads for this speed are used to determine the required propeller self-propulsion advance ratio, Jbs, initial sternplane trim deflection Sst, and weight trim, Wt = Wo — B, where Wo is the initial weight of the submarine. The CG is assumed to be laterally offset from the submarine symmetry plane by yGo to induce an initial roll angle of -2 degrees. The commanded forward speed is set to 6 m/s at t = 0 and the sternplane deflection is given a control sequence to attain and maintain a pitch angle of around 10 degrees during the rise.

5.2. Horizontal Zig-Zag The predictor-corrector method is also tested with a horizontal-plane zig-zag

manoeuvre. In this manoeuvre, the submarine starts with straight-and-level

Sst (deg) 0.813

$o (deg) -2.0

VGo (m) -0.0114

Wt 1.64 x 10-4B

Jbs 0.9783

-9 @ t = 0 s

Ss commanded (deg) 0 @ t = 16 s

3 @ t = 40 s

uc commanded (m/s) 6 @ t = 0 s

Table 2: Simulation parameters for rising scenario S9, 1 efficient-based simulations

(R = 23 X 106) [18]

flight with constant approach speed Uo. At t = 0, the rudder is commanded to deflect to starboard by angle —Sr. This is maintained until the submarine has yawed to starboard by angle ^e relative to its initial course, at which point the rudder angle command is reversed to +5r. The submarine follows a serpentine path as the rudder is repeatedly reversed each time each time a yaw angle of is achieved. A 10°/10° horizontal zig-zag (the first number indicating the magnitude of 5r, the second ^e) with approach speed Uo = 6 m/s is simulated using both Watt's coefficient-based model [18] and 6 DOF URANS. The sternplanes remain at the initial deflection throughout the zig-zag; no compensation is made for depth changes. The commanded forward speed for the propeller is

O maintained at uc = Uo = 6 m/s throughout the manoeuvre. This means that the propeller RPM is held fixed at the speed required to maintain straight and —»speed var.es _ the due

to changing hydrodynamic loads on the submarine.

Following the procedure described in [18], the predicted hydrodynamic loads for u = Uo are used to set propeller Jbs, initial sternplane deflection 5st, weight trim Wt, and determine the initial roll angle that satisfy equilibrium conditions at t = 0. The parameters and initial conditions for coefficient-based and URANS

simulations of the horizontal zig-zag are tabulated in Table 3. The differences in initial conditions and propeller Jbs for the two methods are largely due to the prediction of a lower drag coefficient in the full scale URANS simulation (Reynolds number Re = 350 x 106) than the coefficient-based model, which is based on model-scale data (Re = 23 x 106). Note that these differences were neglected for the simulations of S9 in the work as the coefficient-based initial conditions and Jbs were applied to the URANS simulations.

Re at t = 0 23 x 106 350 x 106

Uo (m/s) 6.0 6.0

5r commanded (deg) ±10 ±10

5st (deg) 0.813 ^ 0.180

$o (deg) -0.56814 -0.35449

VOo (m) 0 0

Wt 1.64 x 10-4B 1.14 x 10~4B

Jbs 0.9783 1.1605

Table 3: Simulation parameters for the horizontal zig-zag manoeuvre.

6. Coeffici

Watt's complete coefficient-based model [18] for hydrodynamic loads was added as a subroutine to the code for simulating unsteady RANS simulations of submarine manoeuvres in CFX. This subroutine can be called to get hydro-dynamic loads from the coefficient-based model instead of solving the RANS equations; Fcfd is completely replaced with Watt's coefficient based model (FH in reference [18]). Since the code for predictor-corrector integration and auxiliary models for propulsion, tailfin loads, and ballast-blowing is the same for the two modes of simulation, it is ensured that any differences between the two prediction methods (i.e. RANS or coefficient-based) are entirely due to the way each method models hydrodynamic loads. A stand-alone Fortran program,

created for rapidly simulating a series of coefficient-based simulations using the same subroutines implemented in CFX, generated the same coefficient-based results as the CFX implementation.

Hundreds of faster than real time coefficient-based simulations were run to quickly test several aspects of the proposed predictor-corrector method and its implementation before performing the 6 DOF URANS simulations, which require on the order of two days to a week to complete with 40 CPUs running in parallel. This testing enabled us to quickly compare fixed-point iteration with Newton iteration, and identify a good model for the Jacobian matrix. The conclusions from the coefficient-based study are then verified with unsteady RANS simulations in Section 7.

6.1. Coefficient-Based Rising Simulations

Several coefficient-based rising simulations were conducted to test the predictor-corrector method. We first verify the order-of-accuracy and temporal discretization errors in the absence of iterative error by converging simulations with different At to machine accuracy using Newton iteration with the full unperturbed Jacobian matrix, J0. These BDF2 predictor-corrector simulations are compared with results from an explicit 4th-5th order Runge-Kutta integration method. We then compare the convergence behaviour for Newton iteration with different models for JF against fixed-point iteration with different under-relaxation parameters. A similar analysis is done for coefficient-based zig-zag simulations in Section 6.2.

CFor coefficient-based rising simulations, Newton iteration with J0 converged to machine accuracy at each time level with around 5 correctors. Fully converged results were obtained for At ranging from 10~5 to 5.0 seconds. The solution with At = 10~5 s is considered to be an "exact" numerical solution for the purposes of evaluating the discretization error of solutions with larger time step sizes. Cumulative discretization errors, EDk, are found for each component k of submarine state vector y at every time step in the solution by calculating the differences with the At = 10~5 solution at corresponding points in time:

3 10-4

I I I 11 ■ 11-1 I I ■ ■ I ■ ■ I-1 I I Mill)-1 I I ■ ■ I ■ ■ I-1 I I Mill)-1 I I ■ ■ I ■ ■ I-1 I I I II

RKF45 ( variable At ) J..

..............................................................

10-4 10-3 10-2 10-1 10° 101 At (s)

Figure 6: Effect of time step size on the temporal discretization error in coefficient-based simulations of rising scenario S9 for the second order BDF2. Results are compared with a 4th/5th order Runge Kutta simulation (RKF45) with variable At [18].

Enk (ti) = y[k](ti) - ye[k](t = ti) (30)

where k = u,v, w,p, q, r, x0, y0, z0, 9, ^ is the component of submarine motion state and ye is the solution with At = 10~5 s. The Lnorm (maximum value, see Eq. (29)) of EDk over the duration of the simulation, \\EDk is shown for each component of y as a function of time step size in Figure 6. The maximum discretization error drops by two orders of magnitude for every of magnitude drop in time step size for 10~3 s < At < 1s, confirming the second order accuracy of the BDF2 integration scheme. The error drops at a slower rate for At < 10~3 s, probably due to round-off error.

The magnitude of the differences between the Runge Kutta and At = 10~5 s BDF2 predictor-corrector simulations are also shown in Figure 6. The discrepancy ranged from 2 x 10~6 for x0 to 7 x 10~5 for r, when normalized by the maximum variation of y[k] during the simulation (i.e., y[k]max — y[k]mi„). This small discrepancy is slightly smaller than the discretization errors for the BDF2

simulations with At = 0.1 s and seems to be a plausible value for cumulative numerical uncertainty in the RKF45 solution, which had an average At around 0.5 s. The close agreement between the integration methods can also be seen in Figure 7, where the At = 10~5 s BDF2 is indistinguishable from the RKF45 solution for all submarine state variables. This close agreement verifies that the implementation of Watt's coefficient-based models, including the auxiliary propulsion, tailfin deflection and ballast-blowing models used in the URANS simulations, is consistent with [18].

Another conclusion from this study is that a stable and fully converged solution can be obtained with the BDF2 predictor-corrector method when using a very large time step size of t = 5 s, which is equivalent to the time that the submarine translates a half body length at the speed attained by the end of the simulation (7.6 m/s). It it is clear from results in Figure 7 that this large At provides only a crude approximation to the exact numerical solution but this is an important result because it indicates that the implicit method is robust enough to allow an appropriate timestep to be selected based on accuracy requirements rather than numerical stability. Although temporal discretization errors can also be seen for the At = 2 s solution in Figure 7, it is considered to be a very good approximation for preliminary computations because modelling errors are expected to be significantly larger. This ability to quickly obtain preliminary results with large time step sizes is even more of a benefit for CFD simulations due to the much larger computational expense.

A second set of coefficient-based simulations of rising scenario S9 were performed to analyse the iterative convergence for Newton iteration with different models for the Jacobian matrix and fixed-point iteration with different values for relaxation parameter 7. Each simulation uses a fixed number of correctors per time step, NC, which was varied from 1 to 10 for At values of 0.1 s, 0.25 s, 0.5 s, 1.0 s, 2.0 s, and 5.0 s. The cumulative iterative error for each component of y, Eik, is evaluated for each simulation by comparing results with a simulation having the same time step size that is converged to machine accuracy, yc:

7.5 7 6.5 6

5 0 -5

^-10 -O

7 ' </

Jf f / -

A*- *■

............. ......, 1 , , , , 1 , , , ,

0 10 20 30 40 50

-1 ' -1.5 -2 -2.5 3

I I I I | I I I I | I I I I | I I I I | I I I I | I I I I

, -15 -20 -25 -30

■ 1 11 11 1 11 11 1 11 11 1 11 11 1 11 11 1 ■ 11 ■

11 10 9 8 7

5 3 z 100 4 ® £

3 2 1 0

0 10 20 30 40 50

■ 1 ■ 11 1 ■ 11 1 11 ■

1 0 -1 -2 -3

-5 -6 -7 -8 -9

20 ■

100 ■

—|-1-1-1-1-1-1-1-

- BDF2, At = 10-5 s

° BDF2, At = 2.0 s. .............. BDF2, At = 5.0 s.

....... RKF45

T-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1—

*0 (m)

Figure 7: Coefficient-based simulations of rising scenario S9, computed using the second order predictor-corrector method (solid line and circles) and a 4th-5th order Runge-Kutta integration scheme with variable time step obtained in [18] (dashed line).

Eik(ti) = y[k](ti) - yc[k](ti)

The Lo norm iterative errors (maximum over time domain, Eq. (29)) for S9 roll angle using Newton iteration with different models for JF (see Table 1) are compared with fixed-point iteration with various relaxation parameters in Figure 8. Results are only shown for roll angle (the primary variable of interest) for simplicity, but the same trends are seen for the other components of submarine motions state. The magnitude of iterative errors are normalizing by \\EDk ||o in these plots. Note that the errors shown are maximums of cumulative errors rather than local errors and that each data point in Figure 8 represents a complete simulation (almost 800 simulations in total) with NC correctors used at each timestep. The main findings from the results in Figure 8 are summarized as follows:

e results

Excellent convergence is achieved with Newton iteration using the full unperturbed Jacobian, J0. For practical At ranging from 0.1 to 2s, only a single corrector is required each timestep for iterative errors to be less than 1% of the discretization error. Iterative errors drop by at a rate of almost 2 orders of magnitude per corrector for At = 2 s, and even faster as At is reduced. Good convergence is even obtained with a very large At of 5 s, with iterative error decreasing at least one order of magnitude per corrector.

The convergence rate for fixed-point iteration increases as less under-relaxation is used (7 is increased), but only up to a stability limit, Ymax, beyond which the iterative solution becomes unstable. As At is increased, this limit becomes more restrictive and more under-relaxation is required for a stable solution. The optimum 7 for convergence is close to the stability limit, as evidenced by the marginally stable At = 2 s, 7 = 0.35 result. Even when a near-optimum 7 is used, fixed-point iteration is much slower than Newton iteration (with J0): it takes 5-6 correctors each timestep to

102 101 100 -

~ 10-1 _s

10-3 10-4 102 101

~ 10-1 _8

W^ 10-2 10-3 10-'

At' = 0.1s.—■

102 101

~ 10-1 _s

wS 10-2 10-3 10-'

J0 -®-JL 1 1 JAMo

J -e-JLp J + AMo

J0 A JlP V JAMo

Fixed : point it. ' ~~'y = 0.4 ]

"""---,/0.45,. ""'0.5,1

~i-1-r

At = 0.5 s.

% 100;

~ 10-1 _8

wS 10-2 10-3

10-40 102 101

W 10-1

1010102 101

">-1-r

~ 10-1 _8

ws10-2 10-3 10-'

Figure 8: Maximum iterative error for roll angle during coefficient-based S9 simulations, normalized by the maximum discretization error, ||-Ed>||~. Solid lines: Newton iteration with models for the Jacobian matrix according to Table 1; dashed lines: fixed-point iteration with an optimum relaxation parameter 7. Superscripts "+" and "-" indicate that all hydrodynamic coefficients in the Jacobian matrix are multiplied by 1.3 and 0.7, respectively.

reduce iterative errors to below 10% of the discretization error for At < 1 s, and more than 8 for At > 1.

Although the rate of convergence of Newton iteration is reduced when linearised models JL and JLp are used in place of Jo, good smooth convergence is still obtained, which significantly outperforms fixed-point iteration for all At. Results for JL and JLp are almost identical, indicating that coefficients included in JL that are neglected in JLp are not important for iterative convergence for this manoeuvre. Only 2-3 correctors per timestep are needed to achieve Ej/ED < 0.1 for At < 1s, increasing to 4 for At = 2 s. This is about half the requirement for fixed-point iteration; this gap widens further for stricter convergence levels.

The performance of Newton iteration with a Jacobian matrix that only contains primary added mass coefficients for the hydrodynamic component, Jamo, is very similar to that of JL and JLp for At = 0.1 s and 0.25 s. However, its performance degrades quickly as At is increased above 0.5 s, to the point of being marginally stable at At = 2 s and divergent for At = 5 seconds. This result is not unexpected because the added mass terms are constant with At while all other terms in the Jacobian grow in proportion to At (see Eq. (28)). The inclusion of the primary trans-lational and rotational coefficients contained in JLp becomes increasingly important as At is increased. Note that an accounting of primary added ass coefficients is essential at all At because simulations diverge when ey are removed from any of the models.

• Newton iteration with J0 and JLp still significantly outperforms fixed-point iteration when all hydrodynamic coefficients are over-predicted by 30% (J+) or under-predicted by 30% (J—) for all At. For this level of uncertainty, both Jacobian matrix models converge at approximately the same rate of about 5 orders in 6-7 correctors when At < 0.25. There is no advantage to including the non-linear hydrodynamic terms for these timestep sizes, which may be used for final results. At larger At, there

appears to be an advantage to including the non-linear terms if errors in all coefficients are at the 30% level. A comparison of J0 and JLp will be revisited in the context of URANS simulations.

ated, es-led mass

The effect of perturbing added mass coefficients in the JaMo model de pends strongly on whether they are overrated or underpin» pecrahy at ,arger At. Convergence is better when the primary added

coefficients are overestimated, and for At > 0.5 s, convergence is actually improved over the unperturbed case. This is believed to be because a result of the over-estimation in added mass compensating for the missing linear translational and rotational coefficients.

This study demonstrates that Newton iteration with a suitable model for the Jacobian matrix is superior to fixed-point iteration for converging the predictor-corrector integration of implicit 6 DOF equations of motion. The simplified hydrodynamic model for the Jacobian matrix, JLp, provides a good trade-off between complexity and performance in the simulation of the emergency rising manoeuvre. There does not appear to be any advantage in including the "secondary" coefficients contained in JL that are not in JLp because there is no difference in convergence between these models. These findings will be tested with a similar study for the horizontal zig-zag manoeuvre. Also, the ability of the Jo and JLp models to handle significant uncertainty in hydrodynamic coefficients suggests that it will also perform well for URANS simulations. The true test of this will be shown by 6 DOF URANS simulations of the same manoeuvres in Se ction 7.

6.2. Coefficient-Based Zig-Zag Simulations

The coefficient-based testing with rising scenario S9 was partially repeated for the 10°/10° horizontal zig-zag manoeuvre. The predicted state of the submarine during the zig-zag is shown in Figure 9.

The yaw angle is predicted to oscillate between approximately ±17° in response to the rudder deflections, giving an overshoot yaw angle of about 7

)s i 0

12 5.6

5.4 5.2 5

-At = 10 s. -At = 1.0 s. At = 2.0 s

uooooouuuo"nnn

uumiooo"""1""'""""

^uuuonDUUOunno

Trajectory. (Ordinaterabscissa scale = 2:1)

x0 (m)

0.4 0.2 0 --0.2 -0.4

Figure 9: Coefficient-based prediction for submarine state during a 10°/10° horizontal zig-zag manoeuvre at 6 m/s.

degrees. The peak yaw angles are achieved between 6.5 s and 7.3 seconds after the rudder angle is reversed. The predicted horizontal motion is sinusoidal, as expected, but there is also some out-of-plane motion caused by the geometrical asymmetry of the submarine due to the sail. The sail is exposed to a local drift angle as the submarine yaws, generating lift. This lift acts above the CG and thus tends to roll the submarine towards the centre of the turn. The coefficient-based simulations also predict that the submarine gradually pitches up to an angle of 2 degrees during the zig-zag. This results from afterbody bound vortic-ity trailing from the sail interacting with the crossflow (the Magnus effect; see for example [41, 42, 43]). A gradual depth decrease is also predicted as a result of the pitch angle.

In terms of the numerical solution, the BDF2 integration scheme was again found to be stable for the zig-zag manoeuvre, even with a very large time step size of 6 s (which is equivalent to the time taken for the submarine to translate approximately 0.5 body lengths at the initial speed for the manoeuvre). As shown in Figure 10, the discretization error follows a similar trend with At as was seen for the rising manoeuvre. The normalized iterative error for the yaw angle is also shown in Figure 10, for the At = 1.0 s case. These results confirm the finding in the rising study that the convergence rate for Newton iteration with J0 and JLp is superior to relaxed fixed-point iteration, even for the perturbed Jacobian matrices. It takes 3-5 correctors for the iterative error to become two orders of magnitude smaller than the temporal discretization error for Newton iteration with perturbed Jacobian matrices while it takes at 8 correctors to reach this level using fixed-point iteration with an optimum relaxation parameter.

As a final note, robust solutions with good convergence were also obtained with the perturbed J0 and JLp models for additional coefficient-based simulations (results not shown) of a hard horizontal turn with a 30 degree rudder deflection and a vertical spiral manoeuvre in which the rudder is deflected to 30 degrees and the sternplanes deflected to -10 degrees.

102p-,-,-,-,-,-,-

: 10°/10° Zig-Zag, - J ■ At = 1 s. J°

101 r J+

N ▼ J0

-e—J,p

> 10° i Q

10-3 r

Fixed pt. ■ 3 = 0.35

" ' - 0.4? „ --0.45 " " 0.5 :

Figure 10: Temporal discretization errors (left) and iterative errors (right) for yaw angle in coefficient-based simulations of a 10° /10° horizontal zig-zag manoeuvre at 6 m/s. Solid lines: Newton iteration; dashed lines: fixed-point iteration.

7. URANS Results and Discussion

The coefficient-based simulations provided a controlled and efficient environment to test the use of approximate Jacobian matrices for 6 DOF simulations of submarine manoeuvres using predictor-corrector integration with Newton iteration. The second phase of this work applies this method to unsteady RANS simulations of the same rising and zig-zag manoeuvres. The unperturbed J0 and JLp models are both tested as approximations to the unknown JF in the URANS simulations. In contrast with coefficient-based simulations, an inner iterative loop with Nf iterations is also needed for converging the URANS equations at of the NC corrector steps.

loop each

7.1. S9 Rising Simulation

The 6 DOF unsteady RANS predictions for submarine state during rising scenario S9 are compared with coefficient-based simulation results in Figure 11. The two methods show a general agreement in submarine motion during the manoeuvre, with both predicting the development of a large roll angle just before

surfacing. However, there are some noticeable differences in the actual quantities due to the differences in modelling assumptions used by the methods. The submarine's motion during the emergency rising manoeuvre has been analysed in more detail elsewhere, with a focus on the roll instability [24, 25, 18] and comparisons between URANS and coefficient-based modelling is discussed in [15] and [12]. For the purposes of the present work, it is important to note that stable URANS solutions were obtained for At ranging from 0.125 s to 2.0 s using the coefficient-based model to estimate JF, despite the modelling differences between URANS and the coefficient-based model seen in Figure 11. As was the case with the coefficient-based model simulations, the error due to temporal discretization in the At = 2 s URANS simulation is considered to be small enough for these results to be valuable for a preliminary analysis. The temporal discretization error becomes too small to be seen in the plots when At is reduced to around 0.5 seconds.

An analysis of the iterative convergence was done by repeating the S9 URANS simulations using different combinations of fluid sub-iterations Nf per corrector and number of correctors NC per time step with At = 2.0 s, for Newton iteration with models J0 and JLp. Figure 12 shows the convergence process during a time step of an unsteady RANS simulation of rising scenario S9 with At = 2 s, Nf = 10, and NC = 4. Only the convergence of roll angle $ and rolling moment K are shown for clarity; similar behaviour is seen for the other components of motion state and hydrodynamic loads. At the start of the first coefficient loop for the time step (cumulative coefficient-loop iteration for the simulation = 1001), the predicted roll angle is extrapolated based on the converged solution for the previous two time levels (t = 50s and t = 48s) using Eqs. (14) and (15). The fluid boundary conditions are set, and held fixed, based on the predicted state for the first 10 coefficient loop iterations. It takes a few coefficient loop iterations for the hydrodynamic loads to converge for these fixed boundary conditions, as can be seen in the convergence of K during iterations 1001-1010 in Figure 12. The converged hydrodynamic loads at the end of the tenth coefficient loop iteration for the time step (cumulative iteration 1010) are

' 1 1 ' 1 ' 1 ' 1 ' / f

1 ....... ■

20 30 40

1 1 // ' jf ' -n i / / / /

OODOUUU«"""

50 60 0

10 20 30 40

10 20 30 40 50

20 30 40

20 30 40

20 30 40

0 10 0 20 40

20 30 40

20 30 40

N? 60 80 100

— At = 0.25 s At = 0.5 s

—- At = 1.0 s

At = 2.0 s (Nf = 5, NC = 4) ° At = 2.0 s (Nf = 5, NC = 7)

— Coefficient-based

0 20 40 60

80 100 120 140 160 180 200 220 240

x„ (m)

Figure 11: Results for URANS simulations of emergency rising manoeuvre S9, compared with simulations using the coefficient-based model. Unless otherwise noted, Nf = 5, Nc = 4, and Newton iteration with model Jo are used for the URANS simulations.

then used to correct the state of the submarine using Newton iteration with JF approximated using coefficient-based model J0. The corrected roll angle is shown by $cl in Figure 12. Once again, the hydrodynamic loads are allowed to converge for 10 coefficient-loop iterations with the fluid boundary conditions fixed for the corrected submarine state. The process continues for a total of 4 correctors. It can be seen that the incremental changes to K and $ become relatively small for the last two correctors; the difference between and is too small to be seen on the scale used for the plot. The residuals for the submarine's position and Euler angles for the same time step are shown Figure 13. Here Ay' is the incremental change in the submarine state during a corrector, divided by the step change in submarine state during the time step (y®+i — y®). The maximum change in any component of submarine state drops to less than 0.5% of the step change for the time step by the 4th corrector.

eg ■6 -20

Figure

I.........I.........I

S9, t = 52 s, At = 2 s, Nf = 10, NC = 4

tnd = iC

-35 1000

l l l l l l l I l l l l l l l l l I l l l l l l l l l I l l l l l l l l l I

1010 1020 1030 Cummulative fluid iterations

e 12: Roll angle and rolling moment on the submarine at the end of each coefficient-loop iteration for the final time step of a URANS simulation of rising scenario S9 with At = 2 s and Newton iteration with Jacobian matrix model Jo. The state of the submarine is corrected after every 10th coefficient-loop iteration.

It can be seen in Figure 12 that the hydrodynamic loads do not change much after the 5th fluid sub-iteration. Considering that the Jacobian matrix is only approximate, it is more efficient to reduce Nf and tolerate a looser convergence

of hydrodynamic loads when correcting the submarine state. This allows more correctors to be performed for a given number of total fluid coefficient-loop iterations per time step, which is an approximate measure of computational cost because the cost of computing the Jacobian matrix and correcting th submarine state is negligible. The simulation was repeated using Nf = 5 and NC = 7. As shown in Figure 13, the normalized residuals for this case were approximately an order of magnitude smaller at the end of the time step than the simulation with Nf = 10 and NC = 4, even though the total number of coefficient-loop iterations was reduced to 35 from 40.

Corrector

2 3 4 5

Corrector

Figure 13: Iterative convergence of the submarine position and Euler angles at time level t = 52 s, for URANS rising simulations S9 with Nf = 10, NC =4 (left) and Nf = 5,NC =7 (right). The residuals are normalized by the change in submarine states over the interval

= 50s <t< ti+1 = 52s.

Additional simulations were performed using Nf = 3, Nf = 4, and Nf = 5 with several values of NC to determine if there is an optimum Nf for iterative performance. For this study, the results for [At = 2.0 s , Nf = 5 , NC = 7] was considered to be the "fully converged" solution, yc, against which all others are compared for estimating cumulative iterative error, Ej, using Eq. (31). For each of the 12 components, k, of submarine motion state, the norm of the

iterative errors k for the simulation are evaluated and normalized by the norm of the discretization errors, ED,k. The temporal discretization errors are estimated by evaluating the difference between the fully converged At = 2.0 s solution and the At = 0.25 s solution. The normalized iterative errors for th 12 components of submarine state are then averaged to determine the average level of iterative error relative to temporal discretization error:

\\EI\L / \\Ed\L = 12 £ \\Ei,k\L / \En,k\L (32)

The average normalized iterative error is shown for all of the At = 2 s results in Figure 14. The total number of fluid coefficient loop iterations per timestep, Nf x NC, is used for the x-axis so that results are compared based on computational effort. First, consider the results for Newton iteration with model J0 for different numbers of fluid sub-iterations (solid lines). Simulations with Nf = 3 and Nf = 5 perform equally well for 15 < Nf x NC < 20 while simulations with Nf = 4 converge slightly faster; Nf = 4 appears to be the optimum. Note that the best number to use for Nf depends on the combination of fluid solver, mesh, time step size, and accuracy of the Jacobian matrix model.

The minimum number of coefficient-loop iterations needed to drop the iterative error to a level that is at least an order of magnitude smaller than the discretization error is 16, which was obtained for Nf = 4 and NC = 4. When a total of 20 - 21 coefficient-loop iterations are used per time step, the iterative error is less than 5% of the temporal discretization error for each case (Nf = 3, 4, 5). This is a relatively small error in terms of the range of values during the simulation, as can be seen by comparing the [Nf = 5, NC = 4] results with the [Nf = 5, NC = 7] results in Figure 11.

Next, notice that simplified model JLp (dashed lines in Figure 14) converges at nearly the same rate as J0. This indicates that JLp provides a good approximation to the "true" URANS Jacobian and that there is little advantage to adding the complexity of a non-linear model. This is an important result because all of the coefficients in JLp are relatively easy to estimate, compute with CFD, or measure in experiments. These coefficients can readily be estimated

i i i i I i i i i I i i i i I i i i i I i i i i I i i i i I i i i i I i i i i I |J

Nf = 3 ■Nf = 4

Nf = 5 Nf = 10

- — "Nf = 4 JLp

- + -N = 5 JLp

ONf = 5 , Y = 0.3

1 11 11 1 11 11 1 11 11 1 11 11 1

10 15 20 25

Figure 14: Iterative error, normalized by discretizatio: and various combinations of Nf and Nc .

y discretization C ^^

tion error,

for S9 simulations with At = 2 s

with negligible computational effort for another submarine geometry using a program such as DSSP50 [44]. Finally, notice the poor convergence of a simulation using fixed-point iteration with a relaxation parameter of 0.3, which was found to be near-optimum for the At = 2 s coefficient-based simulations. Even with 5 fluid sub-iteration and 4 correctors, the iterative errors are 7 times as large as the discretization error. This confirms the finding in the coefficient-based study that under-relaxed fixed-point iteration performs poorly relative to Newton iteration, especially for large At.

7.2. Horizontal Zig-Zag Results

Unsteady RANS simulations of the horizontal zig-zag manoeuvre were conducted using Nf = 4 and NC =4 for At = 1s, 0.5 s, 0.25 s, and 0.125 s. The typical convergence of submarine state during a time step for Newton iteration with models J0 and JLp is shown in Figure 15. The residuals for submarine state drop with corrector iteration in a consistent manner for all time step sizes, being

reduced to the order of 0.5% of (yj+i — yi) by the end of the time step. This indicates that the coefficient based model provides an adequate estimate for the Jacobian matrix for this manoeuvre as well. Very little difference in conver-

gence rate was seen between simulation with J0 and JLp for both timestep sizes, confirming that model JLp contains the important hydrodynamic coefficie:

< 10-3

£ 10-4

10-5 r.

At = 1 s, t = 15 s.

x0~ -O- -

V -o- -

Corrector

Corrector

Figure 15: Residuals versus corrector iteration at t = 15 s for the At = 1 s and At = 0.25 s URANS simulations of the 10°/10° zig-zag manoeuvre, for Newton iteration with Jacobian models Jo and JLp. Residuals are normalized by the step change in values for the time step,

yi+1 - yi.

The URANS predictions for submarine state during the zig-zag manoeuvre are compared with coefficient-based results in Figure 16. There is excellent agreement between the coefficient-based simulations and CFD for predicted yaw le and lateral displacement of the submarine. In contrast, there is substantial discrepancy in the smaller out-of-plane motion; the roll angle, pitch angle, and vertical displacement. Additional analysis has shown that a primary cause of the discrepancy in roll angle is the quasi-steady assumption made in the coefficient-based model. This assumption neglects unsteady effects due to the time lag in vorticity propagation from the sail to the tail, as described in [15]. The time scale for the submarine motion is not long enough relative to the time

t (s) t (s) t (s)

t (s) t (s) t (s)

Figure 16: Results for URANS simulations of a 10°/10° zig-zag manoeuvre, compared with simulations using the coefficient-based model.

for vorticity to propagate from sail to tail (approximately 6 seconds) for the quasi-steady assumption to be valid. For example, during the first 11s when the crossflow is increasing in magnitude, and the lift on the sail (and therefore rolling moment on the submarine) is increasing in proportion to v, the vorticity trailing from the sail has not yet convected aft to the tail where it eventually interacts with the tailplanes to reduce the rolling moment induced at the sail. The quasi-steady coefficient-based model assumes steady state conditions occur instantaneously and does not model the time delay (about 6 s) before the initial rolling moment gets moderated. Therefore it under-predicts the initial rolling moment on the submarine. A simplified analysis estimates that at least 50% of the discrepancy in roll angle between the two methods is due to this unsteady effect, which is accounted for in the URANS solution. Since the pitch rate is a function of the roll angle according to Eq. (3e), the pitch angle and consequently the vertical displacement are also affected by the discrepancy in roll angle.

These discrepancies show the value of comparing coefficient-based simulations with the more comprehensive URANS simulations. They also confirm the point made earlier that only an approximate model is required to estimate JF in the predictor-corrector method. A separate study [12] also shows that stable convergence can be achieved even when the submarine is subject to significant disturbance forces from a passing tanker, which are completely neglected in the model for the Jacobian matrix.

8. Conclusions

CAn implicit predictor-corrector method was developed for integration of 6 DOF URANS simulations of submarine manoeuvres and applied to simulations of emergency rising and a horizontal zig-zag manoeuvres. Two iterative schemes were compared for solving the implicit system each timestep: under-relaxed fixed-point iteration and Newton iteration. It was found that fixed-point iteration suffers from stability constraints on the relaxation parameter that become more restrictive and detrimental to efficiency as timestep size is increased within

a practical range. With suitable approximations to derivatives of hydrodynamic loads, Newton iteration was found to provide significantly better convergence for all At. A simplified linearised coefficient model (JLp) for the hydrodynamic component of the Jacobian matrix gives comparable convergence efficiency with th full non-linear model (Jo), when considering the uncertainty in the coefficients. The JLp model is attractive because all of the coefficients contained within it are relatively easy to estimate semi-empirically, compute using CFD, or measure experimentally. It was found to be stable for a wide range of time step sizes for different submerged submarine manoeuvres, allowing the time step size to be selected based on the desired accuracy rather than stability constraints. An accounting of added mass effects in the Jacobian matrix is essential for all time step sizes, due to the relatively large ratio of added mass to body mass in submarine manoeuvring simulations. The primary linear hydrodynamic derivatives for steady translation and rotation must also be accounted for when timestep sizes are increased towards the higher end of the practical range. The simplest model for the Jacobian matrix, containing only the primary added mass effects, was found to be adequate only for small time steps; it becomes less efficient and unstable as the time step size is increased.

In summary, it is concluded that implicit predictor-corrector time integration using Newton iteration with a simplified linearised Jacobian model for hydro-dynamic loads, JLp, is a good option for integrating the 6 DOF rigid-body equations in unsteady CFD-based simulations of typical full-submerged 6 DOF submarine manoeuvres. Additional work is needed to assess the requirements for the Jacobian model for other scenarios such as submarine motions on the surface in waves. The proposed implicit predictor-corrector method is expected to be useful for other 6 DOF URANS simulations where the added mass/inertia in any direction is comparable to the vehicle mass/inertia (such as surface ships), though some modification may be required to the Jacobian matrix model to suit the problem.

cients.

Acknowledgements

Financial support for this research was provided by Defence Research and Development Canada - Atlantic and by the National Sciences and Engineering Research Council of Canada (NSERC) through an Alexander Graham Bell Canada Graduate Scholarship-Doctoral (CGS D). Some simulations were performed using computational facilities provided by ACEnet, the regional high performance computing consortium for universities in Atlantic Canada. ACEnet is funded by the Canada Foundation for Innovation (CFI), the Atlantic Canada Opportunities Agency (ACOA), and the provinces of Newfoundland and Labrador Nova Scotia, and New Brunswick.

References References

[1] M. Gertler and G.R. Hagen. rd equations of motion for submarine simulation. Technical Report I aval Ship Research and Development Center, June 1967.

[2] J. Feld >RDC revised standard submarine equations of motion. Techni 3TNSRDC/SPD-0393-09, David W. Taylor Naval Ship

Resear dopment Center, June 1979.

oger, F. Davoudzadeh, J.J. Dreyer, H. McDonald, C.G. Schott, W.C. ierke, A. Arabshahi, W.R. Briley, J.A. Busby, J.P. Chen, M.Y. Jiang, . Jonnalagadda, J. McGinley, R. Pankajakshan, C. Sheng, M.L. Stokes, .K. Taylor, and D.S. Whitfield. A physics-based means of computing the flow around a maneuvering underwater vehicle. Technical Report TR 97002, Applied Research Laboratory, Pennsylvania State University, January 1997.

[4] M. Bettle, S.L. Toxopeus, and A. Gerber. Calculation of bottom clearance effects on Walrus submarine hydrodynamics. International Shipbuilding Progress, 57(3-4):101-125, 2010.

[5] G.D. Watt, C.R. Baker, and A.G. Gerber. ANSYS CFX-10 RANS normal force predictions for the Series 58 Model 4621 unappended axisymmetric submarine hull in translation. DRDC Atlantic TM 2006-037, September 2006.

[6] T. L. Jeans, G. D. Watt, A. G. Gerber, A. G. L. Holloway, and

Baker. High-resolution Reynolds-averaged Navier-Stokes flow predictions over axisymmetric bodies with tapered tails. AIAA Journal, 47(1):19-32 2009.

w predi

ce of tu

[7] A.B. Phillips, S.R. Turnock, and M. Furlong. Influence of turbulence closure models on the vortical flow field around a submarine body undergoing steady drift. Journal of Marine Science and Technology, 15(3):201-217, 2010.

and T. Hu;

[8] C.H. Sung, T.-C. Fu, M. Griffin, and T. Huang. Validation of incompressible flow computation of forces and moments on axisymmetric bodies undergoing constant radius turning. In Twenty-First Symposium on Naval Hydrodynamics, Washington, DC, 1997.

[9] S. L. Toxopeus, P. Atsavapranee, E. Wolf, S. Daum, R. Pattenden, R. Wid-jaja, J.T. Zhang, and A.G. Gerber. Collaborative CFD exercise for a submarine in a steady turn. In ASME 31st International Conference on Ocean, Offshore and Artic Engineering, number OMAE2012-83573, Rio de Janeiro, Brazil, July 2012.

10] J.A. Maxwell, J.T. Zhang, G.L. Holloway, G.D. Watt, and A.G. Gerber. Simulation and modelling of the flow over axisymmetric submarine hulls in steady turning. In 28th Symposium on Naval Hydrodynamics, Pasadena, California, September 2010.

[11] S.L. Toxopeus. Practical application of viscous-flow calculations for the simulation of manoeuvring ships. PhD thesis, Delft University of Technology, 2011.

[12] M. Bettle. Unsteady Computational Fluid Dynamics Simulations of Six Degrees-of-Freedom Submarine Manoeuvres. PhD thesis, University of New

Brunswick, 2013.

[13] R. Pankajakshan, M.G. Remotigue, L.K. Taylor, M. Jiang, W.R. Briley, and D.L. Whitfield. Validation of control-surface inducted submarine maneuvering simulations using UNCLE. In Twenty-Fourth Symposium on Naval Hydrodynamics, Fukuoka, Japan, July 2002.

[14] G. Venkatesan and W. Clark. Submarine maneuvering simulations of ONR Body 1. In Proceedings of the International Conference on Offshore Mechanics and Arctic Engineering - OMAE, volume 4, pages 697-705, San Diego, CA, 2007.

[15] M.C. Bettle, A.G. Gerber, and G.D. Watt. Unsteady analysis of the six DOF motion of a buoyantly rising submarine. Computers and Fluids, 38(9):1833-1849, 2009.

[16] J.J. Dreyer and D.A. Boger. Validation of a free-swimming, guided multi-body URANS simulation tool. In 28th Symposium on Naval Hydrodynamics, Pasadena, California, September 2010.

[17] P.M. Carrica, H. Sadat-Hosseini, and F. Stern. CFD analysis of broaching for a model surface combatant with explicit simulation of moving rudders and rotating propellers. Computers and Fluids, 53(1):117-132, 2012.

[18] G.D. Watt. Modelling and simulating unsteady six degrees-of-freedom submarine rising maneuvers. DRDC Atlantic TR 2007-008, February 2007.

[19] Uri M. Asher and Linda R. Petzold. Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. Society for Industrial and Applied Mathematics, Philadelphia, PA, 1998.

[20] K.E. Brenan, S.L. Campbell, and L.R. Petzold. Numerical Solutions of Initial-Value Problem in Differential-Algebraic Equations. Society for Industrial and Applied Mathematics, Philadelphia, PA, 1996.

[21] J. J. Dryer, L. K. Taylor, W. C. Zierke, and F Davoudzadeh. A first-principle approach to the numerical prediction of the maneuvering characteristics of submerged bodies. In FEDSM97-3130, 1997 ASME Fluids Engineering Division Summer Meeting, Vancouver, British Columbia, Canada, June 22-26 1997.

[22] R.V. Wilson, P.M. Carrica, and F. Stern. Unsteady RANS method for ship motion with application to roll for a surface combatant. Computers and

for ship

Fluids, 35:501-524, 2006.

Noack, and

[23] Pablo M. Carrica, Robert V. Wilson, Ralph W. Noack, and F. Stern. Ship motion using single-phase level set with dynamic overset grids. Computers and Fluids, 36:1415-1433, 2007.

[24] G.D. Watt. A quasi-steady evaluation of submarine rising stability: The stability limit. In RTO-AVT Symposium on Advanced Flow Management, June 2001.

[25] G.D. Watt and H.J. Bohlmann. Submarine rising stability: Quasi-steady theory and unsteady effects. In 25th Symposium on Naval Hydrodynamics, St. John's Newfoundland and Labrador, CANADA, August 2004.

[26] M. Mackay. The standard submarine model: A survey of static hydro-dynamic experiments and semiempirical predictions. DRDC Atlantic TR 2003-079, Defence R&D Canada - Atlantic, June 2003.

27] G.D. Watt, B. Tanguay, and K.R. Cooper. Submarine hydrodynamics in the wind tunnel: The DREA Static Test Rig. In RINA Warship 91 International Symposium on Naval Submarines 3, London, May 1991.

[28] G. D. Watt. Re-evaluating the thrust and torque measurement capability of the Mark I propulsion system for the DREA Static Test Rig. DRDC Atlantic TM 95/209, April 1995.

[29] G.D. Watt and E.Y. Fournier. Submarine propulsion testing in the IAR 9 m wind tunnel. In Third Canadian Marine Hydrodynamics and Structures Conference, Halifax, August 1995.

[30] G.D. Watt. Estimates for the added mass of a multi-component, deepl submerged vehicle. DRDC Atlantic TM 88/213, October 1988.

deeply

[31] Albert Strumpf. Equations of motion of a submerged body with va, ying mass. Technical Report SIT-DL-60-9-771, Stevens Institute of Technology,

with vary tute of Technok

_______odeling of

Liild 061102

and Develc

[32] M. Mackay. Semiempirical component based modeling of submarine hydrodynamics and systems: the DSSP21 (build 061102) companion. DRDC Atlantic TR 2007-039, Defence Research and Development Canada - Atlantic, March 2007.

[33] ANSYS Canada Ltd. ANSYS CFX-Solver, Release 12.1: Modeling. Software documentation.

[34] F.R. Menter. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal, 32(8), August 1994.

[35] J.D. Lambert. Numerical Methods for Ordinary Differential Systems: The Initial Value Problem. John Wiley & Sons Ltd., Chichester, 1991.

[36] B.R. Hutchinson, P.F. Galpin, and G.D. Raithby. Application of additive correction multi-grid to the coupled fluid flow equations. Numerical Heat Transfer, 13(2):133-147, 1988.

[37] M.J. Raw. A coupled algebraic multi-grid method for the 3D Navier-Stokes equations. In Proceedings of the Tenth GAMM-Seminar. Notes on Numerical Fluid Mechanics, Keil, Germany, 1994.

[38] B.R. Hutchinson and G.D. Raithby. A multi-grid method based on the additive correction strategy. Numerical Heat Transfer, 9:511-537, 1986.

[39] M.J. Raw. Robustness of coupled algebraic multi-grid for the Navier-Stokes equation. In 34th Aerospace and Sciences Meeting & Exhibit, number AIAA 96-0297, Reno, NV., January 1996.

[40] L. Eca and M. Hoekstra. Evaluation of numerical error estimation based on grid refinement studies with the method of the manufactured solutions. Computers & Fluids, 38:1580-1591, 2009.

[41] M. Mackay and J.T. Conway. Modelling the crossflow body separation on a submarine using a panel method. In International Symposium on Naval Submarines 3, volume 2, London, May 1991. The Royal Institution of Naval Architects.

[42] G.D. Watt, A.G. Gerber, and A.G.L. Holloway. Submarine hydrodynamics studies using computational fluid dynamics. In 8th Canadian Marine Hydrodynamics and Structures Conference, St. John's, October 2007.

[43] D. H. Bridges, J. N. Blanton, W. H. Brewer, and J. T. Park. Experimental investigation of the flow past a submarine at angle of drift. AIAA Journal, 41(1):71-81, 2003.

[44] M. Mackay. DSSP50 Build 090910 Documentation: Part 3: Algorithm Description. DRDC Atlantic TM 2009-228, DRDC Atlantic, Dec 2009.

kay. DS scription. DR

A predictor-corrector method is developed for CFD-based 6-DOF manoeuvring simulations Simplified hydrodynamic models are used to estimate the Jacobian matrix. Submarine zig-zag and emergency rising manoeuvres were simulated with the method. The method was found to be robust and efficient compared with fixed-point iteration. A simple linearized hydrodynamic model was sufficient for the Jacobian matrix.