**3. Numerical formulation**

Equation (9) shows the conservative scheme used to solve the governing equations. The second order Lax-Wendroff method is used for the time integration, and the spatial derivatives are approximated using second order finite differences.

$$\mathbf{Q}\_{i,j}^{n+1} = \mathbf{Q}\_{i,j}^{n} - \Delta t \left[ \boldsymbol{\delta}\_{\xi} \boldsymbol{E}\_{\varepsilon} + \boldsymbol{\delta}\_{\eta} \boldsymbol{F}\_{\varepsilon} + \boldsymbol{\delta}\_{\zeta} \mathbf{G}\_{\varepsilon} - \boldsymbol{\delta}\_{\xi} \mathbf{G}\_{\varepsilon} \right. \tag{9}$$

$$\nabla\_{\xi} \boldsymbol{E}\_{\upsilon} - \nabla\_{\eta} \boldsymbol{F}\_{\upsilon} - \nabla\_{\zeta} \mathbf{G}\_{\upsilon} - \mathbf{}$$

$$\nabla\_{\xi} \boldsymbol{E}\_{RS} - \nabla\_{\eta} \boldsymbol{F}\_{RS} - \nabla\_{\zeta} \mathbf{G}\_{RS} - \boldsymbol{H} \right]^{n}$$

The Roe – Sweby fluxes are responsible for the upwinding and TVD of the scheme, see Eq. (10). For more details about the application of the upwind TVD scheme of Roe –Sweby to the slightly compressible Navier-Stokes equations, refer to Wanderley et al. (2008).

$$\left(E\_{RS}\right)\_{i+1/2,j} = \frac{1}{2}\tilde{A}\_{i+1/2,j}\left(Q\_{i+1,j}^{\eta} - Q\_{i,j}^{\eta}\right) \tag{10}$$

where

$$
\tilde{A} = T\tilde{\Lambda}T^{-1} \tag{11}
$$

A Three-Dimensional Numerical Simulation of the Free Surface Flow Around a Ship Hull 399

The coefficients *a*, *b*, *c* and *d* in Eq. (19) are obtained by the solution of the system of linear

1 11 1 2 22 2 3 33 3 0 00 0

 

A Cartesian computational mesh generator is used to discretize the three-dimensional physical domain. The grid is refined in the region close to the body to allow the capturing of viscous effects close to the body surface. An exponential stretching is used to concentrate grid points close to the body surface and to coarsening the computational mesh in regions away from the body, where the variation of the flow properties is small. The Cartesian computational mesh does not fit the body surface and the no-slip boundary condition is imposed using the Immersed Boundary Method (IBM). Fig. 2 shows a computational mesh plane at the level of the undisturbed free surface. The grid is divided in slices and each slice is superimposed to the neighbors such that the boundary of one side corresponds to interior points of the other. Each slice is solved by different nodes in the computer cluster using the message passing interface (MPI). Between successive iterations, the slices communicate to

*ax by cz d u ax by cz d u ax by cz d u ax by cz d u*

*u ax by cz d R R RR* (19)

(20)

where

Fig. 1. Reconstruction geometry.

where *u0*=0 is the no-slip condition on the body surface.

the neighbors for boundary condition update.

equations shown in Eq. (20).

**4. Grid generation** 

In Eq. (11), the matrix is a diagonal matrix with terms shown in Eq. (12) and the matrix *T* is defined in Eq. (13).

$$\tilde{\mathcal{A}}\_k = \left| \mathcal{J}\_k \right| + \psi\_k \left[ \frac{\Delta t}{\Delta \mathbf{x}} \left| \mathcal{J}\_k \right|^2 - \left| \mathcal{J}\_k \right| \right] \tag{12}$$

$$T = \left[ \vec{\mathbf{x}}\_1, \vec{\mathbf{x}}\_2, \vec{\mathbf{x}}\_3, \dots, \vec{\mathbf{x}}\_k \right] \tag{13}$$

where , *k k x* are the eigenvectors and eigenvalues of the Jacobian matrix *A* defined in Eq. (14).

$$A = \frac{\partial E\_e}{\partial Q} \tag{14}$$

The van Leer [5] flux limiter defined in Eq. (15) is a function of the coefficient *r* defined in Eq. (16) in the wave domain.

$$\nu\_k = \begin{cases} 0 & r\_k \le 0 \\ \frac{2r\_k}{1+r\_k} & r\_k > 0 \end{cases} \tag{15}$$

$$r\_k = \begin{cases} \frac{w\_{i+2,j}^k - w\_{i+1,j}^k}{w\_{i+1,j}^k - w\_{i,j}^k} & \lambda\_k \le 0 \\\\ \frac{w\_{i,j}^k - w\_{i-1,j}^k}{w\_{i+1,j}^k - w\_{i,j}^k} & \lambda\_k > 0 \end{cases} \tag{16}$$

where

$$
\vec{w} = T^{-1}\vec{Q}\tag{17}
$$

The no-slip boundary condition on the body surface is imposed using the Immersed Boundary Method (IBM), Marcelo et al. [4]. In this method, the flow field is informed of the body presence through a force field added to the momentum equations. These forces are computed such that the no-slip condition on the body surface is satisfied. The forces are applied only on the grid points close to the body surface through the reconstruction of the flow properties by linear interpolation. Fig. 1 illustrates the geometry of interpolation.

The properties of the flow field at the pyramid base (grid points 1, 2, and 3) are known because those are computed by the solver. The property at the pyramid pick (point 0 at the body surface) is imposed such that the no-slip condition on the body surface is satisfied. For example, the force added to the x-momentum equation (at point R) is computed as shown in Eq. (18). In this equation, *uR* is the reconstructed velocity and *u* is the velocity obtained by the solver in the iterative process.

$$F\_x = \rho \frac{\mu\_R - \mu}{\Delta t} \tag{18}$$

where

398 Computational Simulations and Applications

In Eq. (11), the matrix is a diagonal matrix with terms shown in Eq. (12) and the matrix *T*

*kkk k k t x*

> *Ee <sup>A</sup> Q*

The van Leer [5] flux limiter defined in Eq. (15) is a function of the coefficient *r* defined in

1

*r*

2, 1, 1, , , 1, 1, ,

*ij i j k k k i j ij*

The no-slip boundary condition on the body surface is imposed using the Immersed Boundary Method (IBM), Marcelo et al. [4]. In this method, the flow field is informed of the body presence through a force field added to the momentum equations. These forces are computed such that the no-slip condition on the body surface is satisfied. The forces are applied only on the grid points close to the body surface through the reconstruction of the flow properties by linear interpolation. Fig. 1 illustrates the geometry of interpolation. The properties of the flow field at the pyramid base (grid points 1, 2, and 3) are known because those are computed by the solver. The property at the pyramid pick (point 0 at the body surface) is imposed such that the no-slip condition on the body surface is satisfied. For example, the force added to the x-momentum equation (at point R) is computed as shown in Eq. (18). In this equation, *uR* is the reconstructed velocity and *u* is the velocity obtained by

> *<sup>R</sup> <sup>x</sup> u u <sup>F</sup> t*

*k k ij ij k k k i j ij*

*w w w w*

*w w w w*

*k k*

*k k k*

*r*

0 0 <sup>2</sup> <sup>0</sup>

*r*

*k*

*k k*

0

0

<sup>1</sup> *w TQ* (17)

(18)

*r r*

 2

are the eigenvectors and eigenvalues of the Jacobian matrix *A* defined in Eq.

 

(12)

*T xxx x* <sup>123</sup> ,,,, *<sup>k</sup>* (13)

(14)

(15)

(16)

 

is defined in Eq. (13).

where , *k k x*

Eq. (16) in the wave domain.

the solver in the iterative process.

(14).

where

$$
\mu\_R = a\chi\_R + by\_R + cz\_R + d\tag{19}
$$

#### Fig. 1. Reconstruction geometry.

The coefficients *a*, *b*, *c* and *d* in Eq. (19) are obtained by the solution of the system of linear equations shown in Eq. (20).

$$\begin{aligned} ax\_1 + by\_1 + cz\_1 + d &= \mu\_1 \\ ax\_2 + by\_2 + cz\_2 + d &= \mu\_2 \\ ax\_3 + by\_3 + cz\_3 + d &= \mu\_3 \\ ax\_0 + by\_0 + cz\_0 + d &= \mu\_0 \end{aligned} \tag{20}$$

where *u0*=0 is the no-slip condition on the body surface.

#### **4. Grid generation**

A Cartesian computational mesh generator is used to discretize the three-dimensional physical domain. The grid is refined in the region close to the body to allow the capturing of viscous effects close to the body surface. An exponential stretching is used to concentrate grid points close to the body surface and to coarsening the computational mesh in regions away from the body, where the variation of the flow properties is small. The Cartesian computational mesh does not fit the body surface and the no-slip boundary condition is imposed using the Immersed Boundary Method (IBM). Fig. 2 shows a computational mesh plane at the level of the undisturbed free surface. The grid is divided in slices and each slice is superimposed to the neighbors such that the boundary of one side corresponds to interior points of the other. Each slice is solved by different nodes in the computer cluster using the message passing interface (MPI). Between successive iterations, the slices communicate to the neighbors for boundary condition update.

A Three-Dimensional Numerical Simulation of the Free Surface Flow Around a Ship Hull 401

Figure 4 presents the temporal series of the total drag (in black), frictional drag (in blue), and pressure drag (in red) coefficients for the sphere for Reynolds number equal to 200. Table 1 presents a comparison between the total drag coefficient obtained in the present work and other experimental and numerical data obtained from the literature. The agreement among

Fig. 4. Frictional, pressure, and total drag coefficients for *Re*=200.

Schlichting [10] <sup>100</sup>

Campregher [2] <sup>100</sup>

Present study <sup>100</sup>

Table 1. Comparison of drag coefficient comparison.

**6. Flow field around a series-60 hull** 

**Reference** *Re Ct* **Comment**

1.100

1.178

1.200

Additional verification of the numerical code is presented for the Series-60 hull for a low Reynolds number of 1.0x103. Figures 5, 6, and 7 show a plane of the computational grid generated around the hull at the level of the undisturbed free surface. The grid has 240x160x160 points and the smallest element is a cube with dimension of 0.0075Lpp. The

0.800 Experimental

0.815 Numerical

0.832 Numerical

200

200

200

grid has 11Lpp in the y direction, and 7.2 Lpp in the *x* and *z* directions.

the three results is remarkable.

Fig. 2. Computational mesh.
