Simulation of incompressible multiphase flows with complex geometry using etching multiblock method

The incompressible two-phase flows are simulated using combination of an etching multiblock method and a diffuse interface (DI) model, particularly in the complex domain that can be decomposed into multiple rectangular subdomains. The etching multiblock method allows natural communications between the connected subdomains and the efficient parallel computation. The DI model can consider two-phase flows with a large density ratio, and simulate the flows with the moving contact line (MCL) when a geometric formulation of the MCL model is included. Therefore, combination of the etching method and the DI model has potential to deal with a variety of two-phase flows in industrial applications. The performance is examined through a series of numerical experiments. The convergence of the etching method is firstly tested by simulating single-phase flows past a square cylinder, and the method for the multiphase flow simulation is validated by investing drops dripping from a pore. The numerical results are compared with either those from other researchers or experimental data. Good agreement is achieved. The method is also used to investigate the impact of a droplet on a grooved substrate and droplet generation in flow focusing devices.


Introduction
Interfacial flows in the domain with the complex geometry are often encountered in nature and industrial applications. These flows are usually difficult for numerical simulations, partially because they may involve with intricate physical phenomena, such as the moving contact line (MCL) at which the gas, liquid, and solid meet. Examples are the striding of insects on the water surface [1] and the floating/sinking of solid objects [2] . It requires appropriate models for the interface and MCLs to perform the accurate computation. Moreover, the complex geometry also demands the generation of appropriate meshes, on which the equations of motion are solved. Traditionally, a body-fitted mesh is generated for this purpose, either structured [3] or unstructured [4] . Nowadays, because of its efficiency and simplicity, immerse boundary (IB) methods become more and more popular in dealing with flow problems in the presence of complex geometry [5][6][7] . In the IB methods, the solid boundaries can be represented by a set of linked Lagrangian points, while the flows are solved on the Eulerian mesh. Since the solid boundaries generally do not coincide with the mesh lines, a feedback force is added in the momentum equation to account for the presence of the embedded solid boundaries. Recently, we combined the IB method [8] with the diffuse interface (DI) model [9] , and successively simulated the flows with MCLs on curved substrates [10] . Despite of the success of IB methods, their accuracy depends on the way of representing the solid boundary, e.g., the use of delta function in the distribution of the feedback force makes this kind of IB methods essentially first-order accurate in modeling the solid boundaries.
Under certain circumstance, the complex geometry can be decomposed into several regular subdomains/blocks, of which each can be discretized by a structured mesh. Therefore, the governing equations can be solved in the individual block, and the communications between the blocks are exchanged by implementing appropriate boundary conditions (BCs). This way of dealing with the complex geometry is often referred to as the multiblock method. In the multiblock methods, we can still enjoy the merits of the structured mesh, such as easiness in constructing numerical schemes of high accuracy and simplicity in the parallel computation. Clearly, the main challenge in the multiblock approach lies in how to properly decompose the blocks and provide the BCs for each block. Wan et al. [11] used block-structured grids with fictitious overlapping interfaces to decompose a complex computational geometry. Raspo [12] solved the incompressible Navier-Stokes (NS) equations by a direct domain-decomposition method coupled with a Chebyshev collocation method. Abide and Viazzo [13] used a direct non-overlapping multidomain method to simulate the incompressible viscous flows, such as von Karman vortex street behind a square cylinder. All these works only consider single-phase flows. This fact motivates us to examine the performance of the multiblock method in simulating multiphase flows, particularly in the presence of MCLs.
In the present study, we use an etching multiblock method coupled with the DI model to investigate incompressible multiphase flows with complex regular geometries. The etching multiblock method defines a regular domain as the global domain, and marks all the embedded solid boundaries. As a result, the blocks can be considered as the patches on a global structured mesh, which naturally establishes the connections between the blocks. The fluid-fluid interface is represented by a DI model [14] . This model replaces the mathematically sharp interface by a DI with the finite thickness, where the fluids are mixed to some extent. The DI model allows for the simulation of two-phase flows with the large density ratio [9] , and has been used to a variety of flow problems, such as vesicle dynamics [15] , Hele-Shaw flows [16] , and head-on droplet collision [17] . We choose the geometric formulation of a DI model to model the motion of contact line on the substrate [18] . This model has been used to many two-phase flows in the presence of MCLs, including droplet spreading [19] , sliding of a three-dimensional droplet on a wall [20][21] , and ejection of satellite droplet in a rapid spreading [22] .
We present the modeling and numerical methods in Section 2, including the governing equations, MCL model, etching multiblock method, and BCs. A variety of numerical experiments are performed to validate the method in Section 3. We first simulate the flow past a square cylinder for the convergence study. Then, we investigate the drops dripping from a pore for code validation, and quantitatively compare the numerical results with experimental observations. We also investigate a drop impacting on a grooved substrate in Subsection 3.3 and droplet generation in flow focusing in Subsection 3.4.

DI model
We investigate here the flows of two incompressible immiscible liquids (A and B) and use a DI method [9] to model the interface between the two liquids. In the DI method, the mathematically Simulation of incompressible multiphase flows with complex geometry 1407 sharp interface is replaced by a DI with the finite thickness of a length scale . The DI is represented by the volume fraction C of the fluid A (0 C 1). Consequently, the time evolution of the interface can be governed by the convective Cahn-Hilliard (CH) equation, where u = (u, v) is the flow velocity, and ψ is the chemical potential defined by in which u and v are the velocity components in the xand y-directions, respectively, φ(C) = C 2 (1 − C) 2 /4 is the bulk energy density, the Cahn number Cn = /L is the dimensionless measure of the thickness of the DI, and L is the characteristic length. P e is the Péclet number that represents the relative importance of the convective fluxes (uC) to the diffusive fluxes (∇ψ). According to Ref. [9], we choose P e = 1/Cn so that the DI would approach the sharp interface limit with the mesh refinement. A geometric formulation of wetting condition [18] is used here to model the MCL, and serves as the BC for C at a solid substrate. In the geometric formulation, the local variation of C along the normal direction to the solid substrate should satisfy where n and τ are the normal vector and the tangential vector of the solid boundary, respectively, and θ s is the microscale contact angle.

NS equations
The equations of motion for the fluids are the dimensionless NS equations and continuity equation, where j denotes the unit vector in the vertical direction, and f s is the surface tension force, The dimensionless averaged density ρ and viscosity μ are defined as where the subscripts A and B denote the liquids A and B, respectively. The dimensionless group includes the Reynolds number Re = ρ A U L/μ A , the Weber number We = ρ A U 2 L/σ, and the Froude number F r = U 2 /(gL), where σ is the surface tension coefficient, g is the gravitational acceleration, and L and U are the characteristic length and velocity, respectively. We shall also refer to the Ohnesorge number Oh = √ We/Re and the Bond number Bo = ρ A L 2 g/σ in the discussion.

Etching method and discretization
We consider here to use an etching multiblock method in the simulation of incompressible viscous two-phase flows in a domain with the complex geometry, which can be decomposed into multiple rectangle subdomains (or blocks). For the convenience of coding, the computational domain is defined on a rectangle that contains the physical domain. In such a way, we can use a uniform Cartesian mesh to discretize the domain, and the boundaries of the physical domain coincide with the mesh lines. As a result, each flow variable such as u, v, p, and C in the domain is stored in one array, and thus can be accessed in the same manner as on a structured mesh, regardless of the complex geometry involved. Therefore, we can adopt the discretization method on the structured mesh to discretize the governing equations.
In the present study, the spatial discretization of the governing equations is performed on a staggered grid, as shown in Fig. 1. All the scalar quantities such as the volume fraction, chemical potentials, and pressure are defined at the cell centers, while the vector quantities such as the velocity components are defined at the cell faces. Unless otherwise mentioned, the central finite difference schemes are used for the spatial discretization of the CH equation and the NS equations. To solve the CH and NS equations in a temporally matched manner, the volume fraction C field is defined at the time levels t n− 1 2 , while the velocity u and the pressure p are defined at t n (n = 1, 2, 3, · · · ). Consequently, the CH equation (1) is discretized at the time level n as follows: where we set the time step Δt ∼ h 2 due to the explicit discretization of the second term in the next equation, and h is the mesh space. We use the Adams-Bashforth and Crank-Nicolson schemes for the discretization of M (u, C) and ∇ 4 C, respectively, where the term M (u, C) is defined as Note that the fluxes at the cell faces (the first term in Eq. (10)) are evaluated with a fifth-order weighted essentially non-oscillatory (WENO) scheme, using the local flow velocity to determine the upwind direction [9] .
A standard projection method is used to couple the momentum equation with the continuity equation. The equations are discretized at the time level (n + 1 2 ) as follows: where we use the Crank-Nicolson scheme in the discretization of the viscous term L = ∇ · (μ n+ 1 2 (∇u + ∇u T )), and the Adams-Bashforth scheme in the discretization of the advective term, H = (∇ · uu). Note that ρ n+ 1 2 The velocity u n+1 is eventually made divergence-free after the pressure correction, An over-relaxation iterative method is used to solve the Poisson equation in Eq. (12).

BCs
How to deal with BCs is the most important issue in the computation with the etching multiblock method. Here, we choose the case of drops dripping from a pore as an example to illustrate the solution procedure. The computational domain is shown in Fig. 2, in which the actual computational domain can be divided into two blocks, i.e., the rectangles AHF E and HBCG, respectively. There are four kinds of flow BCs at the edges of these blocks, i.e., symmetric (AH and HB), inlet (AE), outlet (CB and CG), and solid (EF and F G). There are also BCs particularly related to the etching multiblock method. For example, the edge F H is considered as the block boundary, since it separates the blocks a and b, although it is inside the computational domain. It is referred here to as the connection BC. A mixed BC is enforced at the boundary GH of the block b, which consists of two edges (GF and F H) with different BCs. During the computation, the governing equations are solved in the two blocks one by one, along with the update of the BCs at the boundaries of each block. In the present study, we use one layer of ghost cells around each block to help the implementation of the BCs (see Fig. 3). In particular, the values of the flow variables at the ghost cells are updated according to the local BCs. In the simulation of drops dripping from a pore, the BCs are specified as follows: at the left boundary (inlet), u = u in , v = 0, ∂p ∂x = 0, and C = C in , where u in and C in are the prescribed horizontal velocity and the volume fraction at the inlet, respectively; at the right boundary (outlet), we use a characteristic BC, which is an approximation of the far-field BC, so that all the generated drops can pass smoothly through the boundary. Specially, the characteristic BCs can be expressed as At the bottom boundary, a symmetric BC is enforced, which can be expressed as At the solid wall, the no-slip BC is used for the velocity, and the geometric wetting condition is used for C, i.e., In the implementation of the connection BC, there is no need to update the values at the ghost cells, which are essentially the inner cells. In this sense, we do not require any treatment at the connection BC. For the mixed BC, we only enforce the solid BC at the part of boundary at the solid wall, while leaving alone the part of boundary connected with other blocks. The details of the discretized BCs are listed as follows: (ii) Symmetric BC (iii) Inlet BC where all the flow variables are defined in Fig. 4, and the subscript i = 0 denotes the ghost cells. For simplicity and convenience, we refer the solid BC to as BC 1 , and similarly, the symmetric BC (BC 2 ), the inlet BC (BC 3 ), the characteristic BC (BC 4 ), the connection BC (BC 5 ), and the mixed BC (BC 6 ). In the following, we shall list the BCs for the block in a boundary sequence of left, bottom, right, and top. For example, for the block a in Fig. 2, the BCs can be represented by (BC 3 , BC 2 , BC 5 , BC 1 ).

Results and discussion
In all the tests, the liquid-gas interface is represented by the contour C L = 0.5, and the Cahn number is set to Cn = h.

Flow past square cylinder
The convergence of the etching multiblock method is tested by simulating a two-dimensional single-phase flow past a square cylinder inside a channel. The numerical setup is shown in Fig. 5, in which the computational domain of the size 8D × 50D is split into four blocks (a, b, c, and d), where D is the side length of the square cylinder. The BCs are (BC 3 , BC 1 , BC 6 , BC 1 ) for the block a, (BC 5 , BC 1 , BC 5 , BC 1 ) for the blocks b and c, and (BC 6 , BC 1 , BC 4 , BC 1 ) for the block d. It is assumed that the flow at the inlet is a fully developed laminar channel flow. Therefore, the prescribed velocity at the inlet u in has a parabolic velocity profile, u in = U max (−y 2 /(16D 2 ) + y/(2D)), where U max is the velocity in the middle of the inlet of the channel. The blockage ratio, i.e., the length of the cylinder side over the height of the channel, is 0.125. We choose D and U max as the characteristic length and velocity, respectively. For the flow past a square cylinder, it is known that there exists a steady recirculating flow region at the downstream of the square cylinder at a relatively low Re (0.5 Re 60). In contrast, the vortex shedding occurs at a relatively high Re (60 Re 300) [23] . We select the cases at Re = 40 and 100 to test the etching multiblock method. Figure 6 shows the numerical results in terms of streamlines. It is clear that steady flows are obtained at Re = 40, and unsteady flows with vortex shedding occur at Re = 100 (the time step Δt = 0.003 and the mesh is 256 × 1 600), consistent with the prediction in Ref. [23]. To quantitatively validate the method, we make a comparison of our results with those from other researchers [13,[23][24] in Table 1, where EXP, FVM, LBA, and PDM represent the experiment, the finite-volume method, the lattice-Boltzmann automata, and the projection decomposition method, respectively. In particular, we compare the length of the recirculating flow region (L r ) and the drag coefficient (C d ) at Re = 40, and the averaged drag coefficient (C d ) and the Strouhal number where F x is the drag force along the flow direction exerted by the surrounding fluid on the square cylinder. The present results on the successively refined meshes are listed in Table 1, i.e., 192 × 1 200, 256 × 1 600, and 320 × 2 000, to test the convergence of the method. Clearly, the present results are in agreement with those from other researchers, and they are also converged with the successive mesh refinement.

Drops dripping from pore
Due to the effect of gravity, a drop attached to the ceiling could drip down if its volume is beyond a critical value [25] . This process exhibits many interesting phenomena, such as genera-tion of droplets, formation of a jet, and MCLs. We consider here the axisymmetric simulation of the dripping process of the water that leaks out of a pore. The density and viscosity ratios of the two liquids, water and air, are ρ a /ρ w = 0.001 and μ a /μ w = 0.025, respectively, where the subscripts a and w denote air and water, respectively. Figure 7 shows the numerical setup, in which the domain has a size of 1.5D × 9D, where D is the diameter of the pore. The computational domain can be divided into two blocks (a and b), and the former has a size of 0.5D × D and the latter is 1.5D × 8D. The BCs are (BC 3 , BC 2 , BC 5 , BC 1 ) for the block a and (BC 6 , BC 2 , BC 4 , BC 4 ) for the block b. The simulation is performed on a Cartesian mesh of 1 800 × 300. We choose the liquid velocity U on the entry boundary and the diameter of the pore D as the characteristic velocity and length, respectively. The time step Δt = 0.000 1.  Fig. 7 Sketch of drops dripping from pore, where gray area shows initial space occupied by liquid, blank area is supposed to be filled with air, and direction of gravity is from left to right We simulate two cases with the same Re (= 4.34), Bo (= 0.33), and wettability θ (= 120 • ) but different We: We = 0.16 and 0.2. All the dimensionless numbers are defined by the use of the properties of the water. In the simulation, the liquid that leaks out of the pore gradually accumulates at the end of the pore, and due to the surface tension force, a drop comes into being and grows with the time. When the drop becomes sufficiently large, a droplet breaks up from the drop and drips down under the influence of gravity. During the whole dynamic process, the contact lines are pinned at the corner of the pore, primarily due to the geometry-induced contact-angle hysteresis [26][27] . Figure 8 shows the snapshots of the periodic dripping process from the pore. It is interesting to see that the difference in We leads to different dripping modes, i.e., the single-droplet mode and the double-droplet mode. More details of the dripping process are shown in Fig. 9, in which the time interval between the successive generation of droplets (t d ) is presented as a function of the sequence number of the droplets (N ). Clearly, the dripping process results in a periodic generation of droplets, and the Strouhal number St (= f D d /g, where D d is the diameter of the generated droplet) is 19.2 for the single-droplet mode and 23.6 and 20.1 for the double-droplet mode. Subramani et al. [28] conducted experiments with similar parameters, and they also found the two modes. Moreover, they observed that the volume of the generated droplet was 20.8 for the single-droplet mode, and our numerical result gives nearly the same value (20.7). For the double-droplet mode, the experimental data showed the volumes of the two droplets were 16.6 and 20.8, respectively, and our numerical results predict 18.7 and 21.9, which are in agreement with the experimental results.

Drop impacting on rough substrate
We consider axisymmetric simulations of drop impacting on a rough substrate. The flow phenomena of drop impacting on the smooth substrate have been extensively investigated, and a recent review can be found in Ref. [29]. We represent the rough substrate by a grooved substrate similar to the study of droplet motion inside a grooved channel by Huang et al. [30] . In the present work, we investigate the impact dynamics of a droplet on a grooved substrate (see Fig. 10). The diameter of the droplet is D, and the size of the computational domain is 1.1D × 2.5D (the mesh 440 × 1 000). The grooves are of uniform width and depth, i.e., 0.08D × 0.05D, which also corresponds to the size of the blocks a to e. There are totally ten blocks in the domain (see Fig. 10), and they can be classified into five groups according to their BCs. The BCs are (BC 2 , BC 1 , BC 1 , BC 5 ) for the block a, (BC 1 , BC 1 , BC 1 , BC 5 ) for the blocks b to e, (BC 2 , BC 6 , BC 5 , BC 4 ) for the block f , (BC 5 , BC 6 , BC 5 , BC 4 ) for the blocks g to i, and (BC 5 , BC 6 , BC 4 , BC 4 ) for the block j. The density and viscosity ratios of the air to the liquid are 0.001 and 0.025, respectively. We choose the impact velocity of the droplet V and D as the characteristic velocity and length, respectively, and also use the properties of the liquid to define the dimensionless group. The time step Δt = 0.000 025. We simulate two cases at Oh (= 0.003 5), i.e., Case 1: We = 12.5 and θ = 120 • , and Case 2: We = 32 and θ = 140 • . Figure 11 shows the snapshots of the droplet shapes at different time for the two cases. It is seen that in both cases, the droplet spreads on the substrate horizontally after the impact. However, the droplet in Case 2 also penetrates into the grooves and touches the bottom of some cavities during the spreading. In contrast, the droplet in Case 1 only forms the contact lines pinned at the corners of the grooves. As a result, the droplet in Case 1 rebounds from the substrate at the time t = 3.6, while the droplet in Case 2 generates a satellite droplet ejected into the air, leaving the rest of the liquid attached to the substrate. It is found that the inertia of the droplet plays an important role in the impact dynamics, and may lead to the change of the wetting state on the rough substrate, e.g., partial transition from the Cassie state to the Wenzel state here [31] . It is also clear that the etching multiblock method can successfully capture the contact line pinned at the corners of the grooves, which has been observed in previous experiments [32] .

Flow focusing
The last application is the droplet generation by flow focusing, in which the liquid A in the capillary tube is driven by the liquid B through a pore (see Fig. 12). Flow focusing has been shown to be an effective way to generate uniform and monodispersed droplets [33] . However, because the geometry of the flow focusing device is rather complicated, and the flows involve two-phase flows and MCLs, there are only a few numerical studies of flow focusing phenomena, mainly using commercial softwares [34][35][36] .  Figure 12 shows the numerical setup for axisymmetric simulations. The computational domain is 10D × 2D (the mesh is 1 000 × 200), where D is the inner diameter of the capillary tube. There are five blocks (a, b, c, d, and e), and the BCs of these blocks are (BC 3 , BC 2 , BC 5 , BC 1 ), (BC 3 , BC 1 , BC 5 , BC 1 ), (BC 6 , BC 2 , BC 6 , BC 1 ), (BC 5 , BC 2 , BC 5 , BC 1 ) and (BC 6 , BC 2 , BC 4 , BC 4 ), respectively. We choose the average velocity at the inlet of the capillary tube U A and D as the characteristic velocity and length, respectively. The effect of gravity is considered to be negligible. The properties of the liquid A are used to define the dimensionless group: Re and We. The density and viscosity ratios of the liquid B to the liquid A are 1.03 and 0.2, respectively. The time step Δt = 0.000 025.
We investigate here two cases with Re = 1.7, We = 0.009, and θ = 120 • but different average velocities at the inlet of the fluid B, particularly U B = 0.92 and 1.26, respectively. Figure 13 shows the numerical results of two cases in terms of the interface shape. It is interesting to see that the difference in the flow rate of the liquid B leads to the occurrence of different flow modes. The case with U B = 0.92 appears to enter a dripping mode, in which droplets of uniform size are generated near the pore, and the cone of the liquid A oscillates periodically. The case with U B = 1.26 appears to enter a jetting mode, in which a slender liquid jet is formed at the downstream of the pore, and the shape of the liquid cone is kept unchanged with the time. The second case also continuously generates droplets, which are smaller than those in the first case. We note that the dripping and jetting modes are in consistence with the experimental observations [37] .

Conclusions
We combine the etching multiblock method with the DI method and the geometric formulation of MCL model to simulate incompressible two-phase flows in the presence of a complex geometry. The complex geometry can be decomposed into multiple blocks of regular shapes such as rectangles and squares. One layer of ghost cells is placed around each block, for the convenience of implementing the BCs. The BCs include the physical ones, such as the MCL model for the volume fraction, and those for the connections between the blocks. The latter shall require no treatment at all. As a result, the etching multiblock method can be extended from the single-phase flows to the multiphase flows in a straightforward manner. The etching multiblock method is shown to be able to effectively deal with the incompressible two-phase flows in a domain with multiple blocks. Its performance is examined through a few numerical experiments. The simulation of single-phase flow past a square cylinder shows that the results converge with the mesh refinement. The numerical results of drops dripping from a pore are in good agreement with the experimental results qualitatively and quantitatively. Finally, we use the method to the droplet generation in flow focusing and drop impacting on the rough substrate.