- Email: [email protected]
Adaptive unstructured mesh modelling of multiphase flows
Accepted Manuscript Adaptive Unstructured Mesh Modelling of Multiphase Flows Zhihua Xie, Dimitrios Pavlidis, James R. Percival, Jefferson L.M. A. Gomes, Christopher C. Pain, Omar K. Matar PII: DOI: Reference:
S0301-9322(14)00143-8 http://dx.doi.org/10.1016/j.ijmultiphaseflow.2014.08.002 IJMF 2083
To appear in:
International Journal of Multiphase Flow
Received Date: Revised Date: Accepted Date:
16 April 2014 4 August 2014 7 August 2014
Please cite this article as: Xie, Z., Pavlidis, D., Percival, J.R., Gomes, J.L.M., Pain, C.C., Matar, O.K., Adaptive Unstructured Mesh Modelling of Multiphase Flows, International Journal of Multiphase Flow (2014), doi: http:// dx.doi.org/10.1016/j.ijmultiphaseflow.2014.08.002
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.
Adaptive Unstructured Mesh Modelling of Multiphase Flows Zhihua Xiea,b,∗, Dimitrios Pavlidisa , James R. Percivala , Jefferson L. M. A. Gomesc , Christopher C. Paina , Omar K. Matarb a
Department of Earth Science and Engineering, Imperial College London, UK b Department of Chemical Engineering, Imperial College London, UK c School of Engineering, University of Aberdeen, UK
Abstract Multiphase flows are often found in industrial and practical engineering applications, including bubbles, droplets, liquid film and waves. An adaptive unstructured mesh modelling framework is employed here to study interfacial flow problems, which can modify and adapt anisotropic unstructured meshes to better represent the underlying physics of multiphase problems and reduce computational effort without sacrificing accuracy. The numerical framework consists of a mixed control volume and finite element formulation, a ‘volume of fluid’-type method for the interface capturing based on a compressive control volume advection method and second-order finite element methods. The framework also features a force-balanced algorithm for the surface tension implementation, minimising the spurious velocities often found in such flows. Numerical examples of the Rayleigh-Taylor instability and a rising bubble are presented to show the ability of this adaptive unstructured mesh modelling framework to capture complex interface geometries and also to increase the efficiency in multiphase flow simulations. Keywords: Anisotropic mesh adaptivity, interface-capturing, Rayleigh-Taylor instability, rising bubble, surface tension, two-phase flows
∗
Corresponding author. Email: [email protected]; Tel.: +44(0)20-75943067.
Preprint submitted to International Journal of Multiphase Flow
August 11, 2014
1. Introduction Multiphase flow phenomena with deformable interfaces appear in many engineering applications, e.g. micro-fluidics, oil-and-gas transportation systems, reservoir engineering, porous media flows, nuclear reactors, biomedical applications, geophysical flows, coastal and ocean engineering, and marine renewable energy. These applications typically involve the motion of bubbles, droplets or particles in a continuous phase fluid, featuring tremendous complexity in interfacial topology and a large range of scales: from microns to kilometres. A key requirement for modelling multiphase flows is a method for calculating the shape of the interface. Numerous methods have been proposed and used to simulate interfacial flows on a fixed mesh (Scardovelli and Zaleski, 1999), such as marker-and-cell (Harlow and Welch, 1965), volume-of-fluid (VOF) (Hirt and Nichols, 1981; Rider and Kothe, 1998; Scardovelli and Zaleski, 1999), front-tracking (Unverdi and Tryggvason, 1992), level set (Osher and Sethian, 1988; Sethian and Smereka, 2003), phase field (Anderson et al., 1998; Jacqmin, 1999; Badalassi et al., 2003; Kim et al., 2004) and particle (Monaghan, 1992) methods. In particular, VOF methods are widely used due to the inherent properties of: mass conservation, computational efficiency and easy implementation. From a general point of view, there are two classes of algorithms to solve the transport equation of volume fraction: geometric and algebraic computation (Rider and Kothe, 1998). Most geometric VOF methods are based on a two-stage process (Scardovelli and Zaleski, 1999). Firstly, interfaces are reconstructed from the volume fraction data so that a geometric profile is found which approximates the actual interface location. Then changes in volume fraction are calculated by integrating volume fluxes across cell boundaries, using the geometric profile of the reconstructed interface. In the algebraic computation (Rudman, 1997; Ubbink and Issa, 1999; Chen et al., 2001; So et al., 2011), the interface is captured by solving the transport equation of volume fraction with a differencing scheme. As the dynamics of interfacial flows are highly unsteady and the shape and location of the interface are changing during the simulation, interface calculation methods based on a fixed mesh have to produce a very fine mesh across the whole computational domain in order to capture the details, which will significantly increase computational efforts. The alternative is to consider the use of dynamically adaptive mesh methods, where the mesh resolution can vary in time in response to the evolving solution fields. For example, 2
finer mesh could be placed around the interface during its development while coarser mesh could be used away from the interface while the flow is less dynamic. There are some examples for the use of adaptive mesh refinement for structured meshes with VOF (Popinet, 2009) and hybrid level set/front tracking (Ceniceros et al., 2010) methods, as well as adaptive unstructured meshes with the level set method (Zheng et al., 2005). Unstructured meshes are very attractive to deal with complex geometries in engineering applications. In addition, it is worth mentioning that many computational domains have a high aspect ratio and most multiphase flow phenomena can possess strong anisotropies, therefore mesh resolution may be required to optimally represent the dynamics of the flow. Some examples of anisotropic unstructured mesh adaptivity can be found in Piggott et al. (2009) The motivation for this work is to exploit state-of-the-art methods for multiphase flow modelling with interface-capturing numerical schemes on adaptive unstructured meshes, which can modify and adapt unstructured meshes to better represent the underlying physics of multiphase problems and reduce computational efforts without sacrificing accuracy. The remainder of this paper is organised as follows. Description of the model and numerical methods is given in Section 2. Numerical examples are presented in Section 3. Finally, some concluding remarks and future work are given in Section 4. 2. Mathematical model and numerical methods A multi-fluid modelling framework has been developed based on the multi-component modelling approach with information on interfaces embedded into the continuity equations. Advantages of this approach are: (i) The model is designed for multi-fluid flows with arbitrary number of fluids/materials, and does not rely on priority list of materials (common for traditional multi-material models), which often produces spurious solutions (Wilson, 2009); (ii) The component advection equations are embedded into both pressure and continuity equations resulting in local mass balance (within each control volume).
3
2.1. Governing equations In multi-component flows, a number of components exist in one or more phases (one phase is assumed here but is easily generalised to an arbitrary number of phases or fluids). Let αi be the mass fraction of component i, where i = 1, 2, .., Nc and Nc denotes the number of components. The density and dynamic viscosity of component i are ρi and µi . A constraint on the system is: Nc X αi = 1. (1) i=1
For each fluid component i, the conservation of mass may be defined as, ∂ (αi ) + ∇ · (αi u) = 0, i = 1, 2, .., Nc , ∂t
(2)
and the equations of motion of an incompressible fluid may be written as: ∂(ρu) + ∇ · (ρu ⊗ u) = −∇p + ∇ · [µ(∇u + ∇T u)] + ρg + Fσ , ∂t
(3)
where t P is the time, u is velocity vector, p is the pressure, bulk density PNthe Nc c is ρ = α ρ , the bulk dynamic viscosity is µ = α i=1 i i i=1 i µi , g is the gravitational acceleration vector, and Fσ is the surface tension force. 2.2. Numerical methods The numerical framework consists of a mixed control volume and finite element formulation. In the mixed formulation, the domain is discretised into triangular or tetrahedral elements and in this work, they are P1 DGP2 elements (linear discontinuous velocity between elements and quadratic continuous pressure between elements) (Cotter et al., 2009b). Fig. 1 shows the locations of the degrees of freedom for the P1 DG-P2 element and the boundaries of the control volumes.
4
Pressure Velocity
(a)
(b)
Figure 1: (a) Finite element used to discretise the fluids equations. The central position of key solution variables are indicated here for the P1 DGP2 element pair (Cotter et al., 2009b). (b) Diagram showing the relationship between intersecting control volumes (shaded area) and elements for the P2 element.
When testing Eq. (2) with a control volume (piece-wise constant) basis function Mm and applying integration by parts over each control volume m with the θ-time stepping for the advection terms, and then summing the result over all components, the global continuity equation can be obtained as (Pavlidis et al., 2014): Z
Γm
n+ 12 n+ 12 n n+1 bm n · u n · u dΓ = + (1\ − Θ) Θ m
Z
Vm
1
Mm Sm n+ 2 dV, m = 1, 2, .., M,
(4)
where n is the current time level, n is the outward-pointing unit normal vector to the surface Γm of the control volume Vm , and M is the number of control volumes. n+ 1 The absorption term Sm 2 (from Eq. (4)) is: ! N Nc c n X X 1 αi n+1 1 n+ 2 n+1 m − αi m , (5) 1− αi m − = −wc Sm ∆t ∆t i=1 i=1
5
where the term wc ∈ {0, 1} helps to relax to a solution satisfying the component summation constraint (Eq. (1)) and ∆t is the time step. The space/time flux-limiting functions (from Eq. (4)) are defined as: 1
and:
2 b n+ = Θ m
n+ 12
(1\ − Θ)m
=
Nc X
n+ 1
ˆ i n+1 θi m 2 α m ,
(6)
i=1
Nc X i=1
n+ 12
1 − θi m
α ˆ i nm .
(7)
The initial estimate of these flux-limiting functions in a time step are: n+ 12 n+ 12 \ b = 0. Θm = 1 and (1 − Θ) m
where θ ∈ {0, 1} is the implicitness parameter and the hat decoration (ˆ) represents flux-limited values of the indicated solution variables at the quadrature points on the surface of the control volumes. A transient, mixed, finite element formulation is used to discretise the momentum equations (Eq. (3)). A finite volume discretisation of the continuity equations and a linear discontinuous Galerkin (DG) (Cotter et al., 2009a; Xiao et al., 2013) discretisation of the momentum equations are employed with backward Euler time stepping. Within each time step, the equations are iterated upon using a projection-based pressure determination method until all equations are simultaneously balanced. A non-linear Petrov-Galerkin method is used to stabilise the momentum equations. This involves the addition of a non-linear diffusion term over and above a standard BubnovGalerkin DG discretisation (see Hughes and Mallet, 1986; Pain et al., 2001; Fang et al., 2013, for further details). The main novelty of the framework includes: (a) the new finite element type (P1 DG-P2 ) for multi-fluid flow problems, which ensures exact balance between buoyancy force and pressure gradient; (b) a novel interface capturing scheme based on compressive control volume advection method (Pavlidis et al., 2014), involving a high-order accurate finite element method to obtain fluxes on the control volume boundaries, where these fluxes are subject to flux-limiting using a 6
normalised variable diagram approach (Leonard, 1991) to obtain bounded and compressive solutions for the interface; (c) a forcebalanced algorithm for the surface tension implementation in unstructured mesh methods, minimising the spurious velocities often found in interfacial flows (Francois et al., 2006); (d) anisotropic unstructured mesh adaptivity, which allows the grid resolution to be concentrated in relatively important regions such as the vicinity of interfaces, while lower resolution can be used in other regions. 3. Numerical examples Prior to attempting calculations of interfacial flows, several validation checks of our numerical framework were performed. In this section we demonstrate the convergence properties of the present method and also show the potential of using anisotropic unstructured mesh adaptivity for multiphase flow problems. All the numerical results presented here have been obtained using the P1 DG-P2 finite elements in two spatial dimensions. Here the mesh refinement criterion is based on minimising the interpolation error estimate for the component volume fraction. Between adaptive meshes, Galerkin interpolation is used for the velocity and consistent interpolation is used for the pressure and scalar fields. 3.1. Rayleigh-Taylor instability The Rayleigh-Taylor instability will occur if a heavy fluid initially lies on top of a lighter fluid in a gravitational field. It has been a benchmark test case in the computational fluid dynamics community for numerical methods to study multiphase flow problems. The initial growth and evolution of Rayleigh-Taylor instability have been studied for inviscid (Tryggvason, 1988) and viscous (Guermond and Quartapelle, 2000; Ding et al., 2007) incompressible flows. Here we use the same setup for a rectangular domain [0, d] × [0, 4d] with the density difference, represented by the Atwood ratio At = (ρA − ρB )/(ρA + ρB ), equals to 0.5 and with the Reynolds number ReA = ρA d3/2 g 1/2 /µA = 3000, where ρA and ρB denote the densities of the heavier and light fluids, respectively and µA is the dynamic viscosity of the heavier fluid. The interface between the two fluids is an initially sinusoidal perturbation of amplitude 0.1d. Following Tryggvason (1988), we use the same nondimensional variables, where the length is scaled by d, time is
7
1.5 Tryggvason (1988) Guermond and Quartapelle (2000) Ding et al. (2007) coarse mesh fine mesh
1
Position
0.5
0
−0.5
−1
−1.5
0
0.5
1.5
1
2
2.5
t
Figure 2: Position of the tip of the rising and falling fluid versus time for the RayleighTaylor instability with structured mesh simulations with At = 0.5 and ReA = 3000.
p √ scaled by d/Atg and the velocity is scaled by Atgd. In this case, there is no surface tension, and the ‘interface’ evolves due to gravity. First, computations were carried out using two fixed structured triangular meshes, namely a coarse mesh having 32 elements in the x direction and 128 elements in the y direction with 8192 (32×128×2) elements in total, and a fine mesh (64×256×2=32768 elements) which doubles the resolution of the coarse mesh in each direction. As expected, the heavier fluid falls into the light fluid at the centre while the light fluid rises up at both sides of the domain. Fig. 2 shows the computational results of the evolution of the positions of the tip of the rising and falling fluid, together with corresponding results obtained by the Lagrangian-Eulerian vortex method (Tryggvason, 1988), the variable density FEM method (Guermond and Quartapelle, 2000) and the diffuse interface method (Ding et al., 2007). Inspection of Fig. 2 reveals that convergence has been achieved upon mesh-refinement and that agreement with previous work is also excellent. Numerical simulations were then performed on adaptive, unstructured meshes. Here, the mesh refinement criteria are: the interpola-
8
tion error bound is 0.05, and the specified minimun and maximum length scale are 0.005d and 0.5d, respectively. The evolution of the interface is shown in Fig. 3, where the fixed fine mesh result is shown in the left column in panels (b)-(f) and the adaptive mesh result together with corresponding mesh is shown in the middle and right columns in panels (b)-(f). In order to facilitate comparison, adaptive simulations were started from the same initial profile and mesh for the fixed mesh case at t = 0 in Fig. 3(a). It can be seen that two counter-rotating vortices develop when the heavier fluid falls down. Inspection of Fig. 3 reveals the level of refinement of the unstructured mesh in the interfacial region. This refinement is adaptive and targets the regions of rapid variation in interfacial structure associated with the roll-up phenomena that are known to accompany the development of Rayleigh-Taylor instabilities driven by density contrasts. It can be seen from Fig. 3 that similar results are obtained for both fixed structured and adapative unstructured meshes, which confirms the overall consistency of the adaptive simulation when having varying mesh resolution across the whole flow domain. It is also worth noting that the total number of elements is kept constant as 32768 for the fixed mesh simulation. For the adaptive simulation, the number of elements reduces to 1137 at the first mesh adaption as fine meshes are only needed around the interface while coarse meshes can be placed far away from the interface. The total number of elements gradually increases to 5967 (see Fig. 3(f)) in order to capture the details of the vortex sheet. This demonstrates that in order to achieve the same order of accuracy fewer elements are required using the adaptive mesh, and this significantly reduces the computational efforts. It is also noted that during the adaptive simulation, although the Courant number is larger than unity in some regions with small meshes, the interface is well calculated due to the implicit time stepping for the interface-capturing method. Fig. 4 illustrates the detailed interface shape and mesh distribution at t = 2.0 together with its enlarged version from the left pannel to the right one. It is shown that when the heavier fluid sinks to the lower section of the domain its central region becomes elongated and possesses strong anisotropies. In accordance with this, anisotropic meshes are placed within the vicinity of the interface, thus the method has focused resolution where it is required to optimally represent the deveopment of the interface.
9
(a) t = 0
(b) t = 0.005
(c) t = 1.0
(d) t = 1.25
(e) t = 1.5
(f) t = 1.75
10 Figure 3: Numerical simulation of the Rayleigh-Taylor instability for the same parameters as those used to generate Fig. 2 showing a comparison between structured (left column in panels (b)-(f)) and adaptive unstructured meshing (middle and right columns in panels (b)-(f)). In (a), we show the condition used to initialise all simulations.
Figure 4: Detailed interface shape and mesh profile of the Rayleigh-Taylor instability at t = 2.0. Figures are gradually enlarged from left to right in order to show the anisotropic unstructured mesh around the interface.
3.2. Rising bubble Next we consider the rise of an air bubble in water, as accurate modelling of such flow is challenging due to its deforming interface as well as its surface tension force. Here we follow the same nondimensional setup from previous studies (Francois et al., 2006; Popinet, 2009), where an air bubble with diameter of D = 2/3 is initially placed at (0, 0) in a computational domain [−1, 1] × [−1, 2]. The fluid properties are: ρair = 1.126, ρwater = 1000, µair = 1.137, µwater = 1.78 × 10−5 . The suface tension coefficient σ = 728 2 = 5.983 for water. A force-balanced which gives a Bond number Bo = ρgD σ algorithm was implemented in the model for the surface tension calculation on unstructured meshes, where a diffused interface, reconstructed from the volume fraction on the unstructured meshes, provided good estimation of the interfacial curvature κ in Eq. (3). It is worth noting that these computations are 2D, where curvatures are different compared to 3D cases. Computations were carried out using a fixed structured triangular mesh, having 80 elements in the x direction and 120 elements in the y direction with 19200 (80 × 120 × 2) elements in total, and an adaptive mesh. The mesh 11
refinement criteria are: the interpolation error bound is 0.05, and the specified minimun and maximum length scale are 0.001 and 0.5, respectively. In Fig. 5, it can be seen that the initially spherical bubble rises under the action of buoyancy, and undergoes deformation, resulting in the formation of the well-known cap-shaped bubbles. The results associated with structured, and adaptive unstructured meshes are shown in the left, and middle and right columns of panels (b)-(e); the initial condition for all the simulations is shown in panel (a). It is shown from Fig. 5 that both results closely match with each other until t = 0.35, and they differ slightly after that mainly because the mesh resolution for the fixed mesh simulation is not fine enough to capture the shape of interface when the bubble has large deformation. In contrast, for the adaptive mesh simulation, inspection of Fig. 5 reveals that the unstructured mesh is not dense around the bubble at early times since the latter had not yet been deformed by the flow. As the extent of deformation becomes increasingly pronounced, our numerical method responds by refining the unstructured mesh in the vicinity of the deforming bubble preferentially; an enlarged version of Fig. 5(e) is shown in Fig. 6, which depicts the details of the unstructured mesh around the bubble. In contrast, the regions upstream and downstream of the bubble, which require a lower degree of resolution, have coarser meshes, in order to maximise computational efficiency. It can also be seen from Fig. 6 that the sharp interface is well captured during the adaptive simulation, due to the use of the compressive control volume advection method developed in the framework. In order to demonstrate the efficiency of the adaptive unstructured mesh simulation, Fig. 7 shows the time history of the total number of elements during the simulations. It can be seen that the total number decreases significantly at the first mesh adaption, approximately only 10% of the fixed mesh case. The maximum number of elements used in the adaptive simulation is only 18% of the fixed mesh case, which reduces the computational efforts without sacrificing accuracy. 4. Conclusions In this paper an adaptive unstructured mesh modelling approach, which can reduce computational efforts without sacrifing accuracy, has been presented for interface-capturing in multiphase flow problems with an emphasis on two-phase flows. The interface-capturing method is based on a compres12
(a) t = 0
(b) t = 0.001
(c) t = 0.2
(d) t = 0.35
13 (e) t = 0.5 Figure 5: Numerical simulation of a rising bubble with a Bond number Bowater = 5.983. Here, we show a comparison of the results associated with structured, and unstructured adaptive meshes shown in the left, and middle and right columns of panels (b)-(e), respectively. In (a), we show the condition that was used to initialise all simulations.
Figure 6: An enlarged version of Fig. 5(e) showing the details of the numerical solution obtained using an adaptive, unstructured mesh.
4
number of elements
2
x 10
fixed mesh adaptive mesh
1.5
1
0.5
0
0
0.1
0.3
0.2
0.4
0.5
t
Figure 7: Numbers of total elements for the numerical simulations of a rising bubble with a fixed and adaptive, unstructured mesh. The parameters are the same as those used to generate Fig. 5 and Fig. 6.
14
sive advection method, which uses a novel and mathematically rigorous nonlinear Petrov-Galerkin method that attempts to keep the interface sharp. This was shown to perform well in a standard Rayleigh-Taylor instability benchmark problem as well as bubble rising two-phase flow problems. We also presented results using the new finite element type (P1 DG-P2 ) for multiphase flow problems. The element type ensures that linear buoyancy force and vertical pressure gradient are exactly balanced which eliminates large spurious velocities that might be generated on unstructured meshes if this balance is inaccurate. We also showed the method working with anisotropic mesh adaptivity, which allows the grid resolution to be concentrated in relatively important regions such as the vicinity of interfaces while lower resolution can be used in other regions. In addition, the interface-capturing method is implicit which means that it is stable with relatively large time step sizes exceeding unity Courant numbers. The results presented here established with sufficient confidence that this method can be used to successfully model multiphase flows in a wide range of applications. This approach has the potential to be used for an arbitrary number of components although that has not been demonstrated here. Altough only 2D simulations are considered here, extention to 3D is straightforward. Currently we are working on parallel computing and some fully 3D simulations will be reported in the near future. 5. Acknowledgements We would like to thank the EPSRC MEMPHIS multiphase Programme Grant (number EP/K003976/1), the EPSRC Computational modelling for advanced nuclear power plants project and the EU FP7 projects THINS and GoFastR for helping to fund this work. References Anderson, D.M., McFadden, G.B., Wheeler, A.A., 1998. Diffuse-interface methods in fluid mechanics. Annual Review of Fluid Mechanics 30, 139– 165. Badalassi, V.E., Ceniceros, H.D., Banerjee, S., 2003. Computation of multiphase systems with phase field models. Journal of Computational Physics 190, 371–397. 15
Ceniceros, H.D., Roma, A.M., Silveira-Neto, A., Villar, M.M., 2010. A robust, fully adaptive hybrid level-set/front-tracking method for two-phase flows with an accurate surface tension computation. Communications in Computational Physics 8, 51–94. Chen, Y.G., Price, W.G., Temarel, P., 2001. An anti-diffusive volume of fluid method for interfacial fluid flows. International Journal for Numerical Methods in Fluids 68, 341–359. Cotter, C.J., Ham, D.A., Pain, C.C., 2009a. A mixed discontinuous/continuous finite element pair for shallow-water ocean modelling. Ocean Modelling 26, 86–90. Cotter, C.J., Ham, D.A., Pain, C.C., Reich, S., 2009b. LBB stability of a mixed Galerkin finite element pair for fluid flow simulations. Journal of Computational Physics 228, 336–348. Ding, H., Spelt, P.D.M., Shu, C., 2007. Diffuse interface model for incompressible two-phase flows with large density ratios. Journal of Computational Physics 226, 2078–2095. Fang, F., Pain, C.C., Navon, I.M., ElSheikh, A.H., Du, J., Xiao, D., 2013. Non-linear Petrov-Galerkin methods for reduced order hyperbolic equations and discontinuous finite element methods. Journal of Computational Physics 234, 540–559. Francois, M.M., Cummins, S.J., Dendy, E.D., Kothe, D.B., Sicilian, J.M., Williams, M.W., 2006. A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework. Journal of Computational Physics 213, 141–173. Guermond, J.L., Quartapelle, L., 2000. A projection FEM for variable density incompressible flows. Journal of Computational Physics 165, 167–188. Harlow, F.H., Welch, J.E., 1965. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Physics and Fluids 8, 2182–2189. Hirt, C.W., Nichols, B.D., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics 39, 201– 225. 16
Hughes, T.J.R., Mallet, M., 1986. A new finite element formulation for computational fluid dynamics - IV. A discontinuity-capturing operator for multi-dimensional advection-diffusion systems. Computer Methhods in Applied Mechanics and Engineering 58, 329–336. Jacqmin, D., 1999. Calculation of two-phase Navier-Stokes flows using phasefield modeling. Journal of Computational Physics 155, 96–127. Kim, J.S., Kang, K., Lowengrub, J., 2004. Conservative multigrid methods for Cahn-Hilliard fluids. Journal of Computational Physics 193, 511–543. Leonard, B.P., 1991. The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection 88, 17–74. Monaghan, J.J., 1992. Smoothed particle hydrodynamics. Annual Review of Astronomy and Astrophysics 30, 543–574. Osher, S.J., Sethian, J.A., 1988. Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics 79, 12–49. Pain, C.C., de Oliveira, C.R.E., Goddard, A.J.H., Umpleby, A.P., 2001. Criticality behaviour of dilute plutonium solutions. Nuclear Science and Technology 30, 194–214. Pavlidis, D., Xie, Z., Percival, J.R., Gomes, J.L.M.A., Pain, C.C., Matar, O.K., 2014. Two- and three-phase horizontal slug flow modelling using an interface-capturing compositional approach. International Journal of Multiphase Flow. Accepted. Piggott, M.D., Farrell, P.E., Wilson, C.R., Gorman, G.J., Pain, C.C., 2009. Anisotropic mesh adaptivity for multi-scale ocean modelling. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 367, 4591–4611. Popinet, S., 2009. An accurate adaptive solver for surface-tension-driven interfacial flows. Journal of Computational Physics 228, 5838–5866. Rider, W.J., Kothe, B., 1998. Reconstructing volume tracking. Journal of Computational Physics 141, 112–152.
17
Rudman, M., 1997. Volume-tracking methods for interfacial flow calculations. International Journal for Numerical Methods in Fluids 24, 671–691. Scardovelli, R., Zaleski, S., 1999. Direct numerical simulation of free-surface and interfacial flow. Annual Review of Fluid Mechanics 31, 567–603. Sethian, J.A., Smereka, P., 2003. Level set methods for fluid interfaces. Annual Review of Fluid Mechanics 35, 341–372. So, K.K., Hu, X.Y., Adams, N.A., 2011. Anti-diffusion method for interface steepening in two-phase incompressible flow. Journal of Computational Physics 230, 5155–5177. Tryggvason, G., 1988. Numerical simulations of the Rayleigh-Taylor instability. Journal of Computational Physics 75, 253–282. Ubbink, O., Issa, R.I., 1999. A method for capturing sharp fluid interfaces on arbitrary meshes. Journal of Computational Physics 153, 26–50. Unverdi, S.O., Tryggvason, G., 1992. A front-tracking method for viscous, incompressible, multi-fluid flows. Journal of Computational Physics 100, 25–37. Wilson, C., 2009. Modelling multiple-material flows on adaptive unstructured meshes. Ph.D. thesis. Department of Earth Science and Engineering, Imperial College London. Xiao, D., Fang, F., Du, J., Pain, C.C., Navon, I.M., Buchan, A.G., ElSheikh, A.H., Hu, G., 2013. Non-linear Petrov-Galerkin methods for reduced order modelling of the Navier-Stokes equations using a mixed finite element pair. Computer Methods in Applied Mechanics and Engineering 255, 147–157. Zheng, X., Lowengrub, J., Anderson, A., 2005. Adaptive unstructured volume remeshing ii: Application to two- and three-dimensional level-set simulations of multiphase flow. Journal of Computational Physics 208, 626–650.
18
S0301-9322(14)00143-8 http://dx.doi.org/10.1016/j.ijmultiphaseflow.2014.08.002 IJMF 2083
To appear in:
International Journal of Multiphase Flow
Received Date: Revised Date: Accepted Date:
16 April 2014 4 August 2014 7 August 2014
Please cite this article as: Xie, Z., Pavlidis, D., Percival, J.R., Gomes, J.L.M., Pain, C.C., Matar, O.K., Adaptive Unstructured Mesh Modelling of Multiphase Flows, International Journal of Multiphase Flow (2014), doi: http:// dx.doi.org/10.1016/j.ijmultiphaseflow.2014.08.002
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.
Adaptive Unstructured Mesh Modelling of Multiphase Flows Zhihua Xiea,b,∗, Dimitrios Pavlidisa , James R. Percivala , Jefferson L. M. A. Gomesc , Christopher C. Paina , Omar K. Matarb a
Department of Earth Science and Engineering, Imperial College London, UK b Department of Chemical Engineering, Imperial College London, UK c School of Engineering, University of Aberdeen, UK
Abstract Multiphase flows are often found in industrial and practical engineering applications, including bubbles, droplets, liquid film and waves. An adaptive unstructured mesh modelling framework is employed here to study interfacial flow problems, which can modify and adapt anisotropic unstructured meshes to better represent the underlying physics of multiphase problems and reduce computational effort without sacrificing accuracy. The numerical framework consists of a mixed control volume and finite element formulation, a ‘volume of fluid’-type method for the interface capturing based on a compressive control volume advection method and second-order finite element methods. The framework also features a force-balanced algorithm for the surface tension implementation, minimising the spurious velocities often found in such flows. Numerical examples of the Rayleigh-Taylor instability and a rising bubble are presented to show the ability of this adaptive unstructured mesh modelling framework to capture complex interface geometries and also to increase the efficiency in multiphase flow simulations. Keywords: Anisotropic mesh adaptivity, interface-capturing, Rayleigh-Taylor instability, rising bubble, surface tension, two-phase flows
∗
Corresponding author. Email: [email protected]; Tel.: +44(0)20-75943067.
Preprint submitted to International Journal of Multiphase Flow
August 11, 2014
1. Introduction Multiphase flow phenomena with deformable interfaces appear in many engineering applications, e.g. micro-fluidics, oil-and-gas transportation systems, reservoir engineering, porous media flows, nuclear reactors, biomedical applications, geophysical flows, coastal and ocean engineering, and marine renewable energy. These applications typically involve the motion of bubbles, droplets or particles in a continuous phase fluid, featuring tremendous complexity in interfacial topology and a large range of scales: from microns to kilometres. A key requirement for modelling multiphase flows is a method for calculating the shape of the interface. Numerous methods have been proposed and used to simulate interfacial flows on a fixed mesh (Scardovelli and Zaleski, 1999), such as marker-and-cell (Harlow and Welch, 1965), volume-of-fluid (VOF) (Hirt and Nichols, 1981; Rider and Kothe, 1998; Scardovelli and Zaleski, 1999), front-tracking (Unverdi and Tryggvason, 1992), level set (Osher and Sethian, 1988; Sethian and Smereka, 2003), phase field (Anderson et al., 1998; Jacqmin, 1999; Badalassi et al., 2003; Kim et al., 2004) and particle (Monaghan, 1992) methods. In particular, VOF methods are widely used due to the inherent properties of: mass conservation, computational efficiency and easy implementation. From a general point of view, there are two classes of algorithms to solve the transport equation of volume fraction: geometric and algebraic computation (Rider and Kothe, 1998). Most geometric VOF methods are based on a two-stage process (Scardovelli and Zaleski, 1999). Firstly, interfaces are reconstructed from the volume fraction data so that a geometric profile is found which approximates the actual interface location. Then changes in volume fraction are calculated by integrating volume fluxes across cell boundaries, using the geometric profile of the reconstructed interface. In the algebraic computation (Rudman, 1997; Ubbink and Issa, 1999; Chen et al., 2001; So et al., 2011), the interface is captured by solving the transport equation of volume fraction with a differencing scheme. As the dynamics of interfacial flows are highly unsteady and the shape and location of the interface are changing during the simulation, interface calculation methods based on a fixed mesh have to produce a very fine mesh across the whole computational domain in order to capture the details, which will significantly increase computational efforts. The alternative is to consider the use of dynamically adaptive mesh methods, where the mesh resolution can vary in time in response to the evolving solution fields. For example, 2
finer mesh could be placed around the interface during its development while coarser mesh could be used away from the interface while the flow is less dynamic. There are some examples for the use of adaptive mesh refinement for structured meshes with VOF (Popinet, 2009) and hybrid level set/front tracking (Ceniceros et al., 2010) methods, as well as adaptive unstructured meshes with the level set method (Zheng et al., 2005). Unstructured meshes are very attractive to deal with complex geometries in engineering applications. In addition, it is worth mentioning that many computational domains have a high aspect ratio and most multiphase flow phenomena can possess strong anisotropies, therefore mesh resolution may be required to optimally represent the dynamics of the flow. Some examples of anisotropic unstructured mesh adaptivity can be found in Piggott et al. (2009) The motivation for this work is to exploit state-of-the-art methods for multiphase flow modelling with interface-capturing numerical schemes on adaptive unstructured meshes, which can modify and adapt unstructured meshes to better represent the underlying physics of multiphase problems and reduce computational efforts without sacrificing accuracy. The remainder of this paper is organised as follows. Description of the model and numerical methods is given in Section 2. Numerical examples are presented in Section 3. Finally, some concluding remarks and future work are given in Section 4. 2. Mathematical model and numerical methods A multi-fluid modelling framework has been developed based on the multi-component modelling approach with information on interfaces embedded into the continuity equations. Advantages of this approach are: (i) The model is designed for multi-fluid flows with arbitrary number of fluids/materials, and does not rely on priority list of materials (common for traditional multi-material models), which often produces spurious solutions (Wilson, 2009); (ii) The component advection equations are embedded into both pressure and continuity equations resulting in local mass balance (within each control volume).
3
2.1. Governing equations In multi-component flows, a number of components exist in one or more phases (one phase is assumed here but is easily generalised to an arbitrary number of phases or fluids). Let αi be the mass fraction of component i, where i = 1, 2, .., Nc and Nc denotes the number of components. The density and dynamic viscosity of component i are ρi and µi . A constraint on the system is: Nc X αi = 1. (1) i=1
For each fluid component i, the conservation of mass may be defined as, ∂ (αi ) + ∇ · (αi u) = 0, i = 1, 2, .., Nc , ∂t
(2)
and the equations of motion of an incompressible fluid may be written as: ∂(ρu) + ∇ · (ρu ⊗ u) = −∇p + ∇ · [µ(∇u + ∇T u)] + ρg + Fσ , ∂t
(3)
where t P is the time, u is velocity vector, p is the pressure, bulk density PNthe Nc c is ρ = α ρ , the bulk dynamic viscosity is µ = α i=1 i i i=1 i µi , g is the gravitational acceleration vector, and Fσ is the surface tension force. 2.2. Numerical methods The numerical framework consists of a mixed control volume and finite element formulation. In the mixed formulation, the domain is discretised into triangular or tetrahedral elements and in this work, they are P1 DGP2 elements (linear discontinuous velocity between elements and quadratic continuous pressure between elements) (Cotter et al., 2009b). Fig. 1 shows the locations of the degrees of freedom for the P1 DG-P2 element and the boundaries of the control volumes.
4
Pressure Velocity
(a)
(b)
Figure 1: (a) Finite element used to discretise the fluids equations. The central position of key solution variables are indicated here for the P1 DGP2 element pair (Cotter et al., 2009b). (b) Diagram showing the relationship between intersecting control volumes (shaded area) and elements for the P2 element.
When testing Eq. (2) with a control volume (piece-wise constant) basis function Mm and applying integration by parts over each control volume m with the θ-time stepping for the advection terms, and then summing the result over all components, the global continuity equation can be obtained as (Pavlidis et al., 2014): Z
Γm
n+ 12 n+ 12 n n+1 bm n · u n · u dΓ = + (1\ − Θ) Θ m
Z
Vm
1
Mm Sm n+ 2 dV, m = 1, 2, .., M,
(4)
where n is the current time level, n is the outward-pointing unit normal vector to the surface Γm of the control volume Vm , and M is the number of control volumes. n+ 1 The absorption term Sm 2 (from Eq. (4)) is: ! N Nc c n X X 1 αi n+1 1 n+ 2 n+1 m − αi m , (5) 1− αi m − = −wc Sm ∆t ∆t i=1 i=1
5
where the term wc ∈ {0, 1} helps to relax to a solution satisfying the component summation constraint (Eq. (1)) and ∆t is the time step. The space/time flux-limiting functions (from Eq. (4)) are defined as: 1
and:
2 b n+ = Θ m
n+ 12
(1\ − Θ)m
=
Nc X
n+ 1
ˆ i n+1 θi m 2 α m ,
(6)
i=1
Nc X i=1
n+ 12
1 − θi m
α ˆ i nm .
(7)
The initial estimate of these flux-limiting functions in a time step are: n+ 12 n+ 12 \ b = 0. Θm = 1 and (1 − Θ) m
where θ ∈ {0, 1} is the implicitness parameter and the hat decoration (ˆ) represents flux-limited values of the indicated solution variables at the quadrature points on the surface of the control volumes. A transient, mixed, finite element formulation is used to discretise the momentum equations (Eq. (3)). A finite volume discretisation of the continuity equations and a linear discontinuous Galerkin (DG) (Cotter et al., 2009a; Xiao et al., 2013) discretisation of the momentum equations are employed with backward Euler time stepping. Within each time step, the equations are iterated upon using a projection-based pressure determination method until all equations are simultaneously balanced. A non-linear Petrov-Galerkin method is used to stabilise the momentum equations. This involves the addition of a non-linear diffusion term over and above a standard BubnovGalerkin DG discretisation (see Hughes and Mallet, 1986; Pain et al., 2001; Fang et al., 2013, for further details). The main novelty of the framework includes: (a) the new finite element type (P1 DG-P2 ) for multi-fluid flow problems, which ensures exact balance between buoyancy force and pressure gradient; (b) a novel interface capturing scheme based on compressive control volume advection method (Pavlidis et al., 2014), involving a high-order accurate finite element method to obtain fluxes on the control volume boundaries, where these fluxes are subject to flux-limiting using a 6
normalised variable diagram approach (Leonard, 1991) to obtain bounded and compressive solutions for the interface; (c) a forcebalanced algorithm for the surface tension implementation in unstructured mesh methods, minimising the spurious velocities often found in interfacial flows (Francois et al., 2006); (d) anisotropic unstructured mesh adaptivity, which allows the grid resolution to be concentrated in relatively important regions such as the vicinity of interfaces, while lower resolution can be used in other regions. 3. Numerical examples Prior to attempting calculations of interfacial flows, several validation checks of our numerical framework were performed. In this section we demonstrate the convergence properties of the present method and also show the potential of using anisotropic unstructured mesh adaptivity for multiphase flow problems. All the numerical results presented here have been obtained using the P1 DG-P2 finite elements in two spatial dimensions. Here the mesh refinement criterion is based on minimising the interpolation error estimate for the component volume fraction. Between adaptive meshes, Galerkin interpolation is used for the velocity and consistent interpolation is used for the pressure and scalar fields. 3.1. Rayleigh-Taylor instability The Rayleigh-Taylor instability will occur if a heavy fluid initially lies on top of a lighter fluid in a gravitational field. It has been a benchmark test case in the computational fluid dynamics community for numerical methods to study multiphase flow problems. The initial growth and evolution of Rayleigh-Taylor instability have been studied for inviscid (Tryggvason, 1988) and viscous (Guermond and Quartapelle, 2000; Ding et al., 2007) incompressible flows. Here we use the same setup for a rectangular domain [0, d] × [0, 4d] with the density difference, represented by the Atwood ratio At = (ρA − ρB )/(ρA + ρB ), equals to 0.5 and with the Reynolds number ReA = ρA d3/2 g 1/2 /µA = 3000, where ρA and ρB denote the densities of the heavier and light fluids, respectively and µA is the dynamic viscosity of the heavier fluid. The interface between the two fluids is an initially sinusoidal perturbation of amplitude 0.1d. Following Tryggvason (1988), we use the same nondimensional variables, where the length is scaled by d, time is
7
1.5 Tryggvason (1988) Guermond and Quartapelle (2000) Ding et al. (2007) coarse mesh fine mesh
1
Position
0.5
0
−0.5
−1
−1.5
0
0.5
1.5
1
2
2.5
t
Figure 2: Position of the tip of the rising and falling fluid versus time for the RayleighTaylor instability with structured mesh simulations with At = 0.5 and ReA = 3000.
p √ scaled by d/Atg and the velocity is scaled by Atgd. In this case, there is no surface tension, and the ‘interface’ evolves due to gravity. First, computations were carried out using two fixed structured triangular meshes, namely a coarse mesh having 32 elements in the x direction and 128 elements in the y direction with 8192 (32×128×2) elements in total, and a fine mesh (64×256×2=32768 elements) which doubles the resolution of the coarse mesh in each direction. As expected, the heavier fluid falls into the light fluid at the centre while the light fluid rises up at both sides of the domain. Fig. 2 shows the computational results of the evolution of the positions of the tip of the rising and falling fluid, together with corresponding results obtained by the Lagrangian-Eulerian vortex method (Tryggvason, 1988), the variable density FEM method (Guermond and Quartapelle, 2000) and the diffuse interface method (Ding et al., 2007). Inspection of Fig. 2 reveals that convergence has been achieved upon mesh-refinement and that agreement with previous work is also excellent. Numerical simulations were then performed on adaptive, unstructured meshes. Here, the mesh refinement criteria are: the interpola-
8
tion error bound is 0.05, and the specified minimun and maximum length scale are 0.005d and 0.5d, respectively. The evolution of the interface is shown in Fig. 3, where the fixed fine mesh result is shown in the left column in panels (b)-(f) and the adaptive mesh result together with corresponding mesh is shown in the middle and right columns in panels (b)-(f). In order to facilitate comparison, adaptive simulations were started from the same initial profile and mesh for the fixed mesh case at t = 0 in Fig. 3(a). It can be seen that two counter-rotating vortices develop when the heavier fluid falls down. Inspection of Fig. 3 reveals the level of refinement of the unstructured mesh in the interfacial region. This refinement is adaptive and targets the regions of rapid variation in interfacial structure associated with the roll-up phenomena that are known to accompany the development of Rayleigh-Taylor instabilities driven by density contrasts. It can be seen from Fig. 3 that similar results are obtained for both fixed structured and adapative unstructured meshes, which confirms the overall consistency of the adaptive simulation when having varying mesh resolution across the whole flow domain. It is also worth noting that the total number of elements is kept constant as 32768 for the fixed mesh simulation. For the adaptive simulation, the number of elements reduces to 1137 at the first mesh adaption as fine meshes are only needed around the interface while coarse meshes can be placed far away from the interface. The total number of elements gradually increases to 5967 (see Fig. 3(f)) in order to capture the details of the vortex sheet. This demonstrates that in order to achieve the same order of accuracy fewer elements are required using the adaptive mesh, and this significantly reduces the computational efforts. It is also noted that during the adaptive simulation, although the Courant number is larger than unity in some regions with small meshes, the interface is well calculated due to the implicit time stepping for the interface-capturing method. Fig. 4 illustrates the detailed interface shape and mesh distribution at t = 2.0 together with its enlarged version from the left pannel to the right one. It is shown that when the heavier fluid sinks to the lower section of the domain its central region becomes elongated and possesses strong anisotropies. In accordance with this, anisotropic meshes are placed within the vicinity of the interface, thus the method has focused resolution where it is required to optimally represent the deveopment of the interface.
9
(a) t = 0
(b) t = 0.005
(c) t = 1.0
(d) t = 1.25
(e) t = 1.5
(f) t = 1.75
10 Figure 3: Numerical simulation of the Rayleigh-Taylor instability for the same parameters as those used to generate Fig. 2 showing a comparison between structured (left column in panels (b)-(f)) and adaptive unstructured meshing (middle and right columns in panels (b)-(f)). In (a), we show the condition used to initialise all simulations.
Figure 4: Detailed interface shape and mesh profile of the Rayleigh-Taylor instability at t = 2.0. Figures are gradually enlarged from left to right in order to show the anisotropic unstructured mesh around the interface.
3.2. Rising bubble Next we consider the rise of an air bubble in water, as accurate modelling of such flow is challenging due to its deforming interface as well as its surface tension force. Here we follow the same nondimensional setup from previous studies (Francois et al., 2006; Popinet, 2009), where an air bubble with diameter of D = 2/3 is initially placed at (0, 0) in a computational domain [−1, 1] × [−1, 2]. The fluid properties are: ρair = 1.126, ρwater = 1000, µair = 1.137, µwater = 1.78 × 10−5 . The suface tension coefficient σ = 728 2 = 5.983 for water. A force-balanced which gives a Bond number Bo = ρgD σ algorithm was implemented in the model for the surface tension calculation on unstructured meshes, where a diffused interface, reconstructed from the volume fraction on the unstructured meshes, provided good estimation of the interfacial curvature κ in Eq. (3). It is worth noting that these computations are 2D, where curvatures are different compared to 3D cases. Computations were carried out using a fixed structured triangular mesh, having 80 elements in the x direction and 120 elements in the y direction with 19200 (80 × 120 × 2) elements in total, and an adaptive mesh. The mesh 11
refinement criteria are: the interpolation error bound is 0.05, and the specified minimun and maximum length scale are 0.001 and 0.5, respectively. In Fig. 5, it can be seen that the initially spherical bubble rises under the action of buoyancy, and undergoes deformation, resulting in the formation of the well-known cap-shaped bubbles. The results associated with structured, and adaptive unstructured meshes are shown in the left, and middle and right columns of panels (b)-(e); the initial condition for all the simulations is shown in panel (a). It is shown from Fig. 5 that both results closely match with each other until t = 0.35, and they differ slightly after that mainly because the mesh resolution for the fixed mesh simulation is not fine enough to capture the shape of interface when the bubble has large deformation. In contrast, for the adaptive mesh simulation, inspection of Fig. 5 reveals that the unstructured mesh is not dense around the bubble at early times since the latter had not yet been deformed by the flow. As the extent of deformation becomes increasingly pronounced, our numerical method responds by refining the unstructured mesh in the vicinity of the deforming bubble preferentially; an enlarged version of Fig. 5(e) is shown in Fig. 6, which depicts the details of the unstructured mesh around the bubble. In contrast, the regions upstream and downstream of the bubble, which require a lower degree of resolution, have coarser meshes, in order to maximise computational efficiency. It can also be seen from Fig. 6 that the sharp interface is well captured during the adaptive simulation, due to the use of the compressive control volume advection method developed in the framework. In order to demonstrate the efficiency of the adaptive unstructured mesh simulation, Fig. 7 shows the time history of the total number of elements during the simulations. It can be seen that the total number decreases significantly at the first mesh adaption, approximately only 10% of the fixed mesh case. The maximum number of elements used in the adaptive simulation is only 18% of the fixed mesh case, which reduces the computational efforts without sacrificing accuracy. 4. Conclusions In this paper an adaptive unstructured mesh modelling approach, which can reduce computational efforts without sacrifing accuracy, has been presented for interface-capturing in multiphase flow problems with an emphasis on two-phase flows. The interface-capturing method is based on a compres12
(a) t = 0
(b) t = 0.001
(c) t = 0.2
(d) t = 0.35
13 (e) t = 0.5 Figure 5: Numerical simulation of a rising bubble with a Bond number Bowater = 5.983. Here, we show a comparison of the results associated with structured, and unstructured adaptive meshes shown in the left, and middle and right columns of panels (b)-(e), respectively. In (a), we show the condition that was used to initialise all simulations.
Figure 6: An enlarged version of Fig. 5(e) showing the details of the numerical solution obtained using an adaptive, unstructured mesh.
4
number of elements
2
x 10
fixed mesh adaptive mesh
1.5
1
0.5
0
0
0.1
0.3
0.2
0.4
0.5
t
Figure 7: Numbers of total elements for the numerical simulations of a rising bubble with a fixed and adaptive, unstructured mesh. The parameters are the same as those used to generate Fig. 5 and Fig. 6.
14
sive advection method, which uses a novel and mathematically rigorous nonlinear Petrov-Galerkin method that attempts to keep the interface sharp. This was shown to perform well in a standard Rayleigh-Taylor instability benchmark problem as well as bubble rising two-phase flow problems. We also presented results using the new finite element type (P1 DG-P2 ) for multiphase flow problems. The element type ensures that linear buoyancy force and vertical pressure gradient are exactly balanced which eliminates large spurious velocities that might be generated on unstructured meshes if this balance is inaccurate. We also showed the method working with anisotropic mesh adaptivity, which allows the grid resolution to be concentrated in relatively important regions such as the vicinity of interfaces while lower resolution can be used in other regions. In addition, the interface-capturing method is implicit which means that it is stable with relatively large time step sizes exceeding unity Courant numbers. The results presented here established with sufficient confidence that this method can be used to successfully model multiphase flows in a wide range of applications. This approach has the potential to be used for an arbitrary number of components although that has not been demonstrated here. Altough only 2D simulations are considered here, extention to 3D is straightforward. Currently we are working on parallel computing and some fully 3D simulations will be reported in the near future. 5. Acknowledgements We would like to thank the EPSRC MEMPHIS multiphase Programme Grant (number EP/K003976/1), the EPSRC Computational modelling for advanced nuclear power plants project and the EU FP7 projects THINS and GoFastR for helping to fund this work. References Anderson, D.M., McFadden, G.B., Wheeler, A.A., 1998. Diffuse-interface methods in fluid mechanics. Annual Review of Fluid Mechanics 30, 139– 165. Badalassi, V.E., Ceniceros, H.D., Banerjee, S., 2003. Computation of multiphase systems with phase field models. Journal of Computational Physics 190, 371–397. 15
Ceniceros, H.D., Roma, A.M., Silveira-Neto, A., Villar, M.M., 2010. A robust, fully adaptive hybrid level-set/front-tracking method for two-phase flows with an accurate surface tension computation. Communications in Computational Physics 8, 51–94. Chen, Y.G., Price, W.G., Temarel, P., 2001. An anti-diffusive volume of fluid method for interfacial fluid flows. International Journal for Numerical Methods in Fluids 68, 341–359. Cotter, C.J., Ham, D.A., Pain, C.C., 2009a. A mixed discontinuous/continuous finite element pair for shallow-water ocean modelling. Ocean Modelling 26, 86–90. Cotter, C.J., Ham, D.A., Pain, C.C., Reich, S., 2009b. LBB stability of a mixed Galerkin finite element pair for fluid flow simulations. Journal of Computational Physics 228, 336–348. Ding, H., Spelt, P.D.M., Shu, C., 2007. Diffuse interface model for incompressible two-phase flows with large density ratios. Journal of Computational Physics 226, 2078–2095. Fang, F., Pain, C.C., Navon, I.M., ElSheikh, A.H., Du, J., Xiao, D., 2013. Non-linear Petrov-Galerkin methods for reduced order hyperbolic equations and discontinuous finite element methods. Journal of Computational Physics 234, 540–559. Francois, M.M., Cummins, S.J., Dendy, E.D., Kothe, D.B., Sicilian, J.M., Williams, M.W., 2006. A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework. Journal of Computational Physics 213, 141–173. Guermond, J.L., Quartapelle, L., 2000. A projection FEM for variable density incompressible flows. Journal of Computational Physics 165, 167–188. Harlow, F.H., Welch, J.E., 1965. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Physics and Fluids 8, 2182–2189. Hirt, C.W., Nichols, B.D., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics 39, 201– 225. 16
Hughes, T.J.R., Mallet, M., 1986. A new finite element formulation for computational fluid dynamics - IV. A discontinuity-capturing operator for multi-dimensional advection-diffusion systems. Computer Methhods in Applied Mechanics and Engineering 58, 329–336. Jacqmin, D., 1999. Calculation of two-phase Navier-Stokes flows using phasefield modeling. Journal of Computational Physics 155, 96–127. Kim, J.S., Kang, K., Lowengrub, J., 2004. Conservative multigrid methods for Cahn-Hilliard fluids. Journal of Computational Physics 193, 511–543. Leonard, B.P., 1991. The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection 88, 17–74. Monaghan, J.J., 1992. Smoothed particle hydrodynamics. Annual Review of Astronomy and Astrophysics 30, 543–574. Osher, S.J., Sethian, J.A., 1988. Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics 79, 12–49. Pain, C.C., de Oliveira, C.R.E., Goddard, A.J.H., Umpleby, A.P., 2001. Criticality behaviour of dilute plutonium solutions. Nuclear Science and Technology 30, 194–214. Pavlidis, D., Xie, Z., Percival, J.R., Gomes, J.L.M.A., Pain, C.C., Matar, O.K., 2014. Two- and three-phase horizontal slug flow modelling using an interface-capturing compositional approach. International Journal of Multiphase Flow. Accepted. Piggott, M.D., Farrell, P.E., Wilson, C.R., Gorman, G.J., Pain, C.C., 2009. Anisotropic mesh adaptivity for multi-scale ocean modelling. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 367, 4591–4611. Popinet, S., 2009. An accurate adaptive solver for surface-tension-driven interfacial flows. Journal of Computational Physics 228, 5838–5866. Rider, W.J., Kothe, B., 1998. Reconstructing volume tracking. Journal of Computational Physics 141, 112–152.
17
Rudman, M., 1997. Volume-tracking methods for interfacial flow calculations. International Journal for Numerical Methods in Fluids 24, 671–691. Scardovelli, R., Zaleski, S., 1999. Direct numerical simulation of free-surface and interfacial flow. Annual Review of Fluid Mechanics 31, 567–603. Sethian, J.A., Smereka, P., 2003. Level set methods for fluid interfaces. Annual Review of Fluid Mechanics 35, 341–372. So, K.K., Hu, X.Y., Adams, N.A., 2011. Anti-diffusion method for interface steepening in two-phase incompressible flow. Journal of Computational Physics 230, 5155–5177. Tryggvason, G., 1988. Numerical simulations of the Rayleigh-Taylor instability. Journal of Computational Physics 75, 253–282. Ubbink, O., Issa, R.I., 1999. A method for capturing sharp fluid interfaces on arbitrary meshes. Journal of Computational Physics 153, 26–50. Unverdi, S.O., Tryggvason, G., 1992. A front-tracking method for viscous, incompressible, multi-fluid flows. Journal of Computational Physics 100, 25–37. Wilson, C., 2009. Modelling multiple-material flows on adaptive unstructured meshes. Ph.D. thesis. Department of Earth Science and Engineering, Imperial College London. Xiao, D., Fang, F., Du, J., Pain, C.C., Navon, I.M., Buchan, A.G., ElSheikh, A.H., Hu, G., 2013. Non-linear Petrov-Galerkin methods for reduced order modelling of the Navier-Stokes equations using a mixed finite element pair. Computer Methods in Applied Mechanics and Engineering 255, 147–157. Zheng, X., Lowengrub, J., Anderson, A., 2005. Adaptive unstructured volume remeshing ii: Application to two- and three-dimensional level-set simulations of multiphase flow. Journal of Computational Physics 208, 626–650.
18