A new ghost-cell/level-set method for three-dimensional ﬂows

This study presents the implementation of a tailored ghost cell method in Hydro3D, an open-source large eddy simulation (LES) code for computational ﬂuid dynamics based on the ﬁnite diﬀerence method. The former model for studying the interaction between an immersed object and the ﬂuid ﬂow is the immersed boundary method (IBM) which has been validated for a wide range of Reynolds number ﬂows. However, it is challenging to ensure no-slip and zero gradient boundary conditions on the surface of an immersed body. In order to deal with this, a new sharp-interface ghost-cell method (GCM) is developed for Hydro3D. The code also employs a level-set method to capture the motion of the air-water interface and solves the spatially ﬁltered Navier-Stokes equations in a Cartesian staggered grid with the fractional step method. Both the new GCM and IBM are compared in a single numerical framework. They are applied to simulate benchmark cases in order to validate the numerical results, which mainly comprise single-phase ﬂow over inﬁnite circular and square cylinders for low-and high-Reynolds number ﬂows along with two-phase dam-break ﬂows with a vertical cylinder, in which a good agreement is obtained with other numerical studies and laboratory experiments.


Introduction
There is a need for a reliable method for modelling hydrodynamic forces in many coastal and marine engineering applications, such as coastal erosion, renewable energy, and wave mitigation strategies.The availability of increasing computational resources and the development of Computational Fluid Dynamics (CFD) has made it possible to solve the Navier-Stokes equations with complex geometries together with interface calculation methods for simulating such problems.Traditionally there are four categories of numerical solutions: Direct Numerical Simulations (DNS), Reynolds-Averaged Navier-Stokes (RANS), Large Eddy Simulations (LES), and Detached Eddy Simulations (DES) where a LES model is used in the bulk of the flow and a RANS model is applied near the solid body.DNS produces highly accurate results although at a high computational cost rendering almost impossible the option of applying this solution to practical coastal and ocean engineering studies.In LES simulations only eddies larger than the grid size are resolved while the dissipative effect of smaller eddies is taken into account by the sub-grid scale (SGS) model [1].
S. Martelo Lopez, A. Christou, S. Pan et al. 1.1.Overview of different methods for complex geometries in CFD Simulating fluid flow in the time domain with complex geometries is still a challenge.In the past boundary conforming methods have been adopted widely, such as body-fitted meshes or unstructured mesh methods.In body-fitted methods with overlapping grids, interpolation techniques are used to interpolate the solution from one grid to the other.An overset grid method is an example of a body-fitted method applicable to complex geometries with moving bodies.In an overset grid method, at least two mutually overlapping grids are used [2][3][4].For example, a structured body-fitted grid attached to a moving sphere is overset on a Cartesian grid [5].
Compared to body-fitted grid methods, the Cartesian grid methods typically require less memory and CPU time [6].A number of Cartesian grid methods have been published over the years.One of them is the immersed boundary method (IBM).This family of numerical schemes are an alternative to conventional body-conforming grid methods for modelling the interaction between a viscous fluid and a solid body.These methods permit the use of direct FFT (fast Fourier transform) solvers for the pressure and also have the advantage that they do not require a body-conforming grid of fluid cells at each time step (body-fitted methods).Thus they also decrease the computational time since it is no longer necessary to generate a new grid at every time step.Associated with this, these techniques also eliminate grid-interpolation errors.
There is a major difference between IBMs and body-fitted methods.In the latter, the local orientation of the grid conforms to the direction of the boundary of the fluid domain, which is not the case in IBMs [7].Mainly there are two subcategories for this technique: continuous (or diffuse) and discrete (sharp) methods [8].Additionally, depending on the way boundary conditions are imposed on the surface of the immersed body one can distinguish both continuous (diffuse) and discrete (sharp) immersed boundary methods (IBM), ghost-cell methods (GCM), cut-cell methods (CCM) and hybrid Cartesian/immersed boundary methods (HCIBM).

Continuous/direct forcing (or diffuse-interface) IBM methods
In this approach, a continuous forcing term is added to the Navier-Stokes equations before being discretized.Original formulations of this approach can be found in [9] where the author studied cardiac flows.With this technique, the main advantage is that these schemes are independent of the spatial discretization and can be inserted directly into a Navier-Stokes solver.The disadvantage, though, is that these methodologies create a diffuse boundary between the fluid and the solid.This means that the boundary condition on the immersed surface is not exactly enforced at the location of this interface but within a small neighbouring region.Here the velocity boundary condition was applied through a regularized delta function forcing term in the momentum equations.This interpolation function is a spreading operator that distributes the force from Lagrangian boundary points to Eulerian grid cells.It acts as a filter that smooths out the IBM force defined at the interface over a volume with a thickness of a few grid cells.
Another set of early examples is the ones [10,11] applied to elastic boundaries that are modelled through a smooth external forcing term added to the continuous momentum equation.This term also modifies a certain bounded volume surrounding the fluid-structure interface.Later, this was extended to simulate flow past rigid bodies, as noted, for example, in the works of [12] and [13], in which they modelled the forcing terms of the immersed boundaries by feedback techniques.In [14] rigid boundaries are approximated by highly stiff elastic boundaries.Additional examples can be found in the works of [15][16][17].The main advantage is that it can be easily implemented by adding source terms to the Navier-Stokes equations without major adjustment.On the other hand, the boundary conditions are diffused over the operational area of the forcing terms, which reduces numerical accuracy.
Uhlmann [18] noticed that certain sharp-interface direct forcing IBMs produced large pressure oscillations, and therefore he integrated the direct forcing technique into the regularized delta function formulation of Peskin [19].The resulting diffuse-interface IBM and variants of it have been applied to many computations of turbulent flows with many suspended particles (see for example [20,21]).

Discrete forcing (sharp-interface) IBM methods
The discrete (sharp) method aims to increase accuracy by modifying the discretization scheme of the Navier-Stokes equations or reconstructing the forcing terms in such a way that wall boundary conditions are satisfied creating a sharp interface between the fluid and solid boundary.Nonetheless, in this family of methods, performance is dependent on the spatial discretization scheme.In [22] and [23] it is proposed a reconstruction scheme for the solution at the fluid nodes in close proximity to the surface of the immersed body.In their work, the IBM force is defined as the term that is needed to make the velocity at the forcing points equal to the desired velocity.
This second-order accurate approach works well for bodies that are largely aligned with the grid line and so far has been applied to simulate a number of different flow situations with satisfactory results [6,24].Another example worth mentioning stemming from [23] is the algorithm published in [25] where pressure boundary conditions are included to help enforce mass conservation constraints through a least-square interpolation scheme on a collocated grid.However, it is ambiguous how to select the reconstruction direction, especially in the case of complex geometry.This problem can be tackled by interpolating along the normal line to the fluid-solid interface which was applied in the work of [26][27][28][29].Further bibliography where a forcing term is either explicitly or implicitly inserted into the Navier-Stokes equations can be found in [30][31][32][33][34][35][36].

Other Cartesian methods: ghost-cell, cut-cell and hybrid methods
It has also been experimented with reconstructing the solution on ghost cell nodes (fluid cells inside the immersed body) as can be seen in the work of [23].With this approach, a set of boundary conditions can be totally enforced and create a sharp interface between the fluid cells and the immersed body.Some publications that also employ the ghost-cell approach are [26,33,37,38].In [33] the proposed scheme uses a bilinear reconstruction procedure that is reduced to the one-dimensional linear one when there are no available points in the vicinity of the boundary to support the two-dimensional stencil.
One reference article is the work in [39] for simulating incompressible viscous flow past three-dimensional immersed bodies.It employs ghost cells for enforcing boundary conditions on the surface of the immersed body (immersed boundary).The Navier-Stokes equations are discretized using a cell-centred, collocated (non-staggered) arrangement of the primitive variables.In this scheme, the forcing term is absent in the momentum equations and the pressure and shear stress forces around the immersed body are calculated.The Navier-Stokes equations are advanced in time using the fractional step method.Other recent examples can be found in [40][41][42][43].
Cut-cell methods (CCM) form another sub-class of Cartesian grid methods that improve upon classical stair step methods.In cut-cell methods, Cartesian cells intersecting with the boundary are cut by the boundary, and the discrete mass and momentum conservation laws are also applied to the cut cells, see for example [44][45][46][47][48].
The algorithm presented in [26] is an example of a hybrid Cartesian/immersed boundary method.The approach presented in that article eliminates the previously discussed ambiguities associated with interpolation along grid lines but its applicability is restricted to flows with immersed boundaries that are aligned with one coordinate direction (e.g., two-dimensional or axisymmetric shapes).In such cases, the solution reconstruction is greatly simplified as it needs to be performed in two-dimensional planes.

The contribution and novelty of this study
The objective of the present study is to present a LES-based two-phase flow model with a further enhanced ghost-cell method, which will be referenced later as GCM, for simulating boundary force interactions between fluid and solid surfaces.The scheme for capturing the free surface is a version of the level-set method (LSM) using the re-initialization algorithm and the WENO scheme for computing field derivatives.The main purpose of this work is to present validation results from a range of numerical benchmarks for CFD solvers to assess the reliability of the newly developed GCM and compare it to the direct forcing IBM in the same numerical framework Hydro3D.While the LES, GCM and LSM have been presented separately in other numerical studies for flow simulations, they have rarely been combined.Implementing a GCM model in a finite difference code allows for fast and reliable simulations that evaluate hydrodynamic forces on structures.The current IBM presents a discrepancy between the predicted IBM force and the combined forces of pressure and shear stress.One of the key improvements of the GCM model is that it manages to predict a sharp interface between a solid body and the surrounding fluid.Additionally, other major enhancements are its robustness in high Reynolds flows and the convergence speed toward a steady-state solution.
Another important point to highlight is the capability of defining several layers of ghost points for several kinds of cross-sections.Thus the overall accuracy of pressure and velocity gradients is improved since higher-order differencing and interpolation schemes can be applied.Additionally, when using multiple layers of ghost cells, there is no need to define velocity ghost cell outside the fluid in order to guarantee that at least, the velocity ghost cell is located in the pressure cell face or halfway between two neighbouring pressure ghost cells [7].In contrast, the diffuse-interface IBM blurs the solid-fluid interface and it might suffer from reliability issues at high Reynolds numbers.
The new GCM presented here brings the following advantages to the table compared with diffuse-interface IBM schemes: (i) natural characterization of the boundary layer; (ii) improved accuracy and reliable readings of combined forces of pressure and shear stress compared with the former diffuseinterface IBM scheme; (iii) prediction of a sharp interface between a solid body and the surrounding fluid; (iv) higher convergence speed towards a steady state solution than the former IBM model; In addition, the new ghost-cell method presented here distinguishes itself from other sharp-interface IBM and ghost-cell schemes in: (i) use of a tailored optimization algorithm to find the image point; (ii) combination of a regularized delta function for interpolating velocities and a tailored hp-interpolating scheme for estimating pressure; (iii) several layers of ghost points are defined to increase the overall accuracy of pressure and velocity gradients; The remainder of the paper is organised as follows.The mathematical model and numerical implementation are described in Section 2. Next, benchmark validation results for laminar and turbulent single-and two-phase flows are shown and compared against other numerical simulations and laboratory experimental data in Section 3. Finally in Section 4 the main findings are summarized and a conclusion is included.

Numerical framework
In this section, the numerical implementation of the fluid solver is presented along with the turbulence, free-surface and hydrodynamic force models.

Fluid solver
The governing equations of the fluid model are the spatially-filtered Navier-Stokes equations: where equation ( 1 Eddy-viscosity (WALE) model [49].

Spatial discretization
The filtered Navier-Stokes equations are discretized following a finite difference approach implemented on a three-dimensional uniform Cartesian staggered grid for pressure  and velocity components   .As shown in Fig. 1, on a two-dimensional staggered grid, scalar quantities like pressure and viscosity are evaluated at the black circles, or in other words, the cell centres of the pressure grid depicted by the mesh of blue squares.In this way, their derivatives and velocities are calculated at the face centres of the blue squares.
The convective terms from (2) are approximated with a 4th-order central difference scheme (CDS) instead of using an upwind or blended scheme to avoid numerical dissipation.On the other hand, diffusive terms on the Navier-Stokes equation are discretised with a 2nd-order central difference scheme.Further details are available in [50].

Temporal discretization
The time advancing scheme chosen for Hydro3D-GCM is the fractional step method, first proposed by [51] and further developed by [52], which was implemented by [53] in Hydro3D.
In Fig. 1 the black circles represent pressure cell centres at the centre of the blue cells and the red squares at the centre of the vertical faces represent the cell centres for velocity component  1 , whereas the green triangles at the centre of the horizontal faces represent the cell centres for velocity component  3 .In the first fractional method, the convective and diffusive terms are estimated using the equation below: where  is the intermediate time step of the Runge-Kutta scheme.The time step, Δ is kept fixed and a safety parameter known as the Courant-Friedrichs-Lewy (CFL) coefficient is calculated at the end of each time iteration to check that the time step is small enough to ensure a good numerical solution.The CFL is evaluated with the equation below: )) (5) where   ,   and   are the maximum Cartesian velocities.
Finally, in order to obtain a divergence-free (solenoidal field) velocity vector and comply with the continuity equation, the Pressure Poisson Equation (PPE) is solved through a multi-grid algorithm.In [54] it is explained how such equation stems from the continuity equation and a pseudo-pressure,   , is obtained by solving: (6) which allows to calculate pressure and velocity corrections and update both fields before proceeding to the next time iteration in the fractional step method scheme as: In the present study, a 2:1 reduction in grid cell size between neighbouring sub-domains is imposed on the staggered computational grid to achieve local mesh refinement (LMR) in critical areas.The calculation of ghost cell pressures is achieved by adjusting the coarse pressure gradients at the coarse-fine interface to those computed for neighbouring fine cells, thereby coupling the pressure fields.The calculation of ghost velocities tangential and normal to the non-matching interface is also done in a consistent manner during the prolongation and restriction iteration of the multigrid algorithm for solving the pressure Poisson Equation.Further details can be found in [55].

Turbulence model
Since SGS tensor    depends on the magnitude of velocity gradients, then it is necessary to consider velocity gradients in different directions using the WALE model which captures the underlying physics of shear   = Tensor   denotes a rate-of-strain tensor of the resolved turbulent scale by the mesh and Ω  is the vorticity tensor of the resolved scale of the eddies as well.The proposed turbulent viscosity follows the equation below: where Δ  =   Δℎ is the sub-grid length scale.Here,    is defined as the traceless symmetric part of the square of the velocity gradient tensor: where  2  =     and   is the Kronecker symbol.For a good compromise coefficient   usually takes the value of 0.325 that corresponds to a Smagorinsky coefficient of   = 0.1.
The advantages of the WALE method over Smagorinsky-based methods are the independence from the distance to the wall and the WALE method also takes into account the local shear strain and rotation rates.Additionally, the eddy viscosity goes naturally to zero in the vicinity of the wall so neither an adjustment nor damping function is needed to compute wall-bounded flows.It is also worth noting that this model produces zero eddy viscosity in the case of pure shear.Thus it is possible to reproduce the laminar to turbulent transition through the growth of linear unstable modes.Furthermore, since       = 0 accounts for the rotation of the flow, this model gives more weight to this effect instead of shear since it creates more turbulence.For pure shear flows there is no rotation, so this term       = 0, capturing the behaviour for theoretical cubic decay of the eddy viscosity near the wall.Therefore this model does not generate turbulence in the laminar zone where there is only shear.Another convenient advantage is that WALE is invariant to any coordinate transformation of rotation and translation.Since only local information is needed, it is not required test-filtering operations and knowledge of the closest points in the mesh grid.Thus it is suitable for complex geometries.

Free-surface model
To model the effect of the water surface motion the level-set method (LSM) proposed by [56] and later improved in [57] was chosen for its simplicity.A detailed numerical implementation in Hydro3D of the algorithms is described in [58][59][60].The air-water interface moves with the fluid particles and a pure advection equation is used to model this: where variable  is a signed distance between the fluid particle and free surface or air-water interface. is evaluated at the cell centres of the pressure grid: ∀ = 1, 2, 3, ...,  1 ; ∀ = 1, 2, 3, ...,  2 ; ∀ = 1, 2, 3, ...,  3 where integers  1 ,  2 ,  3 are the number of cells in the direction of the Cartesian axes  1 ,  2 ,  3 .Field variable  is positive for the water phase and negative for the air phase.
The spatial derivatives of field variable  are evaluated through a high-order accurate scheme that avoids spurious oscillations using the weighted non-oscillatory scheme (WENO) [61].This scheme is 5th-order accurate for polynomial functions and employs a 4 cell stencil.If an upwind advection scheme is employed, its cell centred derivatives are: ∀ = 1, 2, 3 where represent the left and right biased stencil, respectively.Additional guidance on how to evaluate these derivatives can be found in [50].A third-order Runge-Kutta scheme is used to march in time advection equation (11): Here, all velocities and derivatives   ,    (∀ = 1, 2, 3) are evaluated at the cell centre coordinates   of the pressure grid.
A Heaviside function is used to avoid discontinuities in fluid properties.This technique allows for a smooth exchange of properties between fluid phases within a transition region of width 2 = 4Δℎ where Δℎ is the maximum grid spacing.The Heaviside function is formulated as: Thus, physical phase properties vary according to the value of the Heaviside function: where subscripts ,  indicate liquid and gas phases.Both phases are assumed to be a continuous medium treated as a fluid.Therefore, pressure and velocities at the interface and transition region are solved with the algorithms discussed in previous sections.In order to ensure convergence the  gradient norm is kept as close as possible to 1 by solving the re-initialization equation within the transition region [62]   * + where  0 = (,  = 0) and  * is an artificial time determined by the grid size multiplied by a safety factor of value less than one, similar to the CFL coefficient.Another term in the previous equation,  (  0 ) , is known as the smoothed signed function, which is defined as: This re-initialization is applied throughout the transition zone within several iterations steps where   represents one grid spacing.This technique is only applied in computational cells lying on the interface, thus there is no need to solve this partial differential equation for the whole domain.

Hydrodynamic forces model
This section explains the principal ideas of a sharp-interface model accounting for the interaction of fluid forces on submerged solid bodies.Such model is based on the ghost-cell method and has been implemented into Hydro3D.In the former IBM scheme, the diffuse immersed boundary forcing term is replaced by an implicit force implied by interpolation of velocity and pressure to Eulerian fluid cells inside the immersed object in order to impose no-slip boundary conditions.These fluid cells inside the solid where the boundary conditions are applied are also known as ghost cells.The resulting effect of this implicit force creates a sharp interface between the fluid and the solid.The total force on the body is evaluated with a surface integral.A pseudo-code describing the steps to be taken at each time step is presented below and a flow chart can be found in Fig. 2. The steps of this approach can be summarized as: (c) Apply equations ( 7) and ( 8), and evaluate new velocity and pressure fields    for the current time iteration  using the updated pressure gradient.(d) Apply no-slip wall boundary conditions at the fluid cells inside the solid body and close to its surface using the newly developed ghost-cell technique mentioned in 2.4.1, similar to [23,27,39].
(e) Calculate force F on the surface  of the solid body with details in 2.4.2 as:

Ghost cell treatment
Among other modifications, an algorithm was coded to automatically select a variable number of ghost layers inside the immersed object sweeping a path defined by the outermost immersed markers of its cross-section.The main idea consists of selecting a set of ghost cells, that is, the Eulerian fluid cells inside the immersed body.After that, the mirror point with respect to the body surface is calculated.These ghost cells will need to be assigned   flow fields that cancel the respective mirror field values if one wants to apply a no-slip boundary condition for the velocities on the surface of a fixed body as described by this relationship: Where upper index ℎ denotes field values in the ghost cell (fluid Eulerian cell inside the immersed body) and  indicates the respective field values for the mirror image point  in Fig. 3.In order to find values of    , fluid velocity at the mirror point location, a delta function interpolation scheme is applied: where  is the delta function operator and Δℎ ,, = 3 √       is the homogeneous grid spacing in the three Cartesian directions.

⃗ 𝑥 𝑖𝑔,𝑗𝑔,𝑘𝑔 𝑖,𝑗,𝑘
and ⃗   ,, are the coordinates of the fluid cells and mirror point  in Fig. 3, respectively.For a certain stencil of a delta function, the sets   ,   ,   of cell centres along the Cartesian axis are defined, with each set having   ,   ,   number of cell centres along the horizontal, transversal and vertical Cartesian directions, respectively: Here, the set of indexes   ,   ,   denote the address of mirror point  in the  − ℎ respective Cartesian grid; while Θ is half of the grid cells in which a support interval is defined for the selected delta function.
In Fig. 3 the black hollow circles depict pressure ghost cell centres at the centre of the blue cells.The red hollow squares at the centre of the vertical faces represent the ghost cell centres for velocity component  1 .The green hollow triangles at the centre of the horizontal faces represent the ghost cell centres for velocity component  3 .
The location of each marker is given by the indexes of the cell faces downstream.This allows for easy implementation of the ghost-cell method since each marker is associated with the closest pressure and velocity cell faces.On the other hand, in order to apply zero-gradient boundary conditions on the body surface, the ghost cell value for pressure will be equal to the pressure magnitude at the mirror point: where index  refers to point b in Fig. 4. In order to estimate the value of pressure at the mirror point a tailored hp interpolation scheme was chosen and it was applied to a stencil  = , , ,  of  = 4 cell centres, equally spaced forming a square, as it is shown in Fig. 4. Thus, the equation below shows the general equation to estimate pressure values at the mirror point.(26) where the coefficients   and   have the following expressions: with pair (  1 ,   3 ) being the location in the -plane of mirror point  while intervals [ 1 ,  2 ] and [ 1 ,  2 ] define the bounding box of the 4 point stencil along the horizontal and vertical directions, as shown in Fig. 4 by the set of  points  = , , ,  .

Force calculation in GCM
Further details of how to evaluate the surface integral in the force calculation ( 21) can be found in Figs. 5 and 6.The pressure force along the  axis acting on an elemental surface spanning markers  −1 and   , separated a distance   , reads: Trapezoidal rule: = 0.5 Simpson's rule: The pressure force along the  axis acting on an elemental surface spanning markers  −1 and   , separated a distance   , reads: Trapezoidal rule: = 0.5  Simpson's rule: In the equations above, Δ is the spacing between the cross-sections along the span-wise direction,  while angle  is the slope between the markers   and  −1 .
In Fig. 6 vectors t and n are the tangent and normal unit vectors to the surface of the immersed body.This work follows the convention of the normal vector pointing outside of the immersed body surface.
The shear stress force along the  axis acting on an elemental surface spanning markers  −1 and   , separated a distance   , reads: Trapezoidal rule: = 0.5 Simpson's rule: The shear stress force along the  axis acting on an elemental surface spanning markers  −1 and   , separated a distance   , reads: Trapezoidal rule: = 0.5 Simpson's rule: Here,  is defined as: where variable  is the dynamic viscosity.Using a second-order approximation for the derivatives we have this relationship: for obtaining the equation above it was assumed that the velocity fields  and  inside the immersed take mirror values from the fluid cells outside the body due to the ghost cell method.This implies: where û and ŵ are the velocity values inside the body that mirror the respective field values outside.

Comparison with other IBMs
The new enhanced GCM presented here brings the following advantages compared with diffusive IBM schemes: • natural characterization of the boundary layer and therefore reliable assessment of combined forces of pressure and shear stress; • prediction of a sharp interface between a solid body and the surrounding fluid; Additionally the new ghost cell method presented here distinguishes itself from other sharp IBM schemes in: • use of a tailored optimization algorithm to find the image point resembling the means scheme [63]; • use of a regularized delta function for interpolating velocities at the image point while a hp interpolating scheme is applied to the pressure field; • use of bilinear interpolation when the previous schemes yield values above (below) the local maxima (minima) of the interpolation stencil; • multiple layers of ghost points are defined increasing the overall accuracy of pressure and velocity gradients;

Numerical validation benchmarks
This section is dedicated to explain in detail how the new GCM model implemented in Hydro3D is validated.Looking at previous works, a few benchmarks are selected to ensure reliable predictions of drag and lift forces for both low-and high-Reynolds number single-and two-phase flows.It is worth mentioning that both the ghost-cell method and the IBM are implemented in the same numerical framework so that we can assess the performance of both methods.The predicted results are compared with available experimental measurements and other numerical results to demonstrate the capability of the proposed method.

Fixed circular cylinder in low Reynolds number flows, RE = 40
For this simulation a circular cylinder is facing a single-phase flow with a Reynolds number  =   = 40, in which  is the inlet velocity,  is the diameter of the cylinder and  is the kinematic viscosity.In Fig. 7 below is the sketch depicting the setup of the computational domain and the location of the cylinder.The finest mesh resolution is  Δℎ = 40, or in other words, there are 40 cells across the cylinder diameter.Additionally, the time step chosen is Δ = 0.005 s to guarantee that the total CFL number will not be above 1 or cause numerical instabilities.The top and bottom surfaces of the computational domain are treated as no-slip walls, which run parallel to the fluid flow and are far away from the body.On the other hand, for the north and south boundaries (frontal and rear sides) of the domain (they are the largest boundaries, which are perpendicular to the cylinder axis) periodic boundary conditions have been applied.The west and east boundaries are modelled as an inlet and outlet, respectively.For the inlet Dirichlet conditions are applied since the velocity is prescribed as  = 0.4 m/s with a kinematic viscosity equal to  = 1 100 m 2 /s.For the outlet, a convective boundary condition is applied as

Table 1
Drag and wake length for a circular cylinder at  = 40.where    is the fluid velocity at the outlet,  ℎ  is the velocity in the first ghost layer of the outlet and  ∞ is the mean stream-wise velocity at the plane where the outlet is located.Fig. 7 shows regions A and B. In region A the grid spacing is Δℎ = 0.05 m and the mesh resolution is  Δℎ = 20 while in region B the grid spacing is Δℎ = 0.025 m and the mesh resolution is  Δℎ = 40.Table 1 summarizes published results by other researchers along with the results obtained using the IBM and ghost-cell method.

𝐶
Mainly, Table 1 compares drag and wake length for a circular cylinder for  = 40 between authors from the literature, the standard IBM and the ghost cell method for a mesh resolution of  Δℎ = 40 and time step Δ = 0.005 s.For such a low value of Reynolds number, there is no vortex shedding downstream, so the lift coefficients are negligible.Also, the prediction of the drag values from both the IBM and GCM agrees with previous numerical results [32,33,38,64].One can also see in Table 1 that the wake length downstream of the cylinder predicted by the GCM also agrees with the data in the literature.
In Fig. 8 the reader can find a comparison between the predicted pressure values around a cylinder with the GCM and the numerical results obtained by [38].Fig. 9 also compares the numerical results of the skin friction coefficient from the same author against the GCM.

Fixed circular cylinder in low Reynolds number flows, RE = 100
For this simulation, a circular cylinder is facing a single-phase flow with a Reynolds number of 100.The same computational setup shown in Fig. 7 is used here, in which region A the grid spacing is Δℎ = 0.05 m and the mesh resolution is  Δℎ = 20 while in   The blue line is the immersed boundary method and the red line indicates the ghost cell method both using the phi2 delta function as an interpolation scheme for velocities and the hp interpolation scheme for pressure.
region B the grid spacing is Δℎ = 0.025 m and the mesh resolution is  Δℎ = 40.In this way, there are 40 cells across the cylinder diameter.Additionally, the time step chosen at first was 0.005 s to guarantee that the total CFL number will not be above 1 and is low enough to not cause numerical instabilities.The boundary conditions are the same as in the previous case.
The results for the drag coefficient value predicted by other authors are around 1.35.In Fig. 10 one can find that the average value for the drag coefficient achieved with the proposed GCM is 1.33 using delta functions for interpolating velocities and a hp scheme for interpolating pressure values.When comparing this against the IBM model with a spreading operator that modifies the boundary layer around the cylinder, the average drag value is   = 1.43, considerably larger and less accurate.One reason for this might be due to the fact that the IBM creates a blurry interface between the fluid and the solid, modifying the boundary layer that ultimately determines the strength of the fluid forces interacting with the circular cylinder.Fig. 11 compares the lift predictions which reasonably agree with data found in published literature.One interesting feature when observing both drag and lift predictions in Figs. 10 and 11 is that the GCM reaches the steady state much faster than the IBM method.Thus, when it comes to moving structures the GCM might be more reliable.Figs. 12 and 13 compare the total hydrodynamic force computed with the GCM and the pressure and skin friction components of these forces.
Table 2 shows the comparison of average drag, root mean square of lift and Strouhal number coefficient predictions for a circular cylinder at  = 100 between different methods available in the literature.The predicted results obtained from the standard IBM and the ghost cell method are for a mesh resolution of  Δℎ = 40 and time step Δ = 0.005 s.Tables 3 and 4 shows the values for drag, lift and vortex-shedding coefficients for several time steps: Δ = 0.005, 0.0025, 0.00125 s and mesh resolutions: ∕Δℎ = 40, 60, 80.With this information, a temporal and spatial convergence study was performed.For the time convergence case, the mesh resolution was kept at ∕Δℎ = 40.For this case, a nearly first-order trend was found for the drag values.It is worth noticing that when the time step is halved by a second time, the change in the value of the average drag coefficient is multiplied by a factor smaller than 1.This fact reflects that as the time step decreases, the error also decreases for the given mesh, The blue line is the immersed boundary method and the red line indicates the ghost cell method both using the phi2 delta function as an interpolation scheme for velocities and the hp interpolation scheme for pressure.Represented in the red line is the total force coefficient from GCM while the green and black lines indicate the pressure and skin friction components of the total force, respectively.

Table 3
Comparison of predictions from the ghost cell method for a circular cylinder and  = 100 with a mesh resolution of which is an indicator that the simulation is converging in time (Fig. 14).On the other hand, when the spatial convergence study was carried out at a constant time-step, Δ = 0.00125 s, looking again at the red line in Fig. 15 the spatial convergence for the GCM was slightly above first-order.This means that numerical results converge faster when decreasing fluid cell size than when decreasing the time step.It was noted that a first-order spatial convergence is observed although a higher-order scheme is used for the advection.This degradation of convergence has also been reported in recent works of IBM for high-order methods [72], which depends on the regularity of the solution across the fluid-solid interface.In order to justify the selection of the domain height, an additional sensitivity analysis was carried out.Table 5 briefly presents the main variations in the force coefficients.
In Fig. 16 the contour plot shows the horizontal velocity field outside and inside the cylinder when using the GCM that employs a delta function interpolation for velocities and hp shape function scheme for pressure interpolation.The mesh resolution is  Δℎ = 40 and the time step used is Δ = 0.005 s.Highlighted in deep blue downstream the cylinder is the wake where negative velocities and re-circulation effects can be observed.Also, depicted in deep blue the frontal cells inside the cylinder show how the ghost layers take the opposite values of the velocities in the fluid outside the cylinder.On the other hand, the rear cells inside the cylinder only have a light blue shade because the fluid velocity outside is close to zero and the interpolation using delta functions predicts conservative values for the mirror velocities.Additionally, in Figs.17a and 17b it can be seen the average horizontal and vertical velocities indicate that the fluid solution is rather symmetric.This can also be seen by looking at the cross-flow shear stress in Fig. 17c.
Figs. 18 and 19 show the comparison of the IBM and GCM employing several interpolations schemes for predictions.Both IBM and GCM use delta functions for interpolating velocities.The most accurate drag prediction with the GCM is the one with the hp shape function for estimating pressure in the ghost cells, shown in red colour, closely followed by bilinear least squares.When using the delta function the pressure force and overall fluid force acting on the surface of the cylinder are significantly underestimated.Both schemes of the GCM for interpolating pressure yield similar maximum lift coefficients.The bilinear scheme tends to slightly    predict higher values than when selecting the hp shape function.On the other hand, when choosing the delta function the maximum lift coefficient is lower than the one predicted by the hp scheme.

Fixed squared cylinder in low Reynolds number flows, RE = 100
A square cylinder in a single-phase flow with a Reynolds number of 100 is considered here.Fig. 20 shows the computational setup along with regions A and B, in which region A the grid spacing is Δℎ = 0.05 m and the mesh resolution is  Δℎ = 20 while in region B the grid spacing is Δℎ = 0.025 m and the mesh resolution is  Δℎ = 40.For this benchmark a time step Δ = 0.005 s was tested using both the hp interpolation scheme for pressure and the delta function for interpolating velocities.Figs.21 and 22 show the predicted time series of the drag and lift coefficient values obtained from both GCM and IBM.The results are quite close to the value of   = 1.5 found in the literature for the proposed GCM while the standard IBM yields much higher values for the drag and lift coefficients.As in the previous case, the GCM also reaches the steady state faster.The results for average drag, root mean square lift and Strouhal number are summarized in Table 6.Fig. 23 presents a snapshot of the streamwise velocity field when using the ghost-cell method that employs a delta function interpolation method for velocities and hp shape function scheme for pressure interpolation.It can be seen that the ghost-cell method is able to deal with the sharp interface of the square cylinder and both the vortex shedding in the wake and the flow field around the square cylinder corner are well captured.Tables 7 and 8 show the values for drag, lift and vortex-shedding coefficients for several time steps: Δ = 0.005, 0.0025, 0.00125 s and mesh resolutions: ∕Δℎ = 40, 60, 80.With this information, a temporal and spatial convergence study was performed.For the time convergence case, the mesh resolution was kept at ∕Δℎ = 40.For this case, a first-order trend was found for the drag values.In Fig. 24 the red line shows a rate of convergence of 1 for the drag coefficient of the square cylinder using the GCM.The black dash and dash-dot lines show theoretical convergence rates of 1st and 2nd order, respectively.In a similar manner, when the spatial convergence study was carried out at a constant time-step, Δ = 0.00125 s, looking again at the red line in Fig. 25 the spatial convergence for the GCM was first-order.26.  curve for a circular cylinder facing a flow of  = 3900 using the ghost-cell method with a delta function interpolation method for velocities and hp shape function scheme for pressure interpolation for time step Δ = 0.005 s and mesh grid resolution  Δℎ = 40.

Fixed circular cylinder in turbulent flows, RE = 3900
In this case, a circular cylinder with a width of  = 4 m is facing a single-phase flow with a Reynolds Number,  = 3900.The size of the computational domain and the location of the cylinder follow the same settings as in Fig. 7.The grid size and time step chosen were  Δℎ = 40 and Δ = 0.005 s.The inlet Dirichlet condition is applied on the west boundary of the computational domain, since the velocity is prescribed as  = 1 m/s with a kinematic viscosity equal to  = 1 3900 m 2 /s and a density  = 1 kg/m 3 .The east boundary is set as an outlet with zero gradient conditions applied to the pressure and velocity fields.Top and bottom boundaries are considered as no-slip walls while in the south and north lateral boundaries, periodic conditions are enforced.The flow at such Reynolds number is categorized under the lower sub-critical range of flow where the bulk of the flow remains laminar beyond separation and the transition to turbulent flow takes place in the free shear layer in the wake of the cylinder where turbulent eddies are shed periodically.As a result, this case presents different scales simultaneously, making a numerical simulation of this flow very challenging.With increasing Reynolds number, the three-dimensional wake behind the immersed cylinder becomes more chaotic.
In order to capture turbulent structures an algorithm using the -criteria [84] was coded and the computational domain was defined with a width  = 4 m along the OY axis in order to capture three-dimensional turbulence effects.Thus the frontal area to consider is  =  = 4 m 2 .To simulate an infinite cylinder, periodic boundary conditions were employed in the spanwise direction.
In Figs. 26 and 27 one can visualize the evolution of drag and lift force coefficients.The drag force calculated in Hydro3D with the ghost-cell method takes into account the pressure and shear stress around the immersed object.Such force is divided by in order to obtain a dimensionless drag coefficient   .Here  denotes density,  is the inlet velocity and  is the cylinder frontal

Table 9
Comparison of drag and lift coefficients predictions for a circular cylinder at  = 3900 between authors from the literature, the standard IBM and the ghost cell method for a mesh resolution of  Δℎ = 40 and time step Δ = 0.005 s.  9 compares results from a comprehensive list of numerical and physical experiments against the IBM and GCM schemes.When it comes to the GCM, the steady-state region compares well with other predicted values found in the literature, such as in [97] where   = 1.1.While the GCM performs reasonably well, the IBM shows a clear deviation from the expected values in the literature for the drag force value.

C𝑑
In Fig. 28 features a snapshot of the horizontal velocity  1 where an alley of vortices is present downstream of the cylinder.In this figure, one can also check that the GCM is working properly inside the cylinder since the first ghost layers the velocity has an opposite sign to the velocity of the fluid cells outside the cylinder.The version of the GCM uses a delta function interpolation scheme for velocities and a tailored hp interpolation method for pressure.On the other hand, the velocity field predicted by the IBM is presented in Fig. 29.This picture shows that inside the cylinder the velocity has a homogeneous pattern close to zero in value.This leads to subsequent errors when evaluating velocity gradients during the stage for predicting the effects of the convective and diffusive terms in the fluid cells close to the solid surface.
Also related to the velocity field predicted by the GCM, Fig. 30 shows coherent turbulent structures using the -criteria for a circular cylinder facing a flow at  = 3900.
The mean and turbulence statistics of the flow obtained from the ghost-cell method are analyzed and presented in Figs.31-36.The results were calculated over a time window spanning  = ∕ ∞ = 50 − 82 s, a period that corresponds with the end of the transient phase, at  = 50 s and subsequent 7 full oscillations of the lift coefficient   .The time step was kept fixed at a value of 0.005 s to ensure a CFL under 0.7.
Mean velocities and turbulent Reynolds stresses from simulations using the GCM are compared against numerical results of [98] and also against experimental results published by [86] and [87].In [86] hot-wire measurements were recorded in the down-stream region spanning from ∕ = 3 to ∕ = 10 while in [86] the near wake was studied.On the other hand, the numerical experiments    from [98] come from a B-spline-based method.Overall, in the figures one can see that the predictions of turbulence statistics from the GCM are in reasonable agreement with both other experimental and numerical studies.In general, the numerical results show a  -shaped patterns for the stream-wise velocity in Fig. 32, while the profile observed in the experiments resembles a  -shape instead inside the region close to the cylinder, at ∕ = 1.06,where flow re-circulation is significant.However, the present simulations do show a  -shape profile in the re-circulation zone edge at ∕ = 1.54.But the size of the re-circulation zone is smaller in the physical experiments.The largest differences between the GCM and the numerical experiments of [98] are in the down-stream region between ∕ = 6 and ∕ = 10, where there is a coarser mesh in the present simulation.

3D dam-break flow with a vertical barrier
The next set of numerical results depicts another case to study the collapse of a three-dimensional column of water and subsequent interaction with a vertical square cylinder with an edge length  = 0.12 m.The computational set-up is defined in Fig. 37 following the experiments and simulations from [99] and all the boundaries are treated as no-slip walls except the top surface that is modelled as an outlet to the atmosphere.Also, the column of water is held in place by a virtual barrier that vanishes during the first time step of the simulation.The mesh resolution is ∕Δℎ = 24, the time step is 5 × 10 −5 s and the  parameter for the level-set method is set to 0.5Δℎ, which means that the transition from the water to the air phase is within one grid spacing Δℎ.
In this study, the force acting on the cylinder is plotted in Fig. 38 and compared with the laboratory experiments conducted at the University of Washington [99].Overall, the simulation results agree with the physical experiments, but the force predicted by the ghost-cell method shows some minor oscillations during the water column impact on the square cylinder.As in the previous turbulent benchmarks, the sub-grid scale model chosen is the Wall Adaptive Local Eddy-Viscosity, aka WALE, with coefficient   = 0.325 and no major change was observed when changing this parameter.However, it was observed that decreasing the  parameter, which    governs the air-water transition in the level-set method, helps to mitigate these force oscillations.Additionally, Fig. 39 shows results of the predicted horizontal velocity for a probe located at a point  = 0.754 m from the west wall and at a height of 0.024 m from the bottom wall.Here the pattern of numerical results reasonably agrees with the physical experiments although there is some offset.The reason for this may be due to the pressure solver and the LSM do not manage to achieve complete convergence of the numerical solution when the maximum number of iterations is reached.Another source of error might be the fact that the GCM is not fully mass-conservative.This last source of error is heavily influenced by the mesh resolution.Thus, it may be possible to obtain better simulations for the velocity field by decreasing the cell size.Unfortunately, by doing this, the time step also needs to be decreased, which leads to an increase in the computational resources needed to perform this kind of simulation with such a small cross-sectional area of the immersed solid.
In Fig. 40 one can see some snapshots of the air-water interface for several time steps.The interface is represented by a contour plot showing an iso-surface where the fluid density is 500 kg/m 3 .As time passes by it is easily noticeable how the column of water collapses and impacts against the vertical square cylinder.

Conclusion
In this work, a new ghost-cell model has been developed and validated for single-and two-phase flows.The novelties introduced in Hydro3D are high-fidelity simulations by means of employing a finite difference open-source code that uses both delta functions and   hp interpolation schemes for enforcing the boundary conditions of the GCM.It also has the capability of handling large geometries of solid bodies defined across several subdomains in parallel computing.
It is worth mentioning that both sharp and diffuse interface methods for complex geometries are compared in a single numerical framework.When comparing the new model implemented in Hydro3D, based on the GCM, against the former IBM model, both methods predict similar drag force coefficient values for low Reynolds numbers.However, for high Reynolds number flow the IBM starts to over-predict the drag force coefficient.Additionally, the GCM tends to reach a steady state much faster than the IBM method.Thus, when it comes to moving structures the GCM may be more reliable.One reason for these effects may be related to the fact that the IBM creates a diffuse interface between the fluid and the solid, modifying the boundary layer that ultimately determines the strength of the fluid forces interacting with the circular cylinder.
The GCM has also been applied for two-phase dam-break flows with a vertical cylinder.It was observed a dependency of small drag force oscillations with the value of the  parameter which governs the air-water transition in the level-set method.When decreasing  the force oscillations also decrease.Another limitation of the proposed GCM is that the velocity field near the immersed solid starts to show discrepancies with physical experiments when the mesh resolution is small.This typically implies that the crosssection of the solid is small in size compared with the mesh cell size.The reason for this source of error may be that the pressure solver and the LSM do not manage to achieve complete convergence of the numerical solution.Another explanation is the fact that the GCM is not fully mass conservative although it is possible to circumvent this issue by decreasing the cell size.It has been shown that Hydro3D with enhanced GCM can be used as a design tool to assess turbulent structures and load forces on both immersed and piercing structures.Future work will be focused on validating the GCM for moving objects in turbulent multiphase flows.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.In order to assemble a shape function, these polynomials need to have certain properties:

Fig. 1 .
Fig. 1.Cartesian grid representing the discretization of a continuous fluid domain into Eulerian Cells.(For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.) ) is the filtered continuity equation and equation(2) depicts the filtered momentum equation, both written in tensor notation.Velocity component   (∀ = 1, 2, 3) denotes the three-dimensional Cartesian velocity components at their respective cell face centre acting on the direction of the th Cartesian axis. is the filtered pressure, and  and  are the filtered fluid density and kinematic viscosity.Finally,   denotes the Cartesian component of gravitational acceleration and the filtered term  -grid scale stresses that have been calculated taking into account the Wall Adaptive Local
(a) Calculate an intermediate velocity  *  for each Eulerian fluid cell using information from the previous time step through equation (3).(b) Ensure divergence-free velocity field by solving the pressure Poisson equation (6).

Fig. 3 .
Fig. 3. Cartesian grid representing the interface (yellow line) between Eulerian cells of a single-phase fluid (filled circles, triangles and squares) and Lagrangian ghost cell of the immersed body (hollow circles, triangles and squares).

Fig. 7 .
Fig. 7. Schematic of the computational setup with local mesh refinement regions A and B. The dimensions are given in terms of the cylinder diameter .

Fig. 8 .
Fig. 8. Pressure coefficient comparison around a circular cylinder facing a single-phase flow with  = 40.

Fig. 9 .
Fig. 9. Skin-friction coefficient comparison around a circular cylinder facing a single-phase flow with  = 40.

Fig. 10 .
Fig. 10.Comparison for prediction of drag coefficient,   , for a circular cylinder and  = 100 with a mesh resolution of  Δℎ = 40 and a time step of Δ = 0.005 s.

Fig. 11 .
Fig. 11.Comparison for prediction of lift coefficient,   , for a circular cylinder and  = 100 with a mesh resolution of  Δℎ = 40 and a time step of Δ = 0.005 s.

Fig. 12 .
Fig. 12.Comparison for prediction of drag coefficient,   , for a circular cylinder and  = 100 with a mesh resolution of  Δℎ = 40 and a time step of Δ = 0.005 s.Represented in the red line is the total force coefficient from GCM while the green and black lines indicate the pressure and skin friction components of the total force, respectively.

Fig. 13 .
Fig. 13.Comparison for prediction of lift coefficient,   , for a circular cylinder and  = 100 with a mesh resolution of  Δℎ = 40 and a time step of Δ = 0.005 s.

Fig. 15 .
Fig. 15.Spatial convergence study for a circular cylinder facing a single-phase flow of  = 100.

Fig. 16 .
Fig. 16. 1 velocity profile for simulation time 200 s of a circular cylinder facing a flow at  = 100.

Fig. 17
Fig. 17.The time-averaged fluid fields of a circular cylinder facing a flow at  = 100.a) Time-averaged horizontal velocity.b) Time-averaged vertical velocity.

.
Fig. 17.The time-averaged fluid fields of a circular cylinder facing a flow at  = 100.a) Time-averaged horizontal velocity.b) Time-averaged vertical velocity.c) Time-averaged cross-flow shear stress.

Fig. 18 .
Fig. 18.Comparison for prediction of drag coefficient,   , for a circular cylinder and  = 100 with a mesh resolution of  Δℎ = 40 and a time step of Δ = 0.005 s.In a blue line, it is depicted the IBM prediction.The other colours indicate the GCM results using a delta function for interpolating velocities and several schemes for interpolating pressure: bilinear-least-squares in green, hp shape function in red and delta function in black.

Fig. 19 .
Fig. 19.Comparison for prediction of lift coefficient,   , for a circular cylinder and  = 100 with a mesh resolution of  Δℎ = 40 and a time step of Δ = 0.005 s.A blue line depicts the IBM prediction.The other colours indicate the GCM results using a delta function for interpolating velocities and several schemes for interpolating pressure: bilinear-least-squares in green, hp shape function in red and delta function in black.

Fig. 20 .
Fig. 20.Schematic of the computational setup with local mesh refinement regions A and B. The dimensions are given in terms of the square cylinder width .

Fig. 21 .
Fig. 21.In blue: drag coefficient curves for the IBM scheme.In red: the ghost cell method using delta functions for velocity interpolation and a hp scheme for pressure.The mesh resolution is  Δℎ = 40, the time step used is Δ = 0.005 s and  = 100.

Fig. 22 . 7
Fig. 22.In blue: lift coefficient curve for the IBM scheme.In red: the ghost cell method using delta functions for velocity interpolation and an hp scheme for pressure.The mesh resolution is  Δℎ = 40, the time step used is Δ = 0.005 s and  = 100.

Fig. 25 .
Fig. 25.Spatial convergence study for a square cylinder facing a single-phase flow of  = 100.

Fig.
Fig.26.  curve for a circular cylinder facing a flow of  = 3900 using the ghost-cell method with a delta function interpolation method for velocities and hp shape function scheme for pressure interpolation for time step Δ = 0.005 s and mesh grid resolution  Δℎ = 40.

Fig. 27 .
Fig.27.  curve for a circular cylinder facing a flow of  = 3900 using the ghost-cell method with delta function interpolation method for velocities and hp shape function scheme for pressure interpolation for time step Δ = 0.005 s and mesh grid resolution  Δℎ = 40.

Fig. 28 .
Fig. 28. 1 velocity profile for a circular cylinder facing a flow at  = 3900.The version of the ghost-cell method uses a delta function interpolation method for velocities and a hp shape function scheme for pressure interpolation.The mesh resolution is  Δℎ = 40 and the time step used is Δ = 0.005 s.

Fig. 29 .Fig. 30 .
Fig. 29. 1 velocity profile for a circular cylinder facing a flow at  = 3900 predicted by a diffuse interface IBM scheme using a delta function interpolation for velocities.The mesh resolution is  Δℎ = 40 and the time step used is Δ = 0.005 s.

Fig. 31 .
Fig. 31.Velocity profile normalized with the inlet velocity on the centre line of the wake region behind a circular cylinder at  = 3900.

Fig. 32 .
Fig. 32.Streamwise velocity profile normalized with the inlet velocity at 3 checkpoints on the re-circulation region behind a circular cylinder at  = 3900.

Fig. 33 .
Fig. 33.Streamwise velocity profile normalized with the inlet velocity at 3 checkpoints on the wake region behind a circular cylinder at  = 3900.

Fig. 34 .Fig. 35 .
Fig. 34.Vertical velocity profile normalized with the inlet velocity at 3 checkpoints on the re-circulation region behind a circular cylinder at  = 3900.

Fig. 36 .
Fig. 36.Reynolds shear stress normalized with the square inlet velocity at 3 checkpoints on the wake region behind a circular cylinder at  = 3900.

Fig. 37 .
Fig. 37. Initial location and dimensions of the water column and vertical barrier as shown in [99].

Fig. 39 .
Fig. 39.Predicted velocities by the GCM near the probe location and comparison with experimental data [99].

Fig. 40 .
Fig. 40.Evolution of the water surface in the dam-break case with a vertical cylinder.

Table 2
Drag, lift and Strouhal coefficients predictions for a circular cylinder at  = 100.

Table 4
Comparison of predictions from the ghost cell method for a circular

Table 5
Comparison of predictions from the ghost cell method for a circular cylinder and  = 100 with a mesh resolution of  Δℎ = 40 and different computational domain heights.

Table 6
Comparison of drag, lift coefficients and Strouhal number predictions for a square cylinder at  = 100 between different methods in the literature.The present results obtained from the standard IBM and the ghost-cell method are for a mesh resolution of  Δℎ = 40 and time step Δ = 0.005 s.