 So again, yeah, so this talk is about the derivation and new set of governing equations for single-phase, multi-species flows, including a full range of diffusive effects, so species diffusion and then viscosity and thermal conductivity as well. So it was initiated by David Young's and my supervisor, Ben Thelumba, and then in the last year or so, I've also come on board and kind of progressing this work a bit further. So I'll be talking about some problems that occur with multi-species finite volume solvers and then the derivation of an existing set of governing equations that eliminates these problems to include the effects of species diffusion. Briefly mention the discretization of these new governing equations and then also give some preliminary simulations and results and then also briefly talk about what we're going to be doing in the future. So our motivation is we want to study the Rikmai-Meshkov instability using high-resolution simulations, so DNS, LES, with shock-capturing finite volume methods and in particular, we want to be looking at the intermediate time, transitional dynamics. So it's our kind of opinion that very early-time dynamics, so the linear regime can be fairly well-described using various theoretical approaches as can the later-time dynamics using something like self-similarity. But then this intermediate time is reasonably hard to describe and so we're initially going to be focusing on the planar Rikmai-Meshkov instability and also hope to transfer that in any insights into the converging cases as well. So a brief review of some fundamentals, so the governing equations for single-phase, compress all multi-species flows, the Navier-Solks equations that are augmented by some additional transport equations with the standard choice just being species mass fractions, so I think the first equation there is just conservation of species mass and then you have conservation of conservation energy and the nice thing about these set of equations is they're fully conservative, so total mass, species mass, momentum and total energy are all conserved. But there's some problems when you actually go to simulate these models and discretize them on a computer. So for finite volume methods, the reconstruction phase, this approach actually leads to oscillations in pressure and temperature across material interfaces and also contact surfaces in a single fluid. It has been shown that for a first-order scheme, pressure equilibrium isn't maintained when the temperatures on either side of a contact surface aren't equal as well as the ratio of specific heats not being equal and also this is further exacerbated when the interface is effected through the computational mesh. So the problem essentially boils down to the fact that because finite volume methods deal with cell averages, the cell average value of the mixture gamma calculated from the species masses is not actually the same value which is required to conserve pressure equilibrium. So, and further problems result when you start using more advanced numerical methods, so more advanced reconstruction approaches, say high-order reconstruction in something like primitive variables or even characteristic variables can actually shift the error, so you can also see these errors occurring in the temperature field as well. In theory, these will become negligible in the DNS limit as the interface between the two fluids is fully resolved and then this cell average value approach is the same value that's required to conserve pressure equilibrium. However, the simulations will be less computationally efficient because you'll require a very large number of cells to actually approach this. So, a layer at our, in this paper down here, showed that you can circumvent the problem by using a non-conservative approach. So, the, based on volume fractions, so this is still the same species mass equation but now written in terms of the partial density and the volume fraction of each species and then these two, these two occasions are also the same and then we have an additional advection equation for the volume fraction of one of the species. And the nice thing about this model is that it gets rid of all those problems that I was just talking about. However, and it's also only non-conservative for this equation, all these equations are still conservative equations but the, unfortunately, it's only really, it was only originally derived for inviscid flows. So, there's no inclusion of effects of viscosity or thermal conductivity or species diffusion. And so, in this work, actually, you take that same model and include these effects and retaining all the favorable properties of the model. So, the fact that it preserves pressure equilibrium. The individual species temperatures are also retained in the mixture, which is quite nice and it's also valid for general equations of state. So, briefly go through the derivation. This is gonna go against what I just said but we assume that the species impression a temperature equilibrium but only at a single point on the face where the diffusive flux acts and then in the rest of the computational cell using isobaric closure, the pressures of, so here's just for two species but yeah, the pressures are all equal but the species temperatures are allowed to vary. From this, we say that using Avogadro's hypothesis because the species at these two points are at the same pressure or temperature, they also have the same number density which then allows us to stay these next two assumptions which is that the volume fractions may be equated with the number density fractions for each species and also for each species, the number-weighted, volumated and mass-weighted velocities are all equivalent but that doesn't hold true for the mixture velocities. And then we also assume the diffusion just to obey a simple fixed law. So, I won't go through the derivation before here but what we do is we start with the evolution equation for total number density, the mixture and then applying the previous assumptions we arrive at this equation here which allows us to see that there are three processes that modify the volume fraction of each species. So, advection with the number-weighted mean velocity which is this guy here, diffusive mixing which is the second term and then also this third term that arises from pressure temperature or the assumption of pressure temperature or collaboration. So, a few more further manipulations gives us our new set of governing equations and this is for binary mixtures but it could also be extended for more than two species. So, these first three equations are all just the same from the previous mass fraction model where so the right-hand side fluxes are the same. Just written the species mass equation in terms of volume fraction instead and then we also have our final volume fraction equation where these terms here basically come from converting from number-weighted to mass-weighted mean velocity. So, now a few notes on how we discretize these equations in our in-house CFD code Flamenco. The existing algorithm itself, we already have a discretization of the original layer volume fraction model using fifth order reconstruction for the inviscid fluxes plus a low mat correction and then second order central differences for the viscous flux and the temporal discretization is second order Runge-Kutta in a method of Lyons approach. And so, we make the following modifications to accommodate the new model. So, the additional terms that modify the upwind direction of the volume fractions. So, all of these terms here are again, they're also computed using second order central differences. The, we modify the intermediate signal speed of the Riemann solver just for the volume fraction to be this instead. And then, we also need to make sure that the discretization is consistent with the assumptions that we made. So, we actually modified these right-hand side terms here to be consistent with, so if under the assumption of pressure and temperature equilibrium, then the mole fraction here will actually be equal to the volume fraction. So, then this ensures that we've conserved, we've ensured consistency with the left-hand side terms because we basically need the diffusion coefficients here to be the same as the ones here as well, which are occurring in the inviscid flux discretization. So, to validate that what we do is actually seems to be correct. We used first year three fundamental cases in 1D. I'll only talk about two of them here because the third one is fairly simple. So, they're both just simple 1D diffusion test cases but simulated with the full set of governing equations. The first case is just, you can see an animation of the first case there, just this diffusion between two initially stationary gases and then case two is the same case but with an additional advection velocity added to the interface. So, the first test case comes from this paper here. It's just in a domain of zero to one which has reflective boundary conditions. We have a 20 to one density ratio and then values of ratio specific heats of each species given there and nice thing about this is that we can actually, there actually exists a analytical solution assuming the solution is fully incompressible and for fairly coarse grid resolutions and this ax also is essentially the same solution of the fully compressible equations. We also fix the diffusion coefficient and the, yeah, we can basically see that both for the mass fraction equations and also for the new volume fraction equations both initially converge at close to exactly second order as we expected, which is a nice validation. What's happening here is the fact that this is an incompressible solution but we're solving compressible equations eventually small compressibility effects become resolved so because this is a reflective boundary conditions basically a very small acoustic waves that are just bouncing back and forward across the domain that we eventually start resolving at higher grid resolutions. You can't quite see here but for the mass fraction equations if you go to even higher grid resolutions the arrows actually jump back up again once these compressible effects are also resolved. The second test case is just the same as the first test case but have this additional velocity interface. We also double the domain and make the boundary conditions periodic so the interface will affect and then to the right and then back again to where it started. So the analytical solution is the same as before which is mirrored about the middle of the domain so actually affecting two interfaces and this test case actually gives us a really nice example of why we set out to derive these new set of equations because basically for cases where the interface is advected at a suitable velocity through the mesh the mass fraction approach will actually generate fairly large errors and so we can see that we're getting equivalent results for around four times less than our points using the volume fraction model which given that the new set of equations are about one and a half times more computationally expensive to solve this equates to roughly a 10 time saving in computational effort just in one day. So the next test case we're looking at is a 2D single mode Rikram-Meshkov instability between air and SF6. The case setup is a Mach 1.5 shock. Pressure and temperature there whereas the ratio of specific heats of both species we make the assumption just for simplicity that the parental estimate numbers are both one. The initial amplitude is a 10th of the wavelength so we're roughly at the end of the linear regime and then we also have initial diffuse interface so that we can fully resolve the interface. The Reynolds number based on the Rikmaya velocity the wavelength and the average viscosity is around 5,000 and we perform simulations in both a stationary and a moving frame of reference so stationary would be with the laboratory frame of reference and then the moving frame of reference is in the post shock reference frame. So we're particularly interested in how the molecular mixing fraction which is given there based on plane averages of volume fractions we're particularly interested in how that converges. So we see for the volume fraction model and the results on the left here the convergence is nice and uniform we actually get to these curves collapse by about 256 points across the initial wavelength whereas for the, so this is in the stationary reference frame and whereas for the mass fraction model we said that the convergence is not uniform at all and actually it doesn't even converge by the largest number of grid points that we considered. And here's an animation it's at a slightly, not at the highest number of grid points at a slightly coarser number to really accentuate the differences between the two approaches. So you see the approach on the left the animation on the left here it's for the new volume fraction equations whereas on the right is for the mass fraction model and store your attention to two things that really stand out the shape of the spike it seems to flatten out a lot more which we deemed to be an error to the true solution and also the interfaces seems to be a lot thinner when you compare there and there it's basically due to the fact that pressure and temperature are actually an error across the interface. So, but in the results for a moving frame of reference so the interface is now stationary with respect to the mesh ignoring the actual evolution of the instability the results are much more consistent and the mass fraction results actually also converge quite nicely and while it's not shown here these two curves actually collapse to the exact, the same single curve. So this is nice but it also so this seems to say that well the mass fraction equation seems to perform all right when if we can perform in this moving frame of reference so that's reasonably stationary with respect to the computational mesh however there are examples such as when you're trying to yep, you're trying to simulate reshock experiments for example that this wouldn't be possible so our new set of equations would seem to be much more beneficial for those kind of situations. And here's just an animation here basically just showing that the two approaches give identical results and finally in the last week or two I've just been looking at running a 3D multi-mode case so this isn't actually to demonstrate the advantage of the new model or the old model, it's more to kind of suss out how the actual algorithm performs what kind of, what's the kind of maximum Reynolds number we're currently capable of simulating as we also have some some ALES results for this case already so it's nice to compare to those so the former DNS of this standard problem from the feeder group collaboration which is recently published on Archive which was a large cross-code comparison of planar Rikmo-Meshkov instability the case setup is given again here it's now a Mach 1.8 shock three to one density ratio the gamma of both cases are the same and the Reynolds number which is based on the initial integral width velocity the mean wavelength and the average viscosity is 455 so the initial perturbation is a narrow band perturbation with a constant power spectrum with the standard deviation is a tenth of the minimum wavelength to ensure linearity of all across all scales and the simulations are performed in a moving frame of reference so looking at convergence of both the integral width which is defined there as well as the molecular mixing fraction again we see that both the original inviscid simulations and also our new viscous simulations all collapse to the same curve for the integral width these are on non-dimensionalizes is non-dimensional time and the integral width is non-dimensionalized by this should be the average wavelength whereas for the molecular mixing fraction we almost achieve a fully converged solution at 512 cubed grid points and also we see that the viscous simulations mix a lot quicker than the original inviscid simulations so we can actually because we're interested in this range of non-dimensional time we can actually use Richardson extrapolation to estimate say the point of minimum mix so this point down here and for the three finest grid resolutions we can there's a grid convergence index for the between 128 and 256 here and then between 256 and 512 and for these three grids we can see that we've pretty much reached the asymptotic range of convergence and then using Richardson extrapolation we extrapolate to find the point of minimum mix and also the non-dimensional time that that occurs at and so here's a nice animation of that case see there's clearly not fully turbulent or you do still see some spikes that get actually advected away from the main mixing layer and also fluid that's initially entrained say right here actually mixes and diffuses quite quickly so I'll let that play through once more one more slide and then just conclusions that's a quick, be real quick so just with thriving set of governing equations demonstrated advantages over traditional approaches have been assessing the efficiency of the new algorithm and basically in future work we plan to extend the color in our algorithm to much more high orders of accuracy so that we can then kind of ideally reach order of magnitude higher Reynolds number with using DNS yeah thank you very much for your time and have you take any questions so for those fundamental test cases it was just fixed but in general we just computed off mixture viscosity and turbulent shipment number for binary diffusion it will be yeah but for more than two species then not necessarily yeah so it's fifth order in one dimension but due to the fact there's only a one dimensional reconstruction stencil in say three days you only a formally second order accurate and also only second or accurate in time and in the viscous fluxes anyway so the plan is to use kind of more state of the art final volume approaches where you can discretize both the left and right hand side of the equations in one kind of unified framework and then so the higher order reconstruction can perform like a multi dimensional reconstruction and also hopefully only have to perform a once a time step so it becomes a lot less expensive yeah well it's based on it's eight of methods so eight of DJ and eight of final volume is what we're looking into here you know it is for the for that last simulation by the time we got down to the 512 degree resolution yeah we were hitting viscous time stepping which we kind of figured was appropriate given that if you're approaching the DNS then we kind of assumed that you would be actually hitting the viscous time step limit we should be at least yeah