A Discontinuous Galerkin Method for Viscous Compressible Multifluids
Abstract
We present a generalized discontinuous Galerkin method for a multicomponent compressible barotropic NavierStokes system of equations. The system presented has a functional viscosity which depends on the pressure of the flow, with the density and the local concentration . High order RungeKutta time discretization techniques are employed, and different methods of dealing with arbitrary coupled boundary conditions are discussed. Analysis of the energy consistency of the scheme is performed in addition to inspection of the relative error of the solution compared to exact analytic test cases. Finally several examples, comparisons, generalizations and physical applications are presented.
Keywords: NavierStokes; discontinuous Galerkin; RungeKutta; mixing; barotropic; compressible; viscous; miscible; multicomponent; multiphase; multifluid; chemical; acoustic.
1§1 Introduction
Much work has been done in the study of the numerics of multicomponent flows. An example of an early yet comprehensive study of computational multiphase mechanics was given by Harlow and Amsden in Ref. [25], where they developed an implicit finite differencing technique for extremely generalized multicomponent settings of both compressible and incompressible flows, including phenomena ranging from bubble formation and cavitation effects, to the formation of atmospheric precipitation and mixing jets. Subsequent and related work in multicomponent flows followed with, for example, the work of J. Dukowicz in Ref. [18] for particlefluid models of incompressible sprays, an approach extended by G. Faeth in Ref. [20] to combustion flows and by D. Youngs in Ref. [56] to interfacial turbulent type flows.
Owing to some of these pioneering works, recent work has demonstrated a resurgence of interest in multicomponent flows, approaches and numerical techniques. The importance of fluidflows comprised of more than one phase, chemical constituent, species or component is represented by a vast array of applications that range across a number of fields. For example, multicomponent flows are essential for any flow demonstrating even rudimentary chemical kinetics; hence, for all (nontrivial) “chemical fluids” [55]. Likewise biological flows often require phase separations, in order to resolve membrane dynamics and interfacial behaviors in cells and cell organelles [46] and medical applications desire estimates in local componentwise variations in blood serosity, which effect the viscosity and flow parameters involved with pulsatile hemodynamics [9]. Likewise we find numerous examples of multicomponent flow applications in the atmospheric [26] and geophysical [52] sciences; as well as in acoustics [39] and astrophysics [42], just to mention a few.
Here we present a new multicomponent numerical scheme based on a mathematically wellposed [38] compressible barotropic system with functional viscosity depending on both the density and the mass fraction of each fluid component. It is wellknown, both experimentally and theoretically, that viscosity has a functional relationship to the density and specie type (for examples see the NIST thermophysical properties server). In addition, these types of mathematical models (with functional transport coefficients) are well understood from the point of view of continuum dynamics, having been extensively studied by Málek, Rajagopal et al. in Ref. [33]. It is further seen in Ref. [38] that the analytic model used in this work a priori satisfies two essential entropy inequalities, much like the shallow water equations [8], which serve as important tools for numerical analysis and implementation.
In this paper we implement a discontinuous Galerkin (DG) finite element method, employing piecewise polynomial approximations which do not enforce or require any type of continuity between the interfaces of “neighboring” elements. This particular implementation is primarily motivated by the works of Cockburn, Shu et al. (see Ref. [12]) and Feistauer, Dolejší et al. (see Ref. [22]). We implement a generalized formulation that is designed to accommodate an arbitrary choice of inviscid, viscous, and supplementary numerical fluxes. We use explicit time discretization methods as described in Ref. [12], which necessitate a conditional stability requirement; namely the time discretization must satisfy the CFL condition. Up to the CFL stability condition we find our method to be very robust and to deal well with arbitrary numbers of fluid components of arbitrary type — up to the additional assumption that a barotropic pressure law is applicable. On the domain boundary data we again strive to generalize our setting. We show two different implementations of boundary conditions, which demonstrate different solvency with respect to interior solutions, initial conditions and phenomenolgically relevant contexts. In both cases arbitrary Robin type BCs may be set.
In we give the general governing system of equations, the mathematical regularity, and the discrete formulation of the problem. In we demonstrate a general way of dealing with boundary conditions by way of the method of characteristics, or alternatively, by way of setting arbitrary data on the boundary. We provide an explicit formulation of the characteristic technique and show the generalized behavior of these types of “characteristic” boundary conditions, while subsequently discussing a number of alternative approaches. In we implement two test cases with exact solutions, which are restrictions placed on the multifluid barotropic governing equations, and show that they are exact up to the possible exception of the boundary data. In we show an example of a bifluid solution using the forward Euler method. We then show the difference between boundary conditions by way of weak entropy solutions versus that of characteristic boundary solutions. The next section, , is used to generalize the setting to fluid components and th order RungeKutta schemes, where the example of an fluid is shown explicitly. Then in we analyse the energy consistency of the modelisation with respect to two entropy inequalities derived in Ref. [38]; one the classical entropy and the second a closely related entropy discovered by Bresch and Desjardins (see Ref. [5]), where it turns out that the numerical scheme from satisfies both energy relations provided the CFL condition is satisfied. Finally in we extend the results to include Fick’s diffusion law, where we inspect the exotic physical setting of a pressure wave traveling through a gas comprised partially of polyynes, and discuss some applications.
2§2 The generalized lfluid
We consider a one dimensional compressible barotropic fluid system governed by the following system of equations:
with initial conditions,
The multicomponent barotropic pressure is chosen to satisfy,
where . The mass conservation (Equation 1), momentum conservation ( ?), species conservation ( ?), and barotropic equation of state (Equation 2) describe the flow of a barotropic compressible viscous fluid defined for . Here the density is given as , the velocity as , the momentum by , and the mass fraction of each component (chemical specie, phase element, etc.) of the fluid is given by , respectively, where corresponds to the emperically determined adiabatic exponent uniquely characterizing each of the species. Furthermore, adopting the notation throughout the paper that , the form of the viscosity functional is fixed to satisfy
for given and as emperically determined constants (see Ref. [33]) and ).
The mathematical wellposedness of such a system (in the case) is given by the following theorem, which was proven by two of the authors in Ref. [38]:
Now notice that for an fluid written with respect to conservation variables, the state vector can be written as the transpose of the row vector
where characterizes the degrees of freedom of our chosen system of equations. Note that we make this choice of a state vector for the sake of flexibility of representation and implementation (see for example ), where the strict degrees of freedom of the system (Equation 1)( ?), due to the multiplicity of (Equation 1) in the conservation form of ( ?), is just . Nevertheless, consistent with our choice of an state vector , we obtain that the inviscid flux vector satisfies
while the viscous flux is given by
In this notation (Equation 1)( ?) can be expressed as
where the notation corresponds to componentwise derivations with respect to the variable .
The Jacobian matrix of the inviscid flux can be written as the matrix:
where is the identity matrix. An important feature of the barotropic pressure law (Equation 2) is that it is not a homogeneous function in , and thus the Jacobian is not formulated to satisfy . This contrasts, for example, with the compressible NavierStokes equations when using the monofluid form of the ideal gas law (see Ref. [22]). It should be noted that some numerical fluxes and schemes are designed or derived by specifically exploiting this homogeneity with respect to the Jacobian matrix of the flux function (for example, see the Vijayasundaram flux as used in Ref. [21]). Nevertheless, our numerical fluxes will be defined independently of the homogeneity property of the corresponding map, where simply satisfies .
For the viscous flux we define the matrix,
where here and below the ’s are zero matrices of the appropriate sizes. Clearly then (Equation 4) satisfies
Further let us introduce the auxilliary function such that we are concerned with solving the coupled system:
The above equations comprise a first order system which can be effectively discretized using the DG method.
Consider the following discretization scheme motivated by Ref. [22] (and illustrated in the one dimensional case in Figure ?). Take an open with boundary , given such that for the closure of . Let denote the partition of the closure , such that taking provides the partition
comprised of elements such that . The mesh diameter is given by such that a discrete approximation to is given by the set . Each element of the partition has a boundary set given by , where elements sharing a boundary point are characterized as neighbors and generate the set of interfaces between neighboring elements. The boundary is characterized in the mesh as and indexed by elements such that . Now for define the indexing set is a neighbor of , and for define contains . Then for , we have and .
We define the broken Sobolev space over the partition as
Further, approximate solutions to (Equation 1)( ?) will exist in the space of discontinuous piecewise polynomial functions over restricted to , given as
for the space of degree polynomials on .
Choosing a degree set of polynomial basis functions for we can denote an approximation to the state vector at the time over , by
where the ’s are the finite element shape functions, and the ’s correspond to the nodal unknowns. Likewise the test functions are characterized by
for and the nodal values of the test function in each .
Letting be a classical solution to (Equation 8) and multiplying through by test functions and and integrating elementwise by parts yields:
Let and denote the values of on considered from the interior and the exterior of , respectivel. It should be noted that for , the restricted functions are determined up to a choice of boundary condition, which we will discuss in more detail in . Then we approximate the first term of (Equation 9) by,
the second term in (Equation 9) by the inviscid numerical flux ,
for the unit outward pointing normal; and the third term on the left in (Equation 9) by,
The viscous term in (Equation 9) integrates by parts to give:
We approximate the first term on the right in (Equation 13) using a generalized viscous flux (see Ref. [1] for a review of choices in the DG framework). We write here for the general viscous flux
while the second term is approximated by:
Finally for the second equation in (Equation 8) we expand it such that the approximate solution satisfies:
where,
given that is the generalized flux term associated with the discontinuous Galerkin method determined up to a congeries of options (please see Ref. [1] for a unified analysis), and using the approximate relations:
We note that these choices of approximations and fluxes define the values of on each element in terms of the values of on that element and adjacent elements. As we shall see later, this indicates that with an explicit timediscretization scheme, computing is a completely local procedure.
Combining (Equation 11), (Equation 12), (Equation 14), (Equation 15) and (Equation 16) and setting,
given the inner product
we define an approximate solution to (Equation 9) as functions and for all satisfying:
We find below that up to a (possibly arbitrary) choice of boundary data, these solutions are quite wellbehaved, extremely robust for arbitrary choice of fluids (we show and here, and have tested up to elsewhere) and readily extended to more complicated systems (e.g. ). The results presented in this paper utilize piecewise linear basis functions, but we have tested quadratic basis functions in our code as well.
3§3 Towards a generalized boundary treatment
Specifying arbitrary boundary data with respect to our approximate solution (Equation 17) is a delicate issue which requires a nuanced understanding of barotropic solutions and the mathematical techniques used to pose them. That is, we wish to determine the nature of an arbitrary boundary state in a way which is wellposed with respect to the system (Equation 1)( ?); which is to say, in such a way that the uniqueness of the solution is maintained.
However, practically speaking, recovering boundary data of an arbitrary nature on poses wellestablished difficulties with respect to the a priori estimates established in Ref. [38], which serve as the cornerstone to the existence and uniqueness result stated in Theorem 2.1. That is, recovering the a priori estimates on the solution is reduced, in the first step, to recovering two entropy inequlities (see for explicit forms) which serve as positive definite functionals over . However, when explicit boundary data is given, these inequalities acquire the addition of the following two unsigned boundary components, respectively (see Ref. [38] for the explicit calculation):
having the consequence of rendering the wellposedness of a formulation which spans any type of boundary data difficult to establish. Instead we offer a number of pragmatic approximate approaches that generalize the solution up to important restrictions, and then discuss some alternative approaches that are aimed at certain specialized types of settings. First we review some known results.
It has been shown by Strikwerda in Ref. [50], and Gustaf’sson and Sundstrom in Ref. [24] that incompletely parabolic systems, such as the shallow water equations and the full NavierStokes equations, may be wellposed with respect to a broad set of initialboundary data. These works additionally demonstrate the appropriate number of boundary conditions expected on incompletely parabolic systems, which differ from completely hyperbolic systems such as Euler’s equations. As the barotropic system (Equation 1)( ?) maintains a formal equivalence to the viscous shallow water equations (see for example Ref. [37]), we might expect (Equation 1)( ?) to behave as an incompletely parabolic system due to Ref. [24]. However, the dependencies of the pressure and viscosity on the mass fractions make showing this nontrivial and require a careful analysis of either incompletely parabolic systems [50], or hyperbolicparabolic systems [28].
The implementation of both incompletely parabolic and hyperbolic systems often rely upon the socalled “characteristic treatment.” In these systems we use characteristic directions to extrapolate values of the system variables on the boundary, while the others become constrained by a set of characteristic relations (see Ref. [29] for the hyperbolic regime). These types of treatments have been extended to treat the full NavierStokes equations [43], the viscous shallow water equations [30], and multifluid systems [51].
We want to consider what we will refer to here and below as characteristic type boundary solutions, which we view as a reduced hyperbolic system (as presented in Ref. [22]). We illustrate the situation for a simple one dimension case, but our analysis easily extends to the multidimensional case. To begin, suppose we have the domain in which to specify a characteristic boundary condition at the boundary point and time . Note that other one dimensional cases can easily be transformed to such a setting with a change of coordinates. We linearize our solution at with respect to a reference solution , which for our purposes represents the numerical solution at taken at a previous timestep. As an approximation, we neglect the viscous terms, resulting in:
where is a linearized approximation to the exact solution. Note that this arrives with a linear hyperbolic system. We consider the initialboundary value problem in the set equipped with the initial condition
and the boundary condition
Our goal is to choose the boundary condition in such a way that the initialboundary value problem is wellposed. To continue, we decompose into characteristic directions. That is, note that since is diagonalizable we have that , where the characteristic directions are the eigenvectors of associated to eigenvalues (see for an example). Then we can formulate the solution in the form
where the initial and boundary data, respectively, satisfy
It then follows (from Ref. [22] chapter 3, for example) that (Equation 18) can be written as initialboundary value scalar problems:
The scalar problems (Equation 21) may be solved via the method of characteristics, from which we obtain the solution,
which provides an explicit form to (Equation 19). From (Equation 22), we obtain the following conditions for the boundary data.
If , then necessarily . This is obtained by extrapolating the solution of to the boundary .
If , then may be freely chosen. However, in some situations it may be useful to choose for this case, such as an impermeable solid wall.
If , then may be freely chosen.
Note that once we have selected wellposed characteristic boundary conditions, we utilize the transformation
to determine the consistent boundary conditions for the conservation variables. It turns out that for (Equation 1)(Equation 3) we can reduce this method to that of the essential choices listed in Table ?. This corresponds with what we know of hyperbolic systems as shown in Ref. [22] and Ref. [29], with respect to the number of free and fixed conditions on the boundaries. In Table ? we also include a number of physically motivated restrictions which should be taken into account when selecting our boundary conditions.
Restrictions  Free  Fixed



Subsonic  inlet 


, 


Supersonic  inlet  and the appropriate  ,  none

supplimentary 


Subsonic  outlet  conditions associated 


with a choice of 


Supersonic  outlet  boundary data,  none 

including: 


Membrane  wall  , etc. 




Membrane  osmotic 




In addition to employing these “characteristic” solutions, we notice that the form of (Equation 18) satisfies the weak entropy solutions of Ref. [12] and Ref. [34] for hyperbolic systems. However, as we show in , even though these two types of solutions are both consistent, they do not display equivalent numerical behavior.
Nevertheless these two choices of boundary data, the characteristic and weak entropy solutions, are not ideal since (Equation 1)( ?) is not a hyperbolic system. We may alternatively consider the route of positing boundary data by a simple extension of the results of Zlotnik (see Ref. [58]) to see that the barotropic system is parabolic in the sense of Petrovskii upon addition of the “quasihydrodynamic” or “quasigasdynamic” auxilliary function (see Ref. [19]) on .
More generally, there exists a large back catalogue of results on compressible barotropic systems, many of which implement differing initialboundary data, and some of which utilize fairly exotic conditions on the boundary. For example, for barotropic inflow problems we can refer to both Ref. [27] and Ref. [36], where in both cases results from Ref. [38] are required and additional extensions are needed to move into the multiphase regime. Likewise solutions exist for free boundary barotropic problems [48], surface tension type boundaries [49], Navier boundary type conditions [7], and various Dirichlet type problems near vacuum states [40]; however, again, all of these results are only strictly satisfied for monofluidic systems, and thus require subtle analysis in order to extend them to the full multifluid regime. In many cases however, such as in Ref. [58], the extension is relatively straightforward.
4§4 Numerical test cases
We inspect two analytic test cases to verify the accuracy of the numerical method presented in and . In both cases we solve a monofluid restriction of (Equation 1)( ?) from the bifluid case (), with and in spatial dimension.
To begin, we specify the DG formulation in the bifluid case. First we define the three vectors , and such that (Equation 1)( ?) are expressed as
whereby setting in (Equation 5) and using (Equation 6) it then follows that
We can thus write a weak form of (Equation 1)( ?) in the same way as (Equation 9).
To solve the system we must first specify the inviscid flux . We test for two choices here. First we implement the local LaxFriedrich’s flux which satisfies
for the outward unit normal and the spectral radius of .
As our second choice of inviscid flux we implement a standard approximate Riemann solver, with flux satisfying:
where and are found from the eigendecomposition given in the appendix, is given by the diagonal matrix of eigenvalues – as also enumerated in the appendix – and the average is given by
Next we specify the viscous flux . Here we use a formulation similar to that presented in Ref. [4], but we adapt it to include the functional dependencies present in the viscosity. We choose
For the numerical flux we use the BassiRebay form, as shown in Ref. [4], which gives
Now we discretize in time, denoting a partition of [0,T] by
for a timestep given as , and implement the forward Euler scheme:
along with a slope limiting scheme in the conservation variables , where the van Leer and Osher MUSCL schemes (as shown in Ref. [53]) have been adopted in this paper.
Now we solve explicitly for (Equation 17). In particular, we show an explicit scheme using the Riemann flux, which is formulated to read: for every find such that
The above formulation lends itself naturally to a staggered scheme. First, given one solves step 3 for . This amounts to a simple, fast, and trivially parallelizable computation as the projection matrix to be inverted is blockdiagonal, with each block corresponding to an individual element. Second, given , one solves step 2 for . This similarly is a trivial computation as the mass matrix to be inverted is blockdiagonal. In fact, with the choice of an orthogonal polynomial basis on each element, the projection and mass matrices become diagonal.
We inspect the first of two numerical test cases. Consider the monofluid steady state case of (Equation 1)( ?), by setting the initial data to and . Clearly here the pressure reduces to unity and the viscosity to a constant . Next we set the periodic boundary condition
The exact solution shows constant solutions in the primitive variables. Our numerical simulations for Equation 25 using both approximate Riemann and LaxFriedrich’s inviscid fluxes have shown that the numerical error in the conservation variables is of the order of machine precision, showing no fluctuation about the steady state in time.
For the second of our test cases, we consider the monofluidic restriction of (Equation 1)( ?) given by taking , and with the additional relations:
Solving this system immediately yields
for , which leads to the ordinary differential equation
Setting the noting that solution is independent of time, we solve the ODE yielding: . Setting the initial data to
with the Dirichlet boundary data provided in the weak entropy sense of via,
we inspect the solution over the domain , with and Note that we enforce the weak entropy boundary conditions by setting in our discontinuous Galerkin formulation. Here we compare the exact solution to the solution obtained using the Riemann flux with the Osher slope limiter (denoted in Figure 1).
In Figure 1 we plot the error of the numerical solution corresponding to a mesh size of and a timestep size of , where it is clear that the relative error over fifty timesteps is of the order of magnitude of the resolution of the mesh. The relative error is zero across the solution at the first timestep, as expected, and remains nearly constant in the interior of the domain in both cases, while the weak entropy implementation displays fluctuations in time of the order of . These boundary fluctuations are neither monotonic nor generally increasing, but show complicated temporal perturbations at the weak entropy boundary points and are seen to weakly propagate into the interior as a function of the timestep. We have obtained similar behavior for the choices of a local LaxFriedrich’s inviscid flux and van Leer’s slope limiter. Further, numerical experiments have revealed that the error of the solution at a fixed time scales like for the choice of a backward Euler scheme, a timestep size of , and piecewise linear basis functions. For a general polynomial order and an explicit time integration scheme of order (see ), we find the error of the solution at a fixed time scales like , as expected, provided the CFL condition is satisfied.
5§5 Example: 2fluid with chemical inlet
Let us show a simple application of the system outlined in and evaluated over two distinct constituents. Consider the bifluid system,
with initial conditions:
The pressure is given by and the viscosity by for and with .
Now as in we easily recover the form
which integrates to (Equation 9). Again we solve for our system in a form equivalent to (Equation 17). We employ the local LaxFriedrich’s inviscid flux , the BassiRebay numerical flux , the usual viscous flux , and the van Leer slope limiter.
All that remains is determining the boundary states . We begin by considering the case of characteristic boundary conditions, and assume that at the boundary we have a subsonic inlet . In our determination of characteristic boundary conditions, we linearize about the state to arrive at an expression for the boundary state at timestep . Then, from Table ?, we see that is fixed by
Now, suppose we want a chemical inlet such that the first chemical constituent is characterized by an influx condition where similarly,
In order to maintain the consistency of our system, we additionally need that and . For , we solve (Equation 23) with the constraint in (Equation 28) to obtain:
where denotes the speed of sound at timestep on the boundary as defined in the appendix. Finally, at the other boundary point we set a transmissive characteristic boundary condition.