level-set method to simulate incompressible two-phase flows with embedded moving solid boundaries

Computers & Fluids 71 (2013) 469–486

Contents lists available at SciVerse ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

An adaptive Cartesian cut-cell/level-set method to simulate incompressible two-phase flows with embedded moving solid boundaries Meng-Hsuan Chung ⇑ Department of Naval Architecture and Ocean Engineering, National Kaohsiung Marine University, No. 142, Haijhuan Rd., Nanzih District, Kaohsiung City 811, Taiwan, ROC

a r t i c l e

i n f o

Article history: Received 18 April 2011 Received in revised form 7 June 2012 Accepted 5 November 2012 Available online 21 November 2012 Keywords: Cartesian grid method Cut cell approach Level set method Adaptive mesh refinement Two-phase flow Projection method

a b s t r a c t In this paper, a Cartesian grid method has been developed to simulate two dimensional unsteady viscous incompressible two-phase flows with embedded moving solid boundaries. The solid boundary is treated by a cut cell approach. The fluid–fluid interface is captured and treated by the level set method. A collocated finite volume method with nominally second-order accurate schemes in space is used for discretization. A pressure-free projection method is used to solve the equations governing incompressible flows. Adaptive mesh refinement is applied to efficiently deploy mesh cells of various sizes. The present method proposes an outer-iteration solution procedure to solve two-phase flows with partially submerged moving bodies and physically treats issues on newly exposed cells. Various test examples were computed and validated against previous analytical solutions, simulations, or experiments to prove the accuracy and effectiveness of the present method. We observed first-order spatial accuracy in a specific validating problem. In brief, the major contribution of the present work is the first-time announcement of a well validated level set method based on the cut-cell approach in moving-body problems. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction The major advantages of the non-body-fitted grid method, mostly the Cartesian grid method, over the body-fitted grid method include the simplicity of grid generation and getting rid of mesh-skewness error. Particularly, a prominent advantage appears when the immersed interface/boundary is moving arbitrarily because, in this case, considerable errors due to grid distortion, re-meshing or interpolation for overlapping grids are almost inevitable for the body-fitted grid method. Further, re-meshing and interpolation are very time-consuming. Large amount of literature can be found on Cartesian grid methods to simulate single-phase flows with embedded solid boundaries. Cut cell approach, as a branch of Cartesian grid method, has been applied to incompressible viscous flows in the literature, e.g. [1–4]. A cut cell is defined as one intersecting with solid boundary. The geometric shape of solid boundary inside this cut cell is explicitly and sharply set as the original or some approximative simple shape. Another branch of Cartesian grid method, immersed boundary method, generally has also succeeded in simulating viscous flows, e.g. [5–8]. The immersed boundary is represented by a discrete set of body or surface forces in these

⇑ Tel.: +886 7 3617141x3076; fax: +886 943 203377. E-mail addresses: [email protected], [email protected] 0045-7930/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2012.11.002

methods, hence programming is easier. However, the cut cell approach could better resolve boundary layer flows because the solid boundary is represented as a sharp interface exactly on which the boundary conditions are imposed. On the other hand, many interface tracking/capturing methods to simulate free-surface/multi-phase flows have been proposed. A qualified numerical method must be able to follow arbitrary motions of irregularly shaped interfaces. The simplest method is the so called height function method [9] which is applicable only in single-valued problems. Harlow and Welch [10] first proposed Marker and Cell (MAC) method. This method tracks a given set of discrete marker fluid particles, the distribution of which marks the region occupied by a certain fluid. Thus the interface can be accurately and precisely captured. Front-tracking methods of [11,12] use a group of discrete points in a one-dimension lower space to represent each interface. The motion of the interface is tracked by updated positions of these points. However, all the above methods can hardly handle topology changes of interface such as merging or splitting. The moving interface problem was also frequently solved by the ALE (Arbitrary Lagrangian Eulerian) method, e.g. [13–15]. Again, because the grid moves and deforms with the interface, numerical difficulties will surface when dealing with violently changing interfaces. Among those capable of treating interface merging and splitting, Volume of Fluid (VOF) method and Level Set (LS) method are most popular. The VOF method was originally proposed by Hirt and Nichols [16] and then improved by others [17–21]. This method

470

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

uses a volume fraction or color function u to identify the phase of a computational cell as pure phase 1 (u = 0), pure phase 2 (u = 1) or mixture (0 < u < 1). The volume fraction is advected by the fluid flow. Though mass conservation and topological change can be handled well, it is difficult to compute accurate local curvatures from volume fractions due to sharp change of volume fraction in interface regions. By the way, in the phase-field simulation of solidification process, e.g. [22–25], the phase-field variable is similarly defined as solid fraction, the space–time variation of which is governed by an equation derived from thermodynamics. The LS method [26,27] represents the interface as the zero-valued surface of a level set function /. The two phases reside in the region of / < 0 and / > 0 respectively. Typically the level set function is maintained as the signed distance to the interface, therefore one can calculate the normal and curvature of the interface by very simple formulas. This method can handle topological merging and breaking naturally. Meanwhile, complexity in 3D implementation is nearly the same as in 2D. Although the level set function follows the same advection equation as for the volume fraction, it suffers from excessive smoothing especially in under-resolved regions, hence severely destroying mass conservation. Attempts to improve mass conservation have led to various reinitialization methods [28,29], higher order ENO/WENO schemes for spatial derivatives in the convection and reinitialization equations [30], combination with Lagrangian front tracking methods [31], use of adaptive mesh refinement (AMR) [32], and coupling with the VOF method [33,34]. In the following we will review in more detail Cartesian grid methods for free-surface/multi-phase flows, particularly with embedded moving solid boundaries. Kleefsman et al. [35] numerically studied some aspects of 2D/ 3D water impact/entry problems. The governing equations of fluid flow were discretized on a staggered fixed Cartesian grid using the finite volume method. The free surface is captured using the VOF method and the solid–liquid interface treated by a cut-cell method. They simulated a dambreak flow with a downstream box as a simple model of green water flow on the deck of a ship. The water entry problem has also been investigated by dropping wedge, cylinder and cone into calm water with a prescribed velocity. Finally, they simulated water impact by a free falling wedge. They showed some success of their numerical method by these test cases. However, some issues remained, including (1) the diffusive term near a solid boundary is discretized as if the boundary has a staircase shape, (2) they did not mention how to treat troubles caused by newly exposed cell centers in moving-body problems, and (3) several spikes appeared in the time signals of pressure. Qian et al. [36] developed a two-fluid solver for 2D flows with moving solid bodies. They solved the incompressible Navier– Stokes equations for a variable density fluid system. The cut cell method is used for treating solid boundaries and the captured contact discontinuity in the density field is regarded as the fluid–fluid interface. The underlying method is a time-accurate artificial compressibility method with a high-resolution Godunov-type scheme. Their test cases include a moving paddle as a wave generator and water impact of a rigid wedge, compared with previous experimental data (transient surface elevation) and theoretical results (pressure distribution and free surface profile) respectively. The whole processes of water entry as well as water exit of a wedge were also simulated. Wang and Wang [37] improved the method of Qian et al. by introducing a body-acceleration term into the pressure boundary condition for moving solid and using the exact Riemann solution to calculate the flux on solid boundary. The improvement on prediction of transient vertical force and vertical velocity in the initial stage for a free-falling wedge is obvious. Further, the predicted slamming coefficient for constant-velocity water entry of a circular cylinder globally agrees with some previous experimental data. However, both works tackled only inviscid

flows, reported intensive consumption of CPU time due to small time step size, and neglected the issue of newly exposed cells. Lin [38] developed a finite difference staggered Cartesian grid method for 2D turbulent free-surface flows with a moving body. The RANS equations are solved by a two-step projection method with a nonlinear eddy viscosity model together with the conventional k–e model to compute the Reynolds stresses. He used the VOF method and cut-cell method to treat the free surface and solid–fluid interface respectively. An irregular body is represented by the volumetric fraction of solid in Cartesian cells (partial cell treatment, PCT). A concept of ‘‘locally relative stationary (LRS)’’ was introduced where a source term is added locally to the conventional continuity equation on body surfaces to take account of body motions. He tested his method by simulating a wave-paddle generated solitary wave, water exit and water impact and entry of a horizontal circular cylinder, fluid sloshing in a horizontally excited tank, and the acceleration/deceleration of an elliptical cylinder near a water surface. However, no comparison was made with other works in literature for the water impact-entry of a circular cylinder and the complicated motion of an elliptical cylinder near a water surface. In other test cases, validation was made only for free-surface shapes. Again, he did not mention how to treat issues due to newly exposed grid points. The CIP (cubic-interpolated propagation/constrained interpolation profile) method [39] tracks and treats the solid–fluid interface and fluid–fluid interface in a unified approach. The solid phase was treated in the same way as the fluid phase, with the solid-phase velocity prescribed or solved from the equations of the motion governing the solid phase. The CIP method has been applied to simulate many 2D and 3D two-phase flows with moving solid boundaries. However, the color function for tracking interfaces in the CIP method is smoothed before advancing in space according to the governing advection equation. Thus interfaces are diffused and solid boundaries possibly distorted. Further, because derivatives of the primitive variables must be stored and solved, a vast additional computational cost is inevitable. Please refer to [40] for some variants of the CIP method and [41] for their applications. Yang and Stern [42] presented a 3D sharp interface staggered finite-difference Cartesian grid method for the large-eddy simulation of two-phase turbulent flows interacting with moving bodies. A four-step fractional-step method is employed for velocity–pressure coupling. They used a second-order direct forcing sharp interface immersed boundary method to treat solid–fluid interface and a level-set/ghost-fluid method for fluid–fluid interface. A field extension procedure was employed to construct previous-step values of the fluid velocity and level-set function at newly exposed mesh points. Validation cases include bubble rising, water entry and exit, landslide-generated waves, and ship hydrodynamics. However, they reported neither pressure distribution nor hydrodynamic force exerted on solid boundaries in all the test cases. Bai et al. [43] applied the SIMPLE algorithm in a collocated finite-volume framework to simulate 2D viscous free surface flows around moving solids. They employed a Cartesian cut cell approach to track and treat the free surface as well as solid boundaries. Apart from some validating cases, their method was applied to simulate radiation waves induced by a forced oscillating submerged circular cylinder. However, a single-valued height function was adopted to evaluate the free surface position, causing their method to fail in simulating breaking or highly distorted non-breaking waves. Meanwhile, treatment for newly exposed cells is unclear. Zhang et al. [44] developed a level set immersed boundary method to simulate 2D constant-velocity water entry/exit of a circular cylinder and free falling of a wedge with a dead-rise angle of 30°. They reported for the circular cylinder a slamming coefficient curve in fair agreement with previous works. However, in the wedge case, they underestimated pressure coefficient, especially

471

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

to a large extent for the peak value, in contrast to previous experimental and numerical works. They also did not state how to overcome issues on newly exposed cells or why these cells do not cause any problem. The present work develops a 2D collocated finite-volume method to simulate unsteady, viscous, incompressible two-phase flows with embedded moving solids. A pressure-free projection method [4] is used as the underlying incompressible flow solver. The solid boundary is treated by the adaptive Cartesian cut cell approach of Chung [45]. The fluid–fluid interface is captured and treated by the level set method. The present method features (1) employing an outer-iteration solution procedure to solve in a coupled way the evolutionary equation of the level set function and the Navier– Stokes equations for flows with partially submerged moving bodies, (2) physical treatment of newly exposed cells, and (3) use of adaptive mesh refinement (AMR). Detailed description of the present method will be given below followed by some verifying and validating examples.

value in phase 1 and negative value in phase 2, then the interface can be represented by the zero level set of /. The initial values of / are prescribed as the signed normal distance to the interface. Consider the following evolutionary equation:

2. Governing equations of fluid flow

The smoothing treatment in Eq. (8) for the discontinuity in the density field is applicable for dynamic viscosity as well. Here a  O(hfs) where hfs is the mesh size near the free surface. We take a = 1.5hfs in this study. For some problems, in order to damp outgoing disturbances generated by solid bodies, numerical sponge layers are introduced close to the upstream and downstream boundaries of the computational domain. If the positive x direction points toward downstream and the still free surface is located on the plane of y = ystill, the level set function obeys the modified evolutionary equation [46]:

In the present work two dimensional flow field of two immiscible incompressible fluids is calculated. The solutions to both phases, denoted as phase 1 and phase 2, are obtained simultaneously. The nondimensional equations of motion are given by the incompressible Navier–Stokes equations:

r  u ¼ 0;

ð1Þ

     1 1 rp þ r  l ru þ ruT ut þ ðu  rÞu ¼ B þ ; Re q

ð2Þ

where u = ui + vj is the fluid velocity, p = p(x, t) the static pressure, q = q(x, t) the fluid density, l = l(x, t) the fluid viscosity, and B = Bx i + By j the body force per unit mass. The superscript ‘‘T’’ represents transpose of matrix. Vectors i(j) are the unit vector directed toward the positive x(y) coordinate axis. For immiscible fluids the density and viscosity are constant following the fluid particle, therefore

qt þ ðu  rÞq ¼ 0;

ð3Þ

lt þ ðu  rÞl ¼ 0:

ð4Þ

The above nondimensional physical quantities are defined as:

u

~ u ; U

p

~ p

q2 U

2

; B

~ B 2

U =L

; q

~t q~ l~ ; l ; t ; q2 l2 L=U

g U 2 =L

j¼

1 Fr2

j:

ð7Þ

It turns out that the zero level set will move in the same way as the interface does. Further, due to smoothness of /, the above equation can be solved without problems caused by numerical instability. Then the fluid density can be obtained algebraically as:

8 if / > a; > < q1 =q2 q¼ 1 if / < a; > : q þ Dq sinðp/=2aÞ if j/j 6 a;

ð8Þ

where



q ¼ ðq1 þ q2 Þ=ð2q2 Þ; Dq ¼ ðq1  q2 Þ=ð2q2 Þ:

ð9Þ

/t þ ðu  rÞ/ ¼ Sd ðx; y; /Þ  md ðxÞðy  ystill  /Þ;

ð10Þ

with







md ðxÞ ¼ C damp Llow  max ðLlow ; 0Þ þ Lhigh  max Lhigh ; 0 ;

ð11Þ

Llow ¼ xlow  x and Lhigh ¼ x  xhigh ;

ð12Þ

where xlow and xhigh denote the x coordinates of the inner end of the upstream and downstream sponge layers respectively. After some preliminary trials on the value of the damping coefficient, Cdamp, we chose Cdamp = 0.04 to obtain satisfactory results in this work. To avoid boundary conditions conflicting with the residual disturbances therein, we further introduced a buffer layer where the level set function is always fixed to the local initial value. Fig. 1 is the schematic presentation of the sponge and buffer layers.

ð5Þ

where U and L are the characteristic velocity and length respectively and ‘‘’’ denotes the dimensional quantity. If the gravitational force pointing to j is the only body force, then

B¼

/t þ ðu  rÞ/ ¼ 0:

ð6Þ

Besides Re and Fr, other nondimensional parameters are q2/q1 and l2/l1. 3. Numerical method 3.1. Level set approach 3.1.1. Surface capturing and property smoothing Due to rapid change of density and viscosity at the interface, there must be severe numerical diffusion when solving Eqs. (3) and (4) by traditional numerical schemes. In the present work the remedy is to capture the interface through the use of level set function [27]. Let / denote the level set function with positive

3.1.2. Reinitialization The smoothing process will properly work only on the premise that the level set function remains a distance function. This important property is inevitably destroyed in the evolution of flow field. Hence / must be reinitialized to recover to a distance function. Among variants of reinitialization, the following method [27] is chosen. Consider a given function /0(x) whose zero level set is the fluid–fluid interface. A signed normal distance function /(x) with the same zero level set as that of /0(x) can be constructed by solving the following problem to steady state:

/t ¼ Sð/0 Þð1  jr/jÞ; /ðx; 0Þ ¼ /0 ðxÞ;

ð13Þ

where S is the smoothed sign function [47]:

/0 : Sð/0 Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 /0 þ jr/0 j2 hfs

ð14Þ

This definition of sign function differs from that of Sussman et al. [27] and exhibits faster rate of convergence in numerical

472

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

Buffer Sponge layer layer

Sponge Buffer layer layer

xlow

xhigh

y

Upstream boundary

Downstream boundary

x

Fig. 1. Schematic demonstration of sponge layer and buffer layer used to damp out outgoing free-surface disturbances.

solution when / changes mildly. Meanwhile, it would not cause flip-flap of sign when / varies drastically. 3.2. Single-grid incompressible flow solver The governing equations, Eqs. (1) and (2) constitute a nonlinear and coupled system of partial differential equations. A cell-centered collocated finite volume Cartesian grid method with the cut cell approach is applied to discretize this set of coupled equations. The coupling between the divergence-free constraint and the pressure gradient in the momentum equation is solved by a pressure-free projection method. In the present Cartesian grid method, we implement boundary conditions on embedded solid-body surface by the cut cell approach. Details of cell merging and interpolation strategies at cut-cell-face centers and the probing point on the body-surface normal can be found in Chung [4]. Particularly, the cell merging technique in the present cut cell approach can better resolve boundary layer flows than other Cartesian grid methods. This implies that hydrodynamic forces can be calculated more accurately. We treat issues on newly exposed cell centers in moving-body problems by the improved technique of Chung [45]. That is, numerical solution of time-dependent equations for these cells is abandoned. Instead, current values of unknown variables (u, v, and /) at these cell centers are linked to those at the neighboring solution nodes by the inverse-distance interpolation scheme. Thus we avoid referring to previous values of unknown variables at newly exposed cell centers which cannot be physically defined. Simulating two-phase flows in this work, the following sub-subsections demonstrate major modifications from the above single-phase Cartesian grid method. Taking a finite volume cell, regular or deformed, as the control volume, cv, with control surface cs, the momentum equation, Eq. (2), is numerically advanced in time in the advection–diffusion step to generate intermediate velocities. Then these velocities are corrected in the pressurecorrection step. 3.2.1. Advection–diffusion step If the control volume does not pertain to any newly exposed cell center in moving-body problems, the following advection–diffusion step is used to obtain intermediate velocities:

Z V cell  ðucc  uncc Þ þ ½uðU  nÞnþ1=2 dS Dt cs Z 1 1 nþ1=2 l runþ1=2  n dS þ nþ1=2 ¼ nþ1=2 qcc Re cs qcc Re Z  ðrlÞnþ1=2 ðruT Þnþ1=2 dV;

ð15Þ

cv

where Vcell is the volume of the control cell, cc denotes value at the cell center, n the time step number, Dt time step size, ucc cell-center intermediate velocity, U cell-face-center velocity.

The terms at the (n + 1/2)th time step are approximated as follows: Z nþ1=2 ½uðU  nÞ dS cs X 8R  n u ðU  nÞdS u ðUn  nÞAfc for fixed-body problems; > > < cs fc X R   > u ðU  nÞAfc for moving-body problems; > : cs u ðU  nÞdS fc

ð16Þ

Z

1 nþ1=2 Re cc

q



cs

1

X

qnþ1 cc Re

fc

Z

1 nþ1=2 Re cc

q



lnþ1=2 runþ1=2  n dS

cv

1

lnþ1 fc

  @u Afc ; @n fc

1

Z

q

nþ1 Re cc

lnþ1 ru  n dS

cs

ð17Þ

ðrlÞnþ1=2 ðruT Þnþ1=2 dV Z

qnþ1 cc Re

ðrlÞnþ1 ðruT Þn dV

cv

V cell

qnþ1 cc Re

T n ðrlÞnþ1 cc ðru Þcc ;

ð18Þ

where fc labels a specific cell-face center, Afc is the cell-face area, P and fc represents summation over all the cell faces on the control surface. Note that ufc in Eq. (16) must be further expressed as a weighting sum of u⁄ values at neighboring solution nodes. The cell-face-center intermediate velocity, U⁄, for a non-boundary cell face is interpolated out from nearby ucc and/or boundary U⁄. However, Un has its own generating mechanism as stated in the next sub-subsection, not through interpolation from neighboring uncc .   The normal gradient of u⁄ at the cell-face center, @u , is approxi@n fc mated by the central difference scheme for regular cell faces, by the first-order one-sided difference scheme for cell faces on the boundary of the computational domain, and by the bilinear interpolation scheme for cut cell faces. In the explicit diffusion term of Eq. (18), ðrlÞnþ1 and ðruT Þncc are approximately calculated using the followcc ing mathematical formula valid for any function F:

ðrFÞcc

1 V cell

Z cv

rF dV ¼

1 V cell

Z cs

Fn dS

1 X F fc Afc nfc : V cell fc

ð19Þ

Note lnþ1 for a non-boundary cell face is interpolated out from fc neighboring lnþ1 and/or boundary lnþ1 cc fc . For newly exposed cell centers in moving body problems, the inverse-distance interpolation scheme [45] is called for to evaluate unþ1 cc . 3.2.2. Pressure-correction step The cell-center velocity at the (n + 1)th time step is calculated by correcting the cell-center intermediate velocity as shown in the following equation:

unþ1  u rpnþ1=2 ¼  nþ1=2 þ B: Dt q

ð20Þ

473

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

If relevant spatial derivatives exist, then we can take divergence of Eq. (20). Spatial integration of this divergence over the control volume gives:

 nþ1   nþ1=2  Z Z u  u rp dV ¼  dV þ r r r Dt qnþ1=2 cv cv cv

Z

 B dV:

ð21Þ

To satisfy the constraint of free divergence of velocity field, Eq. (1), Un+1 must satisfy

Z

nþ1

U

 n dS ¼ 0:

ð22Þ

manner as in the advection–diffusion step except / replacing u and discarding the viscous terms in Eq. (15):

 V cell   /cc  /ncc þ Dt

n+1/2

n+1

¼

Z

1

cs

rpnþ1  n dS ¼ nþ1

q

 Z   U þ B  n dS: cs Dt

ð23Þ

The discretized Poisson equation for pressure is thus

  X Afc @pnþ1 X U A ¼ þ B n : fc @n Dt qnþ1 fc fc fc fc

ð24Þ

Note that B-related terms can be dropped in Eq. (24) if B is not a function of space. The cell-face-center velocity can thus be corrected using the pressure field solved from Eq. (24) as follows:

" nþ1

U



¼ U þ Dt B 

1

qnþ1 fc

# ðrp

nþ1

Þfc :

ð25Þ

That is, Un+1 is obtained by correcting U⁄ as shown in Eq. (25) though U⁄ itself is evaluated by interpolating the neighboring u⁄. The introduction of the cell-face-center velocity not only eliminates the grid-to-grid pressure oscillation in collocated methods, but also can reduce the mass conservation error to machine zero if needed [1,48]. Notably, we employ Eq. (24) for all control volumes and Eq. (25) for all control surfaces, including those associated with newly exposed cell centers in moving-body problems. As for the cell-center velocity, we calculate unþ1 in two different cc ways, depending on the region the cell center is located:

unþ1 cc



8 R nþ1 > ucc þ Dt B  qnþ11V p n dS > cs > cell cc > ! > < X nþ1  1 ¼ pfc Afc nfc for j/cc j > b; > ucc þ Dt B  qnþ1 cc V cell > > fc > > : interpolated out from surrounding Unþ1 for j/cc j 6 b; ð26Þ

where b = 3hfs. That is, we correct ucc to obtain unþ1 in regions some cc distance away from the free surface and interpolate out unþ1 from cc surrounding Un+1 in regions near the free surface to avoid divergent numerical values caused by large variation of fluid density. 3.2.3. Discretization of the evolutionary equation of level set function Because correct solution of / is relevant only near the free surface, we can assign the level set function a constant value wherever distant from the free surface, thus avoiding possibly unwanted numerical solutions in this region. Specifically, / is assigned c + hfs in phase 1 and (c + hfs) in phase 2 whenever |/| > c, where c is taken as 4hfs in this work. Furthermore, the same inverse-distance interpolation scheme as used for the fluid velocity is called for to evaluate /cc at newly exposed cell centers in moving body problems. For cell centers not belonging to the above categories, /cc is obtained by solving the evolutionary equation, Eq. (10), in a similar

½/ðU  nÞnþ1=2 dS ¼ V cell Sd ðxcc ; ycc ; /ncc Þ;

fc

X  nþ1 > > / ðU  nÞAfc > :

for fixed-body or fully submerged moving-body problems ;

ð27Þ

ð28Þ

otherwise :

fc

n+1

Approximating p and q by p and q respectively, applying Gauss’s theorem, and then combining with Eq. (22) eventually produces the Poisson equation for pressure:

cs

where Z ½/ðU  nÞnþ1=2 dS cs 8X > > / ðUn  nÞAfc > <

cs n+1/2

Z

That is, when the body moves and occupies two phases simultaneously, we propose to approximate Un+1/2 by Un+1 instead of Un because the latter is usually unknown at cell faces near the moving-body surface. This implicit approximation calls for an outer iteration loop over three steps: updating the level set function, the advection–diffusion step, and the pressure-correction step. 3.2.4. Discretization and solution of the reinitialization equation of level set function To achieve high-order accuracy in approximating spatial derivatives of the level set function without loss of numerical stability, we perform time integration by the second-order Runge–Kutta method and approximate the first-order derivatives in space by the second-order ENO scheme [26]. Further, according to Peng et al. [47], numerical solution is confined in a region close to the free surface to save on computational cost. The region is chosen as within width of c away from the free surface. The process of reinitialization modifies /⁄ to /n+1, thus obtaining qn+1 and ln+1via Eq. (8). 3.2.5. Solution of other discretized equations In Eqs. (15), (24), and (27), the spatial discretization of advection, diffusion and pressure terms at the cell-face center are all treated by the central difference scheme. For any solution variable (u, v, p, or /), we can establish a system of algebraic equations consisting of discretized equations for all cells. The coefficient matrix of this equation system is sparse and we adopt a free open-source package of libraries, SPARSKIT Version 2, as tool for matrix computation and iterative solution of algebraic-equation system (inner iteration). The source code and documentation can be downloaded from the following web page: http://www-users.cs.umn.edu/ ~saad/software/SPARSKIT/sparskit.html. From this tool kit, the restarted Generalized Minimal RESidual method [49], called GMRES(m), is selected to solve the linear system which is right preconditioned by the incomplete LU factorization with dual truncation strategy (ILUT). Throughout the simulations in this work, the choice of the dimension of Krylov subspace, m, is 3 for the u, v and / equations and 20 for the pressure Poisson equation. The relative dropping threshold in constructing L and U is set at 1  103 for the u, v and Poisson equations and 5  103 for the / equation. The number of fill-in elements in L and U is set at 100 for the u, v and Poisson equations and 60 for the / equation. There are no definite rules for determining or optimizing these numerical values, hence they were settled by trial and error in this work. The convergence criterion for solution of the algebraic-equation system is ðkÞ

ðk1Þ

kZ/  Z/

k1 < TOLðiÞ u;

ð29Þ

where Z/ represents a vector composed of values of any solution variable u at all discrete points, k the inner iteration number and TOLðiÞ u the threshold of inner convergence for variable u. For all

474

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

the simulations in this work, the inner-iteration convergence criteðiÞ ria for p and u are TOLpðiÞ ¼ 103 and TOL/ ¼ 108 . For fixed-body problems, the inner-iteration convergence criteria for u and v are ðiÞ 6 TOLðiÞ u ¼ TOLv ¼ 10 . However, the u and v equations in movingbody problems must be solved in a coupled way because the implicit option is adopted in Eq. (16). That is, one step of such a coupled inner iteration comprises two cycles of m GMRES iterations, the first for u and the second for v equations. The convergence criterion of the coupled inner iteration is exactly the same as Eq. (29) provided that k represents the coupled inner iteration number. In this work, ðiÞ 2 we take TOLðiÞ for moving-body problems. u ¼ TOLv ¼ 10 In partially submerged moving-body problems, Eqs. (15), (24), and (27) must be implicitly solved altogether, as explained in Section 3.2.3. In this case, the convergence criterion for the outer iteration loop is also the same as Eq. (29) with k being the outer ðiÞ iteration number and TOLðoÞ u in place of TOLu the threshold of outer convergence for variable u. Free-surface flow dynamics being the major concern in this work, we use the following thresholds: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðoÞ ðoÞ 2 hfs =Fr2 and TOL/ ¼ 102 hfs . Note that we TOLðoÞ u ¼ TOLv ¼ 10

make judgment only in the region close to the free surface, i.e., with |/| < c, and disregard to monitor the pressure variable. Since the viscous stability limit is removed by advancing the viscous terms implicitly, the time step size must obey the CFL-like condition due to the convective terms and the constraint imposed by gravitation:

Dt ¼ minðDtCFL ; Dt G Þ;

ð30Þ

with



   Dx Dy ; min min ; i2fluid node juj jv j i     Dx Dy min ; min ; i2solid boundary node jun;x j jun;y j i

DtCFL ¼ fCFL  min

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !ffi u u D y D x fs fs ; DtG ¼ fG  t0:5 min ; jg x j jg y j

ð31Þ

ð32Þ

where respectively in the x- and y-directions, Dx and Dy are mesh sizes, un,x and un,y denote the two components of the normal velocity at the solid boundary node, and gx and gy the gravitational accelerations. In this work, Dx = Dy and Dxfs = Dyfs = hfs. Further, for numerical stability, the time step size cannot increase by more than 10% from the current to the next time step. For all the simulations in this work, the ranges of fCFL and fG are 0.8 6 fCFL 6 1 and fG = 1 for stable solutions. 3.2.6. Boundary conditions 3.2.6.1. Intermediate velocity. On the inlet plane and solid boundaries, the intermediate velocity satisfies

U ¼ fðtnþ1 Þ;

ð33Þ

where f(tn+1) denotes the fluid velocity on the inlet plane or the solid-boundary velocity at t = tn+1. On the outlet plane,

 @U  ¼ 0: @n outlet

1

qnþ1

 @pnþ1 1 ¼ ðU  Unþ1 Þ þ B  n: @n Dt

ð35Þ

3.2.6.4. Level set function. The boundary conditions of the level set function are everywhere of zero-gradient type except at the inlet plane where the level set function is prescribed as signed distance from the still-water position. 3.3. AMR The present AMR technique follows largely the work of Chung [45]. In brief, the computational domain is discretized by various hierarchical square-shaped sub-grid blocks which form the nodes of a quad-tree data-structure. All the sub-grid blocks have the same logical structure and cell number in each spatial dimension. They differ only in size, location, parent block, and list of neighboring and child blocks. If the spatial resolution in a sub-grid block is judged to be increased, it spawns four child blocks at one more refinement level, and hides behind as a parent block. The prolongation of bilinear interpolation is used to evaluate the variable values at the cell centers of each spawn child block. This process is called ‘‘refinement’’. If the spatial resolutions in a parent grid block and all of its four child blocks are judged to be decreased, the child blocks merge and disappear with the parent block surfacing. The restriction of arithmetic mean is used to evaluate the variable values at the cell centers of the surfacing grid block. This process is called ‘‘de-refinement’’. Each sub-grid block has (N  N) mesh cells at all refinement levels. For grid blocks cut by or close to (within four times of mesh-cell size away from) the body surface, the refinement level keep unchanged throughout the simulation, called the body refinement level. The refinement level of grid blocks near the fluid–fluid interface, called the interface refinement level, is fixed constant as well. If we know or expect the body surface and fluid–fluid interface will come close together at whatever instant in the simulation process, then the two refinement levels must be the same. For other grid blocks, the basic criteria for refinement and de-refinement are respectively set to

maxjxi;j j > r x xmax i;j

ði ¼ 1; . . . ; N and j ¼ 1; . . . ; NÞ

ð36Þ

and

maxjxi;j j < 0:5r x xmax i;j

ði ¼ 1; . . . ; N and j ¼ 1; . . . ; NÞ;

ð37Þ

where xi,j is the vorticity vector of cell (i, j) and xmax is the maximum magnitude of vorticity vector in the current flow field, and rx = 0.001–0.01 for simulations in this study. Usually, in the bulk region off solid bodies and fluid–fluid interfaces, the prescribed maximum refinement level is smaller than the body or interface refinement level to effectively reduce the total cell numbers. Finally, a second-order accurate scheme with three-point compact stencil is used to approximate the cell-face flux on the co-faces between fine and coarse meshes.

ð34Þ

3.2.6.2. True velocity. The boundary conditions of true velocity are everywhere the same as those of the intermediate velocity except on the outlet plane where Un+1 is first extrapolated from the nearest interior cell center and then added a uniform correction to secure global conservation of mass in the computational domain. 3.2.6.3. Pressure. The boundary conditions of pressure come from projection of Eq. (25) on the direction normal to the boundary:

3.4. Summary and solution algorithm We compare the present method with that of Sussman et al. [32] to summarize and highlight important characteristics of the former different from those of the latter. Sussman et al. presented an adaptive projection method coupled with the level set approach to solve incompressible two-phase flow with surface tension. However, no embedded solid boundary was considered. They employed staggered grid with the pressure variable put at cell corners, in contrast to the present collocated grid. To approximate the

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

nonlinear advection term, their method uses an explicit predictor– corrector scheme. However, explicit schemes have no way to treat newly exposed cells due to previous values being unavailable at these cell centers. The present method thus adopts the fully implicit scheme though order of accuracy could be lowered. To solve the intermediate velocity in the advection–diffusion step, they keep the lagged pressure gradient term which, in the present study, turns to be unclear and would make troubles with newly exposed cell centers. That’s why we call for the pressure-free projection method. As for the level set method, Sussman et al. appealed to an improved reinitialization step [29] for better preserving mass of either fluid phase. This improvement and others alike are not our major concerns because we aim at developing a workable numerical framework to combine the cut-cell Cartesian grid method and the level set method. Successful combination relies on effective treatment of issues on newly exposed cells in moving solid problems. Numerical techniques involved in solving these issues are what we concern only. Other numerical techniques and their various modifications/improvements, e.g. the improved reinitialization step [29], do not affect construction of the framework. In other words, most of them can be straightforwardly implemented within the present framework. On the other hand, Sussman et al. has developed an AMR technique based on the adaptive framework for single-phase flows [50]. Within their framework, the whole grid system is composed of adaptive hierarchy of nested rectangular structured grids, with fine grids inside and overlaid on coarse ones. The present AMR technique utilizes the quad-tree data structure to record and manage relationships between sub-grid blocks. Therefore, the present data structure is much simpler than theirs. Meanwhile, solution by the present technique is required of only leaf sub-grid blocks, all at a time, while Sussman et al. must find solutions on all grids at all refinement levels. Although the present technique causes unstructured solution matrix, we found that the GMRES module in SPARSKIT performed robustly for all the test examples in this work. Finally, we note that Lan et al. [25] developed their adaptive finite volume method using the quad-tree data structure where a node of the tree represents a mesh cell. However, a tree node corresponds to a sub-grid block in the present method. The solution algorithm is summarized and presented in the flowchart of Fig. 2.

4. Test examples 4.1. Fluid sloshing in a fixed rectangular tank We use this test case to demonstrate effectiveness of capturing mildly moving free surface by the present method. Initially (t = 0) a still water partially fills the tank and the water surface is given a half-wavelength sinusoidal profile with an amplitude of 1% average water depth. Then the water starts to slosh due to gravity effect. The computational domain is [1, 1]  [1, 1] with the average water surface being located at y = 0. The no-slip condition is imposed on the four boundaries. Taking pffiffiffiffiffiffithe average water depth, H, as the characteristic length and gH as the characteristic velocity gives Fr = 1 and we performed simulations using Re = 1  104. The ratios of fluid density and dynamic viscosity are q2/q1 = 815.7 and l2/l1 = 63.9 respectively. Three grids with the same maximum refinement level of 9, corresponding to N = 2, 4, and 8 (hfs = 1/256, 1/512, and 1/1024), are used. Fig. 3 depicts the predicted free-surface heights as function of time at the left and right walls. Results for the three grids agree very well. The sloshing period is predicted to be 5.263, corresponding to a relative error of 0.53% against the period derived from the linear theory, 5.235.

475

4.2. Dam breaking To validate the capability of the present method to capture rapidly varying free-surface profiles, this benchmark problem was simulated. A still rectangular water column with height-width ratio of 2 starts to flow at t = 0 due to gravity effect. The computational domain is [0, 5]  [0, 5] with the no-slip condition imposed on the four boundaries. Taking the initial width of water column, pffiffiffiffiffiffiffiffi W, as the characteristic length and gW as the characteristic velocity gives Fr = 1 and we performed simulations using Re = 2.752  106. The ratios of fluid density and dynamic viscosity are q2/q1 = 815.7 and l2/l1 = 63.9 respectively. Four grids corresponding to hfs = 5/128, 5/256, 5/512, and 5/1024 are used. Fig. 4 depicts as function of time the location of water front and the height of wetted wall, compared with experimental data of Martin and Moyce [51]. In brief, the present simulation results agree well with the experimental measurements, especially in terms of the height of wetted wall. Simulations of the four grids show converged location of water front in the early stage, but the prediction exhibits weak convergence as time increases. To explore this phenomenon, we compare predicted water surface profiles at t = 2.5 among the four grids (Fig. 5). The difference in the location of water front is seen to be caused mainly by different water surface profiles near the front. It is suggested that the influence of grid density on the water surface profile near the front strengthens with time because the front thickness decreases with time. Thus the grid density required to achieve grid convergence of this local flow structure increases with development of the dam-breaking flow. Nevertheless, most of the whole water surface profile has grid-converged solution. Fig. 6 presents the finest-grid simulated transient water surface (contour of / = 0) profiles at selected time instants. Also shown is variation of the layout of subgrid blocks in the course of computation. We can observe a kink on the water surface at the upper-right corner of the water body at t = 1. This local fine structure exists in all simulations except the coarse-grid one.

4.3. Uniform flow over a fixed circular cylinder beneath a free surface We use this test case to demonstrate the ability of the present method to accurately simulate the free-surface effect on flow characteristics in the wake of a fixed blunt-body. Fig. 7 illustrates the definition of this problem. The computational domain is [20, 20]  [20, 20] and the cylinder center is located at the origin. We adapt the vertical coordinate of the still-water level to generate various gap ratios, h/D. Simulations were performed with Re = 180 and Fr = 0.2 where D and U are taken as the characteristic length and velocity respectively. On the inlet plane, a uniform horizontal flow is prescribed and the water-surface level keeps unchanged. On the outlet plane, we assume zero-gradient boundary conditions for the two velocity components and the level set function. On the bottom and top boundaries, the horizontal and vertical velocities are set to the inlet velocity and zero respectively and the level set function assumes the zero-gradient boundary conditions. Finally, on the cylinder surface, we impose the no-slip condition for the velocity and the zero-gradient condition for the level set function. Note that the abovementioned sponge and buffer layers are employed to damp out water-surface disturbances generated by the cylinder. The two sponge layers start at xlow = 12 and xhigh = 12 respectively and each of the two buffer layers spans unit length. Further, the outlet velocity must be corrected as stated in Section 3.2.6.2 to meet the global conservation of mass. The ratios of fluid density and dynamic viscosity are q2/q1 = 815.7 and l2/ l1 = 63.9 respectively. We focus on prediction of the Strouhal number, St, based on the oscillation frequency of lift force.

476

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

Fig. 2. Flowchart of solution algorithm.

To examine the mesh-size effect, the case of h/D = 0.4 was simulated using three grids. All grids distribute their smallest mesh size over regions near the cylinder surface and fluid–fluid interface. They are 5/256, 5/384, and 5/512. Predicted force coefficients, CD and CL, as function of time using these three grids are presented in Fig. 8. Here C D  2F D =ðq2 U 2 DÞ and C L  2ðF L  F B Þ=ðq2 U 2 DÞ where FD, FL, and FB denote the drag, lift, and buoyancy respectively. It is observed that the medium and fine grids predict close peak values of either CD or CL. The calculated Strouhal numbers using the lift data during t = 60–100 are 0.1794, 0.1853, and 0.1884 for the coarse, medium, and fine grid respectively. The

method to determine the order of grid convergence follows the work of Roache [52]. If the grid resolution is fine enough, the simulation result would fall into an asymptotic range of convergence and the error reduce in the following manner as grid is refined.

e ¼ f ðDÞ  fexact ¼ C Dq þ H:O:T:;

ð38Þ

where e is the error, f the observed quantity in the simulation, D the grid spacing, C a constant, and q the order of convergence. If simulations using three grids with D = D1, D2, and D3 predict three corresponding quantities f1, f2, and f3, then q can be solved from Eq. (38), neglecting the higher-order terms, as:

477

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

0.015

5

(a)

Free surface height

0.010

Location of Water Front

Coarse Medium Fine

0.005

0.000

Exp. (Martin & Moyce, 1952) Present (coarse) Present (medium) Present (fine) Present (finest)

4

(a) 3

2

-0.005 1

0

0.5

1

1.5

-0.010

2

2.5

3

2.5

3

t 1

-0.015

0

10

20

30

40

50

Height of Wetted Wall

t 0.015

(b)

Coarse Medium Fine

Free surface height

0.010

0.005

0.8

(b)

0.6

0.4

Exp. (Martin & Moyce, 1952) Present (coarse) Present (medium) Present (fine) Present (finest)

0.2

0.000

0

0

0.5

1

1.5

2

t -0.005 Fig. 4. Present simulation results for the dam-breaking problem compared with experimental data of Martin and Moyce [51]: (a) location of water front and (b) height of wetted wall normalized by the initial water column height, both as function of time.

-0.010

-0.015

0

10

20

30

40

50

t Fig. 3. Predicted free-surface heights as function of time at the (a) left and (b) right walls for a fluid sloshing in a fixed rectangular tank.

1 0.5 0 0

ðD3 =D1 Þq  1 f3  f1 ¼ : ðD2 =D1 Þq  1 f2  f1

ð39Þ

Taking f as the Strouhal number, with D3/D1 = 1/2 and D2/ D1 = 2/3, gives the order of convergence by Eq. (39) as q 1.01, indicating a linear convergence. It can also be interpreted as the order of accuracy of the present simulator as a whole in computing real complex flows. Note the relative errors for the three grids are 9%, 6%, and 4% respectively, based on the nominal exact Strouhal number, 0.1972, calculated from this error model. In view of the above analysis, the medium grid is acceptable to simulate cases of other gap ratios. Fig. 9 shows predicted CD and CL as function of time for selected values of h/D using the medium grid. At first, amplitudes of CD and CL are maximized at h/D = 0.4 and 0.6 respectively. Secondly, average level of CL is zero for h/D = 1 and algebraically decreases with decreasing h/D, that is, the cylinder experiences a downward force (not counting buoyancy) which increases with decreasing gap ratio. Thirdly, CD and CL start to oscillate at the very beginning of simulation for finite gap ratios while the oscillation is incepted very lately for h/D = 1. The implication is that the existence of free surface favors triggering the vortex shedding mechanism.

1

2

3

4

5

Fig. 5. Predicted water surface profiles at t = 2.5 for the dam-breaking problem using four grids. Refer to Fig. 4.

Finally, we investigated the normalized Strouhal number, St/St0, as function of the gap ratio and compared our results with those of Reichl et al. [53] as shown in Fig. 10. Reichl et al. numerically simulated this problem using the CFD commercial package, FLUENT. Note that St0 is defined as the Strouhal number for h/D = 1. Reichl et al. did not report the value of St0 they adopted while the present medium-grid simulation predicts St0 = 0.1850. Nevertheless, the two numerical results exhibit very similar trends. 4.4. Water entry of a wedge at constant speed This problem is used to test the capability of the present cut cell approach in treating solid boundaries moving through a free surface. Note the outer iteration loop mentioned in Section 3.2.3 was called for in the following simulations. The computational domain is [7.5, 7.5]  [7.5, 7.5] with the initial water surface and wedge apex being located at y = 0 and 0.2 respectively. To compare with previous results, gravitational force is neglected. The no-slip condition is imposed on the four boundaries. We performed simulations at Re = 5  103 based on the wedge height and entry speed.

478

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

3

2

0.5 2

CD

1

0 0

1

2

3

4

5

3

Force Coefficient

1

0

-1

CL

1.0

Coarse Medium Fine

2 -2

0

20

40

60

80

100

t 1

Fig. 8. Calculated force coefficients as function of time using three grids for the problem of uniform flow over a fixed circular cylinder beneath a free surface with h/ D = 0.4.

0 0

1

2

3

4

5

3

1.5 2

1

0 0

1

2

3

4

5

3

2.0 2

1

0 0

1

2

3

4

5

Fig. 6. Finest-grid predicted transient variation of water surface profile and layout of sub-grid blocks for the dam-breaking problem. Number at the upper-right corner in each plot denotes the time. Each sub-grid block consists of 8  8 mesh cells.

The ratios of fluid density and dynamic viscosity are q2/q1 = 815.7 and l2/l1 = 63.9 respectively. Symmetrical wedges with dead-rise angle a = 30°, 45°, and 60° are considered. If the flow is assumed to be potential, no length scale can be found, implying existence of similarity solution. Therefore, we will present our results mainly in terms of pressure coefficient, Cp, versus normalized vertical coordinate along the edge of wedge, y⁄. Here Cp  2(p  p0)/(qV2) and y⁄  y/(ystill  yapex), where p0 is the pressure at the midpoint on the top surface of the wedge, V the water-entry speed, yapex the vertical coordinate of the wedge apex, and ystill the vertical coordinate of the initial still water surface. To investigate the effect of grid resolution, three grids with the same maximum refinement level of 9, corresponding to N = 4, 6, and 8 (hfs = 15/1024, 15/1536, and 15/2048), are used to simulate the case of a = 30°. Fig. 11 depicts the predicted time-mean pressure coefficient distributions along the edge of wedge using the three grids. Results for the three grids agree very well. Differences in maximum pressure coefficient, Cp,max, among the three grids are graphically indiscernible. The locations where Cp,max occurs, ymax , are 0.466, 0.454, and 0.451 from the coarse to fine grid, indicating a converged solution. Note the time-mean value is obtained by averaging values for a carefully chosen duration as explained below. Before the time when ystill  yapex = 0.2, the simulated flow near the water surface significantly deviates from the similarity solution. We argue that combined effects of air cushion, air–water interaction, and boundary-layer thickness are strong enough to spoil the similarity at this early stage. On the other hand, around at the time when ystill  yapex = 0.6, water has wetted the full edge of wedge, so similarity

Fig. 7. Definition of problems of uniform flow over a fixed circular cylinder beneath a free surface. U: Inflow velocity, g: gravitational acceleration, D: cylinder diameter, h: distance between still fluid surface and cylinder top, qi: mass density of the ith fluid phase, li: dynamic viscosity of the ith fluid phase.

479

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

2

CD

1

Force Coefficient

Force Coefficient

2

0

CL

-1

CD

1

0

-1

CL

h/D =

h/D = 0.5 -2

-2 60

80

100

120

140

160

0

20

40

t

CD

Force Coefficient

1

0

-1

CL

0

-1

CL

h/D = 0.4

-2

-2 0

20

40

60

80

100

0

20

40

t

60

80

100

t 2

2

CD

1

Force Coefficient

Force Coefficient

100

CD

1

h/D = 1

0

-1

CL

1

CD

0

-1

CL

h/D = 0.8

h/D = 0.3

-2

-2 0

20

40

60

80

100

0

20

40

t

60

80

100

t

2

2

CD

1

Force Coefficient

Force Coefficient

80

2

2

Force Coefficient

60

t

0

-1

CL

1

CD

0

-1

CL

h/D = 0.6

h/D = 0.2 -2

-2 0

20

40

60

t

80

100

0

20

40

60

80

100

t

Fig. 9. Calculated force coefficients as function of time using the medium grid for the problem of uniform flow over a fixed circular cylinder beneath a free surface with selected values of h/D.

480

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

1.2

St/St 0

1

Reichl et al.

0.8

Present

0.6 0

1

2

3

h/D Fig. 10. Comparison of normalized Strouhal number as function of the gap ratio between the present results with those of Reichl et al. [53]. St0: Strouhal number for h/D = 1.

Fig. 12. Calculated water surface profiles of the side jet when ystill  yapex = 0.7 using three grids for the water-entry problem of a wedge with a = 30°.

4 3

10 2

9

Coarse Medium

8

1

Fine 7

0

6 5

-1

4 -2 3 -3

2 1 0 -1

-4 -4

-0.5

0

0.5

1

Normalized vertical coordinate y* = y/(ystill-yapex)

-3

-2

-1

0

1

2

3

4

Fig. 13. Fine-grid calculated water surface profile and sub-grid block distribution when ystill  yapex = 0.9 for the water-entry problem of a wedge with a = 30°. Each sub-grid block consists of 8  8 mesh cells.

Fig. 11. Calculated time-mean pressure coefficient distributions along the edge of wedge with a = 30° using three grids for the water-entry problem of a wedge.

cannot sustain in the ensuing flows. Therefore, we choose the duration of ystill  yapex = 0.2–0.6 for taking time-mean values. Fig. 12 presents the calculated water surface profiles of the side jet when ystill  yapex = 0.7 using the three grids for the case of a = 30°. Major differences exist only in the tip area of the side jet. The side jet stretches more and has larger jet angle (measured from the x axis) for higher grid density. The reason for weak grid convergence in predicting this local flow structure is the same as that stated in the dam-breaking simulation. That is, as time advances, finer and finer flow structures make it impossible to achieve grid convergence by employing a single set of grids. However, global flow characteristics has grid-converged solution. Fig. 13 shows the fine-grid calculated water surface profile and sub-grid block distribution when ystill  yapex = 0.9 for the case of a = 30°. Affording computational resources, we used the fine grid to perform all the other simulations. Fig. 14 reports the predicted timemean pressure coefficient distributions along the edge of wedge

for a = 30°, 45°, and 60°, compared with similarity solutions [54,55] and a fully nonlinear numerical simulation [56]. Results agree very well for each dead-rise angle. In particular, for the slamming case of a = 30°, both Cp,max and ymax are highly consistent among the present simulation, the similarity solution [55], and the fully nonlinear numerical simulation [56]. 4.5. Water entry of a circular cylinder at constant speed Water entry of a circular cylinder is a problem stiffer to tackle compared with water entry of a wedge because the dead-rise angle is zero upon the cylinder bottom contacting the initial water surface. The computational domain is [2, 2]  [14, 2] with the initial water surface and cylinder bottom being located at y = 0 and 0.25 respectively. The no-slip condition is imposed on the bottom and top boundaries and the flow is periodic in the horizontal (x) direction. We performed simulations at Fr = 2.5 and Re = 1  104 based on the cylinder diameter and entry speed. The ratios of fluid density and dynamic viscosity are q2/q1 = 815.7 and l2/l1 = 63.9

481

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

7

10

Present (Fine grid) Present (Medium grid) Present (Coarse grid)

Present

9

6

Similarity (Wu et al.) 8

Fully nonlinear (Mei et al.) 5

C S = F/ V 2R

7 6 5 4

4

3

2

3 2

1

1 0 -1

0 -0.5

0

0.5

1

Present Similarity (Hughes) 3

2

1

-0.5

0

0.5

1

2 Present Similarity (Wu et al.)

1

0 -1

0.05

0.1

0.15

0.2

0.25

0.3

Fig. 15. Calculated slamming coefficient as function of normalized submergence depth using three grids for the water-entry problem of a cylinder. F: Total vertical force, excluding the nominal buoyancy force, exerted to the cylinder. V: Waterentry speed. R: Cylinder radius. S: Submergence depth.

Similarity (Wu et al.)

0 -1

0

S/R

4

-0.5

0

0.5

1

Normalized vertical coordinate y* = y/(ystill-yapex) Fig. 14. Pressure coefficient distributions along the edge of wedge for a = 30°, 45°, and 60° for the water-entry problem of a wedge. Comparison between the present numerical simulation and the similarity solution [55]. A fully nonlinear numerical simulation [56] and another similarity solution [54] are referred to additionally for a = 30° and 45° respectively.

respectively. To investigate the effect of grid resolution, three grids are used with hfs = 1/96, 1/128, and 1/256 respectively. Fig. 15 presents calculated slamming coefficient as function of normalized

submergence depth using the three grids. The results differ little among the three grids in the course of slamming except in the very beginning (S/R < 0.04). Fig. 16 illustrates water surface profiles when S/R = 0.5 using the three grid resolutions, showing discernible differences only in the tip region of the splash jet. We compare the fine-grid predicted slamming coefficient with those measured by previous experiments [57,58] in Fig. 17. Note that Cointe and Armand [58] experimented with two entry speeds. In the very early stage (S/R < 0.03), the four results significantly differ from one another. After S/R = 0.03, the prediction lies between the data of Campbell and Weynberg [57] and Cointe and Armand [58], closer to the latter. Finally, although no validation can be made, we provide a snapshot of the fine-grid predicted flow field when S/ R = 2.3 for reference (Fig. 18). We observe that an air cavity zone forms behind the cylinder and the splash jet drastically stretches at this instant. Meanwhile, the distribution of sub-grid blocks demonstrates a pronounced effectiveness of the present AMR technique in clustering mesh cells near the free-surface or highvorticity regions. The fine-grid simulation spent 19.5 h to advance from t = 0 to t = 1.6, totally 1019 time steps consisting of 80 time steps with single outer iteration, two time steps with two outer iterations, 535 time steps with three outer iterations, 356 time steps with four outer iterations, and 46 time steps with P5 outer iterations. The number of mesh cells ranges approximately from 177,000 to 288,000. Accordingly, the CPU time consumption is about 291 ls per cell per time step. Most of simulations in this work were performed on a 2.59 GHz dual core AMD Opteron™ Processor 285 with RAM of 3.25 GB under the operating system of Microsoft Windows XP Professional Version 2002. The code was written in Fortran 90/95 compiled by Compaq Visual Fortran Professional Edition 6.6.0 with all real numbers declared as double precision. 4.6. Rolling motion of a circular-cylinder pair We demonstrate by this test case the capability of the present method to deal with moving multi-body problems. The circularcylinder pair consists of two identical circular cylinders the centers of which are rigidly connected. The two cylinders do not rotate but are permitted to revolve around the fixed midpoint of the centerto-center line, called the revolution center. Taking the cylinder

482

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

7

1

(a) 6

0.6

5

C S = F/ V 2R

0.8

Y

0.4 0.2 0

Campbell and Weynberg Cointe and Armand (V=2.33m/s) Cointe and Armand (V=7.38m/s) Present (Fine grid)

4

3

2

-0.2

1

-0.4 0 0 -0.6 -0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

X

0.05

0.1

0.15

0.2

0.25

0.3

S/R

0.8

Fig. 17. Slamming coefficient as function of normalized submergence depth for the water-entry problem of a cylinder. Comparison between the present numerical simulation and previous experimental measurements of Campbell and Weynberg [57] and Cointe and Armand [58].

1

(b) 0.8 0.6

Y

0.4 0.2 0 -0.2 -0.4 -0.6 -0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

X 1

(c) 0.8 Fig. 18. Fine-grid predicted flow field when S/R = 2.3 for a circular cylinder entering a still water. Thick solid black line: water surface. Flood contour: vorticity strength. Thin solid black line: boundaries of sub-grid blocks. Each sub-grid block has 8  8 mesh cells.

0.6

Y

0.4

pffiffiffiffiffiffi diameter, D, as the characteristic length and gD the characteristic velocity (hence Fr = 1), we prescribe the nondimensional oscillatory revolving motion as

0.2 0

aðtÞ ¼ a0 sinðxtÞ;

-0.2 -0.4 -0.6 -0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

X Fig. 16. Water surface profiles when S/R = 0.5 for a circular cylinder entering a still water using three grid resolutions: (a) coarse, (b) medium and (c) fine.

ð40Þ

where a denotes the angle from the positive x-axis to the center-tocenter line, a0 the amplitude of oscillation, and x the angular frequency of oscillation. The computational domain is [8, 8]  [14, 2] with y = 0 representing the initial water surface and the revolution center being located at (0, 0). The distance between the cylinder centers is two and we simulated the case with a0 = p/12 and x = p/2. The no-slip condition is imposed on the bottom and top boundaries and the flow is periodic in the horizontal

483

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

2

2

0

1

2

1

0

0

-1

-1

-2

-2 -4

-3

-2

-1

0

1

2

3

4

2

-4

-3

-2

-1

0

1

2

3

4

-2

-1

0

1

2

3

4

-2

-1

0

1

2

3

4

-2

-1

0

1

2

3

4

2

0.5

1

2.5

1

0

0

-1

-1

-2

-2 -4

-3

-2

-1

0

1

2

3

4

2

-4

-3

2

1

1

3

1

0

0

-1

-1

-2

-2 -4

-3

-2

-1

0

1

2

3

4

2

-4

-3

2

1.5

1

3.5

1

0

0

-1

-1

-2

-2 -4

-3

-2

-1

0

1

2

3

4

-4

-3

Fig. 19. Water surface profiles (thick black curves) and sub-grid blocks layout at selected times relative to the starting time of the second period for the cylinder-pair problem. Each sub-grid block has 2  2 mesh cells. Time is shown in the upper left corner of each plot.

trates water surface profiles and sub-grid blocks layout at times selected from the second period. Fig. 20 presents transient moments due to pressure exerted on individual cylinders and the cylinder pair as a whole. Note that this problem can be regarded as an elementary model to study the rolling motion of a twin-hull watercraft.

0.7 0.6

Right cylinder

0.5 0.4 0.3

Cylinder pair

0.2

4.7. Rotational motion of a partially submerged flower-shaped cross section

Mp

0.1 0 -0.1 -0.2

Further, in this test case, we will show that the present method can treat complicated geometry as well. The flower-shaped cross section is defined by a trajectory equation in polar coordinates,

Left cylinder

-0.3 -0.4

rðhÞ ¼ 0:5 þ 0:2 sinð5hÞ

-0.5 -0.6 -0.7

0

1

2

3

4

5

6

7

8

t Fig. 20. Transient moments due to pressure (hydrostatic part included) exerted on individual cylinders and the cylinder pair as a whole.

(x) direction. The ratios of fluid density and dynamic viscosity are q2/q1 = 815.7 and l2/l1 = 63.9 respectively. At Re = 1  104, we performed the simulation to t = 8 on a grid with hfs = 1/64. Fig. 19 illus-

ð41Þ

and rotates about its center with a dimensionless angular speed of p/10. The center is fixed at (0, 0.4) and y = 0 depicts the initial stillwater surface. Physical settings and numerical parameters are the same as those for the cylinder-pair case except that we performed the simulation to t = 40 in the domain [2, 2]  [2, 2] using a grid with hfs = 1/128. Fig. 21 illustrates water surface profiles and subgrid blocks layout at selected times. Clearly, the present method successfully handles chopping and merging of fluid–fluid interfaces. Also note the left-going waves generated by petal slamming the water surface at t 34 and 38. Fig. 22 presents transient moments

484

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

1

1

33

37

0

0

-1

-1

-2

-2 -2

-1

0

1

2

-2

-1

0

1

2

-1

0

1

2

-1

0

1

2

-1

0

1

2

1

1

34

38

0

0

-1

-1

-2

-2 -2

-1

0

1

2

-2

1

1

35

39

0

0

-1

-1

-2

-2 -2

-1

0

1

2

-2 1

1

36

40

0

0

-1

-1

-2

-2 -2

-1

0

1

2

-2

Fig. 21. Water surface profiles (thick black curves) and sub-grid blocks layout at selected times for the partially submerged rotational flower-shaped cross section. Each subgrid block has 2  2 mesh cells. Time is shown in the upper left corner of each plot. The arrow identifies a particular petal.

due to pressure exerted on the cross section during the interval t = 32–40. 5. Conclusions The present adaptive Cartesian cut-cell/level-set method has been proved to be capable of efficiently simulating two dimensional unsteady viscous incompressible free-surface/multi-phase flows with moving rigid bodies in laminar regime. For partially submerged moving bodies, we developed an outer-iteration solution procedure which couples solutions of the evolutionary equation of the level set function and the Navier–Stokes equations. This fully implicit solution procedure avoids use of the previoustime-step cell-face fluxes which generally do not exist for some cell

faces close to a moving solid boundary. Meanwhile, a physically reasonable interpolation scheme is used to construct values of variables at newly exposed cell centers. First-order accuracy in the spatial discretization is observed in the validating problem of uniform flow over a fixed circular cylinder beneath a free surface. The simulation results of water entry problems demonstrate the capability of the present method in handling drastically deforming fluid–fluid interfaces which intersect moving bodies, especially in terms of predicting hydrodynamic forces and pressure distributions exerted on solid bodies. The present method can treat multiple moving bodies and complicated geometry as shown in the last two test cases. In summary, the major contribution of this study is the first-time announcement of a well validated level set method based on the cut-cell approach in moving-body problems.

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

0.01

MP

0

-0.01

-0.02 32

33

34

35

36

37

38

39

40

t Fig. 22. Transient moments due to pressure (hydrostatic part included) exerted on the flower-shaped cross section during the interval t = 32–40.

Acknowledgements We would like to thank the National Science Council of Taiwan, ROC for financially supporting this research under Contract No. NSC 99-2221-E-022-021. Mr. Yao-Hsien Chang is appreciated for his assistance in digitizing curve data in literature.

References [1] Ye T, Mittal R, Udaykumar HS, Shyy W. An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries. J Comput Phys 1999;156:209–40. [2] Udaykumar HS, Mittal R, Rampunggoon P, Khanna A. A sharp interface Cartesian grid method for simulating flows with complex moving boundaries. J Comput Phys 2001 11/20;174(1):345–80. [3] Kirkpatrick MP, Armfield SW, Kent JH. A representation of curved boundaries for the solution of the Navier–Stokes equations on a staggered threedimensional Cartesian grid. J Comput Phys 2003;184(1):1–36. [4] Chung MH. Cartesian cut cell approach for simulating incompressible flows with rigid bodies of arbitrary shape. Comput Fluids 2006;35(6):607–23. [5] Peskin CS. Numerical analysis of blood flow in the heart. J Comput Phys 1977;25:220. [6] Goldstein D, Handler R, Sirovich L. Modeling of a no-slip surface with an external force field. J Comput Phys 1993;105:354. [7] Kim J, Kim D, Choi H. An immersed boundary finite volume method for simulations of flow in complex geometries. J Comput Phys 2001;171:132–50. [8] Yang J, Balaras E. An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries. J Comput Phys 2006;215:12–40. [9] Miyata H, Sato T, Baba N. Difference solution of a viscous flow with freesurface wave about an advancing ship. J Comput Phys 1987;72:393–421. [10] Harlow FH, Welch JE. Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface. Phys Fluids 1965;8:2182–9. [11] Unverdi SO, Tryggvason G. A front-tracking method for viscous, incompressible, multi-fluid flows. J Comput Phys 1992;100:25–37. [12] Tryggvason G, Brunnen B, Esmaeli A, Juric D, Al-rawahi N, Han W, et al. A Front-tracking next term method for the computations of multiphase flow. J Comput Phys 2001;169:708–59. [13] Hirt CW, Cook JL, Butler TD. A Lagrangian method for calculating the dynamics of an incompressible fluid with free surface. J Comput Phys 1970;5:103–24. [14] Bach P, Hassager O. An algorithm for the use of the Lagrangian specification in Newtonian fluid mechanics and applications to free-surface flow. J Fluid Mech 1985;152:173–90. [15] Ramaswamy B, Kawahara N. Lagrangian finite element analysis applied to viscous free surface flow. Int J Numer Meth Fluids 1987;7:953–84. [16] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 1981;39:201–25. [17] Kothe DB, Mjolsness RC. Ripple: a new model for incompressible flows with free surfaces. AIAA J 1992;30:2694–700. [18] LaFaurie B, Nardone C, Scardovelli R, Zaleski S. Modelling merging and fragmentation in multiphase flows with SURFER. J Comput Phys 1994;113:134–47.

485

[19] Rider WJ, Kothe DB. Reconstructing volume tracking. Technical report LA-UR96-2375, Los Alamos National Laboratory; 1996. [20] Lemos CM. Higher-order schemes for free surface flows with arbitrary configurations. Int J Numer Meth Fluids 1996;23:545–66. [21] Rudman M. Volume-tracking methods for interfacial flow calculations. Int J Numer Meth Fluids 1997;24:671–91. [22] Beckermann C, Diepers H-J, Steinbach I, Karma A, Tong X. Modeling melt convection in phase-field simulations of solidification. J Comput Phys 1999;154:468. [23] Provatas N, Goldenfeld N, Dantzig J. Adaptive mesh refinement computation of solidification microstructures using dynamic data structures. J Comput Phys 1999;148:265. [24] Tonhardt R, Amberg G. Dendritic growth of randomly oriented nuclei in a shear flow. J Cryst Growth 2000;213:161. [25] Lan CW, Liu CC, Hsu CM. An adaptive finite volume method for incompressible heat flow problems in solidification. J Comput Phys 2002;178:464–97. [26] Osher S, Sethian JA. Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulation. J Comput Phys 1988;79(1):12–49. [27] Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114:146–59. [28] Sussman M, Fatemi E, Smereka P, Osher S. An improved level set method for incompressible two-phase flows. Comput Fluids 1998;27:663–80. [29] Sussman M, Fatemi E. An efficient, interface-preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow. SIAM J Sci Comput 1999;20:1165–91. [30] Jiang G-S, Peng D. Weighted ENO schemes for Hamilton–Jacobi equations. SIAM J Sci Comput 2000;21:2126–43. [31] Enright D, Fedkiw R, Ferziger J, Mitchell I. A hybrid particle level set method for improved interface capturing. J Comput Phys 2002;183(1):83–116. [32] Sussman M, Almgren AS, Bell JB, Colella P, Howell LH, Welcome ML. An adaptive level set approach for incompressible two-phase flows. J Comput Phys 1999;148(1):81–124. [33] Sussman M, Puckett EG. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows. J Comput Phys 2000;162(2):301–37. [34] 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(1):110–36. [35] Kleefsman KMT, Fekken G, Veldman AEP, Iwanowski B, Buchner B. A volumeof-fluid based simulation method for wave impact problems. J Comput Phys 2005;206(1):363–93. [36] Qian L, Causon DM, Mingham CG, Ingram DM. A free-surface capturing method for two fluid flows with moving bodies. Proc Math Phys Eng Sci 2006;462(2065):21–42. [37] Wang W, Wang Y. An improved free surface capturing method based on cartesian cut cell mesh for water-entry and -exit problems. Proc Roy Soc A: Math Phys Eng Sci 2009;465(2106):1843–68. [38] Lin P. A fixed-grid model for simulation of a moving body in free surface flows. Comput Fluids 2007;36(3):549–61. [39] Yabe T, Xiao F, Utsumi T. The constrained interpolation profile method for multiphase analysis. J Comput Phys 2001;169:556–93. [40] Yokoi K. Numerical method for a moving solid object in flows. Phys Rev E 2003;67:045701. [41] Yabe T, Takizawa K, Chino M, Imai M, Chu CC. Challenge of CIP as a universal solver for solid liquid and gas. Int J Numer Meth Fluids 2005;47:655–76. [42] Yang J, Stern F. Sharp interface immersed-boundary/level-set method for wave–body interactions. J Comput Phys 2009;228(17):6590–616. [43] Bai W, Mingham CG, Causon DM, Qian L. Finite volume simulation of viscous free surface waves using the Cartesian cut cell approach. Int J Numer Meth Fluids 2010;63(1):69–95. [44] Zhang YL, Zou QP, Greaves D, Reeve D, Hunt-Raby A, Graham D, et al. A level set immersed boundary method for water entry and exit. Commun Comput Phys 2010;8(2):265–88. [45] Chung MH. Numerical study of rowing hydrofoil performance at low Reynolds numbers. J Fluids Struct 2008;24(3):313–35. [46] Iafrati A, Olivieri A, Pistani F, Campana EF. Numerical and experimental study of the wave breaking generated by a submerged hydrofoil. In: 23rd ONR symposium naval hydrodynamics, Val de Reuil (France); September 2000. [47] Peng D, Merriman B, Osher S, Zhao H, Kang M. A PDE based fast local level set method. Internal report CAM 98-25, Department of Mathematics, UCLA; 1998. [48] Zang Y, Street RL, Koseff JR. A non-staggered grid, fractional step method for time-dependent incompressible Navier–Stokes equations in curvilinear coordinates. J Comput Phys 1994;114:18–33. [49] Saad Y, Schultz MH. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J Sci Stat Comput 1986;7:856–69. [50] Almgren AS, Bell JB, Colella P, Howell LH, Welcome ML. A conservative adaptive projection method for the variable density incompressible Navier– Stokes equations. J Comput Phys 1998;142(1):1–46. [51] Martin JC, Moyce WJ. An experimental study of the collapse of fluid columns on a rigid horizontal plane. Philos Trans Roy Soc Lond: Ser A 1952;244:312–24. [52] Roache PJ. Verification and validation in computational science and engineering. Albuquerque (New Mexico): Hermosa Publishers; 1998. [53] Reichl P, Hourigan K, Thompson MC. Flow past a cylinder close to a free surface. J Fluid Mech 2005;533:269.

486

M.-H. Chung / Computers & Fluids 71 (2013) 469–486

[54] Hughes OF. Solution of the wedge entry problem by numerical conformal mapping. J Fluid Mech 1972;56(01):173. [55] Wu GX, Sun H, He YS. Numerical simulation and experimental study of water entry of a wedge in free fall motion. J Fluids Struct 2004;19(3):277–89. [56] Mei XM, Liu YM, Yue DKP. On the water impact of general two-dimensional sections. Appl Ocean Res 1999;21(1):1–15.

[57] Campbell IMC, Weynberg PA. Measurement of parameters affecting slamming. Report no. 440, Wolfson Unit of Marine Technology, Tech. rep. centre no. OT-R8042, Southampton; 1980. [58] Cointe R, Armand JL. Hydrodynamic impact analysis of a cylinder. ASME J Offshore Mech Arct Eng 1987;109:237–43.