Kinetic formulation of a 2 $\times$ 2 hyperbolic system arising in gas chromatography

TL;DR

This paper reformulates a gas chromatography system as a kinetic equation with an extra variable, develops a numerical scheme, and tests it on Riemann problems to analyze wave behaviors.

Contribution

It introduces a novel kinetic formulation for a hyperbolic system in chromatography and constructs a corresponding numerical scheme for its analysis.

Findings

The kinetic scheme effectively captures wave interactions in the system.

Numerical tests demonstrate the scheme's ability to resolve single and composite waves.

The approach provides new insights into the system's wave dynamics.

Abstract

The PSA system commonly used in the context of gas-solid chromatography is reformulated as a single kinetic equation using an additional kinetic variable. A kinetic numerical scheme is built from this new formulation and its behavior is tested on solving the Riemann problem in different configurations leading to single or composite waves.

Equations

$$\frac{\partial q_i^{\ast }}{\partial c_i}\ge 0,i\in \{1,2\}.$$

$${\partial }_t(c_1+q_1^{\ast }(c_1,c_2))+{\partial }_x(c_{1}u)$$

$${\partial }_t(c_2+q_2^{\ast }(c_1,c_2))+{\partial }_x(c_{2}u)$$

$$c_1+c_2$$

$$\left\{\begin{array}[]{ccl}\partial_{x}(u\,c)+\partial_{t}I(c)&=&0,\\ \partial_{x}u+\partial_{t}h(c)&=&0,\end{array}\right.$$

$$h(c)=q_1(c)+q_2(c)\ge 0,I(c)=c+q_1(c).$$

$$f(c)=c_2{q}_1-c_1{q}_2=q_1(c)-ch(c).$$

$$\left\{\begin{array}[]{ccl}\vspace{2mm}c(t,0)&=&c_{b}(t)\in[0,1],\quad t>0,\\ u(t,0)&=&u_{b}(t)>0,\quad t>0\\ \\ \vspace{2mm}c(0,x)&=&c_{0}(x)\in[0,1],\quad x>0,\end{array}\right.$$

$$0\le c_b,c_0\le 1,$$

$$0<\underset{\mathbb{R}}{\inf }{u}_b\le \underset{\mathbb{R}}{\sup }{u}_b<+\infty .$$

$$Q^{\prime }(c)$$

$${\partial }_x\left(u\psi (c)\right)+{\partial }_{t}Q(c)\le 0,$$

$$\displaystyle\displaystyle\chi(c,\xi)=\mathds{1}_{]0,c[}(\xi)=\left\{\begin{array}[]{cl}1&\mbox{ if }0<\xi<c,\\ 0&\mbox{else}.\end{array}\right.$$

$$\text{mbox}supp(\chi )=[0,c],{\int }_{\mathbb{R}}\chi (c,\xi )d\xi =c$$

$$\forall g\in {\mathcal{C}}^1(\mathbb{R}),{\int }_{\mathbb{R}}{g}^{\prime }(\xi )\chi (c,\xi )d\xi =g(c)-g(0).$$

$$\inf \{{\int }_0^1{S}^{\prime }(\xi )\varphi (c,\xi )d\xi ;\varphi \in L_{\xi }^1(0,1),0\le \varphi \le 1\text{mbox}and{\int }_0^1\varphi (c,\xi )d\xi =c\}.$$

$$<\varphi ,{\partial }_c\chi >=-<{\partial }_c\varphi ,\chi >$$

$$<\varphi (c,\xi ),H(c){\partial }_c\chi >=<H(c)\varphi (c,\xi ),{\partial }_c\chi >={\int }_0^{1}H(\xi )\varphi (\xi ,\xi )d\xi .$$

$$<\varphi (c,\xi ),{\partial }_{c}P>=-<{\partial }_c\varphi (c,\xi ),P>$$

$${\partial }_x(u\chi (c,\xi ))+a(\xi ){\partial }_t\chi (c,\xi )+{\partial }_t\left(h(c)\chi (c,\xi )\right)={\partial }_{\xi }m,$$

$${\partial }_x\left(u\psi (c)\right)+{\partial }_{t}Q(c)\le 0,\text{mbox}where{Q}^{\prime }(c)=H(c){\psi }^{\prime }(c)+h^{\prime }(c)\psi (c).$$

$$-\frac{1}{2}{\partial }_{\xi }\psi (c,\xi )=\chi (c,\xi ),$$

$${\partial }_x\left(u\psi (c,\xi )\right)+{\partial }_{t}Q(c,\xi )=-2m(t,x,\xi ).$$

$${\partial }_x\left(u\chi (c,\xi )\right)+{\partial }_t\left(-\frac{1}{2}{\partial }_{\xi }Q(c,\xi )\right)={\partial }_{\xi }m(t,x,\xi ).$$

$${\partial }_{c}Q(c,\xi )=H(c){\partial }_c\psi (c,\xi )+h^{\prime }(c)\psi (c,\xi ),$$

$$-\frac{1}{2}{\partial }_{\xi }{\partial }_{c}Q(c,\xi )$$

$${\int }_0^c{h}^{\prime }(y)\chi (y,\xi )dy=(h(c)-h(\xi ))\chi (c,\xi ),$$

$$-\frac{1}{2}{\partial }_{\xi }Q(c,\xi )=-\frac{1}{2}{\partial }_{\xi }Q(0,\xi )+H(\xi )\chi (c,\xi )+(h(c)-h(\xi ))\chi (c,\xi ).$$

$${\partial }_x(u\chi (c,\xi ))+H(\xi ){\partial }_t\chi (c,\xi )+{\partial }_t[(h(c)-h(\xi ))\chi (c,\xi )]={\partial }_{\xi }m.$$

$${\partial }_x(u\varphi (t,x,\xi ))+a(\xi ){\partial }_t\varphi (t,x,\xi )+{\partial }_t\left(h(c)\varphi (t,x,\xi )\right)={\partial }_{\xi }m,$$

$${\partial }_x(u\chi (c,\xi ))+{\partial }_t\left((a(\xi )+h(c))\chi (c,\xi )\right)={\partial }_{\xi }m$$

$$\forall T>0,\forall X>0,{\int }_0^T{\int }_0^X{\int }_0^{1}m(t,x,\xi )dtdxd\xi \le \alpha X+T\Vert u_b{\Vert }_{L^{\infty }(0,T)}.$$

Full text

1. 1 Introduction

2. 2 The PSA system

3. 3 Hyperbolicity and entropies

4. 4 Kinetic Formulations of PSA System

5. 4.1 Main kinetic inequality

6. 4.2 Second kinetic inequality

7. 5 A Kinetic scheme for the PSA System

8. 5.1 The Kinetic scheme

9. 5.2 $\displaystyle L^{\infty}$ estimates

10. 5.3 $\displaystyle BV$ estimates

11. 5.4 Discrete entropy inequalities

12. 6 The kinetic scheme and the Riemann Problem

13. 6.1 Exact solution of the Riemann Problem

14. 6.2 One adsorbable component and inert gas

15. 6.2.1 Contact discontinuity

16. 6.2.2 Adsorption step with the BET isotherm: combined waves

17. 6.3 Two adsorbable components with the binary Langmuir isotherm

18. 7 Conclusion

**Key words: ** entropy solution, kinetic formulation, boundary conditions, systems of conservation laws, kinetic schemes.

**MSC numbers: **35L65, 35L67, 35L03, 80M12.

1 Introduction

Since the work of P.-L. Lions, B. Perthame and E. Tadmor ([12, 13]), it is well known that multidimensional scalar conservation laws and some systems can be formulated as a kinetic equation using an additional kinetic variable. The so-called kinetic formulation of nonlinear hyperbolic systems of conservation laws reduces them to a linear equation on a nonlinear quantity related to the conservative unknowns, moreover it allows to recover all the entropy inequalities. It turns out to be a powerful tool to derive mathematical properties such that regularizing effects or compactness results and also efficient numerical schemes. The method was used by several authors who gave further examples of kinetic formulations: the system of chromatography ([10, 11]), the Shallow Water system ([16]) for instance. The objective of the present work is to apply the machinery of entropy to derive a kinetic formulation for the the so-called PSA system introduced in Section 2 and already studied by the authors from various points of view ([3, 4, 6, 7]). As a first application, we construct a kinetic numerical scheme, state some of its properties and test it by solving the Riemann problem.

The paper is organized as follows. In Section 2, we present the PSA system and in Section 3, we recall basics results of hyperbolicity and entropies. These two sections summarize the results of previous work by the same authors essential to understanding the problem addressed and contain no new result.

In Section 4, we build and analyze a kinetic formulation of System (5). In Section 5 a kinetic scheme is built using the preceding formulation. The last section is devoted to the numerical validation of the scheme based on the resolution of the Riemann Problem.

2 The PSA system

Pressure Swing Adsorption (PSA) is a technology that is used to separate some species from a gas under pressure according to these species’ molecular characteristics and affinity for an adsorbent material. PSA is used extensively in the production and purification of oxygen, nitrogen and hydrogen for industrial uses. It can also be used to separate a single gas from a mixture of gases. A typical PSA system involves a cyclic process where a number of connected vessels containing adsorbent material undergo successive pressurization and depressurization steps in order to produce a continuous stream of purified product gas.

As in previous papers by the authors on this subject, we focus on a model describing a step of the cyclic process, restricted to isothermal behavior. As in general fixed bed chromatography, each of the $\displaystyle d$ species ($\displaystyle d\geq 2$) simultaneously exists under two phases, a gaseous and movable one with concentration $\displaystyle c_{i}(t,x)$ or a solid (adsorbed) other with concentration $\displaystyle q_{i}(t,x)$, $\displaystyle 1\leq i\leq d$. Moreover it is assumed that these concentations are at equilibrium, i.e. $\displaystyle q_{i}(t,x)=q_{i}^{*}(c_{1},c_{2})$, where the so-called isotherms $\displaystyle q_{i}^{*}$ satisfy

[TABLE]

In gas chromatography, velocity variations accompany changes in gas composition, especially in the case of high concentration solute: it is known as the sorption effect. This effect is taken into account through a constraint on the pressure, assumed to be constant. The reader can refer for instance to [19]: “Fixed-Bed Adsorption of Gases : Effect of Velocity Variations on Transition Types”.

In [7], the original model, which is nothing else that material balances for two adsorbable components, is written into the dimensionless form:

[TABLE]

In this isothermal model, the constraint of constant pressure is achieved through Eq. (4). In the case where an instantaneous equilibrium is not assumed, the corresponding system was studied from both theoretical and numerical points of view by Bourdarias [1, 2].

Setting $\displaystyle c=c_{1}\in[0,1]$ (then $\displaystyle c_{2}=1-c$), $\displaystyle q_{i}(c)=q_{i}^{*}(c,1-c)$, $\displaystyle i=1,2$, and adding Eqs. (2)-(3), we obtain finally the following system for $\displaystyle x>0$ and $\displaystyle t>0$ which is written in a form ($\displaystyle x$-derivative first) justified in the next section:

[TABLE]

where

[TABLE]

With theses notations, the relations (1) read $\displaystyle q_{1}^{\prime}\geq 0\geq q_{2}^{\prime}$.

Following [19], we introduce a function which will plays a central role in the nonlinear study of the system, namely,

[TABLE]

We will also make use of following functions only depending on the isotherms [4]:

• •

$\displaystyle H(c)=1+q_{1}^{\prime}-ch^{\prime}=1+(1-c)q_{1}^{\prime}-cq_{2}^{\prime}\geq 1$,

• •

$\displaystyle G(c)=\exp g(c)$ where $\displaystyle g^{\prime}=-\displaystyle\frac{h^{\prime}}{H}$.

3 Hyperbolicity and entropies

For self contain, we recall without proofs some results exposed in [4].

As pointed out by Rouchon and al. ([18]), it is possible to analyze System (5) in terms of hyperbolic system of P.D.E. provided the time and space variables are exchanged : this is why System (5) is presented under this unususal form. In this framework, the vector state is \displaystyle U=\left(\begin{array}[]{l}u\\ m\end{array}\right) where $\displaystyle m=u\,c$ is the flow rate of the first species. In this vector, $\displaystyle u$ must be understood as $\displaystyle u\,(c_{1}+c_{2})$, that is the total flow rate. The initial-boundary value problem is then System (5) for $\displaystyle x>0$ and $\displaystyle t>0$ supplemented by the initial ($\displaystyle x=0$) and boundary data ($\displaystyle t=0$):

[TABLE]

where

[TABLE]

For this system, the first two equations of (7) correspond to the initial data and the last one to the boundary data. That is to say that the variable $\displaystyle x$ is progressive, i.e. time-like, and $\displaystyle t$ is a space-like variable. To be clear, we distinguish the physical time $\displaystyle t$ to the mathematical time or hyperbolic time $\displaystyle x$. The mathematical initial value problem is physically relevant for applications because experimenters only control $\displaystyle c_{b},u_{b}$, and $\displaystyle c_{0}$ can be viewed as an equilibrium reached before the beginning of the process.

The eigenvalues of the Jacobian matrix of the flux are $\displaystyle 0$ and $\displaystyle\lambda=\frac{H(c)}{u}$, thus the system is strictly hyperbolic as long as $\displaystyle u>0$: we will show in the last section that it is ensured for the solution of the Riemann Problem thanks to the assumption $\displaystyle\inf_{\mathbb{R}}u_{b}>0$. Moreover $\displaystyle\lambda$ is genuinely nonlinear in each domain where $\displaystyle f"\neq 0$. The Riemann invariants are $\displaystyle c$ and $\displaystyle W=u\,G(c)$ associated to the eigenvalues $\displaystyle 0$ and $\displaystyle\lambda$ respectively.

We have shown in [4] that there are two families of entropies: $\displaystyle u\,\psi(c)$ and $\displaystyle\phi(u\,G(c))$, where $\displaystyle\phi$ and $\displaystyle\psi$ are any real smooth functions. The corresponding entropy flux $\displaystyle Q(c)$ of the first family satisfies

[TABLE]

The first family is degenerate convex (in variables $\displaystyle(u,uc)$) provided $\displaystyle\psi^{\prime\prime}\geq 0$. So we seek entropy solutions which satisfy

[TABLE]

in the distribution sense. The second family is not always convex. There are only two interesting cases where this family is convex, namely $\displaystyle\pm G^{\prime\prime}(c)>0$ for all $\displaystyle c\in[0,1]$. When $\displaystyle G^{\prime\prime}>0$ and $\displaystyle\alpha>1$, we expect to have $\displaystyle\partial_{x}(u\,G(c))^{\alpha}\leq 0$ which reduces to $\displaystyle\partial_{x}(u\,G(c))\leq 0$. In the same way, if $\displaystyle G^{\prime\prime}<0$, we get $\displaystyle\partial_{x}(u\,G(c))\geq 0$.

4 Kinetic Formulations of PSA System

In this section, we consider weak solutions of System (5) and we give two kinetic formulation. This requires the knowledge of a complete family of supplementary conservation laws or more precisely the weak entropy inequalities ([17]). The first general formulation will be used later to build a kinetic scheme. The second formulation restricted to a convex assumption on isotherms is just mentionned in the last subsection.

4.1 Main kinetic inequality

To have a general kinetic formulation we use the family $\displaystyle u\,\psi(c)$, where $\displaystyle\psi^{\prime\prime}\geq 0$. Despite the fact that this family is always degenerate convex, this family has the great advantage to be convex without convex assumption on isotherms.

Let us introduce the classic function $\displaystyle\chi$:

[TABLE]

This function enjoys the following simple properties:

[TABLE]

and

[TABLE]

Since $\displaystyle c\in(0,1)$ we define $\displaystyle\chi$ only for $\displaystyle(c,\xi)\in(0,1)\times(0,1)$ and we have also $\displaystyle\chi(c,\xi)=\mathds{1}_{]\xi,1[}(c)$.

Moreover, this function satisfies the fundamental Gibbs property (also called Brenier’s Lemma): let $\displaystyle S$ be a convex function, $\displaystyle c\in\mathbb{R}$ and consider the minimization problem

[TABLE]

Y. Brenier ([15, 8])) showed that the minimization problem achieves its minimum at $\displaystyle f=\chi$ and that if $\displaystyle S$ is strictly convex the minimizer is unique.

In order to derive the kinetic formulation of (5), we need to state two technical but easy results.

Lemma 4.1

The distribution $\displaystyle\partial_{c}\chi$ satisfies: $\displaystyle\forall\phi\in C^{\infty}_{0}((0,1)^{2})$, $\displaystyle<\phi,\partial_{c}\chi>=\int_{0}^{1}\phi(\xi,\xi)\,d\xi$.

*Proof: *for all $\displaystyle\phi\in C^{\infty}_{0}((0,1)^{2})$ we have, using $\displaystyle\chi(c,\xi)=\mathds{1}_{]\xi,1[}(c)$ :

[TABLE]

$\displaystyle\square$

In the following lemma, $\displaystyle H$ is the function defined at the end of the second section.

Lemma 4.2

The function $\displaystyle P(c,\xi)=H(\xi)\,\chi(c,\xi)$ satisfies $\displaystyle\partial_{c}P=H(c)\,\partial_{c}\chi$.

*Proof: *on one hand, using Lemma 4.1,

[TABLE]

On the other hand,

[TABLE]

$\displaystyle\square$

We are now ready to give our main result:

Theorem 4.1

If $\displaystyle(u,c)$ is a weak entropy solution of System (5), then there exists a nonnegative measure $\displaystyle m(t,x,\xi)$ such that:

[TABLE]

where $\displaystyle a$ is given by $\displaystyle a(\xi)=H(\xi)-h(\xi)=1+f^{\prime}(\xi)$.

*Proof: *we begin to obtain the kinetic formulation (12) by writing entropy inequalities for all $\displaystyle\psi$ such that $\displaystyle\psi^{\prime\prime}\geq 0$:

[TABLE]

With Kruzkhov entropies $\displaystyle\psi(c,\xi)=|c-\xi|-|\xi|$ which satisfy

[TABLE]

we have a nonnegative measure $\displaystyle m$ such that:

[TABLE]

Applying $\displaystyle-\frac{1}{2}\partial_{\xi}$ on the previous equation we get:

[TABLE]

So, we have to compute $\displaystyle-\frac{1}{2}\,\partial_{\xi}\,Q(c,\xi)$. Since

[TABLE]

we have, applying $\displaystyle-\frac{1}{2}\partial_{\xi}$ and using Lemma 4.2:

[TABLE]

Notice that

[TABLE]

thus we have:

[TABLE]

We get finally

[TABLE]

This is valid with the $\displaystyle\chi$ function associated to the special entropy $\displaystyle\psi(c,\xi)=|c-\xi|-|\xi|$. Since Kruzkhov entropies generate all convex functions by convex combinations and density arguments, Theorem 4.1 holds.        $\displaystyle\square$

The conversely of previous theorem is the following one:

Theorem 4.2

If there exist a positive function $\displaystyle u$ such that $\displaystyle\ln u\in L^{\infty}$, a $\displaystyle\chi$-function $\displaystyle\phi(t,x,\xi)=\chi(c(t,x),\xi)$ for some function $\displaystyle c$ and a nonnegative measure $\displaystyle m$ such that

[TABLE]

then $\displaystyle(u,c)$ is a weak entropy solution of System (5).

*Proof: *multiplying Equality (14) by $\displaystyle\psi^{\prime}$ and integrating over $\displaystyle(0,\xi)$ we get (13). With (13) we recover easily System (5): first, using $\displaystyle\psi\equiv\pm 1$ then $\displaystyle Q^{\prime}(c)=\pm h^{\prime}(c)$ and we recover the second equation of (5), next the choice $\displaystyle\psi\equiv\pm c$ gives $\displaystyle Q^{\prime}(c)=\pm I^{\prime}(c)$ and we recover the first equation of (5).        $\displaystyle\square$

Remark 4.1

Writing Eq. (12) under the form

[TABLE]

we highlight an advection velocity $\displaystyle a(\xi)+h(c)$ which is is not purely kinetic, as in [13].

For weak solution we can bound the measure $\displaystyle m$ with respect to $\displaystyle L^{\infty}$ bound of $\displaystyle u$.

Proposition 4.1** (A priori bound for defect measure)**

If $\displaystyle(u,c)$ is a weak entropic solution of System (5) and let $\displaystyle m(t,x,\xi)$ the defect measure satisfying kinetic formulation (12), then there exists a constant $\displaystyle\alpha>0$ depending on $\displaystyle||h||_{\infty}$ such that:

[TABLE]

*Proof: *multiplying (12) by $\displaystyle S^{\prime}(\xi)$ such that $\displaystyle S(0)=0$, we get

[TABLE]

Integrating by parts the previous equality over $\displaystyle(0,T)_{t}\times(0,X)_{x}\times(0,1)_{\xi}$ we have:

[TABLE]

where $\displaystyle A(c)=\int_{0}^{c}a(\xi)\,S^{\prime}(\xi)\,d\xi+h(c)\,S(c).$ Since $\displaystyle 0\leq c\leq 1$, we have only to estimate $\displaystyle u(t,X)$. To control $\displaystyle\int_{0}^{T}u(t,X)\,dt$ we use second equation of System (5): $\displaystyle\partial_{x}u=-\partial_{t}h(c)$ and nonnegativity of $\displaystyle h$, then:

[TABLE]

With $\displaystyle S(c)=c^{2}/2$, a constant $\displaystyle\alpha$ depending only on the supremum of $\displaystyle A$ and $\displaystyle h$ on $\displaystyle(0,1)$ and the last inequality we can conclude the proof.        $\displaystyle\square$

4.2 Second kinetic inequality

A second kinetic formulation not used in this paper is briefly presented.

If $\displaystyle G^{\prime\prime}>0$, for all weak entropy solutions we have, for all $\displaystyle\phi$ such that

[TABLE]

is convex, see Section 3,

[TABLE]

Unfortunately these kind of entropy is related with convexity or not of isotherms. Furthermore we cannot expect to recover all System (5) (but if $\displaystyle G^{\prime\prime}>0$ we must have the Lax entropy condition).

In [4], we prove that is $\displaystyle w\mapsto p(w)$ is nonnegative and nondecreasing function, $\displaystyle\phi(w):=p(w)\exp(w)$ is a convex entropy.

Notice that $\displaystyle w\in\mathbb{R}$ in contrast to $\displaystyle c\in[0,1]$. Let us introduce some notations before exhibiting a new family of entropies:

[TABLE]

[TABLE]

[TABLE]

Furthermore $\displaystyle w(t,x)=\int_{\mathbb{R}}\widetilde{\chi}(w,\zeta)d\zeta.$ So, we deduce easily following second kinetic formulation.

Theorem 4.3** (Second kinetic formulation)**

If $\displaystyle G^{\prime\prime}\geq 0$, $\displaystyle(u,c)$ is a weak entropy solution of System (5), then there exists a nonnegative measure $\displaystyle\mu(t,x,\zeta)$ on $\displaystyle\mathbb{R}^{2}_{+}\times(0,1)$ such that:

[TABLE]

Furthermore we have the a priori bound for all $\displaystyle t>0,X>0$:

[TABLE]

Remark 4.2

In the case of one inert gas with ammoniac or water vapor $\displaystyle G^{\prime\prime}<0$, (17) is valid but with a non positive measure $\displaystyle\mu$.

5 A Kinetic scheme for the PSA System

In this section a kinetic scheme related to the kinetic formulation (14) is proposed . If the velocity $\displaystyle u$ is frozen then the kinetic formulation seems only related to the concentration $\displaystyle c$. Thus a kinetic scheme for scalar conservation laws can be used except that the macroscopic variable $\displaystyle c$ appears in the kinetic velocity. Then, an important step is to update the velocity. To be consistent with the PDE, the second equation of PSA System 5 is used. This scheme, presented in Section 5.1, enjoys some mathematical properties: maximum principle in Section 5.2, $\displaystyle BV$ estimates in Section 5.3 and the scheme satisfies entropy inequality in Section 5.4. Finally the scheme is tested with exact solutions of some Riemann problems and many isotherms in Section 6.

5.1 The Kinetic scheme

Let $\displaystyle T>0$ be the duration of the simulated process. The interval $\displaystyle[0,T]$ is divided in $\displaystyle N$ meshes $\displaystyle]t_{i-1/2},t_{i+1/2}[$ ($\displaystyle t_{1/2}=0$) with same length $\displaystyle\Delta t=T/N$, centered in $\displaystyle t_{i}$. At each time step a new spatial mesh $\displaystyle[x_{n},x_{n+1}]$ ($\displaystyle x_{0}=0$) with lenght $\displaystyle\Delta x^{n}$ is defined, according to some CFL type condition. The discrete unknowns $\displaystyle(u_{i}^{n},c_{i}^{n})$ with $\displaystyle n\in\mathbb{N}^{*}$ and $\displaystyle 1\leq i\leq N$ are the approximations of the velocity and the concentration, respectively, at $\displaystyle x=x_{n}$ in the temporal mesh $\displaystyle Te_{i}=]t_{i-1/2},t_{i+1/2}[$.

The initial and boundary data are taken in account setting:

[TABLE]

Let $\displaystyle I\geq 2$ besome integer. Being given $\displaystyle(u_{i}^{n},c_{i}^{n})$ for $\displaystyle i=1,\cdots,I$, we denote $\displaystyle u^{n}$ (resp. $\displaystyle c^{n}$) the piecewise function with value $\displaystyle u_{i}^{n}$ (resp. $\displaystyle c_{i}^{n}$) on $\displaystyle]t_{i-1/2},t_{i+1/2}[$ and we introduce the following transport equation related to (14) on $\displaystyle[x_{n},x_{n+1}[$ (with $\displaystyle x$ as the evolution variable):

[TABLE]

This equation is solved numerically, using a standard explicit upwind finite volume scheme. More precisely, we set, for each $\displaystyle\xi\in\mathbb{R}$, $\displaystyle\chi_{i}^{n}(\xi)=\chi(c_{i}^{n},\xi)$ and we denote $\displaystyle v_{i}^{n+1}(\xi)$ an approximation of the mean value on $\displaystyle]t_{i-1/2},t_{i+1/2}[$ of the solution $\displaystyle\phi(t,x_{n+1}^{-},\xi)$ of (19)-(20).

Thus $\displaystyle\phi_{i}^{n+1}(\xi)$ is given by the following scheme:

[TABLE]

with

[TABLE]

i.e.

[TABLE]

where $\displaystyle a^{+}=\max(a,0)$ and $\displaystyle a^{-}=-\min(a,0)$.

Notice that $\displaystyle\phi_{i}^{n+1}(\xi)$ is no longer a Gibbs equilibrium, i.e. a $\displaystyle\chi$ function. We recover such an equilibrium setting $\displaystyle c_{i}^{n+1}=\int_{\mathbb{R}}\phi_{i}^{n+1}(\xi)\,d\xi$ and thus getting $\displaystyle\chi_{i}^{n+1}$.

At the macroscopic level, integrating (19) with respect to the kinetic variable $\displaystyle\xi$ we get:

[TABLE]

where we have set, as long as $\displaystyle u_{i}^{n}>0$:

[TABLE]

Finally, we update the velocity $\displaystyle u$ applying a classical finite difference scheme to the second equation of System (5):

[TABLE]

Notice that(18)-(23)-(22)-(24) allow to compute $\displaystyle(u_{i}^{n+1},c_{i}^{n+1})$ for $\displaystyle i=1,\cdots,I-1$ only, because the sign of $\displaystyle a$ is not a priori known, thus the computation will be effective in all the columns for $\displaystyle 0\leq t\leq T=(I-M+1)\,\Delta t$ if $\displaystyle M$ is the number of spatial meshes. We have thus to extend the temporal domain and the functions $\displaystyle u_{b},\,c_{b}$ in a suitable way, $\displaystyle M$ being estimated following Remark 5.1.

In the sequel we note (KS) the kinetic scheme defined by (18)-(23)-(22)-(24).

5.2 $\displaystyle L^{\infty}$ estimates

Proposition 5.1

Assume that $\displaystyle u_{b}\geq\alpha$ for some constant $\displaystyle\alpha>0$. As long as $\displaystyle u_{i}^{n}>0$, if $\displaystyle\Delta x^{n}$ satisfies at each time step the CFL type condition

[TABLE]

thenwe have the following $\displaystyle L^{\infty}$ estimates:

[TABLE]

In (25), the $\displaystyle L^{\infty}$ norms are relative to $\displaystyle[0,1]$.

*Proof: *writing (5.1) under the form

[TABLE]

we obtain $\displaystyle\phi_{i}^{n+1}(\xi)$ as a convex combination of $\displaystyle\chi_{i-1}^{n}(\xi)$, $\displaystyle\chi_{i}^{n}(\xi)$ and $\displaystyle\chi_{i+1}^{n}(\xi)$ as soon as (25) is satisfied. Thus the first inequality of (26) holds by induction. In the same manner we can write (22) as

[TABLE]

where $\displaystyle p_{i-1/2}^{n,+}=\frac{A^{+}(c_{i}^{n})-A^{+}(c_{i-1}^{n})}{c_{i}^{n}-c_{i-1}^{n}}$ and $\displaystyle p_{i+1/2}^{n,-}=\frac{A^{-}(c_{i+1}^{n})-A^{-}(c_{i}^{n})}{c_{i+1}^{n}-c_{i}^{n}}$ satisfy:

[TABLE]

and the second inequality of (26) holds by induction since $\displaystyle\parallel a^{+}\parallel_{\infty}+\parallel a^{-}\parallel_{\infty}\leq 2\,\parallel a\parallel_{\infty}$.        $\displaystyle\square$

Remark 5.1

The previous result is not fully satisfactory because we are not currently able to give a positive lower bound for $\displaystyle u$ as with the Godunov scheme ([4]). With the assumption $\displaystyle\inf u_{b}=\alpha>0$, if we choose to update the velocity using the Riemann invariant $\displaystyle W$, that is setting

[TABLE]

we get immediately, by induction : $\displaystyle u_{i}^{n}\geq\alpha\,\frac{\inf_{[0,1]}G}{\parallel G\parallel_{\infty}}$. Then we can use the uniform CFL condition:

[TABLE]

and of course (26) holds. It is not clear a priori that this variant of the kinetic scheme is able to correctly solve the Riemann problem, but the numerical tests show that the behavior is quite satisfactory from this point of view: see Fig. 6.

5.3 $\displaystyle BV$ estimates

This scheme is a TVD scheme. Let us define the total variation (with respect to the time variable) of $\displaystyle(c_{i}^{n})_{i\in\mathbb{N}}$ by

[TABLE]

Proposition 5.2

The kinetic scheme (KS) is total variation diminishing, that is

[TABLE]

under the (CFL) condition (25).

*Proof: *we have just to show that $\displaystyle c_{i}^{n}$ may be written under an incremental form (see for instance [20, 21]). Now we have:

[TABLE]

with

[TABLE]

So it is easy to verify that under the condition (25) we have

[TABLE]

which ensures the incremental form for $\displaystyle c_{i}^{n}$ and thus concludes the proof.        $\displaystyle\square$

5.4 Discrete entropy inequalities

An important requirement for the scheme is to satisfy some entropy inequalities. This is possible with the choice of updating the velocity $\displaystyle u$.

Proposition 5.3

The discrete unknowns $\displaystyle(u_{i}^{n},c_{i}^{n})$ satisfy the following discrete entropy inequality, where $\displaystyle\psi$ is any real smooth function such that $\displaystyle\psi^{\prime\prime}\geq 0$:

[TABLE]

where $\displaystyle\phi_{i+1/2}^{n}=\phi^{+}(c_{i}^{n})-\phi^{-}(c_{i+1}^{n})$ with $\displaystyle\phi^{\pm}=\int\psi^{\prime}\,a^{\pm}$.

*Proof: *Inequality (28) arises mutiplying (5.1) by $\displaystyle\psi^{\prime}(\xi)$, integrating over $\displaystyle\mathbb{R}$ with respect to $\displaystyle\xi$ and applying Brenier’s lemma (see Section 4) thanks to the $\displaystyle L^{\infty}$ estimate (26) on $\displaystyle f_{i}^{n+1}(\xi)$.        $\displaystyle\square$

6 The kinetic scheme and the Riemann Problem

In this section our aim is to test the ability of the kinetic scheme to select the entropy solution of the Riemann problem for various choices of isotherms following [19]. In particular, the BET isotherm with one inflexion point leads to composite waves which we will see that they are properly calculated by the scheme. For self contain we recall the solution of the Riemann Problem (see [4]), moreover the hyperbolicity condition $\displaystyle u>0$ is ensured for this solution by Proposition 6.4.

6.1 Exact solution of the Riemann Problem

We consider the following Riemann problem:

[TABLE]

and we search a selfsimilar solution, i.e. : $\displaystyle c(t,x)=C(z)$, $\displaystyle u(t,x)=U(z)$ with $\displaystyle z=\frac{t}{x}>0$.

The exact solution is computed using the following results stated for instance in [4].

Proposition 6.1** (Rarefaction waves)**

*Any smooth non-constant self-similar solution $\displaystyle(C(z),U(z))$ of (31) in an open domain

$\displaystyle\Omega=\{0\leq\alpha<z<\beta\}$ where $\displaystyle f^{\prime\prime}(C(z))$ does not vanish, satisfies:*

[TABLE]

[TABLE]

In particular, $\displaystyle\frac{dC}{dz}$ has the same sign as $\displaystyle f^{\prime\prime}(C)$.

Corollary 6.1

Assume for instance that $\displaystyle 0\leq a<c^{-}<c^{+}<b\leq 1$ and $\displaystyle f^{\prime\prime}>0$ in $\displaystyle]a,b[$. Then the only smooth self-similar solution of (31)-(34) is such that :

[TABLE]

where

$\displaystyle z^{+}=\frac{H(c^{+})}{u^{+}}$, $\displaystyle z^{-}=z^{+}\,e^{-\phi(c^{+})}$ with $\displaystyle\phi(c)=\int_{c^{-}}^{c}\frac{f^{\prime\prime}(\xi)}{H(\xi)}\,d\xi$. Moreover $\displaystyle U$ is given by:

[TABLE]

where $\displaystyle u^{0}=\frac{H(c^{-})}{z^{-}}$.

Remark 6.1

It appears that $\displaystyle c$ is always monotone along a rarefaction wave but no longer $\displaystyle u$ because the sign of $\displaystyle h^{\prime}$ may change. Indeed the Riemann invariant $\displaystyle w=\ln u+g(c)$ is constant along such a wave and $\displaystyle g^{\prime}$, $\displaystyle h^{\prime}$ have opposite signs. However notice that in the case where one gas is inert, $\displaystyle u$ is monotone (see also [3]).

We are looking now for admissible shocks in the sense of Liu [14].

Proposition 6.2** ($\displaystyle\lambda-$shock waves)**

If $\displaystyle(c^{-},c^{+})$ satisfies the following admissibility condition equivalent to the Liu entropy-condition:

[TABLE]

then the Riemann problem (31)-(34) is solved by a shock wave defined as

[TABLE]

where $\displaystyle u^{0}$ and the speed $\displaystyle s$ of the shock are obtained through

[TABLE]

with

[TABLE]

Proposition 6.3

Two states $\displaystyle U^{-}$ and $\displaystyle U^{+}$ are connected by a contact discontinuity if and only if $\displaystyle c^{-}=c^{+}$ (with of course $\displaystyle u^{-}\neq u^{+}$), or $\displaystyle c^{-}\neq c^{+}$ and $\displaystyle f$ affine between $\displaystyle c^{-}$ and $\displaystyle c^{+}$.

Finally, concerning the Riemann problem, we make use of the following wave fan admissibility criterion (see [9] for instance): the fan is admissible if each one of its shocks, individually, satisfies the Liu shock admissibility criterion.

Then, in view of the previous results, we get the solution of the Riemann problem (31)-(34) for $\displaystyle c$ in a very simple way, similar to the scalar case with flux $\displaystyle f$.

**Case $\displaystyle c^{-}<c^{+}$: **

we consider the lower convex envelope $\displaystyle f_{c}$ of the function $\displaystyle f$ (see Fig. 2, left). On the subintervals where $\displaystyle f$ is strictly convex (then $\displaystyle f=f_{c}$ ) we get a rarefaction wave according to Corollary 6.1. Elsewhere we get admissible shock waves (or contact discontinuities if $\displaystyle f$ is affine).

**Case $\displaystyle c^{-}>c^{+}$: **

we use the upper convex envelope $\displaystyle f^{c}$ (see Fig. 2, right) and get rarefaction waves where $\displaystyle f$ is strictly concave.

The subsection is concluded by a new general result about the positivity of the velocity witch improves a similar result in [7]. In other words the region $\displaystyle\{0<u,\;0\leq c\leq 1\}$ is an invariant domain for Riemann Problems. Notice that the positivity of $\displaystyle u$ is mandatory to keep the system hyperbolic, and the velocity can blow up [6].

Proposition 6.4

Assume that $\displaystyle f$ has a finite number of inflexion points, then the solution of the Riemann Problem with $\displaystyle u^{+}>0$ involves a positive velocity.

*Proof: *with the previous results, the solution of the Riemann Problem consists in a finite sequence of simple waves and it remains to show that the result holds for a simple wave. In the case of a rarefaction wave we have, using (27), $\displaystyle u^{0}=u^{+}\frac{G(c^{+})}{G(c^{-})}>0$. In the case of a shock wave, we rewrite (35) as follows:

[TABLE]

We have $\displaystyle\frac{[q_{2}]}{[c]}<0<\frac{[q_{1}]}{[c]}$, thanks to (1), and $\displaystyle 0\leq c\leq 1$ thus $\displaystyle A,\,B\geq 1$ and $\displaystyle u^{0}>0$.        $\displaystyle\square$

6.2 One adsorbable component and inert gas

In this subsection we compare the exact solution of a Riemann problem with the approximation given by the kinetic scheme in the case of one active gas and one inert gas, with various isotherms. The following numerical examples show the accuracy of the scheme for contact discontinuities and composite waves.

Assume that $\displaystyle q_{2}^{*}=0$: the first component is the only active gas. The lenght of the column is $\displaystyle L=0.1$ and $\displaystyle 50$ time meshes are used.

6.2.1 Contact discontinuity

In this first test case, we use the following isotherm

[TABLE]

for which the function $\displaystyle f$ is linear: $\displaystyle f(c)=q_{1}^{*}(c)\,(1-c)=K_{1}c$. According to Prop. (6.3), the Riemann Problem in the $\displaystyle(t,x)$ plane is solved by a contact discontinuity connecting $\displaystyle(c^{-},u^{-})=(0.2,0.2)$ to $\displaystyle(c^{-},u^{0})$, with $\displaystyle u^{0}\simeq 0.10385$, followed by a contact discontinuity (due to the linearity of $\displaystyle f$) connecting $\displaystyle(c^{-},u^{0})$ to $\displaystyle(c^{+},u^{+})=(0.7,0.2)$.

In this simulation we have set $\displaystyle K_{1}=1$.

6.2.2 Adsorption step with the BET isotherm: combined waves

The so-called BET isothem, in our adimensional variables, is given by:

[TABLE]

with $\displaystyle 0<c_{s}<1$.

In this simulation we have set $\displaystyle Q=1$, $\displaystyle K=10$ and $\displaystyle\frac{1}{c_{s}}=1.3$. These choices are done in order to obtain a corresponding function $\displaystyle f$ with an inflexion point more easily visible in Fig. 4 below. The Riemann Problem in the $\displaystyle(t,x)$ plane is solved by a contact discontinuity connecting $\displaystyle(c^{-},u^{-})=(0.1,1)$ to $\displaystyle(c^{-},u^{0})$, with $\displaystyle u^{0}\simeq 0.39701$, followed by a shock connecting $\displaystyle(c^{-},u^{0})$ to $\displaystyle(c^{*},u^{*})\simeq(0.41546,0.54985)$ and a rarefaction connecting $\displaystyle(c^{*},u^{*})$ to $\displaystyle(c^{+},u^{+})=(0.7,1)$.

We give below the result obtained by updating $\displaystyle u$ through the relation (27). It turns out that they are quite similar: the shock and the rarefaction are in both cases correctly computed.

6.3 Two adsorbable components with the binary Langmuir isotherm

In this subsection, we assume that the two gases are active and that the process is driven by the binary Langmuir isotherm:

[TABLE]

The following numerical examples show the accuracy of the scheme for contact shock and rarefaction waves.

The lenght of the column is $\displaystyle L=0.1$ and we used $\displaystyle 50$ time meshes. In this simulation, we have set $\displaystyle Q_{1}=Q_{2}=1$, $\displaystyle K_{1}=10$ and $\displaystyle K_{2}=30$ : with these values we get a concave function $\displaystyle f$ (see Fig. 7 below).

The first case, with $\displaystyle c^{-}<c^{+}$ (adsorption step) is solved by a contact discontinuity connecting $\displaystyle(c^{-},u^{-})=(0.2,0.2)$ to $\displaystyle(c^{-},u^{0})$, with $\displaystyle u^{0}\simeq 0.19715$, followed by a shock connecting $\displaystyle(c^{-},u^{0})$ to $\displaystyle(c^{+},u^{+})=(0.7,0.2)$.

The second case, with $\displaystyle c^{-}>c^{+}$ (desorption step) is solved by a contact discontinuity connecting $\displaystyle(c^{-},u^{-})=(0.7,0.2)$ to $\displaystyle(c^{-},u^{0})$, with $\displaystyle u^{0}\simeq 0.0.20289$, followed by a rarefaction connecting $\displaystyle(c^{-},u^{0})$ to $\displaystyle(c^{+},u^{+})=(0.2,0.2)$.

7 Conclusion

We have presented a kinetic formulation of the PSA system, written in an adimensionnal form, which is used in the context of chemical engineering. This formulation, using an additional real variable, consits in a single equation which contains, in some sense, the whole system of two equations and all the entropy inequalities. As a first application, we have built a kinetic scheme, easy to implement and enjoying good properties (positivity and entropy inequality). It has been tested on the resolution of the Riemann problem in various configurations, including the case of an isotherm with at least one inflexion point, as the Langmuir isotherm, leading to composite waves. The good agreement with the analytical solution is an argument for convergence and entropic character of the scheme.