The approach of this physics-based research aimed at maneuvering predictions of self-propelled underwater vehicles is to solve the three-dimensional time-dependent incompressible turbulent Navier-Stokes equations. Because of the magnitude of the problem size with regard to resolving full configurations that include sail, sail planes (or bow planes), stern appendages, and rotating propulsor, the numerical solution of the Navier-Stokes equations is carried out on dynamic relative motion multiblock structured grids. The Reynolds numbers of these flows are extremely large and the viscous regions are resolved to y+ values of near one. This places severe demands on the numerical solution scheme in terms of stability and accuracy. Moreover, the flow is both three dimensional and unsteady in the sense of the Reynolds averaged mean flow. To make the computation of these high Reynolds number unsteady flow problems practical in terms of the total CPU time required, the time step needs to be restricted by the physics of the problem being solved and not the numerics of the scheme used to solve the equations. The equations, numerical flux formulation, and the solution algorithm used to solve the equations and achieve the physics-based time restriction is discussed in this section.
The three-dimensional time-dependent Navier-Stokes equations are first transformed to a time-dependent curvilinear coordinate system. The artificial compressibility idea of Chorin  is then introduced [2, 3, and 4]. The use of artificial compressibility permits the experience gained in the numerical solution of compressible flow problems to be exploited in the numerical solution of incompressible flow problems . The artificial compressibility form of the three-dimensional time-dependent Navier-Stokes equation in general curvilinear coordinates is
K=F, θk=U for k=ξ
K=G, θk=V for k=η
K=H, θk=W for k=ζ
In these equations, β is the artificial compressibility coefficient, with a typical value of 5~10; p is static pressure; u, v, and w are the velocity components in Cartesian coordinates x, y, and z. U, V, and W are the contravariant velocity components in curvilinear coordinate directions ξ, η, and ζ, respectively. Terms where k=ξ, η, and ζ, are the viscous flux components in curvilinear coordinates. J is the Jacobian of the inverse transformation, and kx, ky, kz, and kt with k=ξ, η, and ζ are the transformation metric quantities , where a subscript denotes differentiation.
In this work, the thin-layer approximation is introduced to simplify the full Navier-Stokes equations, an algebraic turbulence model [ 6], a k–ε model and a nonlinear k–ε model  were implemented within the code and used for the turbulent flow computations. The details of treating the viscous terms is explained in the work of Gatlin . In addition, improvements have been made by Chen  with regard to the computation of the wall shear stress by improving the computation of the tangential velocity derivatives normal to a solid surface. This improvement is simple to implement and works extremely well on grids that may be highly skewed .
Equation (1) is discretized into a cell-centered finite-volume form which for one-dimensional flow, for example, can be written as
where the index i corresponds to a cell center and indices i±1/2 correspond to cell faces. In this expression, the dependent variable vector Q is considered to be constant throughout the cell whereas the flux is assumed to be uniform over each surface of the cell. A flux vector is therefore needed at each cell face.
There are numerous ways of developing this flux vector, and the formulation used early on  in this maneuvering underwater vehicle research was the flux difference split scheme of Roe  for the first-order contribution and a hybrid numerical flux vector for the higher-order contribution that was patterned after the flux vector developed for compressible flow . An advantage of this hybrid flux is that the formulation leads more or less naturally to the limiting of characteristic variables which is important for compressible flows