# Asymptotic behavior of a diffusive scheme solving the inviscid one-dimensional pressureless gases system

Laurent Boudin, Julien Mathiaud

<figure><img src="image_1.png" /></figure>

## To cite this version:

Laurent Boudin, Julien Mathiaud. Asymptotic behavior of a diffusive scheme solving the inviscid one-dimensional pressureless gases system. HYP'2012 - Fourteenth International Conference on Hyperbolic Problems, Jun 2012, Padova, Italy. pp.1066. hal-00765620

HAL Id: hal-00765620

https://inria.hal.science/hal-00765620

Submitted on 15 Dec 2012

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# ASYMPTOTIC BEHAVIOR OF A DIFFUSIVE SCHEME
SOLVING THE INVISCID ONE-DIMENSIONAL
PRESSURELESS GASES SYSTEM

LAURENT BOUDIN

UPMC Univ Paris 06, UMR 7598 LJLL,
F-75005 Paris, France
& Inria Paris-Rocquencourt, F-78150 Le Chesnay Cedex, France
F-78150 Le Chesnay Cedex, France

JULIEN MATHIAUD

CEA, DAM, Cesta,
BP2, F-33114 Le Barp, France
& CMLA, ENS Cachan, CNRS, UniverSud, 61 Avenue du Président Wilson,
F-94230 Cachan, France

(Communicated by the associate editor name)

ABSTRACT. In this work, we discuss some numerical properties of the viscous numerical scheme introduced in [8] to solve the one-dimensional pressureless gases system, and study in particular, from a computational viewpoint, its asymptotic behavior when the viscosity parameter $\epsilon > 0$ used in the scheme becomes smaller.

1. **Introduction.** In this work, we focus on the one-dimensional system describing a pressureless gas, which writes, for any $T > 0$,

$$ \partial_t \rho + \partial_x (\rho u) = 0, \qquad (1) $$

$$ \partial_t q + \partial_x (q u) = 0, \qquad (2) $$

in $(0, T] \times \mathbb{R}$, where $\rho(t, x) \ge 0$ is the gas density and $q(t, x) \in \mathbb{R}$ is the momentum at time $t \in [0, T]$ and location $x \in \mathbb{R}$. The gas velocity $u(t, x) \in \mathbb{R}$ must be somehow defined as a quotient of $q$ by $\rho$, but there is trouble when $\rho$ is 0. That is why one needs the notion of duality solutions [5], previously introduced for conservation laws [4]. Clearly, (1) and (2) can be seen as conservation laws, for mass and momentum respectively. System (1)-(2) is supplemented with initial conditions

$$ \rho(0, \cdot) = \rho^{\text{in}}, \quad q(0, \cdot) = q^{\text{in}}. \qquad (3) $$

This system arises from very various physical situations (cold plasmas [3], astro-physics [17, 9], traffic models [1, 13]...) and has been mathematically studied in numerous articles, such as [14, 10, 16, 5].

We here choose a periodic framework: we focus on $[0, 1]$ and impose that mass density, velocity and momentum have the same values at both $x = 0$ and $x = 1$, so that the solutions are 1-periodic.

<footnote n="2000">2000 Mathematics Subject Classification. Primary: 65M12, 58F17; Secondary: 53C35.</footnote>

<footnote>Key words and phrases. Pressureless gases, vanishing viscosity, diffusive scheme, regularity loss, asymptotic behavior.</footnote>
When $\rho$ and $u$ are smooth, and if $\rho$ remains non zero, (2) can be modified into the standard Burgers equation thank to (1):

$$ \partial_t u + \partial_x \left( \frac{u^2}{2} \right) = \partial_t u + u \partial_x u = 0. \qquad (4) $$

The mass density $\rho$ then solves a basic transport equation, where $u$ is given: it does not depend on $\rho$ because of (4). But it is common knowledge that, in finite time, a mass concentration phenomenon can happen, for instance, when $u$ does not increase. Consequently, the regularity properties of $u$ and $\rho$ are lost, and $u$ does not solve (4) anymore.

From the numerical viewpoint, one can think about several ways to discretize (1)–(3). Kinetic schemes [3, 6] or particle methods [12] allow to use the kinetic framework underlying the pressureless gas dynamics. It is also natural to try numerical schemes related to hyperbolic conservation laws [15] or relaxation schemes [2].

In [8], we prove that upwind schemes were not an option, since they failed to ensure the discrete one-sided Lipschitz (OSL) condition introduced [11] for convex scalar conservation laws. Then, following the strategy of [7] at the continuous level, we add an artificial viscosity. The new diffusive scheme we obtain is proven to be $L^\infty$-stable and consistent, hence converging towards the solution of the viscous pressureless gases system. In particular, it satisfies the discrete OSL condition. We here investigate, at the numerical level, the asymptotic behavior of the same numerical scheme when the artificial viscosity vanishes.

## 2. Diffusive numerical scheme.
Let us recall the diffusive scheme we presented in [8]. Consider $\Delta t$, $\Delta x > 0$ such that $N = T/\Delta t \in \mathbb{N}$ and $I = 1/\Delta x \in \mathbb{N}$, and set $\lambda = \Delta t/\Delta x$. Denote $\rho_i^n$, $q_i^n$ and $u_i^n$ the approximate values of $\rho$, $q$ and $u$ at time $n\Delta t \in [0, T]$ and location $(i + 1/2)\Delta x \in [0, 1]$, for $0 \le n \le N$ and $0 \le i < I$. Of course, thanks to the periodicity property, $\rho_i^n$, $q_i^n$ and $u_i^n$ can be extended for any $i \in \mathbb{Z}$. For the sake of readability, in the previous notations, we may drop the time iteration index $n$ and replace $n + 1$ by a prime symbol “'”.

Note that the discrete OSL condition can be written as $n\lambda(u_{i+1}^n - u_i^n) \le 1$, for any $i$ and $n > 0$.

Let us describe step by step the strategy for our scheme.

### 2.1. Periodic initial data.
We choose arbitrary 1-periodic initial data $\rho^{\text{in}} \ge 0$, $u^{\text{in}} \ge 0$. Indeed, the viscous problem in [7] deals with mass density and velocity, and not momentum.

### 2.2. Regularizing initial data.
We take a fixed $\varepsilon > 0$, small enough. We may have to regularize both $u^{\text{in}}$ and $\rho^{\text{in}}$ so that $u^{\text{in}}$, $\rho^{\text{in}} \in C^1(\mathbb{R}; \mathbb{R}_+^*)$ satisfy the assumptions of Theorem 2 in [8], i.e.

$$ \rho^{\text{in}}(x) \ge C\varepsilon^{1/4}, \quad u^{\text{in}}(x) \le C, \quad (u^{\text{in}})'(x) \le \frac{C}{\sqrt{\varepsilon}}, \quad \forall x \in [0, 1], \quad (5) $$

where $C \ge 0$ is a constant not depending on $\varepsilon$. Note that $\rho^{\text{in}}$ must lie in $\mathbb{R}_+^*$, since the continuous viscous model involves a division by $\rho$.
The following quantities

$$
\begin{align*}
U &= \max_{[0,1]} u^{\text{in}} > 0, & V &= \min_{[0,1]} u^{\text{in}} > 0, \\
A &= \max(0, \max_{[0,1]} (u^{\text{in}})') \ge 0, & R &= \min_{[0,1]} \rho^{\text{in}} > 0,
\end{align*}
$$

may depend on $\varepsilon$. More precisely, they must satisfy properties inherited from (5),
i.e.

$$
R \geq C\varepsilon^{1/4}, \quad V \leq U \leq C, \quad A \leq \frac{C}{\sqrt{\varepsilon}}, \qquad (6)
$$

where $C \ge 0$ does not depend on $\varepsilon$.

2.3. Choosing the time and space steps. The steps Δt and Δx > 0 are then chosen such that

$$
0 < \Delta x \le \frac{2V}{1+A}, \qquad (7)
$$

$$
0 < \Delta t \le \min \left( \frac{1}{4A+1}, \frac{1}{4U} \Delta x, \frac{R}{4\varepsilon(1+AT)} \Delta x^2 \right), \quad (8)
$$

where we set $\lambda = \Delta t / \Delta x$ and $\sigma = \Delta t / \Delta x^2$.

2.4. **Writing the scheme.** We eventually write the discretization of the viscous pressureless gases system as

$$
u'_{i} = u_{i} - \lambda \left( \frac{u_{i}^{2}}{2} - \frac{u_{i-1}^{2}}{2} \right) + \frac{\varepsilon\sigma}{\rho_{i}}(u_{i-1} + u_{i+1} - 2u_{i}), \quad (9)
$$

$$
\rho'_{i} = (1 - \lambda u'_{i}) \rho_{i} + \lambda u'_{i-1} \rho_{i-1}. \qquad (10)
$$

2.5. Properties of the scheme. In [8], we prove the following

**Theorem 2.1.** We assume that (7)–(8) hold. Then we have, for any $i$ and $n \ge 0$,

$$
V \le u_i^n \le U, \quad u_i^n - u_{i-1}^n \le \frac{A \Delta x}{1 + A n \Delta t}, \quad \rho_i^n \ge \frac{R}{1 + A n \Delta t} \ge \frac{R}{1 + A T} > 0.
$$

Moreover, the discrete total mass is conserved. Finally, when $\varepsilon > 0$ is fixed, scheme
(9)–(10) is first-order consistent (in both $t$ and $x$) with the viscous pressureless gases
system, and is monotonic.

In other words, both discrete unknown functions satisfy maximum principles, and
the discrete velocity satisfies the OSL condition. Assumptions (8) on Δt are crucial
to ensure the stability of the scheme. They are related to the standard assumptions
to get stability for explicit schemes on transport or diffusion equations.

3. **Numerical study.** The viscosity parameter $\varepsilon > 0$ is chosen in the first place. We use several values of $\varepsilon$ in the following computations: $10^{-4}$, $10^{-5}$, $10^{-6}$, $5 \times 10^{-7}$, $10^{-7}$ and the smallest value $\varepsilon_0 = 10^{-6}/14$, which will be considered as the reference situation. We take $T = 1$ s as the final time. Moreover, we focus on one case with (almost, see the discussion below) smooth initial data. It allows to point out all the problems arising the pressureless gases system, in particular the regularity loss phenomenon, already explained in Section 1.

We set, for any $x \in [0, 1]$,

$$
\rho^{\text{in}}(x) = 1 + \frac{1}{2} \cos(4\pi x), \quad u^{\text{in}}(x) = 2x(1-x) + \frac{1}{2}.
$$
We clearly have $U = 1$, $V = 1/2$, $A = 2$ and $R = 1/2$. When $\epsilon$ belongs to the set of values above, (6) is clearly satisfied. Nevertheless, $u^{\text{in}}$ is not $C^1(\mathbb{R}; \mathbb{R}_+^*)$. Indeed, there is an issue at the bounds of $[0, 1]$. Therefore, we should need to regularize $u^{\text{in}}$. Fortunately, it is not necessary: the regularization near 0 and 1 can happen on intervals whose length may be chosen smaller than $2\Delta x$. We can choose it so that $(u^{\text{in}})'$ remains in $[-2, 2]$, and consequently, $A$ still equals 2.

Let us now choose $\Delta x = 20\epsilon$. Again, if $\epsilon$ belongs to our same set of values, (7) obviously holds. Then $\Delta t$ can be taken as the optimal possible value given by (8), i.e. $\Delta t = \epsilon \min(5, 50/(1+2T)) = 5\epsilon$.

We want to study the behavior of our numerical solution with respect to $\epsilon$. In the following, since $\varrho$ has a measure meaning when $\epsilon$ goes to zero, it is really more convenient to focus on the cumulative mass $M$, defined by

$$M(t, x) = \int_{0}^{x} \varrho(t, y) dy.$$

We first draw the evolution of $M$ and $u$ with respect to $x$, at two different times, before and after approximately $t = 0.5$. Indeed, until $t = 0.5$, at the continuous level, $u$ solves the Burgers equation, and there is no regularity loss. But, near time $t = 0.5$, a regularity loss happens because of mass concentration and subsequent vacuum appearance.

<figure><img src="image_3.png" /><figcaption>FIGURE 1. Behavior, at time 0.4, of $M$ (a) on $[0, 1]$, (b) focused around $x = 0.2$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 2. Behavior, at time 0.4, of $u$ (a) on $[0, 1]$, (b) focused around $x = 0.2$</figcaption></figure>

The solutions are then displayed for various values of $\epsilon$. More precisely, Fig. 1–2 respectively give the behavior of $M$ and $u$ at time $t = 0.4$, i.e. before the regularity loss. Of course, the curves are sharpened with smaller values of $\epsilon$, but there is still no jump at all for either $\rho$ or $u$.
Fig. 3–4 then present the typical profiles of both $M$ and $u$ after the regularity loss, here given at final time $T = 1$. The jump is now visible at $x = 0.59$ on both figures, and $u$ of course still decreases at the jump abscissa.

<figure><img src="image_5.png" /><figcaption>FIGURE 3. Behavior, at final time, of $M$ (a) on $[0, 1]$, (b) focused around $x = 0.59$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Behavior, at final time, of $u$ (a) on $[0, 1]$, (b) focused around $x = 0.59$</figcaption></figure>

Let us now study the evolution of the total momentum

$$Q(t) = \int_{0}^{1} \rho(t, x) u(t, x) dx$$

with respect to time.

As the reader can see on Fig. 5, till the regularity loss, the numerical conservation of $Q$ is quite satisfactorily ensured. The property is not recovered after the loss. This drawback of our scheme (9)–(10) was already pointed out in [8]: we had to choose between the total momentum conservation (if we write a scheme for $\rho$ and $q$) and the OSL condition (for which we wrote a scheme for $\rho$ and $u$). The latter is in fact crucial to ensure that we select the right solution to the pressureless gases system. Nevertheless, we can at least state that the order of magnitude of $Q$ is conserved.

Eventually, we focus on the behavior with respect to $\varepsilon$. We focus on the following quantities

$$\mathcal{E}_M(t) = \max_{x \in [0,1]} |M(t, x) - M_0(t, x)|, \quad (11)$$

$$\mathcal{E}_u(t) = \max_{x \in [0,1]} |u(t, x) - u_0(t, x)|, \quad (12)$$

for the values $\varepsilon$ under study. Quantities indexed by 0 are of course the ones related to the reference value of $\varepsilon$, $\varepsilon_0$. Quantities on coarser grids are projected on the
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Checking the momentum for several values of $\epsilon$</figcaption></figure>

(fine) grid of the reference value ($\epsilon_0$ corresponds to 700,000 space cells) in order to compute the discretized versions of (11)–(12).

Fig. 6–7 respectively show the behavior of $\mathcal{E}_M$ and $\mathcal{E}_u$ with respect to time.

<figure><img src="image_8.png" /><figcaption>FIGURE 6. Behavior of $\mathcal{E}_M$ w.r.t. $t$ for several values of $\epsilon$</figcaption></figure>

These results give us hints of what should be the asymptotic behavior of our solutions, in particular when we loose regularity. The convergence seems very slow in $\epsilon$ but it was expected since the smallest power of $\epsilon$ in the main hypotheses of the scheme is equal to one fourth, see Eqn. 6. That partially explains that the errors are quite large, even for quite small values of $\epsilon$. In fact, we were not able to understand why $\epsilon = 10^{-7}$ gives such a different behavior (in order of magnitude) compared to not so larger values of $\epsilon$, whereas, from $10^{-4}$ to $10^{-6}$, the results are close after the regularity loss.
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Behavior of $\mathcal{E}_u$ w.r.t. $t$ for several values of $\varepsilon$</figcaption></figure>

**4. Perspectives.** As expected, the crucial point in the asymptotic behavior of our scheme is strongly linked to the regularity loss phenomenon which occurs for the pressureless gases. We should be able to prove reasonable convergence estimates with respect to $\varepsilon$ as long as the solutions remain smooth. Beyond the regularity loss, the situation remains quite unclear.

**Acknowledgments.** The authors want to thank Frédéric Lagoutière for the (on-going) discussions they share about this topic.

## REFERENCES

[1] F. Berthelin, P. Degond, M. Delitala, and M. Rascle. A model for the formation and evolution of traffic jams. *Arch. Ration. Mech. Anal.*, 187(2):185–220, 2008.

[2] C. Berthon, M. Breuss, and M.-O. Titeux. A relaxation scheme for the approximation of the pressureless Euler equations. *Numer. Methods Partial Differential Equations*, 22(2):484–505, 2006.

[3] F. Bouchut. On zero pressure gas dynamics. In *Advances in kinetic theory and computing*, volume 22 of *Ser. Adv. Math. Appl. Sci.*, pages 171–190. World Sci. Publ., River Edge, NJ, 1994.

[4] F. Bouchut and F. James. One-dimensional transport equations with discontinuous coefficients. *Nonlinear Anal.*, 32(7):891–933, 1998.

[5] F. Bouchut and F. James. Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness. *Comm. Partial Differential Equations*, 24(11-12):2173–2189, 1999.

[6] F. Bouchut, S. Jin, and X. Li. Numerical approximations of pressureless and isothermal gas dynamics. *SIAM J. Numer. Anal.*, 41(1):135–158 (electronic), 2003.

[7] L. Boudin. A solution with bounded expansion rate to the model of viscous pressureless gases. *SIAM J. Math. Anal.*, 32(1):172–193 (electronic), 2000.

[8] L. Boudin and J. Mathiaud. A numerical scheme for the one-dimensional pressureless gases system. *Numer. Methods Partial Differential Equations*, 28(6):1729–1746, 2012.

[9] Y. Brenier. A modified least action principle allowing mass concentrations for the early universe reconstruction problem. *Confluentes Math.*, 3(3):361–385, 2011.

[10] Y. Brenier and E. Grenier. Sticky particles and scalar conservation laws. *SIAM J. Numer. Anal.*, 35(6):2317–2328 (electronic), 1998.

[11] Y. Brenier and S. Osher. The discrete one-sided Lipschitz condition for convex scalar conservation laws. *SIAM J. Numer. Anal.*, 25(1):8–23, 1988.
[12] A. Chertock, A. Kurganov, and Yu. Rykov. A new sticky particle method for pressureless gas dynamics. *SIAM J. Numer. Anal.*, 45(6):2408–2441 (electronic), 2007.

[13] F. De Vuyst, V. Ricci, and F. Salvarani. Nonlocal second order vehicular traffic flow models and Lagrange-remap finite volumes. In *Finite volumes for complex applications. VI. Problems & perspectives. Volume 1, 2*, volume 4 of Springer Proc. Math., pages 781–789. Springer, Heidelberg, 2011.

[14] W. E. Yu. G. Rykov, and Ya. G. Sinai. Generalized variational principles, global weak solutions and behavior with random initial data for systems of conservation laws arising in adhesion particle dynamics. *Comm. Math. Phys.*, 177(2):349–380, 1996.

[15] L. Gosse and F. James. Numerical approximations of one-dimensional linear conservation equations with discontinuous coefficients. *Math. Comp.*, 69(231):987–1015, 2000.

[16] F. Poupaud and M. Rascle. Measure solutions to the linear multi-dimensional transport equation with non-smooth coefficients. *Comm. Partial Differential Equations*, 22(1-2):337–358, 1997.

[17] Ya.B. Zel'dovich. Gravitational instability: An approximate theory for large density perturbations. *Astron. and Astrophys.*, 5:84–89, 1970.

E-mail address: laurent.boudin@upmc.fr

E-mail address: julien.mathiaud@cea.fr