A 3D conservative sharp interface method for simulation of compressible two-phase flows

Journal Pre-proof A 3D conservative sharp interface method for simulation of compressible two-phase flows

Yi Shen, Yi Ren, Hang Ding

PII:

S0021-9991(19)30812-5

DOI:

https://doi.org/10.1016/j.jcp.2019.109107

Reference:

YJCPH 109107

To appear in:

Journal of Computational Physics

Received date:

15 April 2019

Revised date:

4 November 2019

Accepted date:

5 November 2019

Please cite this article as: Y. Shen et al., A 3D conservative sharp interface method for simulation of compressible two-phase flows, J. Comput. Phys. (2019), 109107, doi: https://doi.org/10.1016/j.jcp.2019.109107.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

Highlights • • • •

3D sharp interface method for compressible flows with topology change of interface. Cell-splitting configurations decrease from 288 configurations to 12 basic ones. Conservation laws are well maintained, even at interfaces with topology changes. Numerical results are validated against various experimental ones.

A 3D conservative sharp interface method for simulation of compressible two-phase flows Yi Shen, Yi Ren, Hang Ding∗ Department of Modern Mechanics, University of Science and Technology of China, Hefei 230027, China

Abstract A three-dimensional (3D) conservative sharp interface (CSI) method, which combines 3D cut-cell methods and finite volume methods in the arbitrary Lagrangian-Eulerian framework, is proposed to simulate compressible inviscid two-phase flows. The 3D cut-cell method allows to reconstruct interfaces on a Cartesian mesh and update unstructured interface cells in the vicinity of the interface, making the reconstructed interface always coinciding with the cell faces of the interface cells. The complexity in generating 3D interface cells is significantly reduced by exploiting the symmetry and rotational symmetry of fluid positions in a cut cell. Compared to our previous two-dimensional work (Lin et al. J. Comput. Phys. 328:140159, 2017), the occurrence of large interface cells is also avoided by using the new strategy of cell assembly. With the help of the 3D cut-cell method, a second-order finite volume method in the arbitrary Lagrangian-Eulerian framework can be used nearly in the whole domain, except for the interface cells with under-resolved interfaces or topology changes, at which a first-order finite volume method is used instead. The convergence of the 3D CSI method is examined by studying the Rayleigh collapse of a gas bubble in water, and secondorder accuracy has been obtained for the interface evolution. Simulation of shock-induced bubble collapse demonstrates that the axisymmetricity of the interface is well preserved in the 3D computation. Moreover, the 3D CSI methods are also given experimental validation, including the interaction of a shock with 3D curved gas-gas interface and droplet deformation after shock impact. In all cases, a good agreement has been achieved with the experimental data, with respect to the interface evolution and flow features. Keywords: compressible two-phase flows; sharp interface; 3D cut-cell method.

1. Introduction Compressible multiphase flows are observed in a variety of industrial applications, such as shockwave lithotripsy [1], ultrasound contrast agents [2], raindrop damage in supersonic flight [3] and fuel injection in internal combustion [4]. The flow features are usually characterized by the presence of interface instabilities, global interface deformation and local ∗

Corresponding author. Email address: [email protected]

Preprint submitted to Journal of Computational Physics

November 11, 2019

topology changes of interfaces (e.g. interface breakup and coalescence). The interaction between the interfaces and the surrounding flows such as shocks has significant effect on the interface dynamics, which in turn affects the flow behaviors. To gain a full understanding of the dynamic process of these flow phenomena, it is necessary to use direct numerical simulations to study time evolution of flow quantities, of which many may not be accessible in experiments. However, accurately resolving interfacial waves of small scale poses big challenges for numerical simulations. In particular, it is desirable to have high resolution in the interface region and high-order approximation of jump conditions across the interface. With this regard, an ideal model to deal with interface is sharp interface method. The sharp interface method is referred here to as a group of numerical models that are able to treat the interfaces in a mathematically sharp manner. More specifically, the interface is assumed to have zero thickness, and the discontinuity in the material properties (e.g. density) is satisfied within one mesh spacing in the vicinity of the interface. We mainly focus here on the conservative sharp interface (CSI) methods, in which the interfaces coincide with the cell faces all the time during the computation, even in the presence of topology changes. To achieve this goal, cut cell methods are often employed to dynamically generate and update the mesh in the vicinity of the interfaces [5]. The main idea of the cut cell methods is: At each time step the interfaces (which is represented here by a level set function φ = 0) are projected onto a fixed structured mesh (e.g. Cartesian mesh) and numerically reconstructed by line segments in the two-dimensional (2D) cases or polygons in the three-dimensional (3D) cases; the regular structured cells are then split by these line segments or polygons, and the unstructured cells that share the cell faces with the interfaces, namely interface cells (ICs), can be generated by assembling the split cells. Cutcell methods are originally designed for simulation of flows involving static rigid bodies with complex geometry on a uniform background mesh, which is split by the boundaries of the rigid bodies instead of the interfaces. Some techniques in the cut cell methods have been well established, such as cell merging [6], flux redistribution [7] and dimension-splitting [8]. With the help of cut cell methods, the interfaces effectively serve as the boundaries that separate the two fluids, with the unstructured ICs in the vicinity of the interface and the regular structured cells in the bulk fluids away from the interface. On the basis of the cut-cell meshing strategy, one can easily implement finite volume methods in the arbitrary Lagrangian-Eulerian (ALE) framework, so that mass, momentum and energy are strictly conserved in a discrete manner, even in the presence of moving boundaries such as interfaces. In addition, the boundary conditions (or jump conditions) at the interface can be resolved by constructing and solving a local Riemann problem at the corresponding cell faces. Note that for moving-mesh ALE approaches, when the mesh undergoes large distortions, a remeshing process is required based on the given boundary segments [9]. Accordingly, a remap technique is used to update the flow variables on the new mesh, preferably in a conservative manner [10]. How to handle the presence of topology changes in geometry (interfaces or structures) is the major challenge for the moving-mesh ALE approach, especially in 3D simulations. By contrast, in the cut-cell ALE algorithm, distorted meshes are circumvented by interface reconstruction and cell assembly (in each time step), even in the presence of topology changes. In this sense, the cut-cell based ALE methods are more general and 2

robust than the moving-mesh ALE approaches. The cut-cell meshing strategy proves to be the most important part of the CSI method. Various cut-cell meshing strategies and corresponding numerical schemes for spatial discretization have been proposed for the simulation of compressible two-phase flows, with particular focus on the treatment of small split cells and the ICs where topology changes occur. Hu et al. [11] proposed to adopt numerical mixing in the ICs, to circumvent the limitation of time step due to small split cells and difficulties in dealing with the topology changes of interfaces. As a result, they essentially used a first-order finite volume scheme for the spatial discretization in the ICs. Lauer et al.[12] extended this method to the three dimension, and used it to study the spherically and non-spherically symmetric bubble collapse. To further improve the robustness of the method, a scale separation approach was proposed to identify the resolved interface structures and the non-resolved ones (likely occur in the presence of topology changes or the interfaces with high curvature), which will be modeled differently in the flow computation [13, 14]. Recently, Hu’s method [11] has also been extended to the 2D simulation of compressible ternary-fluid flows [15]. Cut cell method with the numerical mixing is robust in computation, and simple in implementation. On the other hand, the numerical mixing and first-order spatial discretization in ICs may cause the damage to symmetry preservation of the interface, especially when the local interface resolution is relatively low or a long-time interface evolution is encountered. Chang et al. [16] proposed a second-order 2D CSI method for multiphase flows in the absence of topology changes, in which all the computations on an IC are performed either at the centroids of the cell faces or at its cell center. The second-order CSI method for compressible two-phase flows was validated by simulating the deformation of a liquid drop after shock impact, and the numerical results were shown to agree well with their experiments. Lin et al. [17] proposed to use a delay reinitialization in the 2D second-order CSI methods, in order to suppress unphysical oscillations of flow variables in the presence of topology change. It was shown to be able to capture the complicated interface structures and topology changes in the shock-bubble interaction, in good agreement with the experiments under the same flow conditions. A 3D CSI approach based on cell-merging was proposed by Kim & Liou [18] for compressible multiphase flows. However, the approach cannot deal with interfaces with topology changes. Furthermore, the method has not been validated for compressible multiphase flows, although it was applied to the simulation of rising bubbles due to buoyancy, i.e. essentially incompressible two-phase flows. Despite the success of second-order 2D CSI methods, there are only a few second-order 3D CSI methods available so far for the simulation of compressible two-phase flows, probably due to the significant complexity in the 3D interface reconstruction and the manipulation of cell assembly. In the paper, we propose a 3D CSI method for the simulation of compressible inviscid two-phase flows, which can achieve second-order accuracy at the interfaces. In the 3D cutcell meshing strategy, we exploit the symmetry and rotational symmetry of fluid positions in a single Cartesian cell, and greatly simplify the cell-splitting configurations (decreasing from 284 configurations to 12 basic ones). This also significantly reduces the complexity in the manipulation of cell assembly. Compared to our previous 2D work [17], the cell assembly process is further optimized to reduce the IC size, so as to improve the interface resolution. 3

A second-order finite volume method in the ALE framework is generally used everywhere for the discretization of the Euler equations; the only exception is the IC that contains two or more interfaces (i.e. where the topology changes of interfaces are about to occur or just take place), at which a first-order finite volume method is implemented. The performance of the 3D CSI method is systematically examined. In particular, we study the convergence of the method by simulating the Rayleigh collapse of a bubble in water, and check the symmetry preservation by simulating the axisymmetric bubble collapse induced by a planar shock. More importantly, we numerically investigate the interaction between shock and 3D curved gas-gas interface and droplet deformation after shock impact, and the numerical results are directly compared against the experimental data. To our best knowledge, this is the first time that 3D CSI methods are given experimental validation for compressible two-phase flows. 2. Governing equations We consider here simulations of compressible inviscid two-phase flows, which are governed by 3D Euler equations. The Euler equations in the conservative form can be written as, ∂Q + ∇ · F = 0, ∂t where the conservative variables Q and the flow fluxes F are defined as ⎛ ⎞ ⎛ ⎞ ρ ρu Q = ⎝ ρu ⎠ and F = ⎝ ρuu + pI ⎠ , ρE (ρE + p)u

(1)

(2)

where t is the physical time, ρ the density, u = (u, v, w)T the flow velocity, p the pressure, E (= e + 12 |u|2 ) the specific total energy, e the specific inner energy, and I the identity matrix. The stiffened gas model [19] is chosen as the equation of state to close the Euler equations, p + γP∞ = (γ − 1)ρe,

(3)

where γ is the specific heat ratio and P∞ is a material-dependent constant. In order to track the interface that separates the fluids, a level set method [20] is adopted, ∂φ + vφ · ∇φ = 0, (4) ∂t where vφ is interface velocity and φ is the signed distance function (or level set function). In particular, the contour φ = 0 is used to define the location of the interface, and thus, φ < 0 can represent one bulk fluid and φ > 0 represents the other one. 3. Numerical methods Note that the Euler equations are discretized on a combination of the structured cells (i.e. belonging to the Cartesian mesh) and the unstructured ones (i.e. ICs arising from the cut-cell method; see details in Sec. 4) using a finite volume method, while the level set equation is only discretized on the Cartesian mesh (see Fig. 1). 4

φ>0 φ<0

φ=0 Figure 1: A schematic 2D example of the mesh generated by the cut cell method. The solid line segments denote the reconstructed interface on a Cartesian mesh according to the contour of level set function φ = 0. The square symbols indicate the centers of regular Cartesian cells, where φ is defined, including these in the vicinity of interfaces () and these away from the interface (). The circle symbols denote the centroids of the split cells occupied by the fluid 1 (◦) and the fluid 2 (•) respectively. All the flow variables are defined at the centroids of split cells in the vicinity of the interface and the centers of Cartesian cells away from the interface.

3.1. Finite volume formulation A finite volume method in the ALE framework is used to discretize the Euler equations. Considering a control volume Ω(t) with moving boundary ∂Ω(t), Eq.(1) can be rewritten in an integral form as   ∂ Q dV + F· n dS = 0, (5) ∂t Ω(t) ∂Ω(t) where n is the outward unit vector normal to ∂Ω. Therefore, an explicit time integration of Eq.(5) within one time step Δt yields [17],  Qn+1 Vin+1 = Qni Vin − Δt Fne Sen , (6) i e

where the subscript i (=1 or 2) denotes the respective fluid, e denotes the flow quantities at an arbitrary cell face e, the superscripts n and n + 1 denote the time levels, Se denotes the area of the corresponding cell face, and Fe represents the numerical flux estimated at the centroid of and normal to the cell face e. For the convenience of computation, Fe at either moving or stationary cell faces can be generally defined as ⎞ ⎛ ρ(Un − Ub ) ⎜ ρu(Un − Ub ) + pnx ⎟ ⎟ ⎜ ⎟ Fe = ⎜ (7) ⎜ ρv(Un − Ub ) + pny ⎟ , ⎝ ρw(Un − Ub ) + pnz ⎠ ρE(Un − Ub ) + pUn 5

where Un = u · ne , Ub = ub · ne , and ne = (nx , ny , ny ) is the outward unit vector normal to the cell face and ub is the velocity of the cell face. In the finite volume formulation, all the primitive and conservative flow variables are defined at the centroid of the cells (see Fig. 1). To calculate the numerical fluxes in Eq.(7), the primitive variables at the centroid of the cell faces are first calculated by those at the cell centroid using a second-order approximation scheme, PCF = PC + ψ∇PC · r,

(8)

where P stands for an arbitrary variable and the subscripts CF and C denote the centroid of the cell face and the cell, respectively; r is the displacement vector between the centroids of the cell and cell face. For the approximation of ∇PC , a MUSCL scheme [21] is used for the Cartesian cells away from the interface (i.e. regular cells), along with a slope limiter ψ calculated by the van Leer limiter [22]; by contrast, a second-order least-square method [23] is used for the ICs, along with the Barth limiter [24]. In order to compute the numerical flux Fe , one-dimensional Riemann problem is formulated at the centroid of the cell face along its normal direction, based on the states of flow variables on the two sides of the cell face; note that the velocity components in this Riemann problem are calculated by Un = u · ne . The numerical flux Fe is then obtained by solving the local one-dimensional Riemann problem [25]. In particular, AUSM+ -up scheme [26] is used for the single-phase Riemann problem, while an exact Riemann solver [16] is used for the two-phase Riemann problem. For the two-phase Riemann problem with a stiffened-gas model as equation of state for the liquid, the details on how to implement the exact Riemann solver can be found in the Appendix A of Chang et al. [16]. The time step Δt is chosen adaptively such that the Courant-Friedrichs-Lewy (CFL) condition [27] is strictly satisfied, √ 3 Vi Δt = CFL · min{ }, (9) i∈D |ui | + Ci where the CFL number (= 0.44) is a prescribed constant and is restricted by the presence of irregular ICs, of which the length in one dimension may be much shorter than the other two [6]; the subscript i denotes an arbitrary unstructured IC or Cartesian cell in the computational domain D, C is the local sound speed and V is the volume of the cell. An explicit second-order Runge-Kutta method [28] is used for time integration of the Euler equations. 3.2. Interface advection and reinitialization The level-set function φ is defined at the centers of the Cartesian cells only (Fig. 1), and its advance in time is performed by solving Eq.(4), with a fifth-order weighted non-oscillatory scheme [29] for the spatial discretization and the same scheme for temporal discretization as the Euler equations. In particular, in order to avoid the distortion of the φ field due to the non-uniform flow field, we use the extended interface velocity to advect the φ field [30], rather than the flow velocity. The interface velocity is obtained from the solution of the local Riemann problem at the interface, and is then extended to the Cartesian cells in 6

the vicinity of the interface along the normal direction of the interface, n. This extension process of interface velocity is accomplished by solving ∂vφ ± n · ∇vφ = 0, ∂τ

(10)

where τ is a pseudo time, and the normal vector n is defined as n = ∇φ/ |∇φ|; in particular, +n is used in the region φ > 0 and −n is used in the region φ < 0. In such a way, the φ field is maintained more or less a signed distance function after advection, in the absence of the topology change of interface. For long time simulations or in the presence of topology change or high interfacial curvature, we still need to reinitialize the φ field. In these cases, we solve the reinitialization equation [31], φτ + s(φ0 )(|∇φ| − 1) = 0, (11) where τ is a pseudo time, s(φ0 ) is a sign function, and φ0 corresponds to the φ field before reinitialization, i.e. τ = 0. A second-order ENO finite difference method is used for the spatial discretization of Eq.(11) with subcell fix near the interface [32, 33]. 4. The cut-cell method The main difficulties of the cut-cell method lie in the manipulation of cut cells, i.e. the Cartesian cells contain interfaces. The manipulation includes identification and splitting of the cut cells, assembly of ICs in the vicinity of interfaces, and calculation of the geometry properties of the ICs, to name a few. In the following, details of implementing 3D cut-cell method will be addressed. 4.1. Identification and splitting of 3D cut cells Assume that there is at most one cut point between two adjacent corners of a Cartesian cell. Cut cells can be identified by the local distribution of φ, specifically, at the corners of Cartesian cells. If there is any sign change of φ between two adjacent corners of a Cartesian cell, it implies that an interface intersects with the Cartesian cell at the cell edge between the two corners. As a result, these Cartesian cells are marked as cut cells. Since φ is defined at the centroid of Cartesian cells, one needs to accurately calculate φ at the corners for this process of cell identification. This can be done by averaging the φ at the neighboring Cartesian cells. For example, for an arbitrary corner at (i − 12 , j − 12 , k − 12 ), we can have φi− 1 ,j− 1 ,k− 1 = 2

2

2

1

φi−1,j−1,k−1 + φi,j−1,k−1 + φi−1,j,k−1 + φi,j,k−1 8  + φi−1,j−1,k + φi,j−1,k + φi−1,j,k + φi,j,k .

(12)

In order to keep the interface coinciding with the cell faces, the interface needs to be reconstructed on the Cartesian mesh at each time step during the computation. A piecewise linear interface reconstruction is adopted by connecting the cut points, i.e. the place where an interface intersects with the cell edge. The cut points can be easily determined through 7

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: Basic configurations of single-cut cells, in which the configurations (a)-(f) represent 16, 24, 6, 48, 24, 8 types of single-cut cells, respectively. The black cubes represent the cell corners occupied by the same fluid in the cut cell, and the sphere symbol represents the cut point at which an interface intersects with a cell edge. The connection of the cut points forms the edges of a cut face, i.e. the reconstructed interface in the cut cell. .

the linear interpolation of φ. However, depending on the specific way of interface intersecting, various types of cut-cells can be encountered in the computation, which is crucial for the interface reconstruction. For a 3D cut cell filled with two fluids, the level set function allows to resolve 284 configurations, which is too complicated to consider in practice. We find that the 284 cutcell configurations can be greatly reduced to 12 basic ones, by switching the position of the two fluids and/or rotating the orientation of the cut cell. The 12 basic cut-cell configurations are shown in Figs. 2 and 3, which are further grouped into single-cut and multiple-cut cells. We note that the cut cell is essentially split by the interfaces reconstructed by the cut points, namely cut faces (see Fig. 2). In the single-cut configuration, the number of the cut points ranges from 3 to 6. However, there is only one cut face, which divides the cut cell into two polyhedrons occupied by the two fluids, respectively. The polyhedrons generated by splitting of cut cells are referred to as cut polyhedrons hereafter. In the computation, most cut cells belong to the single-cut configuration, which is therefore needed to be treated regularly. The multiple-cut configuration represents the cut cells where the topological change of interfaces is about to occur or just takes place. The number of the cut faces in the multiple-cut configuration ranges from 2 to 4, resulting in 3 to 5 cut polyhedrons respectively. Interface reconstruction and classification of cut-cell configurations are performed after the level set equation is solved. These processes make it easier to calculate the geometry properties of the cut polyhedrons (e.g. normal vector and area of the cut face; see details in Sec. 4.2) and identify the interface cells with under-resolved interfaces (e.g. the cut 8

(g)

(h)

(i)

(j)

(k)

(l)

Figure 3: Basic configurations of multiple-cut cells, in which the configurations (g)-(l) represent 32, 48, 12, 16, 2, 48 types of multiple-cut cells, respectively. The black cubes represent the cell corners occupied by the same fluid in the cut cell, and the sphere symbol represents the cut point at which the interface intersects with the cell edge.

polyhedron sandwiched by the interface; see details in Sec. 4.5). 4.2. Calculation of geometry properties Finite volume formulation of the Euler equations (i.e. Eq.(6)) requires the geometry properties of the computational cells. However, it is not trivial to calculate the geometry properties, particularly for the cells with the reconstructed interface as its cell face; these cells could be arbitrary polyhedral cells. For the cut cells, these geometry properties include the area, center and normal vector of the cut face and the volume and centroid of the cut polyhedron. It is noteworthy that the cut face consisting of M cut points is probably a non-planar polygon when M > 3. In order to calculate the normal vector and area of the cut face in such cases, we first decompose it into M − 2 triangles. However, there exist 2 types of triangle decomposition for 4 cut points, 5 types for 5 cut points, and 14 types for 6 cut points, characterized by the selection of the target vertex (see Fig. 4). Among all these possibilities, we consider the best decomposition is the one closest to a planar polygon. Suppose that after an arbitrary decomposition we have M − 2 triangles, of which the area and normal vector are denoted by Si and ni , i = 1, ..., M − 2. The average normal vector ˆ = Σi (Si ni )/Σi Si , and the difference in the normal vectors can be defined is calculated by n ˆ |Si )/(Σi Si ). Therefore, the optimal triangle decomposition is the one that as  = (Σi |ni − n has the smallest . With the optimal triangle decomposition, we can calculate the geometry properties of the cut face; see e.g. Fig. 5. For a triangle consisting of three cut points (with the position 9

(a)

(b)

(c)

(d)

Figure 4: Triangle decomposition for cut faces consisting of M > 3 cut points. (a), (b) and (c) show the typical triangle decomposition for cut faces consisting of 4, 5 and 6 cut points, respectively. (d) shows a complete set of triangle decomposition by shifting the target vertex position in (b).

(a)

(b)

CP2

b CP3

xc β

n

a

CP1 CP4 Figure 5: Illustration of interface reconstruction in a cut cell by triangulation. (a) The split cell consisting of four cut points (CP1 , CP2 , CP3 and CP4 ) and two cell corners (CC1 and CC2 ) have six cell faces. (b) The reconstructed interface is triangulated into two triangular cell faces, each of which has its own centroid, area and normal vector.

10

vector xCP1 , xCP2 and xCP3 respectively), its face center can be simply calculated by xc = (xCP1 + xCP2 + xCP3 )/3.

(13)

The area of an arbitrary triangle can be calculated by Si = |a × b|/2, where a = xCP1 − xCP2 and b = xCP2 −xCP3 . The normal vector of the triangle can be calculated by n = (a×b)/|a× b|. Similarly, the geometry properties of other cell faces that enclose the cut polyhedron can also be calculated. Moreover, any triangle decomposition would do since all those cell faces are planar polygons. Note that the interface velocity vφ is obtained by solving the local Riemann problem [34] at the triangles in the decomposition. To extend it to the surrounding Cartesian cells by Eq.(10), the interface velocity for the whole cut face is needed. A area-weighted interface velocity is used for this purpose, vφ = ( i vφ,i Si )/(Σi Si ), (14) where vφ,i corresponds to the Riemann solution of interface velocity at the center of the triangle i. After obtaining the geometry properties of the cell faces that enclose the cut polyhedron, we can calculate the volume (Vcv ) and centroid (rcv ) of the cut polyhedron [35], 1 Vcv = Σj xc,j · nj Sj , 3  1 rcv = Σj (xc,j · nj ) Sj xc,j , 4Vcv

(15)

where xc,j , Sj and nj represent the center, area and normal vector of the jth triangle (or the planar cell face), respectively. 4.3. Formation of interface cells by cell assembly In principle, cut polyhedrons can be directly used as ICs, at which the discretized Euler equations, i.e. Eq.(6), are solved. However, using tiny polyhedrons (e.g. the volume less than 0.1Δx3 ) as the ICs could be troublesome. At the n time level, it would lead to very small time step due to the CFL condition (i.e. Eq.(9)), thereby making the computation very expensive. At the n + 1 time level, tiny ICs may cause unphysical prediction of Qn+1 , i which is computed by dividing by V n+1 (see Eq.(6)). In order to avoid too small ICs, the small polyhedron is merged with a neighboring polyhedron or Cartesian cell, so that an IC of relatively large volume is formed (see also 2D cell merging method [23]). The assembly of the ICs for the fluids 1 and 2 are separately performed, although they may share the same reconstructed interface as their cell faces. It consists of two steps: an initial assembly at the n time level and a correction step at the n + 1 time level. As a result, we can obtain the IC at the n time levle, Ωn , and its counterpart at the n + 1 time level, Ωn+1 . For the convenience of illustration, in the following we only demonstrate the assembly process of 2D ICs. The extension to the assembly of the 3D ICs is straightforward. 11

(a)

(b)

1

1

2

1

2 2 4 4 3

2 2 3 4

φn=0

φn+1=0

Figure 6: Formation of ICs by cell assembly: 2D illustration. (a) The cell assembly for ICs Ωn is firstly predicted according to the interface position φn . (b) The cell assembly for the IC Ωn+1 is obtained based on the initial prediction in (a) and the interface position at the n + 1 time level. Note that the cell assembly for Ωn may need to be corrected, e.g. the 4th IC in (a) will include its right neighboring cell (in grey) in its cell assembly, in order to be consistent with the cell assembly of Ωn+1 .

Fig. 6 shows the examples of assembly processes of the ICs, particularly for the cut cells in the single-cut configuration. At the n time level, we first calculate the volume fractions of the two fluids in the cut cell, according to the splitting by the reconstructed interface. Assume that the fluid 1 occupies a volume fraction of α, and then the fluid 2 would occupy 1 − α. If α (or 1 − α) is smaller than a threshold value (which is a user-defined parameter, and we suggest the threshold value of α ∼ 0.35 − 0.5), the cut polyhedron occupied by the corresponding fluid is considered to be so small that it should be merged with its immediate neighboring cell (namely target cell). The target cell could be either a cut polyhedron or a Cartesian cell. The target cell for merging is chosen according to the normal vector of the interface. For example, for the cut polyhedron in a cut cell (i, j, k) with the normal vector (nx , ny , nz ), its target cell for merging can be determined by the dominant component of the normal vector, ⎧ if |nx | > |ny | and |nx | > |nz |, ⎨ (i ± 1, j, k), (i, j ± 1, k), if |ny | > |nx | and |ny | > |nz |, (16) ⎩ (i, j, k ± 1), if |nz | > |nx | and |nz | > |ny |, where the choice of + or - depends on the condition that the target cell should share the same cell face with the cut polyhedron. The merging of the cut cell with its target cell forms the initial configuration of the IC. After the first-round cell merging, the volumes of all the ICs will be checked again, to see whether they are larger than the volume threshold. If not, the small IC will be further combined with its immediate neighbor (preferably with the least volume), until the volume constraint is satisfied. Exceptions are the under-resolved 12

(a)

(b) n

φ =0 φn+1 = 0

φ

φn+1 = 0 =0

n+2

Figure 7: Illustration of the reassembly processes of ICs at the beginning of each time step in the 3D cases. (a) Cell assembly of ICs at the n time levels. The sphere (cube) symbols denote the cut points at the n (n + 1) time level, and the reconstructed interface is represented by the line segments connecting these cut points. (b) Cell assembly at the beginning of the n + 1 time level. Note that the sub-cells of the IC at the n + 1 time level are different in (a) and (b) because of the reassembly of ICs.

interfaces, such as tiny droplets (or bubbles) and thin filaments, for which special treatments are performed (see Sec. 4.5). The initial assembly process of Ωn is illustrated in Fig. 6(a). At the n + 1 time level, we seek the configuration of Ωn+1 based on the initial assembly of Ωn , which could be corrected by the configuration of Ωn+1 accordingly. When the local interface moves rather slowly, it ˆ n+1 , i.e. the initial guess of may be still contained in the same cut cell; at the same time, if Ω Ωn+1 based on the initial assembly of Ωn , also has a volume larger than the volume threshold, ˆ n+1 can be directly regarded as Ωn+1 . In other words, Ωn+1 has the same configuration as Ω ˆ n+1 is below Ωn in this case; see e.g. the 3rd IC in Fig. 6. In the cases that the volume of Ω n+1 ˆ needs to merge with a neighboring Cartesian cell to form Ωn+1 ; see e.g. the threshold, Ω th the 4 IC in Fig. 6(b). Consequently, the configuration of Ωn is required to be changed accordingly, in order to be consistent with Ωn+1 ; see e.g. the 4th IC in Fig. 6(a). When the interface crosses over the cell corners from n to n + 1 time level, it may lead to a new cut cell at the time level n + 1; see e.g. the 2nd IC in Fig. 6(b). As a result, the new cut cell is ˆ n+1 to form Ωn+1 . Similarly, the passing of the interface added into the configuration of Ω over the cell corners may make a cut cell in the configuration of Ωn become a Cartesian cell ˆ n+1 to generate at the n + 1 time level. Then, this cell is deleted from the cell assembly in Ω n+1 st Ω ; see e.g. the 1 IC in Fig. 6(b). Note that at the beginning of each time step, all the ICs Ωn need to be reassembled, rather than inheriting the configurations from the early time level, e.g. Ωn+1 at the time level one time step earlier; one 3D example is illustrated in Fig. 7. This process ensures that ICs with better resolution are continuously produced to resolve the interface. The dynamic cell assembly gives rise to interface cells with size roughly ranging from 0.65Δx to 1.6Δx. Compared to the strategy of cell combination in the Lin’s method [17], it not only avoids the occurrence of very large interface cells (e.g. ∼ 3Δx), but also restricts the size distribution to a relatively narrow range. Therefore, it significantly improves the interface resolution, and 13

also reduces inconsistency in numerical errors due to the size difference between interface cells. 4.4. Computation and redistribution of conservative variables The computation of the Euler equations is performed at the ICs in the vicinity of the interfaces or at the regular Cartesian cells away from the interfaces (see Fig. 6), by updating the conservative variables. As shown in Sec. 4.3, the ICs are composed of various subcells, such as the cut polyhedrons and Cartesian cells. Due to the fact that the ICs are assembled and disassembled repeatedly during the computation, it is more convenient to store all the conservative variables at the centroids of the subcells than at the centroids of ICs. Therefore, we need to compute Qni at the IC Ωni at the beginning of each time step, and redistribute to the subcells of Ωn+1 at the end of each time step, preferably in a conservative Qn+1 i i manner. Calculation of Qn V n at the IC Ωni is rather straightforward, e.g. by a second-order integration over the subcells of Ωni  Qni Vin = Qni,p Vi,p , (17) p∈Ωn i

where p denotes an arbitrary subcell of Ωni , Qi,p represents the conservative variables defined at the centroid of the subcell p, and Vi,p is the volume of the subcell p. After the conservative variables at the n + 1 time level (i.e. Qn+1 ) are obtained for Ωn+1 , i i n+1 they should be redistributed to the subcells of Ωi , particularly for the ICs. Assuming that q is one of the subcells of Ωn+1 , the conservative variables at its centroid Qn+1 i i,q can be approximated by n+1 n+1 Qn+1 + ∇Qn+1 · (rn+1 ), (18) i,q = Qi i i,q − ri and rn+1 denote the centroid position of Ωn+1 and its subcell q, respectively. where rn+1 i,q i i n+1 n+1 n+1 denotes the gradient of Qi at the centroid of Ωi ; it can be interpolated using ∇Qi the conservative variables at the centroids of Ωn+1 and its neighbors occupied by the same i fluid, similarly to the approximation of ∇P in Eq.(8) (see Sec. 3.1). Note that the stencil for computing ∇Qn+1 is made up of Ωn+1 and its neighbors, i.e. other ICs and Cartesian cells, i i and the gradient approximation is first-order accurate. This ensures that the calculation of n+1 n+1 n+1 Qi,q by Eq. (18) is second-order accurate. Given the fact that Vi = q∈Ωn+1 Vi,q and i n+1 n+1 n+1 n+1 = q∈Ωn+1 (ri,q Vi,q ), it is easy to obtain ri Vi i  n+1 Vin+1 = Qn+1 (19) Qn+1 i i,q Vi,q , q∈Ωn+1 i n+1 where Vin+1 is the volume of Ωn+1 and Vi,q is the volume of the subcell q. In other words, i the redistribution scheme in Eq.(18) is conservative. It is noteworthy that the second-order accurate spatial integration of Euler equations at the ICs is a result from the second-order accuracy of the whole process, i.e., computation and integration of numerical flux Fe , integration of Qn over Ωni (i.e. Qni Vin ) and redistribution of Qn+1 back to the subcells of Ωn+1 . i 14

(a)

(b)

(c)

Figure 8: Interface cells containing under-resolved interfaces. (a) Thin filament, (b) isolated bubble at the cell corner and (c) isolated bubble at the cell center. The solid thick lines represent the reconstructed interface by φn = 0 while the dashed lines represent the one by φn+1 . In (c), the filled circle denotes φ > 0 and the empty squares denote φ < 0.

4.5. Interface cells with under-resolved interfaces Under-resolved interfaces often occur when the interface coalescence/breakup is about to take place. It leads to the cut cells in the multiple-cut configurations (Fig. 3) or thin filaments in the single-cut configuration (Fig. 8(a)). Other under-resolved interfaces include tiny droplets or bubbles arising from the interface breakup (Figs. 8(b) & 8(c)). Specifically, if the diameter of either the filament or the droplet (or bubble) is less than Δx, the local interface related is considered to be the under-resolved one. Accordingly, the relevant computational cells are marked as ICs with under-resolved interfaces (ICUIs). Depending on the specific situations, the ICUIs are formed in different manners. In the multiple-cut configuration, the cut polyhedron sandwiched by the interfaces is related to the ICUI. When the volume of the sandwiched cut polyhedron is larger than the threshold value, the cell is directly marked as the ICUI. Otherwise, it is merged with the largest cell among its neighboring cells that are occupied by the same fluid, in order to form an ICUI. By contrast, the rest of the cut polyhedrons in the multiple-cut configuration is merged with their neighbors to form the ICs, similarly to the single-cut configuration. The accuracy of computation is degraded from the second order to the first order in the ICUIs. Specifically, for the flux approximation, the flow variables at the cell face of the under-resolved cell is directly copied from its centroid, which is equivalent to set ∇P = 0 in Eq.(8). Similarly, ∇Qn+1 = 0 is used for the redistribution of the conservative variables i in Eq.(18). Thus, the variable redistribution is conservative but only first-order accurate. For very tiny droplets (or bubbles) trapped at the vertex (see Fig. 8(b)) or cell center (see Fig. 8(c)), special treatment is usually needed for the update of flow variables to improve the robustness of computation, e.g. through artificial elimination [13]. In the present study, the flow variables at the centroid of these ICUIs at the n + 1 time level are simply kept as the same as those at the n time level, thereby preventing the unphysical violent variation of the flow variables, which could arise from the volume change due to the reinitialization of level set function. Clearly, this treatment would lead to insignificant mass loss. 15

4.6. Solution procedure We summarize the solution procedure of the 3D CSI method with two-stage Runge-Kutta method for time integration as follows: Step 1 Advance the level set equation φn (and conduct the reinitialization if necessary) to obtain φ(1) ; Step 2 Reconstruct the interface according to φ(1) , identify the type of cut cells and calculate the geometry information of the cut polyhedrons; form the ICs by cell assembly for Ωn and Ω(1) respectively, and calculate the conservative variables Qn by integration over the subcells of Ωn by Eq.(17). Step 3 Reconstruct the primitive variables at the centroid of cell faces by Eq.(8), and approximate the fluxes Fne at non-interface cell faces using AUSM+ -up. Compute the fluxes Fne at the cell faces coinciding with the reconstructed interface using Eq.(7) with the knowledge of Un = Ub . Step 4 Integrate Eq.(6) to obtain the conservative variables Q(1) , and redistribute Q(1) to the subcells of the ICs through Eq.(18). Calculate the primitive variables according to Q(1) at the subcells and Cartesian cells. Step 5 Reconstruct the primitive variables at the centroid of the reconstructed interface of Ω(1) by Eq.(8), and perform the exact Riemann solver as discussed in Sec. 3.1. (1) Extend the interface velocity vφ (which is obtained from the solution of the exact Riemann solver) from the reconstructed interfaces to the Cartesian cells by Eq.(10). Step 6 Advance the level set equation φ(1) (and conduct the reinitialization if necessary) to obtain φ(2) , and calculate φn+1 = 12 (φn + φ(2) ). Step 7 Repeat Step 2 using φn+1 , so as to reconstruct the interface, calculate the geometry interface of the cut polyhedrons, form the ICs for Ω(1) and Ωn+1 respectively and calculate the conservative variables Q(1) . Note that the subcells of Ω(1) here may be different from that of Ω(1) formed in Step 2. (1)

Step 8 Approximate the fluxes Fe at the non-interface cell faces and compute the fluxes (1) Fe at the reconstructed interface, in the same way as in Step 3. Step 9 Integrate Eq.(6) to obtain the conservative variables Q(2) , and calculate Qn+1 V n+1 = 12 (Qn V n + Q(2) V (2) ). Redistribute the conservative variables Qn+1 to the subcells of Ωn+1 through Eq.(18), and compute the primitive variables accordingly at the subcells and Cartesian cells. Step 10 Reconstruct the primitive variables at the centroid of the reconstructed interface of Ωn+1 , perform the exact Riemann solver and extend the interface velocity vφn+1 to the Cartesian cells, in the same way as in Step 5. Step 11 Go to Step 1. 16

(a)

(b)

Figure 9: Numerical results of Rayleigh collapse in terms of time variation of (a) bubble volume and (b) 3 /6; in bubble mass. In (a), V ∗ represents the bubble volume normalized by the initial volume Vini = πRini (b) m denotes the instantaneous bubble mass during the computation, and m0 denotes the initial value.

5. Results and discussion In this section, the performance of the 3D CSI method is systematically examined by various numerical experiments, including Rayleigh collapse, shock-induced bubble collapse in water, interaction between shock and curved gas-gas interface, and droplet deformation after shock impact. In particular, numerical results of the last two cases are compared against experimental data, in terms of schlieren at different times. Due to computational efficiency, fine uniform Cartesian mesh is used in the part of the domain where interfaces and their evolution take place, while stretched structured meshes with relatively coarse resolution are used in the rest of the domain. Unless stated otherwise, the resolution is referred to as the fine uniform Cartesian mesh hereafter, and is measured in terms of the number of mesh points per radius (ppr), where the radius could be that of the relevant bubble or droplet at the very beginning of computation. 5.1. Rayleigh collapse Rayleigh collapse refers to the collapse process of an isolated spherical gas bubble in high pressure liquid. The specific case considered here has following initial conditions: the radius of the air bubble Rini is 10 μm, the air density is 1.226 kg/m3 , the pressure inside the bubble is 101.325 kPa, the density of water is 1000 kg/m3 , the pressure in water is 10.1325 MPa, and all the fluids are at rest initially. The flow parameters are: γair =1.4, P∞,air =0 Pa, γwater =4.4 and P∞,water =6 ×108 Pa. The computation is performed in a cubic computational domain of [0, 50] × [0, 50] × [0, 50] μm, which corresponds to one-eighth of the physical domain, because of symmetry of the collapse process. Specifically, a uniform fine mesh is used within [0, 15] × [0, 15] × [0, 15] μm, with the center of the spherical bubble being initially set at (0, 0, 0). Symmetric boundary conditions are used at x = 0, y = 0 and 17

Resolution (ppr) 20 40 80

EV,L2 6.665 × 10−2 2.294 × 10−2 5.576 × 10−3

Order 1.54 2.04

EV,L∞ 1.493 × 10−3 5.478 × 10−3 1.474 × 10−3

Order 1.45 1.89

Table 1: Convergence study of the Rayleigh collapse of a bubble in water, with respect to the relative error of the bubble volume at t = 0.085 μs.

z = 0, while an extrapolation boundary condition is used for other boundaries. The flows are resolved on four meshes with successive mesh refinement, i.e. 20, 40, 80 and 160 ppr, respectively. In particular, the results with the mesh resolution of 160 ppr are chosen as the solution for comparison. The result comparison is shown in Fig. 9(a), with respect to the time evolution of the bubble volume. At the early stage of Rayleigh collapse, the significant pressure difference between the gas bubble and its surrounding liquid induces the inward motion of the interface. As a result, the bubble shrinks with time, accompanied by rapid increase of pressure inside the bubble. When the pressure inside the bubble is sufficiently high to prevail over the inward motion of the liquid, the bubble would start to expand (at t ≈ 0.09μs). After then, the bubble volume grows with time until the pressure in the bubble becomes too low to sustain the expansion. It is clear that the results of this process converge with mesh refinement. In order to assess the convergence rate, we calculate the volume difference between the 3D results roughly at the time with smallest bubble volume, i.e. t = 0.085 μs. We define the L2 -norm and L∞ -norm errors of the instantaneous bubble volume as follows,  2 i,j,k |Vlow,i,j,k − Ve,i,j,k | max|Vlow,i,j,k − Ve,i,j,k | (20) EV,L2 = , EV,L∞ = , V V e,i,j,k e,i,j,k i,j,k i,j,k where Ve and Vlow represent the volumes obtained with a mesh resolution of 160 ppr and lower, respectively. Numerical results of EV,L2 and EV,L∞ at t = 0.085 μs are shown in Table 1, and we can see that second-order accuracy has been achieved. Fig. 9(b) shows the time evolution of mass loss of the bubble on the different meshes. The mass loss is primarily due to the accumulation of machine errors, and is generally negligible in all cases (less than 10−10 ). It is also noteworthy that the mass loss does not increase with time on the fine meshes, i.e. 80 ppr and 160 ppr. These demonstrate the conservative properties of the method. 5.2. Shock-induced bubble collapse in water Shock-induced bubble collapse in water can be seen in applications such as shockwave lithotripsy [1] and ultrasound guided injection [2]. It has been investigated by Hawker and Ventikos[36] by 3D simulations with a front-tracking scheme. To validate our CSI method, their results are selected as the benchmark solutions, along with the results from axisymmetric simulations. Due to the symmetry of the problem, a quarter of the physical 18

Collision time (μs) Jet speed at collision (m/s) Water-hammer shock pressure (GPa)

Present Hawker & Ventikos [36] 0.541 0.540 2747 2754 5.20 4.04

Table 2: Quantitative comparison of numerical results with Hawker & Ventikos [36] for the bubble collapse induced by shock.

domain is simulated. The computation is performed in a cubic domain of [−1, 1] × [0, 1] × [0, 1] mm, with an air bubble of 0.5 mm in radius resting at (0,0,0). Initially, the bubble has a density of 1.184 kg/m3 and pressure of 101.325 kPa; the density of the water before the shock is 995.972 kg/m3 and the water pressure is 101.325 kPa. The pressure after the shock is 1 GPa. The shock is initially at plane x = 0.55mm, and then propagates to the x+ direction. A uniform fine mesh is used in [−0.7, 0.6] × [0, 0.6] × [0, 0.6] mm, which corresponds to a resolution of 200 ppr. Symmetric boundary condition is used at y = 0 and z = 0 sides while an extrapolation boundary condition is used for the other boundaries. Fig. 10 shows the snapshots of the bubble during its collapse process, along with the corresponding flow features in terms of pressure contours and numerical schlieren. Once the shock hits the air bubble, a re-entrant liquid jet occurs at the upstream of the bubble and moves in the same direction as the incident shock due to the baroclinic effect [36]. When the jet collides with the downstream interface of the bubble, a water-hammer shock emits from this collision process (see e.g. t = 0.545μs in Fig. 10(b)). The collision process is characterized by the time span from the shock-bubble contact to the occurrence of collision (namely collision time), the jet speed and the pressure of the water hammer shock at the occurrence of collision. After then, the jet penetrates through the bubble and merges with the surrounding liquid. Fig. 11 shows the comparison of the bubble shapes (at 0.415 μs and 0.585 μs) with the axisymetric results on the mesh with the same resolution (200 ppr). The two results are virtually overlapped. It is clear that the 3D CSI method can provide high-fidelity simulations, with the symmetricity of the solution being well maintained. A quantitative comparison with the benchmark solution [36] is provided in Table 2, with respect to the collision time, the jet speed at the collision and the pressure of the water hammer shock. A good agreement has been obtained regarding the collision time and the corresponding jet speed. As the strength of the water hammer shock is an instantaneous value and more sensitive to the mesh resolution, it has relatively bigger difference compared to the prediction of Hawker & Ventikos [36]. 5.3. Interaction between shock and curved gas-gas interface We consider here the interaction of a planar incident shock with curved gas-gas interfaces, consisting the same 3D convex and concave interfaces as in the experiments of Ding et al. [37] (see Fig. 12). The interfaces initially rest between two parallel plates separated by a distance of h = 20 mm, and the radius of the interface in contact with the plate r1 is fixed to 17.5 mm. The 3D convex interface (see Fig. 12(a)) can be described in a cylindrical coordinate 19

(a)

(b)

(c)

(d)

Figure 10: Numerical results of bubble collapse induced by a planar shock with respect to 3D bubble shape and contours of pressure and numerical schlieren in the plane of z=0, at times: (a)0.415μs, (b)0.545μs, (c)0.585μs and (d)0.630μs, respectively.

20

(a)

(b)

Figure 11: Comparison of interface shapes in the shock-induced bubble collapse between the 3D simulation (specifically, in the plane of z=y; the axis orientation can be seen in Fig. 10(a)) and an axisymmetric one. (a) t = 0.415μs and (b) t = 0.585μs.

(a)

(b)

Figure 12: Initial shape of the 3D curved interfaces: (a) Convex interface and (b) concave interface.

21

Figure 13: Dynamics of interaction between the shock and 3D convex interface, in terms of experimental schlieren [37], numerical schlieren and snapshot of 3D interface in simulations, respectively. The number at the right-top corner of each panel in the experimental schlieren denotes the time in μs.

22

system (r, z) as

−rzz 1 Δp + = , (21) 2 3/2 2 1/2 (1 + rz ) r(1 + rz ) σ where σ is the surface tension coefficient of the gas-gas interface, Δp is the pressure difference across the interface, and rz = dr/dz and rzz = d2 r/dz 2 . In the experiments of Ding et al. [37], Δp is set to 4.3 Pa and σ = 3.5×10−2 N m−1 . The 3D concave interface is obtained by setting Δp = 0 in the experiment, thereby yielding a ‘minimal surface’, of which the shape can be obtained analytically r z = cosh( ), (22) r0 r0 where r0 = 13.7 mm is the radius of the concave surface at the symmetry plane; see Fig. 12(b). Initially, the pressure before the shock is the ambient pressure, i.e. 101.325 kPa, and the ambient gas is SF6 (with density 6.05 kg/m3 and γ = 1.094). In the case of interaction between a planar shock and the 3D convex interface, the gas enclosed by the interface is N2 (with density 1.16 kg/m3 ). The initial shock Mach number is 1.28, and γ is set to 1.4 for N2 . In the case of interaction between a planar shock and the 3D concave interface, the gas enclosed by the interface is a mixed gas, i.e. 95% N2 + 5% SF6 (of density 1.41 kg/m3 ). The Mach number of the incident shock is 1.29, and γ is set to 1.34 for the mixed gas. The shock is initially at plane z = 10mm, and then propagates to the z+ direction. Due to the symmetry of the flow problem, the computation is only performed in one-fourth of physical domain, i.e. [0, 70]× [0, 10]× [0, 200]mm, and a uniform mesh of 350 × 50 × 1000 is used. The center of the 3D curved interfaces is located at (0, 0, 40) mm initially. The symmetric boundary conditions are enforced at the boundaries of x = 0, x = 70 mm, y = 0 and y = 10 mm, respectively. An extrapolation boundary condition is implemented at z = 0 and z = 200 mm. Since the schlieren is only two-dimensional in experiments, we need to appropriately give an integrated view of numerical schlieren along the y direction similar to the experimental counterparts [37] for visualization and comparison. The numerical schlieren in the (x, z) plane is integrated as follows: We firstly calculate the equivalent density ρei,k , jR j=jL (ρ1,i,j,k V1,i,j,k + ρ2,i,j,k V2,i,j,k ) . (23) ρei,k = j R j=jL (V1,i,j,k + V2,i,j,k ) where jL and jR denote the leftmost and rightmost index in the y direction respectively, and V1 and V2 represents the volume fraction of the fluids 1 and 2 respectively. Then we numerically compute ∇ρe by a second-order central difference scheme. The numerical schlieren function ϕ is defined as ϕ = exp(−

|∇ρe | ). |∇ρe |max

(24)

Fig. 13 shows the evolution of the 3D convex interface after being impacted by the planar shock, with respect to the experimental schlieren, numerical schlieren and 3D snapshots of 23

(a)

(b) 60

(mm)

50 40 30 20 10

0

200

400

600

800

1000

t (μs) Figure 14: (a) Quantitative comparison between the experimental data (symbols) [37] and numerical results (curves), with respect to the bubble length in the streamwise direction in the symmetry plane () and plate plane (), and the bubble width in the symmetry plane ( ) and plate plane (), respectively. (b) Variation of the bubble mass m with time. m0 is the initial mass of the bubble.

the interface, respectively. We can see that after the incident shock hits the convex interface (by firstly contacting the interface in the symmetry plane), a transmitted shock and reflected rarefaction waves occur simultaneously; the former moves towards the downstream while the latter moves towards the upstream (see e.g. t = 45 μm in the first row and second column). A second transmitted shock forms at the downstream of the interface roughly at t = 165 μs, with notable 3D shock structures. Moreover, the interface on the plate moves faster than the part in the symmetry plane; consequently, the interface collision takes place on the plate plane roughly at t = 665 μs, while it occurs roughly at t = 845 μs. At the late stage (e.g. t > 1000 μs), the interface keeps evolving and exhibits complex topology changes, although the incident shock is away from the interface. In general, the numerical schlieren compare well with the experimental observations during the whole process, and the 3D interface snapshots provide a new view of complicated interface evolution in the shockinterface interaction. We also perform a quantitative comparison with the experimental data in Fig. 14(a), with respect to the length and width of the convex interface in symmetry and plate planes. A good agreement is also obtained, except that there is a slight difference in the length and width of the interface in the plate plane at the late stage of interaction. However, we note that the experimental image of interface on the plate is blurred at this time due to the complicated interface structures. Fig. 14(b) shows the variation of mass during the computation. Mass loss of the bubble appears to coincide with the occurrence of topology change (or the ICUIs), e.g. interface breakup on the plate at t ≈ 600μs and in the middle of the curved surface at t ≈ 900μs (see Fig. 13). Overall, the total mass reduces roughly 1%, suggesting that mass conservation is well maintained in the CSI method. Fig. 15 shows the evolution of the 3D concave interface after being impacted by the planar shock, with respect to the experimental schlieren, numerical schlieren and 3D snapshots of 24

Figure 15: Dynamics of interaction between the shock and 3D concave interface, in terms of experimental schlieren [37], numerical schlieren and snapshot of 3D interface in simulations, respectively. The number at the right-top corner of each panel in the experimental schlieren denotes the time in μs.

25

(a)

(b)

Figure 16: Comparison of droplet shapes with the experiments [38] at t = 7 μs (a) and 10μs (b). The numerical results are obtained with a resolution of 100 ppr.

the interface, respectively. Compared to the convex interface, it is the part of interface close to the plate that protrudes out of the concave interface. Especially, the interface is significantly stretched at the late stage of shock-interface interaction. From the comparison of schlieren between the numerical simulations and experiments, it is clear that both the flow features such as the shocks and rarefaction waves and the interface structures have been accurately resolved in the computation. Also, the CSI method shows to be able to capture the dynamic 3D interface structures in the presence of topology changes, which may not be accessible in experiments. 5.4. Droplet deformation after shock impact We consider here the simulation of droplet deformation in air after being exposed to a shock at a Mach number of 2.67. The same flow problem has been experimentally studied by Theofanous et al. [38]. They observed that the shock-induced interface waves initially occurred at the upstream of the droplet; after then the interface waves propagated towards the downstream of the droplet, accompanied by the deformation at the rear of the droplet. However, it is difficult in experiments to quantitatively assess the growth of the interface waves and droplet deformation in the shock-droplet interaction. Therefore, it is necessary to resort to the numerical simulation for this purpose. We perform the 3D simulations under the same flow conditions as in [38]. Because of the high flow speed, the surface tension and gravitational force are considered to be negligible. The initial radius of a glycerol droplet is 0.895 mm, and the density of the glycerol droplet is 1263 kg/m3 with ambient pressure 101.325 kPa. The air density before the shock is 1.25 kg/m3 , and the shock Mach number 26

(a)

(b)

Figure 17: Flow features at times t = 7 μs (a) and 10 μs (b), in terms of the rescaled pressure p∗ = p/p0 with p0 = 2.5 × 106 P a in the cross section of z = 0 and the rescaled vorticity (in y-direction) field ω ∗ = ω/ω0 with ω0 = 107 m2 s−1 in the cross section of y = 0. Solid curves represent the interface in the cross section of z = 0.

is 2.67. The planar shock is initially located at x = −1 mm, propagating towards the x+ direction. The stiffened gas model is used, and γG = 1.5048 and P∞,G = 2.99 × 108 Pa for glycerol. Due to the symmetry of the flows, only one-fourth of the physical domain is considered. Therefore, the computation is performed in a domain of [−3, 3] × [0, 3] × [0, 3] mm, in which a uniform fine mesh of 100 ppr resolution is used in the sub-domain of [−1.3425, 1.3425] × [0, 1.3425] × [0, 1.3425] mm. Symmetric boundary conditions are enforced at y = 0 and z = 0, while the extrapolation boundary condition is implemented at the other boundaries. The droplet initially rests at (0, 0, 0), and the time t = 0 corresponds to the moment when the shock comes into contact with the droplet. Fig. 16 shows the snapshots of the droplet at t = 7 μs and 10 μs. It is clear that fish-scalelike 3D interface waves are developing at the droplet surface. In particular, the interface waves are notable around θ ≈ 110 − 135◦ at t = 7 μs (see the definition of θ in Fig. 16(a)), and gradually propagate towards the downstream with the increase of time. The wave front reaches the peak of the droplet (θ = 90◦ ) roughly at 10 μs, accompanied by the occurrence of a groove-like interface structure at the upstream of the droplet around θ = 135◦ and the surface deformation at the downstream of the droplet. It is worthy mentioning that these aforementioned droplet morphology agree well with the experimental observations of Theofanous et al. [38]. Fig. 17 shows the corresponding flow features at these two times, with respect to the pressure distribution in the cross-section of y = 0 and the vorticity distribution in the cross-section of z = 0. We can see that the two vorticity streams with 27

opposite signs are generated due to the baroclinic effect, and detach from the droplet because of flow convection. The relatively low pressure in the recirculation region leads to the global deformation of the droplet [39]. It is also interesting to see that there is a complex propagation of transmitted waves inside the droplet, which may have effect on the occurrence of the surface groove. 6. Conclusions A second-order 3D CSI method has been proposed for the simulation of compressible inviscid two-phase flows, by a combination of 3D cut-cell methods and finite volume methods in the ALE framework. The interfaces are reconstructed on a Cartesian mesh, and effectively split the Cartesian cells that contain interfaces. Utilizing the symmetry and rotational symmetry of fluid positions in a cut cell allows us to significantly reduces the complexity in the subsequent cell assembly. As a result, unstructured interface cells are generated in the vicinity of the interface, while the regular cells are maintained in the bulk fluid away from the interface. Therefore, a second-order finite volume methods in the ALE framework are used for the discretization of the Euler equations in the whole domain, except for the interface cells with under-resolved interfaces or topology changes, at which a first-order finite volume method is used. The performance of the 3D CSI method has been systematically examined. The convergence study of the Rayleigh collapse of a gas bubble in water shows that secondorder accuracy has been achieved. Numerical results of axisymmetric shock-induced bubble collapse demonstrate that the axial symmetry of the 3D interface is well preserved in the 3D computation. Moreover, the interaction between shock and 3D curved gas-gas interface and droplet deformation after shock impact are numerically investigated, and the numerical results are in good agreement with the experimental data, quantitatively and qualitatively. Therefore, this second-order 3D CSI method is able to deal with compressible two-phase flows with large density ratio and topology changes. Also, this is the first time that 3D CSI methods are given experimental validation for compressible two-phase flows. In this work, we only consider 3D CSI simulations of compressible inviscid flows, in which fluid viscosity and surface tension are neglected. It is known that sharp treatment of interface is particularly important for the prediction of interfacial instability for viscous flows [23], and the extension of the present algorithm incorporating fluid viscosity and surface tension is straightforward. This is beyond the scope of this work, and will be considered in our future work. Acknowledgements We are grateful for the support of the National Nature Science Foundation of China (grant nos 11425210, 11672288, 11702288, 11932019), the Chinese Academy of Sciences (grant no. XDB22040103), the Science Challenge Project (Grant No. TZ2016001) and the Fundamental Research Funds for the Central Universities. [1] T. Johnsen, E. and. Colonius, Shock-induced collapse of a gas bubble in shockwave lithotripsy, J. Acoust. Soc. Am. 124 (4) (2008) 2011–2020.

28

[2] Q. Wang, K. Manmi, M. L. Calvisi, Numerical modeling of the 3d dynamics of ultrasound constrast agent microbubbles using the boundary integral method, Phys. Fluids 27 (2) (2015) 022104. [3] T. Theofanous, Aerobreakup of newtonian and viscoelastic liquids, Annu. Rev. Fluid Mech. 43 (2011) 661–690. [4] B. E. Gelfand, Droplet breakup phenomena in flows with velocity lag, Prog. Energy Combust. Sci. 22 (1996) 201–265. [5] J. Glimm, X. Li, Y. Liu, Z. Xu, N. Zhao, Conservative front tracking with improved accuracy, SIAM J. Numer. Anal. 41 (5) (2003) 1926–1947. [6] P. Barton, B. Obadia, D. Drikakis, A conservative level-set based method for compressible solid/fluid problems on fixed grids, J. Comput. Phys. 230 (2011) 7867–7890. [7] L. Schneiders, D. Hartmann, M. Meinke, W. Schr¨ oder, An accurate moving boundary formulation in cut-cell methods, J. Comput. Phys. 235 (2013) 786–809. [8] N. Gokhale, N. Nikiforakis, R. Klein, A dimensionally split cartesian cut cell method for hyperbolic conservation laws, J. Comput. Phys. 364 (2018) 186–208. [9] M. Berndt, J. Breil, S. Galera, M. Kucharik, P.-H. Maire, M. Shashkov, Two-step hybrid conservative remapping for multimaterial arbitrary lagrangian-eulerian methods, J. Comput. Phys. 230 (2011) 6664– 6687. [10] C. W. Hirt, A. A. Amsden, J. L. Cook, An arbitrary Lagrangian-Eulerian computing method for all flow speeds, J. Comput. Phys. 14 (1974) 227–253. [11] X. Y. Hu, B. C. Khoo, N. A. Adams, F. L. Huang, A conservative interface method for compressible flows, J. Comput. Phys. 219 (2) (2006) 553–578. [12] E. Lauer, X. Y. Hu, S. Hickel, N. A. Adams, Numerical modelling and investigation of symmetric and asymmetric cavitation bubble dynamics, Comput. Fluids 69 (2012) 1–19. [13] L. Han, X. Hu, N. Adams, Scale separation for multi-scale modeling of free-surface and two-phase flows with the conservative sharp interface method, J. Comput. Phys. 280 (2015) 387–403. [14] J. Luo, X. Hu, N. Adams, Efficient formulation of scale separation for multi-scale modeling of interfacial flows, J. Comput. Phys. 308 (2016) 411–420. [15] S. Pan, L. Han, X. Hu, N. A. Adams, A conservative interface-interaction method for compressible multi-material flows, J. Comput. Phys. 371 (2018) 870–895. [16] C.-H. Chang, X. Deng, T. G. Theofanous, Direct numerical simulation of interfacial instabilities: A consistent, conservative, all-speed, sharp-interface method, J. Comput. Phys. 242 (2013) 946–990. [17] J.-Y. Lin, Y. Shen, H. Ding, N.-S. Liu, X.-Y. Lu, Simulation of compressible two-phase flows with topology change of fluid-fluid interface by a robust cut-cell method, J. Comput. Phys. 328 (2017) 140–159. [18] H. Kim, M.-S. Liou, Adaptive cartesian cut-cell sharp interface method (aC3 SIM) for three-dimensional multi-phase flows, Shock Waves (2019) https://doi.org/10.1007/s00193–019–00902–6. [19] R. R. Nourgaliev, T. N. Dinh, Adaptive characteristics-based matching for compressible multifluid dynamics, J. Comput. Phys. 213 (2) (2006) 500–529. [20] S. Osher, J. A. Sethian, Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79 (1) (1988) 12–49. [21] B. van Leer, Towards the ultimate conservative difference scheme V. A second-order seqeual to Godunov’s method, J. Comput. Phys. 32 (1) (1979) 101–136. [22] B. van Leer, Towards the ultimate conservative difference scheme ii: Monotonicity and conservation combined in a second order scheme, J. Comput. Phys. 14 (4) (1974) 361–370. [23] R. R. Nourgaliev, M.-S. Liou, T. G. Theofanous, Numerical prediction of interfacial instabilities: Sharp interface method (SIM), J. Comput. Phys. 227 (8) (2008) 3940–3970. [24] T. J. Barth, D. C. Jespersen, The design and application of upwind schemes on unstructured meshes, AIAA Paper (1989) 89–0366. [25] E. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, 2009. [26] M.-S. Liou, A sequel to AUSM, part II : AUSM+ -up for all speeds, J. Comput. Phys. 214 (2006) 137–170.

29

[27] R. Courant, K. Friedrichs, H. Lewy, On the Partial Difference Equations of Mathematical Physics, IBM J. Res. Dev 11 (2) (1967) 215–234. [28] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys. 77 (2) (1988) 439–471. [29] S. Osher, R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Springer, 2003. [30] D. Adalsteinsson, J. A. Sethian, The fast construction of extension velocities in level set methods, J. Comput. Phys. 148 (1999) 2–22. [31] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1) (1994) 146–159. [32] G. Russo, P. Smereka, A remark on computing distance functions, J. Comput. Phys. 163 (2000) 51–67. [33] C. Min, On reinitializing level set finctions, J. Comput. Phys. 229 (2010) 2764–2772. [34] M. J. Ivings, D. M. Causon, E. F. Toro, On Riemann solvers for compressible liquids, Int. J. Numer. Meth. Fluids 28 (3) (1998) 395–418. [35] Z.-J. Wang, Improved formulation for geometric properties of arbitrary polyhedra, AIAA J. 37 (10) (1999) 1326–1327. [36] N. A. Hawker, Y. Ventikos, Interaction of a strong shockwave with a gas bubble in a liquid medium: A numerical study, J. Fluid Mech. 701 (2012) 59–97. [37] J. Ding, T. Si, M. Chen, Z. Zhai, X. Lu, X. Luo, On the interaction of a planar shock with a threedimensional light gas cylinder, J. Fluid Mech. 828 (2017) 289–317. [38] T. G. Theofanous, V. V. Mitkin, C. L. Ng, C.-H. Chang, X. Deng, S. Sushchikh, The physics of aerobreakup. II. viscous liquids, Phys. Fluids 24 (2012) 022104. [39] J. C. Meng, T. Colonius, Numerical simulation of the aerobreakup of a water droplet, J. Fluid Mech. 835 (2018) 1108–1135.

30

Declaration of interests ‫ ܈‬The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ‫܆‬The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Conflict of interest Statement None