A Discontinuous Galerkin Method for Viscous Compressible Multifluids

A Discontinuous Galerkin Method for Viscous Compressible Multifluids


We present a generalized discontinuous Galerkin method for a multicomponent compressible barotropic Navier-Stokes 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 Runge-Kutta 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: Navier-Stokes; discontinuous Galerkin; Runge-Kutta; 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 particle-fluid 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 fluid-flows 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 well-posed [38] compressible barotropic system with functional viscosity depending on both the density and the mass fraction of each fluid component. It is well-known, 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 Runge-Kutta 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 l-fluid

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 well-posedness 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 component-wise 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 Navier-Stokes 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:


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 time-discretization 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 well-behaved, 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 well-posed 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 well-established 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 well-posedness 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 Navier-Stokes equations, may be well-posed with respect to a broad set of initial-boundary 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 hyperbolic-parabolic systems [28].

The implementation of both incompletely parabolic and hyperbolic systems often rely upon the so-called “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 Navier-Stokes 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 initial-boundary 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 initial-boundary value problem is well-posed. 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 initial-boundary 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 well-posed 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.

Choice of boundary conditions
Restrictions Free Fixed

Subsonic inlet


Supersonic inlet and the appropriate , none


Subsonic outlet conditions associated

with a choice of

Supersonic outlet boundary data, none


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 initial-boundary 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 Lax-Friedrich’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 Bassi-Rebay 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 block-diagonal, 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 block-diagonal. In fact, with the choice of an orthogonal polynomial basis on each element, the -projection and mass matrices become diagonal.

Figure 1:  The two graphs show the solution to () in terms of the linear Riemann flux and the Osher limiter, denoted u_{r} and \rho_{r}, versus the exact solution.
 The two graphs show the solution to () in terms of the linear Riemann flux and the Osher limiter, denoted u_{r} and \rho_{r}, versus the exact solution.
Figure 1: The two graphs show the solution to () in terms of the linear Riemann flux and the Osher limiter, denoted and , versus the exact solution.

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 Lax-Friedrich’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 Lax-Friedrich’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: 2-fluid 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 .

Figure 2: The left plot shows miscible species at t=12 given the characteristic chemical inlet conditions from () with \mathscr{C}=0.9 and on the boundary a = 0, with the first order transmissive conditions on b=x_{ne} (see Figure ). The right plot shows the same solution using the weak entropy formulation. Here we have a miscible solution of methanol and water at \vartheta=500K and initial \rho_{0} =5, u_{0}=0, and \mu_{1,0}=\mu_{2,0}=0.5.
The left plot shows miscible species at t=12 given the characteristic chemical inlet conditions from () with \mathscr{C}=0.9 and on the boundary a = 0, with the first order transmissive conditions on b=x_{ne} (see Figure ). The right plot shows the same solution using the weak entropy formulation. Here we have a miscible solution of methanol and water at \vartheta=500K and initial \rho_{0} =5, u_{0}=0, and \mu_{1,0}=\mu_{2,0}=0.5.
Figure 2: The left plot shows miscible species at given the characteristic chemical inlet conditions from () with and on the boundary , with the first order transmissive conditions on (see Figure ). The right plot shows the same solution using the weak entropy formulation. Here we have a miscible solution of methanol and water at K and initial , , and .

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 Lax-Friedrich’s inviscid flux , the Bassi-Rebay numerical flux , the usual viscous flux , and the van Leer slope limiter.

Figure 3: Here we show the difference between the first chemical constituent of the weak entropy \mu_{1,w} and characteristic \mu_{1,c} solutions shown in Figure , where \xi = \mu_{1,w}-\mu_{1,c}. In this figure, for emphasis, we show only the reduced spatial interval (0,9).
Figure 3: Here we show the difference between the first chemical constituent of the weak entropy and characteristic solutions shown in Figure , where . In this figure, for emphasis, we show only the reduced spatial interval .

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.

Figure 4: Here we have the complementary difference between species two of the weak entropy \mu_{2,w} and characteristic \mu_{2,c} solutions, where \eta = \mu_{2,w}-\mu_{2,c}
Figure 4: Here we have the complementary difference between species two of the weak entropy and characteristic solutions, where