Georgina Reece, Benedict D. Rogers, Steven Lind, Georgios Fourtakas

Department of Mechanical, Aerospace and Civil Engineering, University of Manchester, Manchester, UK

Abstract: This paper assesses the ability of smoothed particle hydrodynamics (SPH) to simulate mixing of two-phase flows and their transition to instabilities under different flow regimes. A new measure for quantification of the degree of mixing between phases in a Lagrangian framework is also developed. The method is validated using the lid-driven cavity and two-phase Poiseuille flow cases. The velocity along the centre of the cavity is compared with results from the literature, whilst commercial volume-of-fluid code STAR-CCM+ provides a benchmark for the mixing and different mixing measures are considered. The velocity of two-phase Poiseuille flow along the channel is compared to the analytical solution, and the appearance of interfacial instabilities with perturbation theory.This is the first time SPH has been used to investigate the onset and development of these instabilities. In particular, it is able to model the deforming shape of the interface, which is not given by analytical studies, while also offering improved predictions over conventional mesh-based computational fluid dynamics simulations.

Key words: Smoothed particle hydrodynamics (SPH), multi-phase, mixing, instability

Introduction

Numerous industrial applications involve mixing of two fluid phases. In almost all cases, it is important to know how well-mixed the product is. Yet taking a measurement at a point within the mixture can range from difficult to dangerous, as well as it being expensive to run repeated experiments until the desired confidence in the result is found. Therefore, a computational model where different experimental configurations can be tested to find out which results in the most homogeneous mixture is valuable. A computational model also gives a view of the whole domain at any time, giving an increased understanding of the physical processes which produce the resulting mixture.

Mesh-based methods are regularly used to model mixing, such as the viscous multi-phase mixing models of Chang et al.[1]and Kunz et al.[2]. However,all mesh-based multi-phase methods are subject to the extra computational cost of tracking the interface.

The meshless method smoothed particle hydrodynamics (SPH) offers a different approach with distinct advantages since its Lagrangian mesh-free nature allows it to naturally follow a complex,deforming interface. It can also deal with interpenetration of phases without any additional expense while also conserving important physical properties, such as mass. SPH is a Lagrangian discretisation scheme,originally developed by Gingold and Monaghan[3]and Lucy[4]for astrophysical modelling. It has since been applied to many fields, including free-surface flows[5],multi-phase flows[6]and fluid-structure interactions[7].The fluid is discretised as a set of points, each representing a volume of the fluid with fixed mass,which are free to move and carry material properties.

Simulating mixing is a natural application of multi-phase SPH models which were first introduced by Monaghan in 1997[8], before being developed further in a number of works including[9-11]. Mixing of two phases using SPH has also been investigated many times[12-18], whilst other papers concentrate on viscous mixing[19-21].

To report mixing, a number of aspects may be considered, depending on the application and desired result. The distribution of each phase, for example using a volume fraction, gives information about how well-mixed the fluid is at any point and can be quantified globally using a norm over the domain.Robinson[22]introduced a method to compare the local proportion of phases with the global proportion, which is especially useful when there are more than two phases. However, to understand why some regions are more thoroughly mixed than others, there must be further insight into the physical processes involved.The finite time Lyapunov exponent (FTLE) has been used a number of times in conjunction with SPH[22-25]as a means of identifying Lagrangian coherent structures, and therefore regions within the flow where mixing is promoted or inhibited. However,Lagrangian methods such as SPH offer unique insights that are further explored herein.

In this work, the open-source (weakly compressible) SPH solver DualSPHysics[26]is modified so that two phases of different viscosity may be modelled. The DualSPHysics code is rigorously validated and widely used in both academia and industry, providing a highly optimised code that runs on different hardware including central processing units (CPUs) or graphics processing units (GPUs)[26].This modified DualSPHysics formulation is validated using results from the literature[27]for the lid driven cavity case, as well as the analytical velocity of the two-phase Poiseuille flow case, which exhibits instabilities at the interface under certain conditions.

Prior to mixing, under specific conditions, the interface between two phases is known to become unstable. Numerous papers in the literature consider the interface instability which appears in two-phase Poiseuille flow. Herein, the particular case of two phases of different viscosity but constant density is the focus. The first person to characterise this instability was Yih in 1967[28], who derived an Orr-Sommerfeld equation which governs stability. This leads to an eigenvalue problem with respect to the wave speed.Assuming constant density and equal widths of phases,a relationship was found in terms of viscosity ratio( m) which governs when the flow becomes unstable.

A relationship between viscosity ratio, m , and phase width ratio, n, for instability development from long-wavelength disturbances was also found by Yiantsios and Higgins[29]using perturbation analysis.Gravity is included but is insignificant for the case of constant density. Conditions were found for interfacial wave growth in terms of the relationship between phase width ratio and viscosity ratio.

Experiments have also been performed to investigate the instability. Charles and Lilleleht[30]used oil and water in a rectangular channel to find that interfacial waves occurred at the turbulent transition in the less viscous (water) phase. These waves increased when the flow rate of water and ratio of water to oil were both higher. Kao and Park[31]used a similar set up to find a critical Reynolds number, below which instability growth at the interface is fully damped.

The challenge when investigating instabilities is to find the resolution required to accurately reproduce their onset. In addition, the interpolation and discretisation errors introduced through the SPH formulation and investigated by Quinlan et al.[32]must be quantified and controlled in order to be confident of the accuracy of instability prediction.

The appearance of these instabilities is often compared with analytical results. However, since almost all literature on this case is analytical, there is little to compare the development of the instability once initiated. Most of the analytical investigations(including Pinarbasi and Liakopoulos[33], Valluri et al.[34], Hooper[35-36]and Anturkar et al.[37]) use perturbation theory to find a relationship between wavenumber and formation of an interfacial instability.In a computational setting, identifying the formation and growth of an interfacial instability is less obvious,and here a method to measure the numerical growth rate of the instability from the simulation is also investigated.

It is demonstrated that SPH is an ideal method to investigate and understand both the growth rate of the instability and development of the interface shape.Mixing measures to quantify the degree of mixing taking place are compared with a new measure (that combines volume fraction with the FTLE) for both lid driven cavity mixing and an unstable case of two-phase Poiseuille flow.

1. Numerical methodology

The open-source SPH solver DualSPHysics[26]is used as a base for modelling these cases. Therefore,the formulation follows that detailed by Crespo et al.[26], with changes described later in this section.

Governing equations-Navier-Stokes equations:The Navier-Stokes equations for fluid flow are given by

Continuity equation

Momentum equation

Here D/Dt is the material derivative with respect to time, ρ is the density, I is the identity matrix, v is the velocity, p is the pressure, τ is the stress tensor and F is the body force acting on the fluid.

Assuming the fluid is weakly compressible, we use an equation of state to calculate pressure from density. This is described later in this section.

SPH method: SPH particles interact via an interpolation over neighbouring particles, weighted by a kernel function W . For an arbitrary fluid property Φ, its value at a point with neighbourhood Ω is given by

where r is the particle position, h is the smoothing length of the kernel and - denotes the approximate integral interpolation.

In the SPH method, this integral is discretised as a sum over all particles within the kernel radius,shown in Fig. 1. For a particle i, the value of property Φ is given by the following sum over neighbouring particles j SPH formulation: The weakly compressible SPH formulation is used, with governing equations

where f is the body force, ν is the kinematic viscosity. The subscripts refer to individual SPH particles and Wijis the kernel value of particles i and j. To calculate pressure, the Tait equation of state[38]

The Wendland kernel[40], a fifth-order approximation of the Gaussian kernel, is used due to its high accuracy[41]. The 2-D kernel with smoothing length h=2 is shown in Fig. 1. The standard boundary condition in DualSPHysics is the dynamic boundary condition (DBC) model[42]and is ideal for acceleration offered by a GPU. However, it can sometimes suffer from excessive particle separation from the boundary that can only be mitigated by increasing the resolution.In Section 2, the modified dynamic boundary conditions (MDBC) presented by English et al.[43]are used to obtain the correct velocity and pressure profiles for this case, which is relevant for the mixing investigated therein. The Fickian based shifting[44]and density diffusion[45]formulations built into DualSPHysics are employed to mitigate against the numerical errors introduced by the SPH discretisation.

The viscosity term in the conservation of momentum equation (Eq. (5)) is usually modelled with the operator of Lo and Shao[46]

whereir is the position of a particle i, rijis the displacement between particles i and j and ε is a small constant introduced to avoid singularities.

In order to model cases with two or more phases of different viscosity, a new variable is associated to each particle in the form of a viscosity array within DualSPHysics. Each particle is initialised with a viscosity which remains constant throughout the simulation. The viscous operator must also be adapted to account for two interacting particles with different viscosities. This is done by replacing the constant viscosity in Eq. (7) with the average of the two viscosities of the particles involved. Thus, the original operator is recovered when the interacting particles have the same viscosity and the operator becomes

Fig. 1(a) (Color online) Kernel around particle i of radius 2h,showing interacting particles in blue

Fig. 1(b) (Color online) Wendland kernel with smoothing length h=1, showing the weighting on particle contributions at different distances from particle i at r=0

In this work, it is assumed that density （ρ） is the same in each phase and viscosities are constant with time.

Mixing measures: In order to compare the quality of mixing cases, the degree of mixing at any point within the domain must be quantified. Four measures are considered here, including a new measure, with their features and relative strengths for this application assessed in the Numerical Results sections.

(1) The volume fraction gives the concentration of one of the phases at a point in the domain. When there are two phases, one minus the concentration of the first phase naturally gives the concentration of the second phase. For phase 1, its volume fraction at point i is calculated as the quotient of two SPH summations

where Vjis the volume of particle j, J is the set of all particles j within a kernel radius of point i and J～?J is the subset of these particles which are also of phase 1. For comparison with the benchmark,an L2error norm is calculated.

(2) The Lyapunov exponent (λ) is useful to characterise chaotic flow or mixing, since if all λ<0 then all solutions are asymptotically stable. Lyapunov exponents can be used to locate manifolds in the flow,which show regions where mixing is promoted or inhibited (manifolds act as separatrices of the flow).From these, the quality of mixing can be inferred.However, the definition takes a limit as time tends to infinity and is therefore not practical to apply to simulation data.

The maximum is taken over all particles in J, which are within a kernel radius of i at time t=0. Sun et al.[47]extended the use of FTLEs in SPH to 3-D.Grid-based methods use tracer particles to calculate FTLE.

(3) Robinson’s[22]local measure of mixing was introduced in order to quantify the mixing between phases at a given point. For this work, the number of particles is replaced with the local volume fraction,with the added assumption that there are two phases whose total volume fraction at any point (and globally)is equal to one. The measure at a point i can then be calculated using only the local and global volume fraction of phase 1 as

where F1is the global volume fraction of phase 1.

This mixing measure will be used to quantify local mixing between phases in the two-phase lid driven cavity validation case. It will also be compared to the new measure developed to quantify the degree of mixing between phases.

(4) New mixing measure: Although each of the previous measures gives a valuable insight into the mixing achieved for a given case, they do not individually give all the information required to assess the quality of mixing. A new measure, taking into account both the degree and composition of mixing throughout the domain, is introduced with the aim of giving the full picture. It is defined as a combination of the FTLE and the volume fraction

where the bars denote that the FTLE and volume fraction are both normalised with respect to the maximum FTLE at that time and the global volume fraction of phase 1 respectively. Firstly, both components are normalised so they are 1 if most mixed, 0 if least, such that

where σmaxis the maximum σ for T and F is the global volume fraction. Therefore, the measure has values from 0 to 1, with high values achieved only when both phases are present and also particles have diverged from those nearby whom they were initially close. This measure, defined in Eq. (13), holds when there are two phases whose total volume fraction sums to one. For more phases, the definition will need to be adapted.

2. Numerical results: Lid driven cavity

A square cavity of width L=1m, surrounded by boundary walls, is filled with a fluid which is driven by moving the top wall at horizontal velocity Vlid=1m/s . A diagram of the case can be seen in Fig.2. The flow behaviour is characterised by the Reynolds number

Fig. 2 Set-up of lid driven cavity case

The shifting configuration used in the following results is not applied to particles within 2h of the boundaries and does cause some clumping at the right hand wall. However, particles were lost through the top right corner of the cavity when shifting is applied to all particles. The Lind et al.[48]shifting formulation with tensile correction reduces the number of particles lost but does not prevent it completely.

Velocity profile–single phase: The case when the fluid within the cavity is of a single phase has been comprehensively studied in the literature. Ghia et al.[27]provide a large set of data in their paper, which is recommended by the SPHERIC organisation as a benchmark[49]. The case is determined to have reached a steady state according to

where E=0.5∑i mivi·viis the total kinetic energy of the system and the dot denotes the derivative with respect to time. The constant ε must not be chosen too small since the energy measure displays fluctuations due to its discrete nature, yet small enough to be a realistic steady state. This was added to the main code to be output at each time step. Figure 3 shows how both E and E˙ change with time for a given case, as well as highlighting the point of steady state. An approximately steady state (determined by setting ε=0.1) is reached for this case at 11.25 s.

Fig. 3 (Color online) Energy (a) and change in energy (b) over time to determine steady state of lid driven cavity case.The point at which steady state is reached is marked with a black circle. Re=100, =0.01 dp

Comparisons with this data for Re=100, at steady state are shown in Fig. 4. It can be seen that SPH produces a velocity close to the literature data.The velocity profile given by STAR-CCM+ also comes close to the benchmark results.

It can be seen in Fig. 5 that increasing the resolution generates a velocity curve closer to the data of Ghia et al.[27]. The2L error norm of the velocity in Fig. 6 shows convergence with increasing resolution.These results confirm the ability of this SPH formulation to reproduce lid driven cavity flow, which shall be further used for mixing of two phases in the cavity.

Fig. 4 (Color online) Horizontal component of velocity along a vertical line through the centre of the cavity. Re=100,dp=0.02

Fig. 5 (Color online) Convergence of velocity profile towards Ghia et al. data with reducing particle spacing

Fig. 6 (Color online) L2 error norm of centre line velocity with 1st and 2nd order convergence lines

Streamlines–single phase: Streamlines of the same Re=100 case are plotted in Fig. 7 for SPH.The main features of this flow are captured by the SPH model. The main large vortex has formed at the top of the cavity, slightly to the right, and the formation of the smaller anticlockwise vortices is captured at the bottom corners. Therefore, the model is capable of capturing flow behaviour near the boundaries. An increase in particle number could also improve the resolution of these smaller vortices.

Fig. 7 Streamlines of lid driven cavity flow at steady state for Re=100

Two-phase lid-driven cavity flow: In the absence of an analytical solution, to validate the 2-phase mixing flow, the results from the SPH model are compared with a reference numerical solution from the commercial code STAR-CCM+. For the two-phase case, the cavity is split horizontally in half into two different fluid phases, as shown in Fig. 8.Each fluid phase has a different constant viscosity, but the same density. Here the upper phase has viscosity ν1=0.01m2/s (Re=100), and the bottom phase ν2=0.02m2/s (Re=50).

Fig. 8 (Color online) Set-up of lid driven cavity case with two phases

Mixing measures: A natural comparison with a VOF method is to plot the volume fraction, calculated using Eq. (10). The results for the volume fraction are shown in Fig. 9, coloured by the volume fraction of phase 1,1f. Since there are two phases in this case,the volume fraction of phase 2 is21 f=1-f at any point. The third column shows the absolute difference in volume fraction between the two methods, where 0 means they are identical at that point. The two approaches produce many qualitative similarities,such as the clear line of phase 1 towards the bottom of the mixing and the shape of unmixed phase 2 after 60 s. Due to the volume fraction approach that does not produce a sharp interface between phases, it is not possible to get complete convergence of the STARCCM+ results at the interface, even with a high resolution. The break-up and diffusion of the interface can be seen here around (0.9, 0.8) at 5 s.

The volume fraction of the SPH method is compared to STAR-CCM+ quantitatively using the L2error norm in Fig. 10. The difference initially increases as movement begins, then decreases more gradually after 30 s. This shows the methods tending towards convergence as time increase to 60 s.

Figure 11 shows the SPH mixing measure results at three different times for the finite time Lyapunov exponent (FTLE), Robinson’s local mixing measure and the new mixing measure. The left-hand column in Fig. 11 shows the FTLE where the scale is different for each time in order to make the features visible.The Lagrangian approximation for the forward-in-timeFTLE in Eq. (10) is used in post-processing. A greater value of FTLE indicates that SPH particles starting near this point at t=0 have moved furthest away from each other at the presented time. This is useful for understanding the movement of particles but does not show the mixing between phases relevant for this work.

Fig. 9 (Color online) Volume fraction for STAR-CCM+ and SPH, with difference between methods in right column. Rows at 5 s, 20 s, 60 s Upper phase ν1=0.02 m2/s , lower phase ν2=0.01m2/s

Fig. 10 2L1 norm of volume fraction with time, =0.01mdp

Robinson’s mixing measure[22](Eq. (11)) in the centre of Fig. 11 was designed to quantify mixing locally at a point in the domain. The areas of mixing,as well as the regions of unmixed phases, are clear to see. The lid driven cavity clearly does not produce a homogeneous mixture, at least with phase viscosities as here. This demonstrates that the lid-driven cavity is a poor mixing case, but shows the value of a mixing measure in determining the quality of mixing.

Fig. 11 (Color online) Mixing measures at times 5 s, 20 s, 60 s. Upper phase ν1=0.02m2/s , lower phase ν2=0.01m2/s

The right hand column of Fig. 11 shows the new combined VF-FTLE measure (Eq. (13)), which illustrates both the mixing of phases and the movement of SPH particles. Here the volume fraction component must be calculated at SPH particle locations at t=0 so that these values are at the same points as the FTLE values. This new measure appears similar in distribution to Robinson's mixing measure,but appears to be smoother. The definition has a clear effect towards the bottom right of the cavity at 60 s,where f is higher but FTLE is small. Therefore,both phases are present but there is little divergence between particles which begin within each-other’s neighbourhood and the resulting mixing measure is low.

Whereas the volume fraction gives information about the distribution of phases at a given time and the FTLE shows how much particles starting near a point have diverged from each other, both the Robinson and VF-FTLE measures display the amount of movement between particles of different phases. In particular, the new VF-FTLE measure can exhibit regions where either both phases are present but particles have remained close to those they were near at t=0, or particles have moved considerably relative to each other but there has been no mixing between phases.For example, around (x, y)=(0.9,0.2) at 20 s.

3. Numerical results: Two-phase Poiseuille flow

Velocity profile: Unlike the single-phase case, the introduction of two phases of different viscosities means the velocity profile is not necessarily symmetric about the centre of the channel and may not even be smooth at the interface between phases. An analytical solution for the velocity profile of two-phase Poiseuille flow is given by Bird et al.[50], but assumes that the account for cases wherethephases are not ofequal interface is at the centre of th e channel (d 1=d 2). To width, the analytical solution derived by Rowlatt and Lind[51]for three-phase Poiseuille flow is modified for only two phases. The analytical solution for velocity in each phase v～ is then

with constants

The fluid phases start from stationary and are put into motion by an acceleration along the channel of

2(0.8,0)m/s.

Herein, five cases are presented for different ratios of the viscosities m =ν2/ν1and width n =d2/d1as given in Table 1.

Figure 13 shows that the method detailed in Section 1 reproduces the shape of the velocity profile well for the unstable Case 0. This simulation gives

Fig. 12 (Color online) Two-phase Poiseuille flow set-up, showing kinematic viscosities 1ν and 2ν and width of each fluid phase. Motion is from left to right

Table 1 Test points for two-phase Poiseuille flow instability

close results to the analytical solution, capturing the qualitative behaviour of the velocity well–the shapes of the velocity profile in each phase and the discontinuity at the interface. The maximum velocity and velocity at the interface are both less than 3% higher than the analytical value. The differences are explained by the discrete approximation of the numerical scheme and error from the interpolation kernel, as well as the assumption by Bird et al.[50]of incompressibility. The MDBC used for the lid driven cavity case are found to have negligible difference here.

Fig. 13 (Color online) Velocity profile along channel for case 0,compared with analytical solution

For case v, ν1=0.10m2/s , ν2=0.20m2/s ,d1=0.3333m, d2=1.6667m , which is stable (see instability subsection), a particle convergence study with particle spacing, dp, has been undertaken to demonstrate that this method converges to the analytical solution as resolution is increased.

Equation (17) is used to determine when the simulation has reached a steady state. The points at which steady state has been reached are circled at approximately 13 s in Fig. 14, which shows how Eq.(17) changes with time for each particle spacing. Note that the difference between all curves remains bounded.The velocity profile for each resolution at steady state are shown in Fig. 15 and can be seen to get closer to the analytical solution as dp is decreased.

To quantify the convergence, the2L norm of the velocity along the channel in comparison to the analytical solution, v～x, is taken over all N SPH particles i

Fig. 14 (Color online) Steady state measure with time for various particle spacings, for case v

Fig. 15 (Color online) Velocity profiles of convergence study cases for various particle spacings, for case v

The resulting norms are plotted against particle spacing on logarithmic axes in Fig. 16. A line through these points gives the rate of convergence to be of O(1). This figure also confirms the results from Fig.15 that the simulations are converging.

Fig. 16 (Color online) 2L norm of velocity for each resolution,showing first order rate of convergence

Fig. 17 (Color online) Regions of stability for two-phase Poiseuille flow, with points chosen to test agreement of stability with perturbation theory marked. Green regions are stable and white unstable. Case v not plotted for clarity

Of the cases already considered, case 0 in Fig. 13 lies in an unstable region and case v in Fig. 15 is stable.Further points are chosen in each of the regions, shown in Fig. 17.

Each of the test points marked in Fig. 17 agrees with the stability of the region in which it is located.For example, Case iii producestheparticle distribution shown in Fig. 18 after the interfacial instability has started to grow. An investigation revealed that the use of the particle shifting technique does not prevent the onset of instability and only influences the shape of the growing instability at later time, when perturbation theory is no longer valid.

Fig. 18 (Color online) Distribution of phases of case iii for SPH after 21 s for n=2,=0.5m

Growth rate of instability: The growth rate of the instability is estimated by following the position across the channel of particles which are initially located at the interface. A case investigated by Valluri et al.[34]with theoretical growth rate 0.34 is reproduced using the SPH model described in Section 1. A maximum vertical position (rz) is taken over M particles,equally spaced along the line of the initial interface, at all output times. The onset of instability is defined to be at the first time max rzexceeds dp above the interface. To focus on linear growth only, a period of 0.2 s after instability onset is plotted in Fig. 19 and an exponential line fitted.

The results for different numbers of M points are given in Table 2. The points are regularly spaced throughout the domain in order to capture all amplitudes. It is found that 2 points are required per wavelength (≈1 mhere) in order to accurately capture the growth rate. This gives a growth rate of 0.35 for 2 or 4 particles, demonstrating the ability of this SPH method to model the early stages of instability growth.

Fig. 19 (Color online) Growth rate of Poiseuille flow instability for case defined by Valluri et al.[34]

Table 2 Growth rate of interfacial instability for different numbers of sample points

The commercial code STAR-CCM+ gives results in agreement for velocity along the channel, as shown in Fig. 21. However, it is unable to produce the deforming shape of the interface in Fig. 20 as SPH does in Fig. 20. For any instability that does develop at the interface using STAR-CCM+, growth is in the form of mixing of the two phases and as the instability grows it becomes more difficult to capture the position of the interface. This is due to the diffusive nature of the VOF method used for multi-phase models. Hence, it is also not possible to capture the growth rate and evolution of the interfacial instability accurately.

Fig. 20 (Color online) Two-phase Poiseuille flow for case iii simulated using STAR-CCM+ (same case as Fig. 18)

Since VOF is a mesh-based method, the exact position and shape of the interface must be approximated using reconstruction techniques which can be computationally expensive[52]and the accuracy depends on the resolution of the mesh. This is particularly difficult for a complex dynamic interface[53]and may require VOF to be coupled with another method[54]. Since exact information for the interface is lost and the new volume fluxes for time integration depend on the inferred interface, this error grows with time[54].

Fig. 21 (Color online) Velocity profiles for different methods,point (i)

Mixing measures performance: In order to compare results with other methods, it is valuable to be able to quantify the degree of mixing achieved throughout the domain. The mixing measures described in Section 1 are applied to the two-phase Poiseuille flow case. All figures shown here are for the same conditions as Fig. 18 after 45 s.

For periodic simulations such as two-phase Poiseuille flow, special consideration must be taken into account for particles near the periodic boundaries since they are both close together and far apart. For the volume fraction, particles within 2h of the periodic boundaries must consider particles at the opposite boundary to maintain a full kernel support. When calculating the FTLE for a case with periodic boundaries such as this, fluid particles lying within copies of the domain at either boundary are also considered and, for any particle within the original domain, the distance to another particle is taken to be the minimum between all pairs of copies.

The volume fraction plotted in Fig. 22(a) and defined by Eq. (9) clearly shows where the two phases are separate and where there is mixing between the two.However, it is not clear how much the particles within a phase have moved relative to each other and therefore there is no insight into the physical processes driving the mixing.

Conversely, the FTLE defined in Eq. (10) shows particle divergence relative to the initial configuration.Higher values indicate that particles that started close together near this point are relatively far from each other at the time of measurement. However, in Fig.22(b) it can be seen that there is no way of knowing the distribution of each phase. Therefore, the degree of mixing achieved cannot be gathered.

Fig. 22 (Color online) Mixing measures for two-phase Poiseuille flow (n=2,m=0.5 at 45 s)

Robinson designed the measure plotted in Fig.22(c) to quantify the degree of mixing locally at any point in the domain. This is done by comparing the local volume fraction to the global volume fraction for both phases present. The areas of mixing and regions where no mixing has taken place are clear to see in the figure. However, this again does not capture the relative particle movement.

The new measure introduced in Eq. (13) aims to possess the benefits of both volume fraction and FTLE by combining the measures. In this way, the VF-FTLE measure shows both the distribution of phases and the divergence from the initial conditions. The results in Fig. 22(d) are generally in agreement with Robinson's measure, but appear to be smoother with regions of high value more clearly defined.

4. Conclusions

This paper investigated two-phase mixing,interfacial instability onset and growth using the open-source code DualSPHysics. The results for the two-phase lid-driven cavity case show qualitative similarities between SPH and the commercial VOF code STAR-CCM+ but differences can be seen in the mixing measures. The results from the STAR-CCM+VOF model show capability for modelling mixing, but the diffusive nature means a deforming interface, and subsequent interfacial instabilities, cannot be tracked accurately.

Modifying the weakly compressible SPH scheme of DualSPHysics to include two phases with different viscosities has been shown to accurately predict the velocity in a case with two phases of different viscosity. Also, it is able to replicate the onset of interfacial instabilities under conditions in agreement with analytical results, with the growth rate in agreement with theory.

Several mixing measures for are compared and their value for these applications shown. In particular,a new measure which can show both the degree of mixing at a particular time and the relative movement of phases is presented. This new mixing measure has been demonstrated to provide useful insight into the mixing of particles of different phases at any given time.

Acknowledgement

The work was supported by the EPSRC and National Nuclear Laboratory (Grant No. 1961431).