Discontinuous Galerkin Methods for the Vlasov–Stokes System

1 Introduction

In this paper, we study a mathematical model describing sprays and aerosols, bubbly flows or suspension and sedimentation phenomena in two dimensions. It is the following coupled system of 2D incompressible steady generalized Stokes equation and a Vlasov type equation:

The Vlasov–Stokes system has been rigorously derived from the dynamics of a system of particles in a fluid by Jabin and Perthame [23]. The existence of solution with its large time behaviour of a similar system has been studied by Hamdache [18]. In this work, the authors take the unsteady incompressible Stokes equation as the fluid equation for the background medium. The main idea of the existence proof in [18] is to first regularize the main problem, followed by an application of the Schauder fixed point theorem, to arrive at a solution of the regularized problem. Finally, the regularized solution in the limiting case yields the existence of a weak solution to the original problem. One of the crucial steps in the proof is the use of non-negativity of 𝜌, which makes the generalized Stokes system into a damped Stokes problem, leading to an application of the inf-sup condition. Based on the Banach fixed point argument and using ρ ≥ 0 in a crucial way, Höfer [20] has proved the existence of a unique weak solution to problem (1.1)–(1.2). In [22], the authors proved the existence of a unique strong solution in 3D using the Banach fixed point argument. In [22], the authors consider the unsteady Stokes equation for the background medium. For some regularity result, see [13].

In this paper, we develop discontinuous Galerkin methods for both the equations (dG-dG). Our work is inspired by the works [5, 11] which dealt with the Vlasov–Poisson system. Recently, in [21], a discontinuous Galerkin method for the Vlasov–viscous Burgers system is developed. In those three papers, the authors have analysed conservation and stability properties of the discrete system and have also derived optimal error estimates with the help of a special projection operator. The Vlasov equation is a form of transport problem which conserves total mass. Out of the several numerical methods, the discontinuous Galerkin method has the property to preserve the total mass of the discrete system, and hence we propose dG for the Vlasov equation. One can also choose local discontinuous Galerkin or any conforming numerical method for the Stokes system. One of the main difficulties in proving results like existence and uniqueness as well as optimal error estimates for the semi-discrete system is the currently unavailable results demonstrating the non-negativity of the discrete local density ρ h , and hence the corresponding generalized discrete Stokes problem becomes non-coercive and indefinite. This creates serious difficulties in the analysis, and therefore, a more refined analysis is needed to derive results under smallness condition on the discretizing parameter ℎ. In the literature, there have been some numerical works for a similar type of models in [16, 17].

Our major contributions in this paper are as follows.

• Conservation properties and stability of the discrete system are analysed.

• Since the generalized discrete Stokes operator fails to satisfy the inf-sup condition, a more in-depth analysis helps to achieve the result under a smallness assumption on the discretizing parameter ℎ.

• Based on a modified generalized Stokes projection and a special projection, optimal error estimates are derived for the semi-discrete system.

• Existence of a unique pair of solutions for the discrete system is established.

• Finally, using a splitting algorithm, we perform some numerical experiments.

2 Qualitative and Quantitative Aspects of the Model Problem

Through out this paper, we use standard notation for Sobolev spaces. Denote by W m , p the L p -Sobolev space of order m ≥ 0 . Set

We further denote a special class of divergence-free (in the sense of distribution) vector fields by

Set the local macroscopic velocity 𝑉 as

We define some physical quantities, mass and momentum, respectively as

the energy and the velocity moments m k ⁢ f for k ≥ 0 , respectively, as

Throughout this manuscript, any function defined on Ω x is assumed to be periodic in the 𝑥-variable.

The following properties hold for the solution ( f , u , p ) of the Vlasov–Stokes equations (1.1)–(1.2), whose proof can be found in [18].

• (Positivity preserving) For any given non-negative initial datum f 0 , the solution 𝑓 is also non-negative, i.e., if f 0 ⁢ ( x , v ) ≥ 0 , then f ⁢ ( t , x , v ) ≥ 0 . Further, the local density is ρ ⁢ ( t , x ) ≥ 0 .

• (Mass conservation) The distribution function 𝑓 conserves mass in the sense that

  (2.1) ∫ Ω x ∫ R 2 f ⁢ ( t , x , v ) ⁢ d v ⁢ d x = ∫ Ω x ∫ R 2 f 0 ⁢ ( x , v ) ⁢ d v ⁢ d x for all ⁢ t ∈ [ 0 , T ] .

• (Total momentum conservation) The distribution function 𝑓 conserves momentum in the sense that

  (2.2) ∫ Ω x ∫ R 2 v ⁢ f ⁢ ( t , x , v ) ⁢ d v ⁢ d x = ∫ Ω x ∫ R 2 v ⁢ f 0 ⁢ ( x , v ) ⁢ d v ⁢ d x for all ⁢ t ∈ [ 0 , T ] .

• (Energy dissipation) The energy of the system dissipates, i.e.,

  (2.3) d d ⁢ t ⁢ ∫ R 2 ∫ Ω x | v | 2 ⁢ f ⁢ ( t , x , v ) ⁢ d x ⁢ d v ≤ 0 for any ⁢ f 0 ⁢ ( x , v ) ≥ 0 .

This is referred to as the energy identity. As f ≥ 0 , the third term on the left-hand side becomes non-negative, and hence

Moreover, if ∫ Ω x ∫ R 2 | v | 2 ⁢ f 0 ⁢ ( x , v ) ⁢ d v ⁢ d x < ∞ , then the above inequality yields u ∈ L 2 ⁢ ( 0 , T ; J 1 ) . A use of the Sobolev inequality yields u ∈ L 2 ⁢ ( 0 , T ; L p ) for all 2 ≤ p < ∞ .

The following result shows integrability estimates on the local density and on the momentum. As these appear as source terms in the Stokes equation, these estimates are crucial in deducing the regularity of a solution to the Stokes problem. The proof of the following result can be found in [18, Lemma 2.2, p. 56].

The following lemma gives the L ∞ estimate in time and space variable for the local density. The proof of this can be found in [19, Proposition 4.6, p. 44].

3 Discontinuous Galerkin Approximations

In this section, we design and analyse a discontinuous Galerkin scheme to approximate solutions to the continuum model (1.1)–(1.2). Note that a compactly supported (in the 𝑣 variable) initial datum results in a compactly supported (in the 𝑣 variable) solution f ⁢ ( t , x , v ) to the Vlasov equation (1.1), provided t ∈ [ 0 , T ] , where 𝑇 depends on the size of f 0 and the magnitude of the background fluid velocity 𝒖. So, if we restrict ourselves to compactly supported (in the 𝑣 variable) data, we can assume without loss of generality that there exist constants L > 0 and T > 0 such that supp f ( t , x , ⋅ ) ⊂ Ω v = : [ − L , L ] 2 for all t ∈ [ 0 , T ] and x ∈ Ω x .

Let T h x and T h v be two shape regular and quasi-uniform families of Cartesian partitions of Ω x and Ω v , respectively. Let { T h } be defined as the Cartesian product of these two partitions, i.e.,

The mesh sizes ℎ, h x and h v relative to these partitions are defined as follows:

We denote by Γ x and Γ v the set of all edges of the partitions T h x and T h v , respectively, and we set Γ = Γ x × Γ v . The sets of interior and boundary edges of the partition T h x (respectively, T h v ) are denoted by Γ x 0 (respectively, Γ v 0 ) and Γ x ∂ (respectively, Γ v ∂ ) so that Γ x = Γ x 0 ∪ Γ x ∂ (respectively, Γ v = Γ v 0 ∪ Γ v ∂ ).

Our objective is to develop a discontinuous Galerkin scheme to approximate solutions to the continuum model. We will consider the following broken polynomial spaces:

where P k ⁢ ( T ) is the space of scalar polynomials of degree at most 𝑘 in each variable. More precisely, we will approximate the background velocity 𝒖 by u h ∈ U h , the fluid pressure 𝑝 by p h ∈ P h , the distribution function 𝑓 by f h ∈ Z h .

Let n − and n + be the outward and inward unit normal vectors on the element 𝑇. The interior trace ψ − of 𝜓 and the outer trace ψ + of 𝜓 on Γ x are defined by

Following [4, 28], we set the average and jump of a scalar function 𝜓 and a vector-valued function Ψ at the edges as follows:

For 0 ≤ δ ≤ 1 , the weighted average of a vector-valued function Ψ is defined by

Note that, for a fixed edge, n − = − n + . For boundary edges, the jump and average are taken to be ⟦ ψ ⟧ = ψ n and { ψ } = ψ .

We will be considering the following discrete spaces:

We further denote W m , 2 ⁢ ( T h ) by H m ⁢ ( T h ) for m ≥ 1 . We define the following discrete norms:

for all 1 ≤ p < ∞ and ∥ w ∥ L ∞ ⁢ ( T h ) = ess ⁢ sup w ∈ T h ⁡ | w | . For F ∈ Γ x , w h ∈ X h ,

For all ( w h , q h ) ∈ U h × P h ,

Next, we recall certain function inequalities in the aforementioned discrete function spaces. These inequalities will be used in our subsequent analysis.

• Inverse inequality (see [12, Lemma 1.44, p. 26]): For any w h ∈ X h , there exists a constant C > 0 such that

  ∥ ∇ x w h ∥ 0 , T x ≤ C ⁢ h x − 1 ⁢ ∥ w h ∥ 0 , T x for all ⁢ T x ∈ T h x .

• Trace inequality (see [12, Lemma 1.46, p. 27]): If w h ∈ X h , then for all F ∈ Γ x and T x ∈ T h x , we have

  ∥ w h ∥ 0 , F ≤ C ⁢ h x − 1 2 ⁢ ∥ w h ∥ 0 , T x .

• Norm comparison (see [12, Lemma 1.50, p. 29]): Let 1 ≤ p , q ≤ ∞ and w h ∈ X h . Then, for all T x ∈ T h x ,

  (3.1) ∥ w h ∥ L p ⁢ ( T x ) ≤ C ⁢ h x 2 p − 2 q ⁢ ∥ w h ∥ L q ⁢ ( T x ) .

• A Poincaré–Friedrichs type inequality is proved in [3, Lemma 2.1, p. 744], which says that

  ∥ v h ∥ 0 , Ω x ≤ C ⁢ | | | v h | | | for all ⁢ v h ∈ H 1 ⁢ ( Ω x ) .

• Projection operators: Let k ≥ 0 . Let P x : L 2 ⁢ ( Ω x ) → X h , P x : L 2 → U h , P v : L 2 ⁢ ( Ω v ) → V h and P h : L 2 ⁢ ( Ω ) → Z h be the standard L 2 -projections onto the spaces X h , U h , V h and Z h , respectively. Note that P h = P x ⊗ P v (see [8, 1]).

  By definition, P h is stable in L 2 , and it can be further shown to be stable in all L p -norms (see [10] for details). Let 1 ≤ p ≤ ∞ . Then, for any w ∈ L p ⁢ ( Ω ) , ∥ P h ⁢ ( w ) ∥ L p ⁢ ( T h ) ≤ C ⁢ ∥ w ∥ L p ⁢ ( Ω ) .