On the application of the single-phase level set method to naval hydrodynamic flows

Computers & Fluids 36 (2007) 868–886 www.elsevier.com/locate/compfluid

On the application of the single-phase level set method to naval hydrodynamic flows A. Di Mascio, R. Broglia *, R. Muscari INSEAN – Istituto Nazionale Studi ed Esperienze di Architettura Navale, Via di Vallerano 139, 00128 Rome, Italy Received 20 December 2005; received in revised form 1 August 2006; accepted 13 August 2006 Available online 17 November 2006

Abstract The application of the single-phase level set approach to the numerical simulations of three-dimensional free surface flows around complex geometries, at both non-breaking and breaking regimes is presented. In this approach only the liquid phase is simulated and the level set function is used as tracking device to locate the free surface position. The extrapolation of the solution in the dummy points in the gaseous phase is such that second-order accuracy is maintained also in the points adjacent to the free surface; the time evolution of the level set function and the re-initialization step have been merged so to get a function which is a distance function everywhere, and satisfies, at the same time, the kinematic condition on the free surface. The implementation of this technique into a general purpose Reynolds averaged Navier–Stokes (RANS) equations solver developed at INSEAN [Di Mascio A, Broglia R, Favini B. A Second Order Godunov-type Scheme for Naval Hydrodynamics. Kluwer Academic/Plenum Publishers; 2001, p. 253–61], is described in details; capabilities of the algorithm in dealing with non-breaking and breaking flows in the naval hydrodynamic context will be demonstrated by using a submerged hydrofoil and two different ship hulls in straight course as test cases. Comparisons with both experimental data and numerical surface fitting computations are presented; convergence properties of the algorithm, as well as validation and verification assessment will be also discussed.  2006 Elsevier Ltd. All rights reserved.

1. Introduction The numerical simulation of free surface turbulent flows is a goal of paramount importance for the design and the optimization of ship hulls. Peculiar to this field of research is the presence of a material interface, whose location is one of the unknowns of the problem; as a matter of fact, the simulation of this interface is still a numerical research issue. The results of many efforts dedicated to this difficult topic are summarized, for instance, in the proceedings of the various workshops on CFD in ship hydrodynamics; in the workshop held in Gothenburg [2] in 2000, it can be seen that most of the computations were performed by means of surface fitting approaches, i.e. with algorithms where the grid moves to fit the actual position of the free *

Corresponding author. Tel.: +39 06 50299297; fax: +39 06 5070619. E-mail address: [email protected] (R. Broglia).

0045-7930/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2006.08.001

surface. Although excellent results can be obtained by means of these techniques (see for example [3,4]), they turn out to be quite complicated when dealing with practical cases, because of the geometrical shape of the domain boundary (hull), and/or as a consequence of the cell aspect ratio; moreover, due to the need to regenerate the mesh at each time step, they could also be computationally expensive. Furthermore, when dealing with unsteady problems of practical interest, like a ship maneuvering in waves, the deformation of the free boundary can be so large that it would require either a change in grid topology (for breaking waves, for instance) or strong cell deformation. In both cases, it would be almost impossible to handle these problems by a general simulation algorithm. A remedy to these problems is the use of methods in which the grid is fixed, and the position of the free surface is defined by either massless particles (front-tracking methods) or through some kind of scalar function (front-

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

capturing methods). One of the earliest attempt to use massless particles was the Marker-and-Cell method proposed by Harlow and Welch [5]; in this method, massless particles, which move with the underlying flow, are distributed all over the computational domain and are used to identify each fluid or each phase of a single fluid; the interface is reconstructed by following these particles. More recently, Unverdi and Tryggvason [6] have proposed a front-tracking method where massless particles, distributed over the interface, explicitly track the surface of discontinuity, and the multi-phase fluid is modelled as a single fluid with variable properties (density and molecular viscosity). Front-tracking methods have proved to be both accurate and robust; however, in order to prevent the clustering or the rarefaction of marker particles, regridding algorithms must be employed; moreover, special care is required when dealing with an interface that changes in topology, as, for instance, in the case of a free surface that breaks and/or merges [6,7]. In front-capturing approaches, such as Volume-of-Fluid [8] (VOF) or Level Set [9,10] (LS) methods, no special procedure is required to model interfaces that change in topology. In the former, a volume fraction function, convected by the underlying fluid flow, is defined in the whole domain; this function is set to unity everywhere in one fluid (or phase) and zero in the other. In those regions where both fluids or phases are present, the value of this function represents the fraction of the volume occupied by the first fluid or phase; density and viscosity are then defined as a smooth function of the volume fraction. The major advantages of this technique lie in its capability in treating complex interface evolutions, even breaking and reconnecting surfaces, and in its mass conservation properties; the major disadvantages of the VOF method lies in the difficulty to compute the geometric properties of the interface as, for instance, the curvature or the normal vector [11], which are needed, for example, when surface tension effects are taken into account. For VOF methods, the interface reconstruction is still a research topic and moreover modelling errors stem from the difficulties to control the interface thickness, because of artificial spreading due to numerical diffusion. In the LS approach, the interface is simulated by the evolution of the zero level of a smooth function convected by the underlying flow. This approach allows the simulation of complex interface evolutions including breaking or merging phenomena; moreover, the computation of geometrical parameters of the interface is straightforward. Nevertheless, there are still some issues to be addressed. In fact, in classical LS approaches, both air and water are simulated as a single fluid whose properties vary continuously across the free surface; in order to get a stable algorithm, the change of fluid properties must be spread over several grid cells. This can be a source of significant loss of accuracy, especially in three-dimensional computations, because of the limited amount of cells at disposal in the free surface region. Moreover, it has been shown by Sussman et al. in [10] that, even if a conservative algorithm is

869

employed for the evolution of the level set function, this does not imply conservation of volume of each fluid (or phase); this problem can be partially cured by keeping the level set function as the signed distance from the interface. This procedure, known as re-initialization step has been widely investigated in the literature [12–15] and it has to be used to avoid excessive deformation of the isosurfaces of the level set function around the zero level and to ensure constant thickness of the transition region during the evolution of the interface. However, it has to be taken into account that the re-initialization step can be computationally expensive. It is worth to mention that, in order to get the best from both VOF (mass conservation properties) and LS (easy computation of geometric properties of the interface) some authors have developed combined LS/VOF methods (see, for instance, the works done by Sussman et al. [10], Sussman and Puckett [16], Sussman and Dommermuth [17] and Sussman [18]). In spite of the above-mentioned problems, LS technique has been successful applied to a wide variety of problems, ranging from fluid mechanics, dendritic growth and solidification problems, combustion and image processing (for an accurate review of the method and an overview of the latest applications, see [19,20]), as well as in the computation of the flow around surface piercing ship hulls (see [21–23], among others); these approaches have also been applied for the simulation of complex three-dimensional flows with breaking waves [24,25] and for the simulation of flows around floating bodies with large amplitude motion [26]. Moreover, the increase in the use of front-capturing methods such as LS or VOF approach in the field of naval hydrodynamics is demonstrated by the number of authors that have used these methodologies to perform the steady and unsteady test cases proposed by the last workshop on CFD in ship hydrodynamics held in Tokyo in 2005 [27]. Numerical solutions obtained with a LS approach for the steady free surface flow around practical ship hulls have been presented by Broglia et al. [28], Wilson et al. [29], Miller et al. [30], Kim et al. [31], Tahara et al. [32] and Hino and Sato [33]; whereas Rhee and Skinner [34], HSU et al. [35], Deng et al. [36] and Chao [37] used a VOF technique. Front-capturing methods have been also used to perform the unsteady diffraction around the naval combatant DDG51; Carrica et al. [38] and Cura–Hachbaum et al. [39] have made use of a single-phase and two-phase LS approach respectively, whereas Deng et al. [36] presented results obtained by using a VOF technique, demonstrating the capabilities of these approaches in dealing with seakeeping and maneuvering problems. Nevertheless, the need to spread the transition of the fluid properties on several grid cells (required for numerical stability), can be a source of inaccuracy when simulating the flows around complex geometries in three spatial dimensions. In order to maintain the interface sharp and to improve mass conservation properties, the so called Ghost Fluid Method [40–42], or hybrid LS/Marker particles methods (see [43,44]), or a single-phase method in

870

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

the framework of both VOF [45] and LS approaches [46– 48,23], can be used. In the single-phase LS approach, the solution is computed only in the liquid phase, and the level set function is used as a tracking device to locate the intersection of the free surface with the lines of the fixed underlying grid. Moreover, a single fluid being computed, there is no transition zone across which the fluid properties vary, and consequently no uncertainty is related to interface smoothing. On the other side, the solution in the gaseous phase being computed only for numerical reasons, the single-phase LS approach is not suitable neither for problems with air entrapment phenomena or bubbles, nor for problems where stresses on the free surface due to air play an important role, as in wind–induced waves problems. Some interesting applications of single-phase LS can be found in [48] for the simulation of the 3D free surface flow around simple geometries, in [47] for 2D free surface flow with breaking around submerged profiles and in [46] for the simulation of 3D steady flows around ship hulls without breaking waves. The computation of 3D steady flows with breaking waves around a wedge and at the bow of a practical ship hull appeared in [24,25]. Carrica et al. [49] applied unsteady single-phase LS to the study of the forward speed diffraction problem around a combatant ship. The single-phase LS algorithm adopted here is peculiar because the level set function satisfies the free surface kinematic boundary condition only on the free surface itself, whereas it can be shown that it turns out to be a distance function at any other point, which enhances accuracy. Moreover, the formulation is such that it can be easily implemented in conjunction with any RANS solver developed for the simulation of constant density fluid flows. The formulation described in the following is applicable in the present form only to steady state problems, because no attempt was done to address the issue of solution prolongation in those points where the state changes from ‘‘air’’ to ‘‘water’’. This topic is fundamental for the extension of the single-phase LS to unsteady simulations. An original solution to this problem is proposed by Carrica et al. in [50]. The objectives of this work are: (i) to describe in details the single-phase LS approach as it is applied to a general purpose in-house solver for incompressible turbulent flow; (ii) the verification of the method and the validation of the numerical results obtained with this approach in both two- and three-dimensional test cases; (iii) to assess the capabilities of this approach in dealing with complex three-dimensional free surface flows in the presence of strong breaking waves. The paper is organized as follows: the physical problem is defined first, in terms of conservation equations governing the flow of an incompressible viscous fluid. In the following section, the numerical algorithm used for the simulation of the bulk flow and the free surface is summarized and the single-phase LS approach for the treatment of the free surface will be described in details. In Section 4 results obtained with the single-phase LS method will be presented for the two-

dimensional flow past a submerged hydrofoil and the three-dimensional flow around a Series-60 ship hull; both cases will be considered as tests for the verification and validation of the methodology for non-breaking free surface flows. Finally, as an example of practical application, the flow around a naval combatant ship hull (characterized by a bulbous bow and a transom stern), in both non-breaking and breaking flow regimes, will be presented; verification and validation in breaking regime will be also carried out. Wherever possible, the numerical solutions provided will be compared with available experimental and other numerical data. In the last section conclusions will be drawn. 2. Governing equations The turbulent motion of an incompressible (constant density) viscous fluid can be described by the Reynolds averaged Navier–Stokes (RANS) equations, that read, in integral form I ui ni dS ¼ 0; SðVÞ Z I o ui dV þ ½ui uj nj þ pni  sij nj  dS ¼ 0 i ¼ 1; 2; 3; ot V SðVÞ ð1Þ where V is a control volume, SðVÞ its boundary, and ~ n the outward unit normal. A reference length l and velocity U1 have been chosen to make the previous equations nondimensional; ui is the ith Cartesian component of the velocity vector (in the following, the Cartesian components of the velocity will be also denoted with u, v, and w); p is a new variable related to the pressure P and the acceleration of gravity g (parallel pffiffiffiffiffi to the vertical axis z) by p = P + z/Fr2, Fr ¼ U 1 = gl being the Froude number. Finally, sij = mt(ui,j + uj,i) is the stress tensor, mt = 1/Re + mT is the global kinematic viscosity, with Re = U1l/m the Reynolds number, m the kinematic viscosity and mT the turbulent viscosity. In the present work, the turbulent viscosity has always been calculated by means of the Spalart and Allmaras one-equation model [51]. The problem is closed by enforcing appropriate conditions at physical and computational boundaries. On solid walls, velocity is set to zero (whereas no condition on the pressure is required); at the (fictitious) inflow boundary, velocity is set to the undisturbed flow value, and the pressure is extrapolated from inside; on the contrary, the pressure is set to zero at the outflow, whereas velocity is extrapolated from inner points. At the free surface, whose location is one of the unknowns of the problem, the dynamic boundary condition requires continuity of stresses across the surface; if the presence of the air is neglected, the dynamic boundary condition reads: j z p ¼ sij ni nj þ þ ; sij ni t1j ¼ 0; sij ni t2j ¼ 0; ð2Þ We Fr2

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

where We ¼ qU 21 l=r is the Weber number (q being the density of the fluid and r the surface tension coefficient) and j the surface curvature; ~ n; ~ t1 and ~ t2 are the surface normal and two tangential unit vectors, respectively. The kinematic boundary condition states that the free surface F(x, y, z, t) = 0 is a material surface, and allows to determine its unknown configuration

2

ð3Þ

Initial conditions have to be specified for the velocity field and for the free surface configuration F ðx; y; z; 0Þ ¼ F ðx; y; zÞ:

ð4Þ

3. Numerical model The major difficulty in solving the incompressible RANS equations lies in the solenoidal constraint on the velocity vector field. When only the average steady state has to be computed, this system of equations can be conveniently replaced by the pseudo-compressible formulation [52], that reads, in integral form Z Z o ~ ð~ F cs  ~ F ds Þ dS ¼ 0; ð5Þ q dV þ ot V SðVÞ where ~ q ¼ ðp; u; v; wÞT is the state variables vector for pseudo-compressible flows and ~ F cs ¼ ðbul nl ; u1 ul nl þ pn1 ; u2 ul nl þ pn2 ; u3 ul nl þ pn3 ÞT ; ~ F ds ¼ ð0; s1l nl ; s2l nl ; s3l nl ÞT ð6Þ are the convective (inviscid and pressure) and diffusive normal fluxes through the surface S of the finite volume V, and b is the pseudo-compressibility factor. This system of equations (that include an evolution equation for the pressure) is hyperbolic in its Eulerian part; therefore, all the numerical methods developed for compressible flow simulation can be adapted and applied to the above system. In the scheme adopted here, these equations are approximated by a finite volume technique, with pressure and velocity co-located at the cell center; the algorithm for the discretization of the RANS equations is identical to the one described in [1] and it will be only briefly recalled here; the reader is referred to the cited paper for details. The residual on each control volume is computed as an interface flux balance Z 6 Z X o ~ ð~ F cs  ~ F ds Þ dS ¼ 0; ð7Þ q dV þ ot Vijk S s s¼1 where the surface integrals are evaluated by means of the (second order) trapezoidal rule. The computation of the viscous fluxes requires the value of the stress tensor at cell interfaces; for instance, at the interface ði þ 12 ; j; kÞ the velocity gradients in

  oum oul þ oxl oxm iþ1;j;k

ð8Þ

2

are computed by means of a standard second-order centered finite volume approximation  Z oum  1 ¼ um nl dS þ OðD2 Þ; ð9Þ oxl iþ1;j;k V iþ12;j;k R 1 2

DF ðx; y; z; tÞ ¼ 0: Dt

ui ðx; y; z; 0Þ ¼  ui ðx; y; zÞ;

sml jiþ1;j;k ¼ miþ12;j;k

871

iþ ;j;k 2

where the integral is extended to the surface Riþ12;j;k of the control volume V iþ12;j;k that includes the cell face S iþ12;j;k and is overlapped to half the cell (i, j, k) and half the cell (i + 1, j, k). For the inviscid part, a second-order Essentially NonOscillatory (ENO) scheme [53] has been adopted. The fluxes are computed as the solution of a Riemann problem ~ F ciþ1;j;k ¼ ~ F c ð~ qiþ12;j;k Þ ¼ ~ F c ð~ ql ;~ qr Þ, whose right and left states 2 are given by 1 ~ qji1=2;j;k ; D~ ql ¼ ~ qi;j;k þ minmod ðD~ qjiþ1=2;j;k Þ; 2 1 ~ qjiþ1=2;j;k ; D~ qr ¼ ~ qiþ1;j;k  minmod ðD~ qjiþ3=2;j;k Þ; 2 where D~ qjiþ1=2;j;k ¼ ~ qiþ1;j;k  ~ qi;j;k and function (to be applied to each vector as  0 minmod ðx; yÞ ¼ signðxÞ minðjxj; jyjÞ

ð10Þ

minmod(x, y) is a component) defined if xy 6 0; if xy > 0:

ð11Þ

In order to simplify the algorithm, a second-order accurate solution of the Riemann problem [1] is used in place of the exact one, which should be computed iteratively, given the nonlinearity of the problem. It is easy to prove that the resulting scheme is second-order accurate, and yields oscillation-free discrete solutions, also when the exact solutions are discontinuous (see [1] and the literature cited there). A two-levels multi-stages second-order Runge–Kutta scheme [54], is used for the integration in the pseudo-time of the discrete model; convergence to steady state is accelerated by local time stepping and a multi-grid algorithm [55]. 3.1. Single-phase level set approach The single-phase level set approach is described here. For the sake of completeness some details about the twophases level set method are also recalled. The readers interested in this topic are referred to the cited papers. In level set approaches [9,10], a smooth function /(x, y, z, t), whose zero level coincides for t = 0 with the free surface, is defined in the whole physical domain (i.e. in both liquid and air phase) as the signed distance from the free surface. This function is enforced to be a material surface for t > 0 o/ðx; y; z; tÞ þ~ u  r/ðx; y; z; tÞ ¼ 0; ot /ðx; y; z; 0Þ ¼ dðx; y; zÞ:

ð12Þ

872

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

Here ~ u is the velocity of the underlying flow, and d(x, y, z) the signed distance from the free surface at t = 0. Of course, the zero level of /(x, y, z, t), being convected by the flow (Eq. 3), will always represent the free surface location. Moreover, as /(x, y, z, t) is the signed distance from the surface of discontinuity, the sign of the level set function remains unchanged for t > 0 at material points. In two-phase level set approaches, the density and the molecular viscosity of the fluid are then assumed to depend on the sign of the level set function. In order to avoid sharp discontinuities, that give rise to numerical difficulties, density and viscosity are assumed to have a smooth transition around the zero level of /(x, y, z, t) [10]. Unfortunately, Eq. (12) does not enforce the function /(x, y, z, t) to remain a distance function for t > 0, whereas for the two-phase approach, this requirement is crucial in order to maintain the thickness of the interface constant in time and avoid mass loss [10]. Therefore, in the so-called re-initialization step, /(x, y, z, t) is replaced, at each time step, by a new ~ y; z; tÞ with the same zero level, that correfunction /ðx; ~ y; z; tÞ can sponds to the distance from the interface. /ðx; be computed as the steady state solution (with respect to the pseudo-time s) of the Eikonal equation ~ o/ ~  1Þ ¼ 0 þ signð/Þðjr/j ð13Þ os It is important to highlight that, due to the hyperbolic nature of Eq. (13), the solution evolves as a distance function close to the zero level set first, then it propagates outward in the whole domain (see, for instance, [10,47]); since in two-phase level set approaches the density and the viscosity are smoothed around the zero iso-surface of /(x, y, z, t), it is not necessary that the level set function be a distance function throughout the whole computational domain, but only close to its zero level, and therefore it is not necessary to reach the complete steady state solution of Eq. (13) [56]. In single-phase level set methods only the liquid phase of the fluid is computed and the level set function is used as a tracking device to locate the actual position of the free surface on a fixed grid. The reason for doing so is that, when attention is focused on high Reynolds number flows around complex geometries, the need to limit the number of grid points is of paramount importance. In fact, when using a two-phase approach, the size of the transition zone between water and air has to be O(6–8) grid cells in order to have a stable integration. This boundary can result excessively thick, given the limited number of points available in this region when simulating flows around complex geometries. It is worth to note that, by definition, the single-phase level set method maintains the interface sharp. In the approach described here, the computational domain is formally decomposed in three regions: (i) cells close to the surface of discontinuity (full circles and full squares in Fig. 1): these points are defined

Air i+1/2

i1/2

FS

i,j,k j+1/2

i+1,j,k

η

Water j1/2

Fig. 1. Computational domain: water and air regions; free surface detection and extrapolation of pressure.

as the nodes where at least one of the following six inequalities holds: /i;j;k /i1;j;k < 0;

/i;j;k /i;j1;k < 0;

/i;j;k /i;j;k1 < 0; ð14Þ

i.e. where the cell lies in one phase and at least one of its six neighborhoods is in the other. In these points, the level set function is computed by means of the evolution equation (12); velocity and pressure in the liquid region (full squares) are computed by means of the Navier–Stokes equations (7), whereas, in the air phase (full circles), pressure is evaluated by using the dynamic boundary condition and the velocity is extrapolated, as explained in the following; (ii) cells in the liquid phase region (empty squares in Fig. 1): these points are simply defined as those nodes where the level set function is negative and none of the inequalities (14) holds. Here the solution is computed by the numerical solution of the governing equation (7), and the level set function is enforced to be a distance function by means of Eq. (13); (iii) cells in the air region (empty circles in Fig. 1): these are the points where the level set function is positive and again none of the inequalities (14) holds. The level set function is computed from Eq. (13) in order to get a distance function, and an extension velocity is computed as done by Adalsteinsson and Sethian [57]. It has to be noted that in the present single-phase formulation, the solution outside the water region is not required; however, the extension of the velocity field outside the water region ensures second-order accuracy also close to the interface of discontinuity. 3.1.1. Computation of the level set function For the nodes of the computational domain close to the surface of discontinuity, i.e. where one of the conditions

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

(14) holds (full symbols in Fig. 1), the level set function is computed by the numerical solution of Eq. (12). When computing free surface flows around ship hulls with curvilinear grids, it is convenient to split the level set function as /(x, y, z, t) = u(x, y, z, t) + z, where u(x, y, z, t) is the solution (from Eq. (12)) of ou þ~ u  ru þ w ¼ 0: ot

ð15Þ

By doing so, it is easier to assign the boundary condition for the level set function at inflow, that reduces to u(x, y, z, t) = 0. Moreover, on any curvilinear mesh, a uniform flow is described exactly by the null function. The previous equation is discretized by using an ENO technique similar to the one used for the bulk flow; to this aim, the equation is first rewritten in terms of curvilinear coordinates ou ou þ Um þ w ¼ 0; ot onm

ð16Þ

m U m ¼ ui on being the contra-variant components of the oxi velocity vector. The derivatives of the function u(x, y, z, t) at cell center are approximated by an upwind based second-order finite difference formula; for instance, at the cell center (i, j, k), the derivative along n1 is

ou ¼ uiþ12;j;k  ui12;j;k ; on1

ð17Þ

where the interface value at ði þ 12 ; j; kÞ is computed as 1 uiþ12;j;k ¼ ui;j;k þ minmod ðDujiþ1;j;k ; Duji1;j;k Þ; 2 2 2

ð18Þ

if U 1i;j;k P 0 or 1 uiþ12;j;k ¼ uiþ1;j;k  minmod ðDujiþ3;j;k ; Dujiþ1;j;k Þ; 2 2 2

ð19Þ

if U 1i;j;k < 0, with Dujiþ1;j;k ¼ uiþ1;j;k  ui;j;k . It has been 2 proved in [53] that this procedure yields a second-order approximation to (15) in space. Time integration of Eq. (16) is achieved by using a standard two-level multi-stage second-order Runge–Kutta scheme. For the nodes which are not close to the surface of discontinuity (empty symbols in Fig. 1), the level set function is enforced to represent the signed distance from the interface (when asymptotic solution is attained). It has to be highlighted that, in two-phases level set approaches, the level set function has to be reinitialized as a distance function because of the need to keep the thickness of the transition zone between air and water constant (otherwise it would depend on the local behavior of the flow). In the present approach, this necessity is not felt, because there is no transition zone at the free surface. Nevertheless, it is useful to have this requirement satisfied, because the local structure of the function /(x, y, z, t) can spoil the accuracy of the interpolation required to locate the free surface; moreover, it can deeply affect the computation of the surface curvature, required when surface tension effects are taken into account. To this aim, the constraint j$/j = 1

873

is enforced by the computation of the asymptotic solution of the evolution equation   o/ r/ þ signð/Þ  r/  1 ¼ 0; ð20Þ ot jr/j then, by using the definition for u(x, y, z, t), the previous equation can be rewritten in term of this function as   ou r/ / þ sign ð/Þ  ru þ z  1 ¼ 0: ð21Þ ot jr/j jr/j Finally, by defining ~ w ¼ sign ð/Þ

r/ jr/j

 and

b ¼ sign ð/Þ

/z 1 jr/j

 ð22Þ

Eq. (21) can be written as ou þ~ w  ru þ b ¼ 0: ot

ð23Þ

This equation is solved with respect to u(x, y, z, t) with an ENO scheme analogous to the one used to solve the kinematic equation (15). It can be seen that the discretization of this equation yields a marching scheme that updates u(x, y, z, t) moving away from the free surface [10,47]. Therefore, if the solution converges toward a steady state, it provides a function /(x, y, z, t) that satisfies Eq. (12) on the surface /(x, y, z, t) = 0 and it is a distance function at any other points. Note that /(x, y, z, t) is a distance function also at points close to the free surface (full symbols in Fig. 1) because of Eq. (27) in the next section. 3.1.2. Free surface detection and solution prolongation Once the function u(x, y, z, t) is known throughout the whole domain, the free surface is located by the surface /(x, y, z, t) = u(x, y, z, t) + z = 0. The intersection of this surface with the underlying grid line is computed as follows (see Fig. 1). Consider, for instance, the coordinate line n1; when in two adjacent points (i, j, k) and (i + 1, j, k) the condition /i,j,k/i+1,j,k 6 0 holds, it means that the free surface cuts the segment ~ P iþ1;j;k  ~ P i;j;k at some point ~ P FS (~ P i;j;k is the position vector for the point (i, j, k)). The segment fraction below the free surface is g¼

j/i;j;k j j~ P FS  ~ P i;j;k j ¼ ~ ~ j/ jP iþ1;j;k  P i;j;k j iþ1;j;k  /i;j;k j

if /i;j;k < 0

ð24Þ

the level set function being the distance from the interface. A similar relation holds for /i+1,j,k < 0, with i and i + 1 interchanged. The computation of the residuals for the RANS equations at those points whose neighboring cells are not into the water region needs some attention. In fact, in these cases the numerical convective and viscous fluxes at interfaces that separate two cells, one of

which is in the air region (like the interface i þ 12 ; j; k in Fig. 1), must be evaluated; in these points the proper information to compute the correct flux are needed to retain second-order accuracy. To circumvent this difficulty, the procedure

874

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

described in what follows is applied. The pressure at the

interface i þ 12 ; j; k is extrapolated as 1 ðp  pi12;j;k Þ; þ g FS 2

piþ12;j;k ¼ pi12;j;k þ 1

ð25Þ

with pi12;j;k ¼ ðpi;j;k þ pi1;j;k Þ=2. Note that the use of pi12;j;k instead of pi,j,k to calculate piþ12;j;k provides a well behaved extrapolation also when g ! 0. In the above formula, pFS is computed from the dynamic boundary condition (2) pFS ¼

z j ; þ sij ni nj þ 2 We Fr

ð26Þ

with the derivatives in the second and third term on the right hand side computed approximation. by centered

The normal velocity at i þ 12 ; j; k is then computed by solving a Riemann problem with a constant pressure boundary condition, and the tangential velocity is simply extrapolated along the normal ~ n ¼ r/=jr/j to the free surface, as in the following Eq. (27). The remaining dynamic boundary conditions for the tangential stresses in (2) are explicitly enforced when

computing the viscous fluxes at the interface i þ 12 ; j; k . Outside the water region, extension velocities are computed as r/  rui ¼ 0;

i ¼ 1; 2; 3;

ð27Þ

which guarantees that /(x, y, z, t) evolves as a distance function also at the points adjacent to the free surface (full circles in Fig. 1; see, for the proof [57]). Note that, as the mesh size tends to zero, points in the water region for which Eq. (14) hold, approach the zero level of / (x, y, z, t), i.e. the free surface: this is the reference value from which the distance function is computed. The steady state solution being the goal of the computation, the previous relation is substituted by an evolution equation for the velocity components ui oui þ r/  rui ¼ 0; ot

ð28Þ

which is solved by a second-order ENO scheme analogous to the one used to solve the kinematic equations (15) and (23), with characteristic speed $/. All the other quantities, like turbulent viscosity, dissipation rate, or turbulent kinetic energy, are computed in the air region as in (28). Note that the values of the velocity, pressure and turbulent viscosity in the air phase are not used for the computation of residuals, as the values at the interface are computed from the dynamic free surface boundary conditions and the values of the solutions only in the water region. Therefore they do not affect neither the solution itself, nor the formal accuracy of the numerical scheme. Nevertheless, their estimation is of great importance during the iterative procedure at those points that change their physical state from air to water, for which an initial estimate is needed.

4. Results Three test cases have been considered in order to assess the performances of the single-phase LS algorithm: the flow past a submerged NACA-0012 profile is considered first, then the three-dimensional flow around a ship hull, namely the Series-60 Cb = 0.6 hull, is analyzed; finally, the three-dimensional flow, in both non-breaking and breaking regimes, around a ship hull with a complex geometry (characterized by a bulbous bow and a transom stern) is taken as an example of practical application. Verification and validation has been carried out in order to assess the reliability of the computed results for each test case. To this aim, each numerical solution has been computed at least on three grid levels, each obtained by removing every other point from the previous finer one, i.e. with a grid refinement ratio equal to 2. Verification and validation is carried out by using the methodology proposed by the ITTC Quality Manual [58] (see also [59]). For all the simulations performed, the iterative convergence error and uncertainty (hereafter dI and UI) have been found negligible with respect to the grid convergence components (dG and UG), and therefore, they have been neglected in the estimation of the numerical errors and uncertainties (i.e. dI  dG and UI  UG, and therefore, dSN  dG and USN  UG). Moreover, when non-monotonic convergence is observed, numerical uncertainty is evaluated by the formula proposed in [60], i.e. USN = 1.25jdREj, dRE being the correction factor computed from Richardson’s extrapolation. The solutions obtained with the single-phase LS algorithm are compared with available experimental data and, for the non-breaking flows, with other numerical solutions provided by surface fitting solvers. 4.1. Submerged NACA-0012 hydrofoil The wavy flow generated by an hydrofoil moving beneath the free surface has been the subject of extensive experimental studies (see, for instance, [61]) and thus it is often used as test case for the validation of numerical methods developed to compute free surface flows [47,62,63]. In this work, it has been chosen as a validation test case for the proposed algorithm because it is amenable for a throughout calculation of convergence properties and numerical uncertainties. In the experimental set-up, a NACA-0012 hydrofoil, whose chord is c = 20.3 cm, is towed in a tank with speed U1 = 0.8 m/s with an angle of attack a = 5 (see Fig. 2); the depth of submergence s is varied by varying the water level, whereas the distance between the profile and the bottom of the basin is kept fixed (H = 17.5 cm from the mid-chord of the profile). The relevant non-dimensional parameters based on the chord length and the free-stream velocity are Fr = 0.567 and Re = 1.42 · 105. Simulations will be reported for a non-dimensional depth of submergence equal to s/c = 1.034. In this case the leading wave of the wave train does not break [61] and, therefore, a classical surface fitting

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

875

1. 5

h(x,y,t) 1. 0

s

z /c

0. 5

U∞ 0. 0 -0.5 -1.0 -1.0

H

α

c

-0.5

0.0

0.5

1.0

1.5

-1.0

2.0

-0.5

0.0

0.5

1.0

1.5

2.0

x/c

x/c

Fig. 2. Flow around a submerged NACA-0012 hydrofoil: left, sketch of the problem; right, computational grid (every other point is shown).

method [1] can also be applied; the results obtained with the two numerical methods are compared with the experimental data. The two-dimensional computational domain has been discretized by means of a non-orthogonal non-conformal block structured grid; 19 blocks have been employed. The total length of the domain is about 19 chord lengths, with the inflow boundary placed five chord lengths upstream from the leading edge of the hydrofoil. In Fig. 2, a portion of the mesh is shown, where for the sake of clarity, every other point has been removed, and thick lines represent block boundaries. In order to study the convergence properties and to speed up the convergence rate, the solution has been computed by a multi-grid algorithm with four grid levels, each one obtained from the previous finer level by removing every other point. The finest computational grid contains about 100,000 control volumes, with 144 cells distributed on the hydrofoil. The height of the volumes on the surface of the hydrofoil is about 5 · 105 (’0.5 wall units), and approximately 30 cells are within the boundary layer. About 9400 cells lay in a layer above the unperturbed free surface level (the gaseous phase at the rest). A uniform velocity field with the upstream value, zero pressure and a flat free surface have been considered as initial conditions. Free-stream velocity is enforced at the inflow and at the

bottom boundaries, whereas zero pressure is assigned at the outflow boundary. The L2 norm of x-momentum residual and the values of drag and lift force coefficients are presented in Fig. 3. Good iterative convergence properties of the algorithm are clearly appreciable; indeed, time history for L2 norm of the x-momentum residual shows regular convergence toward steady state, with an almost constant rate for all grid levels. Moreover, from the time history of the hydrodynamic forces, it has to be noted that iterative convergence (in terms of constant values of horizontal and vertical forces) is attained on all grid levels after a few thousand multi-grid cycles (i.e. when the residual level is around 103). The solution at steady state on the finest grid level is presented in Fig. 4, where pressure and level set function contours are shown; as it can be seen, the level set function behaves as a distance function from the free surface in both water and air regions. The grid convergence study for the hydrodynamic force coefficients (lift CL and drag CD) is summarized in Table 1, where the computed values on the three finest grids (S3 denotes coarse, S2 medium, and S1 fine calculations), the measured orders of accuracy r, the correction factors CG and the numerical uncertainties USN are shown. For the drag component, the measured rate of convergence is equal to 1.75, very close to the

100

0.6

0.10

-1

10

0.08 0.4

-3

10

0.04

10-4

0.2

0.02

-5

10

-6

10

CL

0.06

CD

L2(u)

10-2

0

10000

20000

Cycle

30000

40000

0.00

D L 0

10000

20000

30000

0.0 40000

Cycle

Fig. 3. Flow around a submerged NACA-0012 hydrofoil: convergence histories of, left, L2 norm of the x-momentum residual and, right, lift and drag coefficients.

876

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886 1.5 0.2

1.0 -0.05

0.5

0.15

-0.20

0.05

z /c

-0.20

-0.15

0.15

-0.05

0.05

-0.4 -0.8

0.0

0.1 0

-1.2

0.00

-0.5

-1.6

-1.0 -0.5

0.5

1.5

2.5

3.5

4.5 -0.5

0.5

1.5

2.5

3.5

4.5

x /c

x /c

Fig. 4. Flow around a submerged NACA-0012 hydrofoil: left, pressure contours, pmin = 0.9, pmax = 0.5, Dp = 0.05; right, level set function contours, /min = 1.9, /max = 0.3, D/ = 0.1. The thick line is /(x, y, z, t) = 0.

Table 1 Flow around a submerged NACA-0012 hydrofoil: grid convergence and numerical uncertainty studies S3

S2

S1

r

CG

USN

CD

0.0265

0.78

0.5237

0.0191 (8.2%) 0.5373 (0.7%)

1.75

CL

0.0208 (21.5%) 0.5413 (3.4%)

–

–

1.56 · 104 (0.8%) 1.67 · 103 (0.3%)

r

CG

USN

2.47 1.72 1.49

1.52 0.76 0.60

2.81 · 103 4.79 · 104 3.92 · 103

D4h 2h

D2h h 2

3

3.49 · 10 1.50 · 102 2.00 · 102

u w p

6.29 · 10 4.57 · 103 7.10 · 103

theoretical value of 2; small values for the numerical uncertainty (about 2 · 104, i.e. about 1%), due to small variations with grid refinement, indicate that grid convergence is almost attained. For the lift component, instead, although convergence is not monotonic, the difference between the computed force on the medium and fine meshes, which is less than 1% (and around 3% between coarse and medium calculations), and the numerical uncertainty of about 2 · 103 (less than 1% of the finest value), clearly indicate that grid converged results have been attained. In Table 1, the grid convergence study for field variables, namely the horizontal and vertical velocity components, and the pressure, is also provided. In this table D2h h

indicates the L2 norm of the differences between the solution computed on two successive grids (fine and medium in this case); moreover, since the correction factors for u and p are considered far from one, the so-called ‘‘lacking confidence’’ formula is used for estimating the numerical uncertainties [58]. As it can be seen, fairly good rates of convergence have been obtained also for field variables; numerical uncertainties are of the same order of magnitude as those of global coefficients. Fig. 5 shows the wave elevations computed on three grid levels: the reduction of the difference between the solutions on subsequent grids clearly indicates the convergence properties of the algorithm when refining the grid. The convergence rate for the wave elevation is 1.88 (measured within

0.12

Wave Elevation

0.08 0.04 0.00 -0.0 4 -0.0 8 -0.1 2 -2.0

0.0

2.0

x /c

4.0

6.0

-2.0

0.0

2.0

4.0

6.0

x /c

Fig. 5. Flow around a submerged NACA-0012 hydrofoil: left, wave profile computed on (—) fine, (- - -) medium and (    ) coarse mesh; right, comparison between (—) level set solution, (- - -) surface fitting solution [63], and (d) experimental data [61].

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

the x-range show in Fig. 5), and hence very close to the expected theoretical value. Moreover, also for the wave profile the very small difference between the solutions provided by the medium and the fine calculations clearly indicates that grid converged results have almost been reached. In Fig. 5, the wave elevation computed with the singlephase level set approach is compared with the experimental data by Duncan [61] and with the solution from a surface fitting calculation [63] obtained by using the same grid. It can be seen that the two numerical solutions are very close to each other, even though the level set solution seems to under predict slightly the wave height with respect to the surface fitting solution. However, it has to be noted that, in the surface fitting calculation (the grid being adapted to the actual free surface position), all the control volumes are in the water, and as a results, the grid at the free surface is more refined. Nevertheless, the agreement with the experimental data is pretty good almost everywhere. 4.2. Flow around Series-60 Cb = 0.6 ship hull In this section, the described single-phase algorithm is applied to the simulation of the three-dimensional steady flow around a surface piercing ship hull, namely the Series-60 Cb = 0.6. Results are compared with both experimental data, measured in the towing tank at the Iowa Institute of Hydraulic Research (IIHR) [64], and numerical results obtained with a surface fitting algorithm. Grid convergence analysis is carried out for the total resistance coefficient, as well as for its friction and pressure contributions. The values of the non-dimensional parameters, based on the ship length between perpendiculars Lpp and free stream velocity U1, are Fr = 0.316 and Re = 4 · 106. The physical domain is discretized by means of a two blocks grid with C–O topology (see Fig. 6) with a total of 144 · 96 · 48 volumes along the stream-wise, normal and girth-wise direc-

877

tions, respectively; points are clustered toward the solid wall in order to resolve the boundary layer (no wall function was used). Inlet, outlet and side boundaries are located at 1.5Lpp far from the hull surface. Five grid levels have been employed for multi-grid acceleration, whereas only the three finest ones are used for convergence assessment. Convergence histories for the L2 norm of x-momentum residual and for the total resistance coefficient C T ¼ R= 12 qU 21 S 0 (S0 being the wetted surface at rest, and R the total dimensional force along longitudinal axis) are provided in Fig. 7; it can be seen that iterative convergence and steady value for the force are attained on all grid levels, after a small number of multi-grid cycles. Moreover, from the history of the resistance, it can be observed that the uncertainty due to incomplete iterative convergence is negligible. In the same figure the convergence of the computed numerical values toward experimental data (dashed line [64]) can be clearly inferred, as the value provided by the level set calculation on the finest grid almost recovers the measurement. Grid convergence study for the total resistance coefficient CT, and for its pressure RP and friction RF contributions is summarized in Table 2. The values in parentheses represent the variation of the solution with respect to the previous coarser mesh computation; the experimental data D for the total resistance coefficient [64] is also reported. The actual convergence rate is shown in the second last column, and, finally, the numerical uncertainties are reported in the last column; for this quantity, it has to be noted that, for CT and CP, the correction factors obtained could be considered close to one, and, therefore, both the extrapolated values and the numerical uncertainties could be computed (as well as the uncertainties in percent of the extrapolated value); on the other side, for the friction contribution the ‘‘not having confidence’’ criteria has been employed for the estimation of the numerical uncertainty, CG being far from one.

Fig. 6. Flow around Series-60 ship hull: particular of the mesh (every other point is shown) in the hull region. Free-stream flow is from right to left.

878

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886 101

15

10

10-1

3

10 CT

L2(u)

100

10-2

5

10-3 10-4

0

500

1500

1000

2000

2500

0

0

500

Cycle

1000

1500

2000

2500

Cycle

Fig. 7. Flow around Series-60 ship hull: convergence histories of, left, L2 norm of the x-momentum residual, right, total resistance coefficient (dashed line, experimental value from [64]). Table 2 Flow around Series-60 ship hull: grid convergence and numerical uncertainty studies for the resistance coefficients 103S3

103S2

103S1

103SC

103D

r

103USN

CT

7.92

5.96

1.61

3.38

1.65

N/A

1.74

CF

4.54

6.06 (7.0%) 1.77 (17.3%) 4.29 (2.0%)

5.91

CP

6.52 (17.7%) 2.14 (36.7%) 4.38 (3.5%)

–

N/A

0.83

0.07 (1.2%) 0.03 (2.1%) 0.21 (2.7%)

It can be seen from Table 2 that the agreement between the computed value of the total resistance coefficient on the finest grid level and the value provided by the experiments is very good (the difference is less than 2%), and its extrapolated value almost recovers the experimental data (with an error equal to 0.8% of D). The variations of the computed value when refining the grid are very small and consequently equally small is the numerical uncertainty for CT (around 1%). Since the error with respect to the experimental data is less than the numerical uncertainty, the corrected value of CT is validated (even not considering the uncertainty on the experimental measurements) within 1.2% of D. Similar considerations can be done for the pressure contribution, even if nothing can be said about validation because of the obvious lack in the experimental data; moreover, in spite of the large variations of the computed data when refining the grid, the measured order of accuracy is very close to the theoretical value. On the other side, the

convergence order for the friction resistance coefficient is far from the expected value; however, it has to be noted that the very small grid dependence of this coefficient (less the 3.5% for the medium mesh, and even less for the finest one) suggests that a grid converged result is attained also for CF. In Fig. 8, the computed wave profiles and the comparison with both experimental data and the solution obtained by a surface fitting approach on a similar grid are presented. The wave profiles on the three finest levels clearly show grid convergence, with remarkably small variations for the bow crest between the fine and the medium grid calculations. Finally, good agreement with both experiments and surface fitting simulation can be clearly inferred from the figure. The computed wave pattern is presented in Fig. 9, again in comparison with both experimental measurements and surface fitting results; in Fig. 10 some wave cuts at different

0.015

Wave Elevation

0.010 0.005 0.000 -0.005 -0.010 -0.015 -0.5

-0.3

0.0

x/Lpp

0.2

0.5 -0.5

-0.3

0.0

0.2

0.5

x/Lpp

Fig. 8. Flow around Series-60 ship hull: left, wave profile computed on (—) fine, (- - -) medium and (    ) coarse mesh; right, comparison between level set solution (—), surface fitting solution (- - -), and experimental data (j) from [64].

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

879

0.6

y/Lpp

0.4

0.2

0.0 0.6

y/Lpp

0.4

0.2

0.0 -0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

x/Lpp Fig. 9. Flow around Series-60 ship hull: wave pattern; top, comparisons of the level set solution (—) with the surface fitting solution (- - -); bottom, comparison with the experiments [64] (- - -).

0.012

y/Lpp = 0.1

y/Lpp = 0.2

y/Lpp = 0.3

Wave Elevation

0.006

0.000

-0.006

-0.012 -0.8

-0.4

0.0

0.4

0.8 -0.8

-0.4

x/Lpp

0.0

x/Lpp

0.4

0.8 -0.8

-0.4

0.0

0.4

0.8

x/Lpp

Fig. 10. Flow around Series-60 ship hull: wave cuts at different distances from longitudinal axis, (—) level set, (- - -) surface fitting and (j) experimental data [64].

distances from the longitudinal axis are shown. Very good agreement with both surface fitting solution and experimental data can be clearly observed. In particular, for both the single-phase LS and the surface fitting solution, the differences with experimental data are quite small in the region close to the hull surface, as it can be clearly seen from both the wave pattern figure and from the longitudi-

nal wave cuts at y/Lpp = 0.1, 0.2; the departure of the simulation from data becomes relevant only in the far field (see, for example, the wave cut at y/Lpp = 0.3), where the grid is too coarse to get a good prediction (see Fig. 6). It is worth to note that, despite of the error in wave amplitude caused by numerical dissipation, the agreement in wave phase is still remarkably good.

880

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

4.3. Flow around a surface combatant From the previous test cases, it has been demonstrated that the present single-phase LS algorithm is robust and accurate in the simulation of both two and three-dimensional free surface flows, and it is at least as accurate as a similar surface fitting solver. In this section, as an example of practical application to flows with breaking waves, the numerical simulations of the flow around a realistic ship hull are presented. The ship under consideration is the naval combatant DDG-51, whose geometry is characterized by both a sonar dome and a transom stern (see Fig. 11). This ship has been taken as benchmark for an international collaboration between three institutes (namely the David Taylor Model Basin, the Iowa Institute of Hydraulic Research and the Italian Ship Model Basin, INSEAN), and therefore, extensive experimental results are available [65]; moreover, in the last two workshops on CFD for ship hydrodynamics held in Gothenburg in 2000 [2] and in Tokyo in 2005 [27], it has been used as representative of modern ship hull forms for a surface naval combatant. The test analyzed here corresponds to the towing condition in still water at fixed trim and sinkage; two different flow regimes are considered: the first one is the same proposed for the workshops, with a Froude number equal to 0.28 and a Reynolds number of 1.26 · 107 at model scale; the second case is characterized by the breaking of the bow wave, with Froude and Reynolds numbers equal to 0.41 and 1.85 · 107, respectively. The geometry parameters correspond to a ship model with a length between perpendiculars equal to 5.72 m (scale factor k = 1/24.8) that corresponds to the model used by both the David Taylor Model Basin and INSEAN for measurements in the towing tank. The physical domain is discretized by a block structured non-conformal grid with 19 blocks and C–O topology; the total number of volumes is

about 2.5 · 106; points are clustered toward both solid walls and the free surface region. Part of the grid is shown in Fig. 11. Five grid levels have been used for multi-grid acceleration, and, as usual, the finest three are used for both grid convergence study and uncertainty assessment. The grid used for the two test cases were similar as total number of volumes and block topology, but the higher Froude number simulation the volumes were clustered in the breaking region. An overview of the flow field is given in Fig. 12, where the free surface shape at steady state for the higher Froude number is shown. This test case is very severe to check the capabilities of free surface simulation algorithms in dealing with complex three-dimensional wave breaking phenomena, because of the formation of a water sheet on the hull at the bow, the formation of multiple jets, impingements and splash-ups (see Fig. 12 and the close-up view of the bow region in Fig. 13). The wave patterns obtained by the numerical simulations are presented in Fig. 14 for both Froude numbers, where the contour lines of the wave elevation g are plotted. As expected, at Fr = 0.41 larger wave elevations can be observed, minimum and maximum wave heights being g = 0.013, +0.026; for the lower Froude number, extreme values are almost halved (g = 0.007, +0.014). For the higher Froude number, it is interesting to note the divergent wave system downstream the bow breaking wave, whose crests are almost aligned with the free stream, and the so-called rooster tail system at the transom stern. In Figs. 15 and 16, it can be noticed that some interesting features of the breaking waves are well resolved by the numerical simulation (see the works by Dong et al. [66] and Waniewski et al. [67]): as the flow impinges on the bow, a strong vertical velocity component induces a motion which creates a liquid sheet along the wall (see Fig. 16a and the enlarged view of the bow region in Fig. 13); the formation of a jet that impinges on the free

Fig. 11. Flow around a surface combatant: particular of the multi-block mesh (every other point is shown). Free-stream flow from right to left.

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

881

Fig. 12. Flow around a surface combatant: perspective view of the wave pattern at Fr = 0.41.

Fig. 13. Flow around a surface combatant: close-up view of wave pattern at Fr = 0.41 around the bow region.

surface follows (Fig. 16b and c). A second jet due to the splash-up and its plunging can be clearly observed in Fig. 16d, followed by the reconnection phase (Fig. 16e and f). The effects of the breaking bow wave on the boundary layer and the free surface vertical flow is beyond the scope of the paper. A top view of the breaking wave pattern is shown in Fig. 14 (the dashed lines correspond to the cuts in Fig. 16), where the formation of the jet, the point where impact occurs for the first time, as well as lines of the first and second impacts are highlighted by thick

lines. In these figures, the prediction of the shoulder wave, and the wave system at the dry transom, composed by both the transverse and the divergent front waves can be clearly observed. Like for the bow breaking wave, the flow around a transom stern is a typical flow which would be extremely difficult to compute with a surface fitting approach, because the severe grid distortion around the sharp edge and the possible topological change that would happen if the surface deformation is such that an entire block moves completely outside the water region. Differently from

882

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

0.50

Fr=0.28

y/Lpp

0.25

0.00

-0.25

Fr=0.41 -0.50 -0.6 0.00

-0.4

-0.2

0.0

0.2

y/Lpp

Jet Profile

0.4

0.6

0.8

1.0

Impact Line

First Impact

-0.10

Second Impact

-0.20

-0.5

-0.4

-0.3

-0.2

-0.1

x/Lpp Fig. 14. Flow around a surface combatant: top, wave pattern at Fr = 0.28 and Fr = 0.41; bottom, close-up view of the bow region at Fr = 0.41. gmin = 0.013, gmax = 0.026, Dg = 0.001.

Fig. 15. Flow around a surface combatant: left, transverse wave cuts at the bow region, right, longitudinal wave cuts at the transom region (Fr = 0.41).

experimental observations, breaking waves at the transom stern did not appear in the simulation, because of lack in the grid resolution. Convergence histories for the L2 norm of x-momentum residual and for the total resistance coefficient are reported in Fig. 17; iterative convergence and steady value for the

force are attained on all grid levels, for both Froude numbers. Convergence of the computed numerical values toward experimental data can be clearly inferred. In Fig. 18, the computed wave profiles on the hull when convergence is attained are compared with the results provided by the surface fitting calculation on a similar grid

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

a

b

c x=-0.43

x=-0.49

e

d

883

x=-0.39

f x=-0.30

x=-0.31

x=-0.33

10 0

20

10 -1

15

10 3CT

L2(u)

Fig. 16. Flow around a surface combatant: transversal wave cuts at different x stations (Fr = 0.41).

10 -2 10 -3 10 -4

10 5

0

500

1000

0 0

1500

500

Cycle

1000

1500

Cycle

10 0

15

10 -1 10

10 3CT

L2 (u)

10 -2

5 10 -3 10 -4

0

2000

4000

6000

8000

Cycle

0

0

2000

4000

6000

8000

Cycle

Fig. 17. Flow around a surface combatant: convergence histories of, left, L2 norm of the x-momentum residual, right, total resistance coefficient; top, Fr = 0.28, bottom, Fr = 0.41; dashed line, experimental value from [68].

(which was possible only for the case at Fr = 0.28), and with the experimental data measured by Olivieri et al. [68] at both Froude numbers. At the lower Froude number the agreement with the experimental data is satisfactory, in spite of the large deviation at the first crest; such disagreement has been already observed in all the calculations in [2], and it should be due to the way the experimental wave profile is taken. In fact, a photographic and digitizing sys-

tem was used; therefore, if the maximum wave height is not attained on the hull surface (as highlighted in Fig. 16), the result in the photo does not show the profile on the hull but the envelope of the maximum wave heights. This conjecture is supported by the DTMB measurements of the wave profile on a similar model; in this institute a waterproof pencil is used and the first crest of the wave profile is somewhat lower (see Fig. 5 in [65]). However, the agreement

884

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

Wave Elevation

0.020

0.03

0.015

0.02

0.010 0.01 0.005 0.00

0.000 -0.005 -0.5

-0.3

0.0

0.2

0.5

-0.01 -0.5

-0.3

0.0

0.2

0.5

x/Lpp

x/Lpp

Fig. 18. Flow around a surface combatant: wave profiles for, left, Fr = 0.28 and, right, Fr = 0.41. (—) level set, (- - -) surface fitting [3] and (j) experimental data [68].

with experimental data for the crest at the bow is well within the experimental uncertainties [69]; moreover, for the shoulder and stern waves, the agreement is very satisfactory, and even better than the surface fitting computation. For the higher speed, the agreement with the experimental data is very good along the whole hull surface; it is needless to say that numerical algorithms based on surface fitting technique fail for the computation of this flow, because of the overturning waves. Finally, in Table 3 the computed extrapolated values for the total resistance coefficient together with experimental data D taken from [68] and two other numerical solutions (computed with surface fitting approaches) are reported for comparison purposes (SI from [3] by using a similar grid and SII from [4]). As it can be seen, extrapolated values are very close to the measurements, the error being less than 4% and 1% of D for low and high Froude number case, respectively. Moreover, for the low Froude number case, the value provided by the LS calculation is in good agreement with the surface fitting solutions. Verification and validation for the total resistance coefficient is given in Table 4; experimental uncertainties UD in this table are computed following [69], in order to take into account also facility uncertainty. It can be clearly noticed that small numerical Table 3 Flow around a surface combatant: comparisons with both experimental and numerical data for 103CT Fr

SC

SI

SII

D

0.28

4.39 (3.8% D) 6.72 (0.6% D)

4.27 (1.0% D) –

4.36 (3.1% D) –

4.23

0.41

6.76

Table 4 Flow around a surface combatant: verification and validation of total resistance coefficient CT Fr

r

E

UD

USN

UV

0.28

1.68

0.41

1.72

0.16 · 103 (3.8% D) 0.04 · 103 (0.7% D)

0.27 · 104 (0.7% D) 0.26 · 103 (3.8% D)

0.92 · 104 (2.1% SC) 0.76 · 104 (1.1% SC)

0.96 · 104 (2.3% D) 0.27 · 103 (4.0% D)

uncertainties USN are attained for both non-breaking and breaking computations and the observed convergence rate is close to the theoretical value of two. It can be also noticed that the error E with respect to the experimental value is larger for the low Froude number case than for the higher velocity; this could be attributed to the fact that for Fr = 0.41 the wave pattern is more resolved, the wave being longer. In spite of the low numerical uncertainty (2.1% SC), the calculation is not validated for Fr = 0.28, the error (3.8% D) being larger than the validation uncertainty (2.3% D). This is probably due to modelling errors, like turbulence model and lack of the surface tension effects, that become relevant when wave resistance does not contain strong dissipation effects due to breaking. On the contrary, for the higher Froude number case, these modelling errors are overshadowed by the breaking wave; consequently the comparison error is much lower than the validation uncertainty and, therefore, the total resistance is validated within 4% of D. 5. Conclusions The single-phase Level Set approach has been applied to the simulations of free surface viscous flows at high Reynolds number. The capabilities of the algorithm in dealing with problems in the naval hydrodynamics context have been proved by means of a number of test cases, ranging from simple two-dimensional flow around a submerged hydrofoil to the three-dimensional flow in breaking regime around a surface piercing ship hull. In the single-phase LS method, only the liquid phase is simulated and the level set function is used as a tracking device to locate the actual position of the free surface. As a consequence, the interface between water and air remains sharp, and, therefore, differently to two phases LS approaches, there is no transition zone across which the fluid properties have to be smoothed; this also allows to increase the resolution in the free surface region without a large increase in both CPU time and memory requirements. On the other side, this approach is not suitable neither for problems where air entrapment phenomena are present, like bubbly flows, nor for problems where stresses on the free surface from the air region

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886

play an important role, like the simulation of wind-induced waves. The features of the single-phase LS algorithm in dealing with complex two-dimensional and three-dimensional free surface flows in both non-breaking and breaking regimes have been assessed in terms of numerical uncertainty and comparison with experiments. No problems have been observed when simulating complex free surface phenomena, like breaking and reconnection of the interface, as well as no numerical difficulties emerged when dealing with the flow around complex naval configurations such as a ship hull with a transom stern and, in general, a satisfactory agreement with available experimental and numerical data have been observed. Moreover, numerical experience showed that this approach is robust and provides good results in terms of both grid convergence properties and reliability of the results; second-order convergence rate and low validation uncertainties are achieved in both two-dimensional and three-dimensional simulations, even in the presence of complex geometries and/or breaking flow phenomena. The long term goal being the development of a general tool for the analysis of practical naval problem like seakeeping and maneuverability by means of unsteady RANS simulations, subjects of future research will be the extension of the single phase-level set method to unsteady problems, in the framework of dynamical overlapping grids discretization. In particular, the analysis of the flow field around ships with prescribed or predicted motion, ship maneuvering in confined waters, hulls with moving appendages, and the analysis of the hull/propeller interaction with the flow around the operating propeller simulated in real geometry, will be topics of future research. References [1] Di Mascio A, Broglia R, Favini B. A Second Order Godunov-type Scheme for Naval Hydrodynamics. Kluwer Academic/Plenum Publishers; 2001. p. 253–61. [2] Larsson L, Bertram V, Stern F. Proceedings of GOTHENBURG 2000. In: A workshop on CFD in ship hydrodynamics, Gothenburg, Sweden, 2000. [3] Di Mascio A, Muscari R, Broglia R. Computational of the flow past the US Navy combatant DTMB 5415 by a Godunov-type scheme. In: Proceedings of GOTHENBURG 2000, A workshop on CFD in ship hydrodynamics, Gothenburg, Sweden, 2000. [4] Wilson R, Paterson E, Stern F. Verification and validation for RANS simulation of a naval combatant. In: Proceedings of GOTHENBURG 2000, A workshop on CFD in ship hydrodynamics, Gothenburg, Sweden, 2000. [5] Harlow F, Welch J. Numerical calculation of time-dependent viscous incompressible flow of fluids with free surface. Phys Fluids 1965;8:2182–9. [6] Unverdi SO, Tryggvason G. A front-tracking method for viscous, incompressible, multi-fluid flows. J Comput Phys 1992;100:25–37. [7] Tryggvason G, Bunner B, Esmaeeli A, Juric D, Al-Rawahi D, Tauber W, et al. A front-tracking method for computations of multiphase flow. J Comput Phys 2001;169:708–59. [8] Hirt CW, Nichols BD. Volume-of-fluid (VOF) method for dynamics of free boundaries. J Comput Phys 1981;39:201–21.

885

[9] Osher S, Sethian JA. Fronts propagating with curvature-dependant speed: algorithms based on Hamilton–Jacobi formulations. J Comput Phys 1988;79:12–40. [10] Sussman M, Smekerda P, Osher SJ. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114:146–59. [11] Scardovelli R, Zaleski S. Direct numerical simulation of free-surface and interfacial flow. Ann Rev Fluid Mech 1999;31:567–603. [12] Chang YC, Hou TY, Merriman B, Osher S. A level set formulation of Eulerian interface capturing methods for incompressible fluid flows. J Comput Phys 1996;124:449–64. [13] Russo G, Smereka P. A remark on computing distance functions. J Comput Phys 2000;163:51–67. [14] Sussman M, Fatemi E, Smereka P, Osher S. An improved level set method for incompressible two-phase flows. Comput Fluids 1998;27:663–80. [15] Sussman M, Fatemi E. An efficient interface preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow. J Sci Comput 1999;20:1165–91. [16] Sussman M, Puckett G. A coupled level-set and volume-of-fluid method for computing 3D and axisymmetric incompressible twophase flow. J Comput Phys 2000;162:301–37. [17] Sussman M, Dommermuth D. The numerical simulation of ship waves using cartesian grid methods. In: Proceedings of 23th Symposium on Naval Hydrodynamics, Val de Reuil, France, 2000, p. 46–61. [18] Sussman M. A second order coupled level set and volume-of-fluid method for computing growth and collapse of vapor bubbles. J Comput Phys 2003;187:110–36. [19] Sethian JA, Smereka P. Level set methods for fluid interfaces. Ann Rev Fluid Mech 2003;35:341–72. [20] Osher S, Fedkiw RP. Level set methods: an overview and some recent results. J Comput Phys 2001;169:463–502. [21] Rhee S, Hino T. Unstructured grid flow solver for practical ship hull. In: Proceedings of GOTHENBURG 2000. A workshop on CFD in ship hydrodynamics, Gothenburg, Sweden, 2000. [22] Cura Hochbaum A, Vogt M. Toward the simulation of sea-keeping and maneuvering based on the computation of the free surface viscous flow. In: Proceedings of 24th Symposium on Naval Hydrodynamics, Fukuoka, Japan, 2002, p. 141–55. [23] Di Mascio A, Muscari R, Broglia R. A level set approach for naval applications. In: Proceedings of 13th international offshore and polar engineering conference, vol. III, Honolulu, Hawaii, USA, 2003, p. 220–6. [24] Di Mascio A, Muscari R, Broglia R. Computation of free surface flows around ship hulls by a level-set approach. In: Proceedings of the 8th international conference on numerical ship hydrodynamics, Busan, Korea, 2003. [25] Broglia R, Di Mascio A, Muscari R. Numerical simulations of breaking wave around a wedge. In: Proceedings of 25th symposium on naval hydrodynamics, St. John’s Newfoundland and Labrador, Canada, 2004. [26] Wilson R, Carrica P, Hyman M, Stern F. A steady and unsteady single phase level set method for large amplitude ship motions and maneuvering. In: Proceedings of 25th symposium on naval hydrodynamics, vol. IV, St. John’s Newfoundland and Labrador, Canada, 2004, p. 88–104. [27] Hino T. Proceedings of Tokyo 2005 Workshop. In: A workshop on CFD in ship hydrodynamics, Tokyo, Japan, 2005. [28] Broglia R, Di Mascio A, Muscari R. Computation of free surface turbulent flows around ship hulls by a RANS solver. In: Proceedings of the workshop on CFD ship hydrodynamics, March 09–11, Tokyo, Japan, 2005, p. 498–503. [29] Wilson R, Carrica P, Stern F. RANS simulation of a container ship using a single-phase level set method with overset grids. In: Proceedings of the workshop on CFD ship hydrodynamics, March 09–11, Tokyo, Japan, 2005, p. 546–51. [30] Miller R, Gorski J, Wilson R, Carrica P. RANS simulation of a naval combatant using a single-phase level set method with overset grids.

886

[31]

[32]

[33]

[34]

[35]

[36]

[37]

[38]

[39]

[40]

[41] [42]

[43]

[44]

[45]

[46]

[47]

A. Di Mascio et al. / Computers & Fluids 36 (2007) 868–886 In: Proceedings of the workshop on CFD ship hydrodynamics, March 09–11, Tokyo, Japan, 2005, p. 575–80. Kim J, Park I, Van S. RANS computations for KRISO container ship and VLCC tanker using the WAVIS code. In: Proceedings of the workshop on CFD ship hydrodynamics, March 09–11, Tokyo, Japan, 2005, p. 598–603. Tahara Y, Wilson R, Carrica P. Comparison of free-surface capturing and tracking approaches in application to modern container ship and prognosis for extension to self-propulsion simulator. In: Proceedings of the workshop on CFD ship hydrodynamics, March 09–11, Tokyo, Japan, 2005, p. 604–11. Hino T, Sato Y. Ship flow computations by an unstructured Navier– Stokes solver. In: Proceedings of the workshop on CFD ship hydrodynamics, March 09–11, Tokyo, Japan, 2005, p. 641–5. Rhee S, Skinner C. Unstructured grid based Navier–Stokes solver for free-surface flow around surface ship. In: Proceedings of the workshop on CFD ship hydrodynamics, March 09–11, Tokyo, Japan, 2005, p. 504–9. Hsu K, Chen Y, Chau S, Chien H, Kouh J. Ship flow computation of DTMB 5415. In: Proceedings of the workshop on CFD ship hydrodynamics, March 09–11, Tokyo, Japan, 2005, p. 516–9. Deng G, Guilmineau E, Queutey P, Visonneau M. Ship flow simulations with ISIS CFD code. In: Proceedings of the workshop on CFD ship hydrodynamics, March 09–11, Tokyo, Japan, 2005, p. 530–8. Chao K. Numerical propulsion simulation for the KCS container ship. In: Proceedings of the workshop on CFD ship hydrodynamics, March 09–11, Tokyo, Japan, 2005, p. 539–45. Carrica P, Wilson R, Stern F. Linear and nonlinear response of forward speed diffraction for a surface combatant. In: Proceedings of the workshop on CFD ship hydrodynamics, March 09–11, Tokyo, Japan, 2005, p. 552–7. Cura-Hachbaum A, Pierzynski M. Flow simulation for a combatant in head waves. In: Proceedings of the workshop on CFD ship hydrodynamics, March 09–11, Tokyo, Japan, 2005, p. 593–7. Fedkiw R, Aslam T, Merriman B, Osher S. A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the Ghost fluid methods). J Comput Phys 1999;154:383–427. Kang M, Fedkiw R, Liu L. A boundary condition capturing method for multiphase incompressible flow. J Sci Comput 2000;15:323–60. Nguyen D, Fedkiw R, Kang M. A boundary condition capturing method for incompressible flame discontinuities. J Comput Phys 2001;172:71–98. Enright D, Fedkiw R, Ferziger J, Mitchell I. A hybrid particle level set method for improved interface capturing. J Comput Phys 2002;183:83–116. Enright D, Nguyen D, Gibou F, Fedkiw R. Using the particle level set method and a second order accurate pressure boundary condition for free surface flows. In: Proceedings of 4th ASME JSME joint fluids engineering conference, Honolulu, Hawaii, USA, 2003. Kim M, Lee W. A new VOF-based numerical scheme for the simulation of fluid flow with free surface. Part I: New free surface tracking algorithm and its verification. Int J Numer Methods Fluids 2003;42:765–90. Bet F. Ha¨nel D, Sharma S. Numerical simulation of ship flow by a method of artificial compressibility. In: Proceedings of 22th symposium on naval hydrodynamics, vol. II, Washington, DC, USA, 1998, p. 173–82. Vogt M, Larsson L. Level set methods for predicting viscous free surface flows. In: Proceedings of the 7th international confer-

[48]

[49]

[50]

[51] [52] [53]

[54]

[55]

[56]

[57] [58] [59] [60] [61] [62] [63] [64]

[65]

[66] [67] [68]

[69]

ence on numerical ship hydrodynamics, Nantes, France, 1999, p. 2.1.4–19. Foster N, Fedkiw R. Practical animation of liquids. In: Proceedings of SIGGRAPH 2001 conference, Los Angeles, CA, USA, 2001, p. 15–22. Carrica P, Wilson R, Stern F. Unsteady RANS simulation of the ship forward speed diffraction problem. Comput Fluids 2006;35: 545–70. Carrica P, Wilson R, Stern F. An unsteady single-phase level set method for viscous free surface flows, Int J Numer Methods Fluids, in press. Spalart PR, Allmaras SR. A one-equation turbulence model for aerodynamic flows. La Recherche Ae´rospatiale 1994;1:5–21. Chorin A. A numerical method for solving incompressible viscous flow problems. J Comput Phys 1967;2:12–26. Harten A, Engquist B, Osher S, Chakravarthy SR. Uniformly high order accurate essentially non-oscillatory schemes. J Comput Phys 1987;71:231–303. Jameson A, Schmidt W, Turkel E. Numerical solutions of the Euler equations by finite volume methods using Runge–Kutta timestepping schemes, AIAA paper 81-1259. Favini B, Broglia R, Di Mascio A. Multi-grid acceleration of second order ENO schemes from low subsonic to high supersonic flows. Int J Numer Methods Fluids 1996;23:589–606. Sethian JA. Evolution, implementation, and application of level set and fast marching methods for advancing fronts. J Comput Phys 2001;169:503–55. Adalsteinsson D, Sethian JA. The fast construction of extension velocities in level set methods. J Comput Phys 1999;148:2–22. ITTC Quality Manual, Resistance Committee of 23th, ITTC 2001. AIAA, Guide for the verification and validation of computational fluid dynamics simulations, G-077-1998, 1998. Roache PJ. Quantification of uncertainty in computational fluid dynamics. Ann Rev Fluid Mech 1997;29:123–60. Duncan JH. The breaking and non-breaking wave resistance of a two-dimensional hydrofoil. J Fluid Mech 1983;126:507–20. Rhee SH, Stern F. RANS model for spilling breaking waves. J Fluid Eng 2002;124:1–9. Muscari R, Di Mascio A. A model for the simulation of steady spilling breaking waves. J Ship Res 2003;47:13–23. Toda Y, Stern F, Longo J. Mean-flow measurements in the boundary layer and wake and wave field of a serie 60 CB = 0.6 ship model for Froude numbers .16 and .316, IIHR Report, No. 352, 1991. Stern F, Longo J, Penna R, Olivieri A, Ratcliffe T, Coleman H. International collaboration on benchmark CFD validation data for surface combatant DTMB model 5415. In: Proceedings of 23th symposium on naval hydrodynamics, Val de Reuil, France, 2000. Dong RR, Katz J, Huang TT. On the structure of Bow waves on a ship model. J Fluid Mech 1997;347:77–115. Waniewski TA, Brennen CE, Raichlen F. Bow wave dynamics. J Ship Res 2002;46(1):1–15. Olivieri A, Pistani F, Avanzini G, Stern F, Penna R. Tank experiments of resistance, sinkage and trim, boundary layer, wake and free surface flow around a naval combatant INSEAN 2340 model, IIHR Report, No. 421, 2001. Stern F, Olivieri A, Shao J, Longo J, Ratcliffe T. Statistical approach for estimating intervals of certification or biases of facilities or measurement systems including uncertainties. J Fluid Eng 2005;127(3):604–10.