{'fraction_non_alphanumeric': 0.026257414283647588, 'fraction_numerical': 0.011187732073106043, 'mean_word_length': 4.767626199415937, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A five-equation model for the simulation of miscible and viscous compressible fluids

Typical multispecies compressible Navier-Stokes computations employ conservative equations for mass fraction transport. Upwind discretisations of these governing equations produce spurious pressure oscillations at diffuse contact surfaces between gases of differing ratio of specific heat capacities which degrade the convergence rate of the algorithm. Adding quasi-conservative equations for volume fraction can solve this error, however this approach has been derived only for immiscible fluids. Here, a five-equation quasi-conservative model is proposed that includes the effects of species diffusion, viscosity and thermal conductivity. The derivation of the model is presented, along with a numerical method to solve the governing equations at second order accuracy in space and time. Formal convergence studies demonstrate the expected order of accuracy is achieved for three benchmark problems, cross-validated against two standard mass fraction models. In these test cases, the new model has between 2 and 10 times lower error for a given grid size. Simulations of a two-dimensional air-SF_6 Richtmyer-Meshkov instability demonstrate that the new model converges to the solution with four times fewer points in each direction when compared to the mass fraction model in an identical numerical framework. This represents an ~40 times lower computational cost for an equivalent error in two-dimensional computations. The proposed model is thus very suitable for Direct Numerical Simulation and Large Eddy Simulation of compressible mixing.


Introduction
This paper is concerned with the simulation of miscible, multicomponent fluids with single momentum equation for the mass-weighted mean velocity of the mixture, including the effects of viscous, thermal and molecular diffusion which may be solved using shock-capturing numerical methods, also known as diffuse interface methods. By making standard simplifying assumptions for flows of this type (see e.g. [1]), the governing equations can be shown to be the compressible Navier-Stokes equations, complemented by an additional transport equation(s) that is used to determine the composition of the mixture. One possible choice is to add one or more additional equation(s) for the transport of (N-1) species mass fractions, resulting in a system of governing equations that is fully conservative, commonly referred to as the mass fraction or fourequation model. However, it has been well documented that when using the standard four-equation model to simulate flows where the ratio of specific heats varies with mixture composition, spurious pressure oscillations are generated when material interfaces are advected through the computational mesh [2].
Other types of conservative discretisations, such as the level-set model, also suffer from this error under the same conditions [3].
The primary motivation of this current research is to enable accurate, highresolution viscous computations of compressible turbulent mixing problems.
These types of problems are significant in applications such as inertial confinement fusion [4], supersonic combustion [5] and astrophysics [6]. Of crucial importance is the understanding of fundamental flow instabilities that trigger transition to turbulence, such as the Richtmyer-Meshkov instability [7,8]. An example application is in the computation of the implosion of inertial confinement fusion capsules [4,9]. Initially there are sharp interfaces at solid/solid or solid/gas interfaces. The very high temperatures and pressures generated by the implosion process convert all the materials to dense gaseous plasmas and at various stages of the implosion there may be both sharp and diffuse gas/gas interfaces. The need to model both sharp contact surfaces and diffuse boundaries is also a requirement in shock tube turbulent mixing experiments (see, e.g. [10]).
For these applications the numerical technique must be able to compute flows involving multiple gaseous species which are initially separate but become mixed as time proceeds due to various hydrodynamic instabilities. The diffuse boundary layer between components may initially be small compared to the mesh size i.e. there is a well-defined contact surface. As the flow evolves, nearhomogeneous mixing due to species diffusion may occur at a scale which is at sub-cell level at early times, then becomes of the order the mesh size at time evolves. When simulating these transitional, non-equilibrium flows numerically, stability and high order of accuracy are required in order to properly capture all of the necessary physics. The purpose of the current paper is to derive an advection equation which may be added to the standard system of Navier-Stokes plus mass fraction equations to accurately model both of these limits, 4 including species diffusion, without spurious numerically-generated pressure oscillations.
For multiphase single velocity models, this problem has been thoroughly studied and has prompted a number of non-conservative or quasi-conservative approaches to be proposed. Karni showed that the pressure oscillations could be eliminated, either by using a primitive variable formulation (with corrections for leading order conservation errors) [11], or a conservative formulation augmented with a non-conservative pressure evolution equation [12]. Jenny et al. [13] proposed a correction for the energy equation, rendering the computation of the conserved variables a single-fluid computation and hence reducing (but not eliminating) pressure oscillations. The same paper also gave an expression for the relative error in pressure generated across a contact discontinuity, propagating with constant velocity u, that is initially aligned with the cell interface x i− 1 2 and has the CFL criterion 0 < σ = u∆t/∆x < 1, where T k are the temperatures for species k on either side of the contact surface, and γ k the ratio of specific heats for each species.
This applies to a conservative scheme that is first order accurate in time and space, includes a complete Riemann solver, with the ideal gas equation of state.
Hence differences in temperature and specific heat ratio across the interface, as 5 well as the convective velocity, all contribute to the size of the pressure error.
Abgrall [14] showed that such pressure errors could be eliminated by using an additional transport equation (in advection form) for a given function of the ratio of specific heats. This approach was dubbed quasi-conservative as it produces results with extremely small conservation errors. This class of algorithm was extended to a wider range of equations of state and multiple dimensions by Shyue [15] and Saurel and Abgrall [16]. A key idea behind these extensions was to have as many additional transport equations as there are parameters in the equation of state (e.g. γ for an ideal gas), while not requiring any additional transport equations to simulate mixtures of more than two fluids. However this results in increasing complexity of the algorithm with an increasingly complex equation of state. Other approaches for eliminating pressure oscillations have also been proposed that involve non-conservative modifications to the solution procedure rather than the addition of transport equations (see e.g. [3], [17], [18]) however they will not be discussed in further detail here.
Shyue [15] also gave a reformulation of the γ-based model of Abgrall in terms of the four-equation volume fraction ,which consists of (in 1D) the three conservative equations for density, momentum and energy plus a non-conservative advection equation for the volume fraction of one of the species. Allaire et al.
[19] and Massoni et al. [20] showed that by extending this to a five-equation It was shown by Johnsen and Ham [24] that although the models of [14,15] conserve total mass, momentum and energy, they do not discretely conserve the mass of each species and also generate temperature errors at material interfaces. These temperature errors are irrelevant in inviscid flows, however 7 they become important once physical diffusion effects are added. The authors proposed adding a transport equation for species mass fraction to the model of Abgrall [14] (in conservative form, thus enforcing species mass conservation) and showed that temperature errors were prevented if the reconstruction and upwinding for both the continuity and mass fraction equations is consistent.
A potential downside of this approach is that the location of the interface is not uniquely defined, although this difference is very small. Finally, Coralic and Colonius [25] applied the numerical framework of [23] to the five-equation model of Allaire et al., with the inclusion of the effects of viscosity as in [22].
Species diffusion was not included however.
At very large scales it may be sufficient to assume that the dissipative effects of physical mechanisms are approximated by those of the numerics (referred to as ILES, see e.g. [26,27]), however at smaller scales the effects of viscosity, thermal conduction and species diffusion become important [28] and require explicit modeling so as to gain confidence in the results. Previously, any computation of miscible compressible turbulent mixing with species diffusion has had to use the fully-conservative mass fraction model, even when the ratio of specific heats is not constant [29,30,31,32]. As outlined above, this approach is likely to suffer from errors in pressure and/or temperature generated at material interfaces. These errors are important even in relatively well mixed flows as species gradients are always present in the flow and thus spurious numerically generated pressure fluctuations may be large compared to physically generated 8 pressure fluctuations.

Standard Four-Equation System
The standard four-equation model describing advection and diffusion of multiple differing species may be written as follows: where, Note that q = −κ∇T where T = pm/Rρ is the mixture temperature with mixture molecular weight m = ( k Y k /m k ) −1 , R the universal gas constant, m K the species molecular weights. The enthalpy diffusion term where species sensible enthalpy h s,k = c p,k T . The specific total energy e = ǫ + |u| 2 /2, the internal energy ρǫ = k ρY k ǫ k , with ǫ K = c v,k T K and each fluid satisfies the perfect gas equation of state To close the above system without supplementary equations, it is assumed that the species are intimately mixed and are thus in pressure and temperature equilibrium. Given pressure and temperature equilibrium then p = (γ − 1)ρǫ, As stated in the introduction, it has been demonstrated that this system of equations suffers from spurious pressure oscillations in the inviscid limit when used to model the transportation of species of differing thermodynamic properties in a finite volume Godunov-type approach [2,14,3]. This has a substantial practical impact on the accuracy of computations undertaken at high Reynolds numbers, for cases evolving from an initially sharp interface, or with strongly varying temperatures and fluid properties. It is essential that a governing model for diffusive problems be able to represent gradients which may be steep compared to the local grid size without suffering from these numericallygenerated errors. Previous papers have added a extra equation for this purpose but only for the case when species diffusion is absent [20,19].

The Five-Equation Model
In order to model species diffusion in the presence of steep gradients in mixture composition, an additional equation is now derived which enables the preservation of pressure-constancy at contact surfaces in the sharp-interface limit, and is valid for problems with viscous, thermal and species diffusion. It should be emphasised that the basic physics is the same as for the standard four-equation model. However, the integral form of the equations used for numerical integration is different, as explained in Section 3, and it is this which eliminates spurious numerical pressure oscillations.
During diffusion, as mixing occurs at a molecular level, both pressure and temperature equilibrium is expected during the mixing process. Avogadro's hypothesis, ([33] 2.42), is assumed to hold: gases at the same pressure and temperature have the same molecular number density. As each molecule for a given species has the same mass, the number weighted mean velocity, the volumeweighted mean velocity and the mass-weighted mean velocity for species k are all the same and denoted by u k . The equation for evolution of species k particle number density n k : where u k is now the mean velocity of each component. By summation, the equation for the evolution of total number density is: where u = k n k u k N is the number-weighted mean velocity of the mixture.
Next, the individual number density equations are rewritten using the number fraction X k = n k N . Employing ∇ · (φF) = ∇ · φF + φ∇ · F where φ is a scalar, leads to: Thus there are three terms which modify X k ; advection with the numberweighted mean velocity, diffusive mixing, and pressure-temperature equilibriation. According to the last term, if a parcel of fluid k moves into a region with different pressure and temperature the number fraction adjusts to the local value.

Binary Mixtures
For a binary mixture the terms in the number fraction equation can be specified in a simple form. If Fickian diffusion is assumed dominant [33] then Substituting into Eq. (14) implies X k (u k − u) = −D 12 ∇X k and the number fraction equation becomes: The number-weighted mean velocity of the fluid u is related to the mass weighted velocity u as follows: where m k denotes the molecular mass for species k. Thus for two species, This form is useful as the further manipulations may result in numerical difficulties associated with dividing by number fractions or mass fractions which could be zero.
For binary diffusion, D kj are equal and gradient diffusion is dominant, then diffusive fluxes can be related to mass fraction gradients [34], and Eq. (14) can also be written as This is employed as the diffusion term for the evolution of ρY k . The full set of equations for compressible binary mixture advection diffusion problems is then: where, and and number density N = p/k b T where k b is the Boltzmann constant.

Multiple Species
The testcases presented in this paper will focus exclusively on binary mixtures, however the model can be extended to multiple species. In the case where diffusion is dominated by the ∇(X k ) terms [33], The number weighted velocity is then given by Here a simplified form of the number fraction equation is derived for the case where each species has an effective mixture diffusion coefficient D k , such as the well known approximation of Hirschfelder and Curtiss [35]. Giovangigli [36] demonstrated that the Hirschfelder and Curtiss approximated diffusion velocities can be defined as The velocity correction is employed to ensure that diffusive mass fluxes sum From Eq. (25) it follows that the difference between the number weighted mean and the mass weighted mean is given by It also follows that which leads to the following equation under the assumption of pointwise temperature and pressure equilibrium: Note that the number fraction diffusion fluxes sum to zero, as required. For numerical discretisation, it may be preferable to set Y j /X j = m j /m which would be less likely to suffer from inaccuracy at very low number or mass fractions due to machine error. Note that for fundamental studies of multispecies mixing using Direct Numerical Simulation it may be useful to assume constant diffusivity , D kj = D [37]. In that case Eq. (24) can be solved directly to give

Integral Form
The new set of governing equations for n fluid species can be written in integral form as: Note that the final equation need only be solved for n−1 components. Writing ∇N N and using the identity ∇·(ψV) = V·∇ψ+ψ∇·V, the number fraction equation may be written as a conservative equation minus a correction: For the applications considered here, interfaces may be initially sharp (subcell), then become diffuse as time proceeds. Thus it is important to clarify the definition of the volume averaged X k -variable. For a computational cell where X K is thus the volume averaged number fraction for species k. This differs from the number fraction based on the number of molecules of species k in the computational cell:- where N is the total number density. If there is a sharp interface within a cell, X k is equal to the cell averaged volume fraction z k commonly employed in multiphase systems [19,20,21]. If there is homogeneous mixing within a cell, Moreover, if pressure and temperature are uniform within a cell, then so is N (according to Avogadro's hypothesis: equal volumes of gas at the same temperature and pressure contain the same number of molecules) and in this more general case X k =X k . However, in typical cases X k andX k will not be equivalent, as volume averaged species temperatures can vary. It should be emphasised that the purpose of using the X k -equations is to give improved numerical accuracy; the basic physics is the same as that used in the mass-fraction based model.

Equation of State for a Finite Volume
The governing equations must be closed by an equation of state for the mixture within a finite volume. For cases where diffusion is negligible, or the interface is sharp relative to the mesh size, then the algorithm must robustly advect this interface without spurious numerically generated pressure waves. Here it is important to consider the 'volume-averaged' equation of state, rather than a 'mixture' equation of state, since the latter implies a mixture in the sense of Avogadro's hypothesis, whereas the former does not require equal species pressures and temperatures since gases may occupy separate regions of space.
In this paper, the system is closed by assuming that each species has the same volume-averaged pressure, however, volume-averaged species temperatures are different. At first glance this appears contrary to the application of Avogadro's hypothesis during the derivation of the system of equations, but 20 is a consequence of solving the integral form of the governing equations in a finite-volume framework. The fact that a flow where all species are in local pressure and temperature equilibrium can exhibit different volume-averaged species temperatures is entirely physical as will now be outlined.
The classical approach is to consider the advection of an isolated contact surface at very high Reynolds number where the effects of diffusion, viscosity and conduction may be neglected. The Riemann problem is shown schematically in Fig. 1. In this case, the sharp interface is perfectly aligned to an interface between two finite volumes. The exact solution of this simple advection problem is a constant pressure contact surface moving at velocity u, and represents the simplest case of inviscid advection. This sub-section will demonstrate the form of the volume-averaged equation of state which does not produce spurious numerically generated pressure oscillations in this limit.
The contact surface separates a left gas mixture L and right gas mixture R with properties: where the velocity u > 0, the mixture pressures are equal but have different temperatures and mixture compositions. All species are in local pressure and temperature equilibrium.

21
At time t + ∆t the contact surface has swept out a volume V L = u∆t of the downstream cell. This volume now contains mixture L at a temperature T L , and the remaining volume of the cell V R = V − V L contains mixture R at a temperature T R . The fraction of the cell volume which is occupied by mixture Volume-averaged number fractions X k can be defined in the downstream cell at t + ∆t as: and internal energy is given by, Defining a species partial density and note that as each individual mixture is in local pressure and temperature equilibrium then Given ǫ k,(L,R) = c v,k T (L,R) then internal energy may be written as: The analytical solution to this problem is a constant pressure contact surface.
To satisfy this condition, it is clear that the appropriate volume-averaged effective ratio of specific heats is: With this condition, the downstream internal energy ρǫ = p/(γ − 1), and pressure in the mixed cell is maintained at the original pressure p following the advection step as required by the analytical solution. This simple analysis can be extended to more complex equations of state following analagous derivations for immiscible fluids [19,21,20], and illustrates the key useful property of the governing model which will hold when the equations are discretised with a contact-resolving Riemann solver. Subsequent time steps may be treated in a similar fashion by considering the advection of sub-components of mixtures each at their respective pressure and temperature equilibria.
Note that the above derivation is equivalent to specifying that the appropriate volume-averaged temperature for each species is T K = pX k /R k ρY k . In the limit of an infinitely small cell then T k = T as expected from the assumptions underlying the governing equations. Assuming constant specific heats for each species such that T k = ǫ k /c vk , the volume-averaged species temperatures T k are mass-weighted averages of the initial temperatures T L,R . Enforcing temperature equilibrium would change the number fractions, X k , and this will give a different mean γ for the cell and give rise to a change in pressure if γ 1 and γ 2 are different. Hence homogenizing a mixed cell will, in general, give a non-physical change in pressure at the contact (unless γ 1 = γ 2 ).
The isobaric closure proposed above is designed to correct this problem, which manifests itself clearly in the standard mass fraction model, where even if the species mixtures are initially in local pressure and temperature equilibrium (as in the current problem), the analytical cell averaged species temperatures can be expected to differ following an advection step. The above closure is only possible with the addition of the equation for transport of X k .

Numerical Methods
The governing equations are implemented within the University of Sydney code Flamenco [38]. This code has an existing inviscid implementation of the volume fraction model of Allaire et al. which has been well documented in previous publications [19,39]. Given that the inviscid part of the governing model proposed here is in the same mathematical form as the model of Allaire et al., the discretisation may be undertaken in an analogous fashion.
The existing scheme is based on a Godunov-type method of lines approach in a structured multiblock framework, which has been modified to provide reduced dissipation in low Mach number regions of the flow [40,41,42]. This modification is low cost but restores the accuracy of the compressible algorithm demonstrated for M > 10 −4 . Most computations with this algorithm utilise a nominally fifth order reconstruction for the inviscid component [43], and second order central differences for the viscous and diffusive terms. Temporal integration is achieved via a second order accurate in time TVD Runge-Kutta method [44]. Note that in more than one dimension the formal order of accuracy would be second, given the one dimensional reconstruction stencils.
This approach is applied in a straightforward manner to all of the conservative equations, however a different approach must be chosen for the number fraction model.
The Jacobian for the two-species system written in primitive variables [ρY 1 , ρY 2 , u, P, X 1 ] T is as follows: where for an ideal gas the speed of sound c 2 = k ρY k h k /ρξ, ξ = 1/(γ − 1) = k X k /(γ k − 1). The system has a full set of eigenvectors, with the 25 corresponding eigenvalues: The form of the Jacobian given above is not conventional as it depends on derivatives of the solution variables. These terms are included here on the left hand side of the system of equations for numerical reasons (discussed in detail in Section 3.3), namely that they modify the upwind direction of X k .
A justification for this is given by applying Cattaneo's relaxation approach [45,46], whereby a new variable ψ is introduced in place of ∂X 1 ∂x along with a transport equation that relaxes ψ towards ∂X 1 ∂x by means of a stiff source term, thus enabling the Jacobian to be written in terms of solely the solution variables. Further exploration of using this type of relaxation approach to allow for a more efficient discretisation will be the subject of future work.
Here, the X k equation is discretised following Abgrall [14], where the cell interface flux F i+1/2 and F i−1/2 in physical space for cell i is computed as where (.) RS indicates a term arising from the solution of the Riemann problem at the cell interface, which in the current algorithm employs the HLLC approximate solver.

26
The where u RS is gained from the solution of the classical Reimann problem. The upwind number fraction (X RS k ) required to compute the number fraction fluxes is determined via a modified HLLC approach where the signal speeds incorporate the additional diffusion velocities, i.e. the contact surface in the Riemann problem for the number fraction is assumed to advect at a velocity U, with the gradients given by the aforementioned second order central difference approximation. It is worth noting that recent advances could permit an arbitrarily higher order discretisation, which will be the subject of future work [47].
Next, the diffusion terms must be discretised in a manner which is consistent with the underlying physics. As diffusive fluxes are computed at the cell interfaces, Avogadro's hypothesis must be applied. This is a consequence of the hypothesis of local pressure and temperature equilibrium invoked in the derivation of the governing equations, i.e. the species diffusion terms require pressure and temperature equilibrium at the cell interface. If pressure and temperature equilibrium is not enforced during the computation of the diffusive fluxes, then the resultant species densities and temperatures will be in error. Within the current algorithm, this is addressed by rewriting the mass fraction gradient in the diffusive flux as: where W k is the species k molecular weight and W is the mixture molecular weight.
All cases use a fifth-order accurate limited scheme to interpolate from cell averaged quantities to cell interface values in the computation of the inviscid terms [43], however it must be noted that the current discretisation is formally second order accurate due to the discretisation of the quasi-conservative number fraction equation, viscous and diffusive terms.
For explicit time integration, as employed here, the stable time step is chosen based on the minimum in the whole domain of the following: For example, in the highly compressible two dimensional Richtmyer-Meshkov problem detailed in Section 5, a diffuse boundary forms after shock passage, and there is little pressure variation across the this boundary layer. Hence 1D test cases based on diffusion at a boundary between two gases at uniform pressure provide highly relevant test cases of the new model. If diffusive velocities are small compared the sound speed, the mixing process is near incompressible and results can be checked against analytic solutions for the incompressible limit. In applications, the diffuse boundary layer will be advected through the mesh and this is a more challenging situation for the numerical method.
Hence, results are also given for a case when a uniform velocity is added to the flow. Note that the numerical results are expected to agree with the analytical results only up to the point where the difference between the numerical and analytical results is not dominated by compressibility effects, which are absent from the analytical solution. This occurs at the finest grid resolutions where the error due to the discretisation (truncation error) becomes lower than the difference due to compressibility and thus the convergence rate measured relative to the analytical solution stalls. It occurs at the same error magnitude regardless of whether the mass or number fraction models are employed.
The first three cases are one-dimensional, where Case 3 is the same as Case 1 but with a net mean velocity added. The one-dimensional problem of diffusion at a plane boundary between two gases is a useful test case for numerical implementations. For mixing between two gases at the same pressure to be an incompressible process, there should be no change in pressure when the two gases mix. This will be true for two different gases at the same temperature or for two identical gases at different temperatures. In Case 1 and Case 3 heat conduction is unimportant. In Case 2 heat conduction is a dominant process.
Case 1 and 2 are used to verify the algorithms for near-pure diffusion cases, whereas Case 3 tackles a more realistic case with both advection and diffusion.
In these three one-dimensional cases viscosity is assumed to be zero, and equal diffusivities are used (D kj = D and κ = ρc p D).
Note that the variable reconstruction employed within Flamenco includes limiting at maxima and minima for all test cases presented in this Section.

Species
This test case has been employed in Kokkinakis et al. [50], however here the initial conditions are described in full, along with a rigorous initialisation process which enables a formal convergence study, up to the limit of assumed incompressibility. If the initial uniform pressure is sufficiently high, the mixing process is quasiincompressible and the analytic solution for the number fraction distribution is given by where x 0 = 0.5, h 0 = 0.02 and the diffusion coefficient D = 0.01.
In the initialisation of the problem, compressibility effects are minimised by assuming that the volume weighted mean velocity is initially zero. Hence the initial mass weighted mean velocity is given by u = − D ρ ∂ρ ∂x .
Noting that a finite volume algorithm requires the cell averaged quantities for the initialisation, the following exact initialisations are given for completeness: For the standard test case, the initial uniform pressure is p 0 = 10000. The For p 0 = 10000, the flow is near-incompressible with pressure and temperature remaining approximately constant, thus the incompressible analytical solution may be employed as a reference solution to demonstrate the scheme's order of accuracy. Note that for Case 1, the (u − u) · ∇X 1 term in the number fraction equation is essential, but the ∇N/N term vanishes in the incompressible limit (Case 2 will validate this term).

Quasi-Incompressible Diffusion
This section focusses on the results for the quasi-incompressible case where p 0 = 10000. Tables 1 and 2 document the L 1 , L 2 and L ∞ errors in the simulated number fraction profile at t = 0.5 for the number fraction and mass fraction        feature of the governing equations.

Compressible Diffusion
This section focusses on the results for the compressible case where p 0 = 10.
In this case there is no analytical solution as the density distribution varies slightly from the analytical, hence Figure 5 plots the comparison between the mass fraction, number fraction and Lagrange-remap results. The key result here is that all models converge to the same solution in the fully compressible limit using 256 cells.

Case 2: Diffusion of a Contact Surface with a Two Identical Species at two Temperatures
The computation domain is 0 ≤ x ≤ 1, and grid sizes from 32 → 2048 cells have been employed. Reflective boundary conditions are used at the left and right hand boundaries. The two fluids have the properties ρ 1 = 20, ρ 2 = 1, γ 1 = γ 2 = 5/3. The specific heats satisfy the requirement that c v1 = c v2 , and D = 0.01. A uniform initial pressure is achieved by specifying that T 2 = 20T 1 .
As previously, p 0 = 10000 D = 0.01 is used for both the species diffusvity and the heat diffusivity, and the solution is run to t = 0.5. As the heat diffusivity is constant, the error function solution for the density distribution is again applicable for the quasi-incompressible case, The initial density and mass fraction distributions are the same Case 1, implying that the initial number fractions are different. Thus the initial cell averages ρ, ρY 1 and ρu are the same as Case 1. The initial distribution (and incompressible analytical solution) of X 1 is equal now to Y 1 and is given by: This function does not have an exact antiderivative, so the initial conditions for X 1 are computed by a five-point Gaussian quadrature scheme.
In the incompressible limit, the density and mass fraction distributions are the same as for the uniform temperature case. However, number fractions are Table 3 Case 2 convergence rates for the number fraction model.  Table 4 Case 2 convergence rates for the mass fraction model.  observed convergence rates. Figure 6 plots the error norms as a function of number of points N x for the mass fraction and number fraction models.
As this case includes two species with the same ratio of specific heats, then the differences between the mass and number fraction formulations should be lower. This is the case, and both algorithms converge at the expected second order of accuracy up to the point where compressibility impacts the agreement with the incompressible solution. As with Case 1, the errors with the mass fraction model increase at 4096 points which confirms this hypothesis.
The number fraction model again has lower error at a given grid resolution, 41 however the differences are not as substantial as the previous case. This is to be expected as the principle errors in the mass fraction approach are generated when the ratios of specific heat vary. For x ≤ 1, the analytical solution for the number fraction distribution is given by where x 0 = 0.5. The number fraction profile is then mirrored about x = 1, thus for 1 ≤ x the number fraction is given by where x 1 = 1.5, h 0 = 0.02 and the diffusion coefficient D = 0.01. The initial cell averaged quantities are gained by simply mirroring the expressions detailed Table 5 Case 3 convergence rates for the number fraction model.  Table 6 Case 3 convergence rates for the mass fraction model.    To demonstrate the applicability of the new algorithm to compressible mixing problems, Case 4 undertakes direct numerical simulations of a two dimensional single mode Richtmyer-Meshkov instability triggered by a shock wave passing from air to SF 6 . The setup for this case is adapted from that proposed by Tritschler et al. [29]. The purpose of this case is to illustrate the advantage of using the proposed number fraction formulation in a more practical setting, since there is widespread experimental use of air-SF 6 combinations in shockinduced turbulence (see, for example, [51,52,53,54]).
To conduct an accurate simulation, the initial conditions must be very well resolved. Here, a single mode perturbation of wavelength 0.5mm and initial amplitude 0.025mm is computed. The initial diffuse layer thickness is 0.1mm, suf-49 ficiently low to prevent substantial damping of the initially imparted impulse.
The initial state of the unshocked gases are p 0 = 23000Pa and T = 298K.
For simplicity a fixed viscosity µ = 2.243 × 10 −5 Pa-s is specified [55], along with Sc = 1 and P r = 1. Note that inviscid computations have also been run to represent the very high Reynolds number limit. Air is assumed to have γ air = 1.4, and γ SF 6 = 1.1, and the molecular weights are M air = 28.964 and M SF 6 = 146.057 [29]. The shock Mach number is 1.5, and the shock has an initial offset of 1mm from the mean interface position. The interface diffuses slightly prior to shock interaction, increasing the integral width by 2.7% at the time of shock interaction in the cases incorporating viscosity, diffusion and conduction. All results are scaled by the modal wavelength λ and Richtmyer's velocity.
The details of the initial conditions for this two-dimensional problem are derived as follows. The incompressible limit of mixing by species diffusion gives the following non-zero divergence of velocity [56]: Where the diffusion velocity is given by: In 2D this becomes: Then integrating gives: which is applied only at the interface (not at the shock). Finally, taking cell averages: Where: With F (x, y) being the minimum distance from the point (x, y) to the perturbed interface. Equations (60) and (61) are calculated numerically using a 5-point Gaussian quadrature rule, while the minimum distance F (x, y) is also calculated numerically using a bisection method iteration. example [57,49,58,59,60]) are defined as: These quantities are plotted in Figure 13 for the mass fraction and num-   Thus the observations of errors in the prediction of Θ and Ξ are as expected from the previous cases, however W is not. Figure 15 plots the x variation of the planar averages X 1 X 2 and X 1 X 2 for the mass fraction model and 54    Figure 15. Plot showing the specific terms employed to compute the integral properties for the mass fraction (MF) and number fraction (NF) models at converged (512 points) and non-converged (128 cells) grid resolutions.
The integral width is the area under the curve X 1 X 2 . From this Figure   it can be seen that the mass fraction model has larger errors in X 1 X 2 than the number fraction model, however those errors are approximately equal on both sides of the mixing layer. Thus although the integral width appears to converge faster for the mass fraction model, it is simply that the errors are cancelling each other.
Turning to X 1 X 2 , the integral of which is the numerator in the Θ equation, there is a very large error when using the mass fraction model compared to the number fraction model at 128 ppw.
Overall, the solution of the number fraction model at 64 ppw has an equivalent error to the mass fraction model 256 ppw computed using identical variable reconstruction, Riemann solver and discretisation of the viscous terms. This represents an enormous computational saving for equivalent error, a factor of 16 fewer points in 2D along with a factor of 4 fewer time steps. Even 56 with an observed (unoptimised) increase in computational cost of ≈ 50%, the computations are on the order of 40 times faster for equivalent accuracy.
For three-dimensional flows the gain would be expected to be greater than 100. This clearly demonstrates the superiority of the newly proposed number fraction based governing equations.

Conclusions
This paper has presented a new number fraction based five-equation model for miscible fluids incorporating viscosity, species diffusivity and thermal conduction. An equation is solved for the variable X which reduces to the usual definition of volume fraction when there is a sharp interface in a computational cell and to the species number fraction if there is homogeneous mixing within the cell. A numerical scheme has been derived to solve these equations which is second order accurate in space and time, and respects the underlying physics.
Four test cases have been proposed, and computations have been run using the standard mass fraction model and the newly proposed number fraction model.These test cases were designed such that they validate and verify all key terms in the governing equations. Three of the cases have analytical solutions in the incompressible limit which were used to compute the observed order of accuracy.

57
The discretisation of the number fraction and mass fraction models were both demonstrated to be second order accurate in space. This was shown to be valid up to the point at which compressibility effects became important. The first two cases are relatively straightforward for both models, having either constant temperature or constant ratio of specific heats. However, the errors of the number fraction model were consistently two or three times lower than that of the mass fraction model for the same grid.
The third test case was designed to be more challenging. The third test case consists of a diffusing contact surface between fluids of different thermodynamic properties moving with a mean velocity. Although both models converge at second order accuracy, the number fraction model is one order of magnitude more accurate for a given grid resolution.
The final test case is of a commonly employed gas combination in experiments of shock-induced turbulent mixing. A shock wave impinges on to a twodimensional perturbed interface between air and SF 6 , triggering the growth of a Richtmyer-Meshkov instability. In this test case, the number fraction model is approximately converged on all measures at 128 cross-sectional grid resolution, whereas the mass fraction model is not converged until 512 points. This represents a computational saving of approximately 40 times for the equivalent accuracy.
Based on these results, the number fraction model is clearly superior to the standard mass fraction approach for the computation of compressible turbu-lent mixing problems of miscible fluids with distinct thermodynamics properties. This formulation is directly applicable to Direct Numerical Simulations (as undertaken here) and Large-Eddy-Simulations.