Implicit-Explicit WENO schemes for the equilibrium dispersive model of   chromatography

Rosa Donat; Francisco Guerrero anad Pep Mulet

arXiv:1705.00220·math.NA·May 2, 2017

Implicit-Explicit WENO schemes for the equilibrium dispersive model of chromatography

Rosa Donat, Francisco Guerrero anad Pep Mulet

PDF

Open Access

TL;DR

This paper develops efficient, fully conservative implicit-explicit numerical schemes for the nonlinear equilibrium dispersive model of chromatography, effectively capturing sharp chromatographic fronts in convection-dominated regimes.

Contribution

It introduces a novel IMEX scheme combining explicit and implicit discretizations for the ED model, leveraging the one-to-one relation between variables for conservation.

Findings

01

The scheme accurately captures sharp chromatographic fronts.

02

It maintains stability similar to hyperbolic problems.

03

Numerical experiments demonstrate high accuracy and efficiency.

Abstract

Chromatographic processes can be modeled by nonlinear, convection-dominated partial differential equations, together with nonlinear relations: the adsorption isotherms. In this paper we consider the nonlinear equilibrium dispersive (ED) model with adsorption isotherms of Langmuir type. We show that very efficient, fully conservative, numerical schemes can be designed for this mode by exploiting the relation between the conserved variables of the model and the physical concentrations of the multi-component mixtures. We show that this relation is one to one and admits a smooth global inverse, which cannot be given explicitly but can be easily computed by using a convenient root finder. These results provide the necessary ingredients to implement fully conservative numerical schemes for the model considered. Implicit-Explicit (IMEX) techniques can be used in the convection-dominated…

Figures16

Click any figure to enlarge with its caption.

Tables2

Table 1. Table 1: D a = 0 subscript 𝐷 𝑎 0 D_{a}=0 . Conservation of total mass for the single-component elution test.

	NCS		CS		NCS		CS
m	100	500	100	500	100	500	100	500
$T = 0.5$	0.184	0.186	0.207	0.202	0.203	0.198	0.207	0.202
$T = 1.0$	0.171	0.173	0.207	0.202	0.200	0.194	0.207	0.202
$T = 1.4$	0.166	0.166	0.207	0.202	0.197	0.192	0.207	0.202
	First order, $Δ t / Δ z = 0.9$				Second Order, $Δ t / Δ z = 0.9$

Table 2. Table 2: Observed maximum values of k / h 𝑘 ℎ k/h for different number of cells m 𝑚 m .

m	100	500	1000	2000
${(Δ t / Δ x)}_{m a x}$	0.9	0.7	0.5	0.4
$u + 2 D_{a} / Δ z$	1.1	1.5	2	3
$Δ t / (u / Δ z + 2 D_{a} / Δ z^{2})$	0.99	1.05	1	1.2

Equations172

D_{a} = \frac{Lu}{2 N _{t}},

D_{a} = \frac{Lu}{2 N _{t}},

\frac{\partial c}{\partial t} + \frac{1 - ϵ}{ϵ} \frac{\partial q}{\partial t} + u \frac{\partial c}{\partial z} = D_{a} \frac{\partial ^{2} c}{\partial ^{2} z}

\frac{\partial c}{\partial t} + \frac{1 - ϵ}{ϵ} \frac{\partial q}{\partial t} + u \frac{\partial c}{\partial z} = D_{a} \frac{\partial ^{2} c}{\partial ^{2} z}

u c - D_{a} \frac{\partial c}{\partial z}_{z = 0} = u c_{inj} (t), \frac{\partial c}{\partial z}_{z = 1} = 0,

u c - D_{a} \frac{\partial c}{\partial z}_{z = 0} = u c_{inj} (t), \frac{\partial c}{\partial z}_{z = 1} = 0,

q_{i} = \frac{a _{i} c _{i}}{1 + \sum _{i = 1}^{N} b _{i} c _{i}}, i = 1, 2, \dots, N .

q_{i} = \frac{a _{i} c _{i}}{1 + \sum _{i = 1}^{N} b _{i} c _{i}}, i = 1, 2, \dots, N .

\frac{\partial w}{\partial t} + \frac{\partial ( u c )}{\partial z} = D_{a} \frac{\partial ^{2} c}{\partial ^{2} z}, w = W (c) = c + \frac{1 - ϵ}{ϵ} q (c)

\frac{\partial w}{\partial t} + \frac{\partial ( u c )}{\partial z} = D_{a} \frac{\partial ^{2} c}{\partial ^{2} z}, w = W (c) = c + \frac{1 - ϵ}{ϵ} q (c)

(I + \frac{1 - ϵ}{ϵ} \frac{\partial q}{\partial c}) \frac{\partial c}{\partial t} + u \frac{\partial c}{\partial z} = D_{a} \frac{\partial ^{2} c}{\partial z ^{2}}

(I + \frac{1 - ϵ}{ϵ} \frac{\partial q}{\partial c}) \frac{\partial c}{\partial t} + u \frac{\partial c}{\partial z} = D_{a} \frac{\partial ^{2} c}{\partial z ^{2}}

\frac{\partial w}{\partial t} + \frac{\partial f ( w )}{\partial z} = D_{a} \frac{\partial ^{2} C ( w )}{\partial ^{2} z}, f (w) = u C (w),

\frac{\partial w}{\partial t} + \frac{\partial f ( w )}{\partial z} = D_{a} \frac{\partial ^{2} C ( w )}{\partial ^{2} z}, f (w) = u C (w),

W_{i} (c) = c_{i} (1 + \frac{η _{i}}{1 + \sum _{j = 1}^{N} b _{j} c _{j}}), η_{i} = \frac{1 - ϵ}{ϵ} a_{i}, 1 \leq i \leq N

W_{i} (c) = c_{i} (1 + \frac{η _{i}}{1 + \sum _{j = 1}^{N} b _{j} c _{j}}), η_{i} = \frac{1 - ϵ}{ϵ} a_{i}, 1 \leq i \leq N

1 < λ_{1} < d_{1} < λ_{2} < \dots < d_{N - 1} < λ_{N} < d_{N}, d_{i} = 1 + η_{i} (1 + j \sum b_{j} c_{j})^{- 1} .

1 < λ_{1} < d_{1} < λ_{2} < \dots < d_{N - 1} < λ_{N} < d_{N}, d_{i} = 1 + η_{i} (1 + j \sum b_{j} c_{j})^{- 1} .

J_{ij} (c) = \partial_{j} W_{i} (c) = δ_{i, j} (1 + \frac{η _{i}}{p ( c )}) - \frac{η _{i} b _{j} c _{i}}{p ( c ) ^{2}} .

J_{ij} (c) = \partial_{j} W_{i} (c) = δ_{i, j} (1 + \frac{η _{i}}{p ( c )}) - \frac{η _{i} b _{j} c _{i}}{p ( c ) ^{2}} .

d_{i} := 1 + η_{i} / p (c), τ_{i} := - \frac{η _{i} c _{i}}{p ( c ) ^{2}}, i = 1, \dots, N

d_{i} := 1 + η_{i} / p (c), τ_{i} := - \frac{η _{i} c _{i}}{p ( c ) ^{2}}, i = 1, \dots, N

J = D + τ b^{T}, D := diag (d_{i})_{i = 1}^{N} .

J = D + τ b^{T}, D := diag (d_{i})_{i = 1}^{N} .

S_{c} (λ) := 1 + b^{T} (D - λ I)^{- 1} τ = 1 + \frac{1}{p ( c ) ^{2}} i = 1 \sum N \frac{η _{i} b _{i} c _{i}}{λ - d _{i}} .

S_{c} (λ) := 1 + b^{T} (D - λ I)^{- 1} τ = 1 + \frac{1}{p ( c ) ^{2}} i = 1 \sum N \frac{η _{i} b _{i} c _{i}}{λ - d _{i}} .

S_{c} (\pm \infty) = 1, λ \to d_{i}^{\pm} lim S_{c} (λ) = \pm \infty, S_{c} (1) = \frac{1}{p ( c )} > 0.

S_{c} (\pm \infty) = 1, λ \to d_{i}^{\pm} lim S_{c} (λ) = \pm \infty, S_{c} (1) = \frac{1}{p ( c )} > 0.

1 < λ_{1} < d_{1} < λ_{2} < d_{2} < \cdot \cdot \cdot < d_{N - 1} < λ_{N} < d_{N} .

1 < λ_{1} < d_{1} < λ_{2} < d_{2} < \cdot \cdot \cdot < d_{N - 1} < λ_{N} < d_{N} .

Q (λ) = S_{c} (λ) Π_{i = 1}^{N} (λ - d_{i}),

Q (λ) = S_{c} (λ) Π_{i = 1}^{N} (λ - d_{i}),

x^{i} := - (D - λ_{i} I)^{- 1} τ \neq = 0.

x^{i} := - (D - λ_{i} I)^{- 1} τ \neq = 0.

0 = S_{c} (λ_{i}) = 1 + b^{T} (D - λ_{i} I)^{- 1} τ = 1 - b^{T} x^{i},

0 = S_{c} (λ_{i}) = 1 + b^{T} (D - λ_{i} I)^{- 1} τ = 1 - b^{T} x^{i},

J x^{i} = D x^{i} + τ (b^{T} x^{i}) = (D - λ_{i} I) x^{i} + λ_{i} x^{i} + τ = λ_{i} x^{i} .

J x^{i} = D x^{i} + τ (b^{T} x^{i}) = (D - λ_{i} I) x^{i} + λ_{i} x^{i} + τ = λ_{i} x^{i} .

R_{w} (y) = 1 - y + i = 1 \sum N \frac{y}{y + η _{i}} b_{i} w_{i} .

R_{w} (y) = 1 - y + i = 1 \sum N \frac{y}{y + η _{i}} b_{i} w_{i} .

y \to - η_{i}^{\pm} lim R_{w} (y) = \mp \infty, y \to \pm \infty lim R_{w} (y) = \mp \infty, R_{w} (0) = 1

y \to - η_{i}^{\pm} lim R_{w} (y) = \mp \infty, y \to \pm \infty lim R_{w} (y) = \mp \infty, R_{w} (0) = 1

- η_{N} < ρ_{N} < \cdot \cdot \cdot < - η_{1} < ρ_{1} < 0 < ρ_{0}

- η_{N} < ρ_{N} < \cdot \cdot \cdot < - η_{1} < ρ_{1} < 0 < ρ_{0}

R_{w} (p (w)) = - i \sum b_{i} w_{i} + i \sum b_{i} w_{i} \frac{p ( w )}{p ( w ) + η _{i}} = - i \sum b_{i} w_{i} \frac{η _{i}}{p ( w ) + η _{i}} < 0

R_{w} (p (w)) = - i \sum b_{i} w_{i} + i \sum b_{i} w_{i} \frac{p ( w )}{p ( w ) + η _{i}} = - i \sum b_{i} w_{i} \frac{η _{i}}{p ( w ) + η _{i}} < 0

Q (y) = R_{w} (y) Π_{i = 1}^{N} (y + η_{i}),

Q (y) = R_{w} (y) Π_{i = 1}^{N} (y + η_{i}),

C_{i} (w) := \frac{w _{i}}{1 + η _{i} / ρ _{0} ( w )}

C_{i} (w) := \frac{w _{i}}{1 + η _{i} / ρ _{0} ( w )}

W_{i} (c) = c_{i} (1 + \frac{η _{i}}{p ( c )}) \equiv c_{i} = \frac{W _{i} ( c )}{1 + η _{i} / p ( c )},

W_{i} (c) = c_{i} (1 + \frac{η _{i}}{p ( c )}) \equiv c_{i} = \frac{W _{i} ( c )}{1 + η _{i} / p ( c )},

p (c) - 1 = j = 1 \sum N b_{j} c_{j} = j = 1 \sum N \frac{b _{j} W _{j} ( c )}{1 + η _{j} / p ( c )}

p (c) - 1 = j = 1 \sum N b_{j} c_{j} = j = 1 \sum N \frac{b _{j} W _{j} ( c )}{1 + η _{j} / p ( c )}

0 = 1 - p (c) + p (c) j = 1 \sum N \frac{b _{j} W _{j} ( c )}{p ( c ) + η _{j}} = R_{W (c)} (p (c)),

0 = 1 - p (c) + p (c) j = 1 \sum N \frac{b _{j} W _{j} ( c )}{p ( c ) + η _{j}} = R_{W (c)} (p (c)),

C_{i} (W (c)) = \frac{W _{i} ( c )}{1 + η _{i} / ρ _{0} ( W ( c ))} = \frac{W _{i} ( c )}{1 + η _{i} / p ( c )} = c_{i},

C_{i} (W (c)) = \frac{W _{i} ( c )}{1 + η _{i} / ρ _{0} ( W ( c ))} = \frac{W _{i} ( c )}{1 + η _{i} / p ( c )} = c_{i},

p (C (w)) = 1 + i = 1 \sum N b_{i} C_{i} (w) = 1 + i = 1 \sum N b_{i} \frac{w _{i}}{1 + η _{i} / ρ _{0} ( w )}

p (C (w)) = 1 + i = 1 \sum N b_{i} C_{i} (w) = 1 + i = 1 \sum N b_{i} \frac{w _{i}}{1 + η _{i} / ρ _{0} ( w )}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsComputational Fluid Dynamics and Aerodynamics · Numerical methods for differential equations · Fractional Differential Equations Solutions

Full text

Implicit-Explicit WENO schemes for the equilibrium dispersive model of chromatography.

Rosa Donat 111Email: [email protected], Francisco Guerrero 222Corresponding author. Email: [email protected] Tel: +34963543232 Fax:+34963543922 and Pep Mulet 333Email: [email protected] Department of Mathematics. Universitat de València. Av. Dr. Moliner, 50; 46100-Burjassot-Valencia. Spain

Abstract

Chromatographic processes can be modeled by nonlinear, convection-dominated partial differential equations, together with nonlinear relations: the adsorption isotherms. In this paper we consider the nonlinear equilibrium dispersive (ED) model with adsorption isotherms of Langmuir type. We show that very efficient, fully conservative, numerical schemes can be designed for this mode by exploiting the relation between the conserved variables of the model and the physical concentrations of the multi-component mixtures. We show that this relation is one to one and admits a smooth global inverse, which cannot be given explicitly but can be easily computed by using a convenient root finder. These results provide the necessary ingredients to implement fully conservative numerical schemes for the model considered.

Implicit-Explicit (IMEX) techniques can be used in the convection-dominated regime in order to increase the efficiency of the numerical scheme. We propose a second order IMEX scheme, combining an explicit Weighted-Essentially-non-Oscillatory discretization of the convective fluxes with an implicit treatment of the diffusive term, in order to ilustrate the numerical issues involved in the application of IMEX techniques to this model. Through a series of numerical experiments, we show that the scheme provides accurate numerical solutions which capture the sharp discontinuities present in the chromatographic fronts, with the same stability restrictions as in the purely hyperbolic case.

keywords:

Numerical methods , WENO schemes , Chromatography , Implicit-Explicit methods , Conservation laws.

1 Introduction

Chromatography is a powerful tool for the separation of complex mixtures. In liquid batch chromatography, a pulse of fluid mixture (the solute) is injected at one end of a long cylindrical column filled with a porous medium (the stationary phase), followed by a continuous flow of liquid (the mobile phase) along the column. The solute interacts with the porous medium and is distributed between the liquid and solid phases, and the components of the mixture begin to separate according to the strength of their interaction with the stationary phase. For a sufficiently long column, band profiles of single component-concentration travel along the column and it is possible to collect pure fractions of components at the outlet of the device. These tools are used for difficult separation tasks when a high purity of the product is demanded, as it is often the case in the pharmaceutical industry.

It has been long recognized that chromatographic processes can be modeled by considering non-linear, convection-dominated partial differential equations [1, 2], coupled with some algebraic relations between the concentrations of the components of the mixture in the mobile and solid phases. Under reasonable assumptions, such as negligible dispersion effects and transport resistances, these equations become systems of first order non-linear conservation laws. Understanding the mathematical theory of these systems can enlighten many of the engineering aspects [1], in particular the formation and evolution of shock waves, which are an essential ingredient in the formation of band profiles of pure components. In addition, and since analytical solutions can only be obtained in very simple situations, it is important to develop tools that are able to perform accurate numerical simulations using these models. As observed in [3], robust and reliable numerical techniques can help practitioners to reduce the need for costly trial-and-error empirical experimentation.

In this paper we concentrate on the equilibrium-dispersive (ED henceforth) model. This is an ideal model based on the following assumptions (see e.g. [4, 2])

There is a permanent equilibrium between the solid and mobile phases at all positions in the column. 2. 2.

The compressibility of the mobile phase is negligible and there is no interaction between the solvent (carrier) and the solid phase. 3. 3.

The porous medium in the column is homogeneous. Then, the adsorption process is uniform in time and axial direction. 4. 4.

There are no radial concentration gradients in the column. 5. 5.

Only axial dispersion causes band broadening. The column efficiency is characterized by an apparent axial dispersion coefficient $D_{a}$ , related to the height of the column, $L$ , the (constant) velocity of the mobile phase, $u$ , and the number of theoretical plates $N_{t}$ , see [2], through the following relation

[TABLE] 6. 6.

Any additional factor that could influence the adsorption behavior (such as the temperature) is neglected.

The mass balance equation of the ED model involves the concentrations of the $N$ components of the mixture in the mobile phase, ${\boldsymbol{c}}=(c_{1},\ldots,c_{N})^{T}$ , and the solid phase, ${\boldsymbol{q}}=(q_{1},\ldots,q_{N})^{T}$ , and takes the following form

[TABLE]

where $\epsilon$ is the total porosity of the solid phase, $t$ is the time and $z$ the axial coordinate along the column, that is normalized to have unit height, so that the top is at $z=0$ and the bottom at $z=1$ . Under the assumptions listed above, the equilibrium relationship between the solid phase and liquid phase concentrations is given by the adsorption isotherm $\boldsymbol{q}=\boldsymbol{q}(\boldsymbol{c})$ , which is usually a non-linear function [2]. Appropriate boundary conditions for this model are proposed in [2]:

[TABLE]

for a known function $c_{inj}(t)$ .

The form of the adsorption isotherm determines the mathematical structure of the solutions to the ED model. When dispersion is negligible, the model equations (1) and the algebraic relation $\boldsymbol{q}=\boldsymbol{q}(\boldsymbol{c})$ form a system of nonlinear, first order partial differential equations. The mathematical structure of the model for $N=1$ , i.e. single-component chromatographic elution, has been described in [1] for various types of adsorption isotherms.

In this paper we consider multi-component mixtures for which the adsorption isotherms are of Langmuir type, that is

[TABLE]

where $a_{i}>0$ are the Henry coefficients, and the coefficients $b_{i}>0$ quantify the nonlinearity of the isotherm. For $N=1$ and $D_{a}=0$ , the analysis of the resulting hyperbolic conservation law carried out in [1] shows that the solutions are characterized by continuous or discontinuous composition fronts that propagate along the separation unit. For $0<D_{a}<<1$ , (1) becomes a parabolic, convection dominated PDE whose solutions may display very sharp fronts. The mathematical theory for the multi-component case seems to be much less developed.

Numerical simulations involving the nonlinear system (1) require efficient numerical techniques that can accurately describe discontinuous fronts. As reported by various authors (see e.g. [4] and references therein), Finite element (FE) methods, normally used for diffusion dominated problems, often lead to numerical oscillations in convection dominated problems whose solutions display sharp gradients, and it is also well known that spurious numerical oscillations are also observed when classical finite difference schemes (FD) are used for such problems.

In [4], the ED model (1)-(3) is rewritten as

[TABLE]

and the authors propose to use a conservative discretization of the convective terms, $\partial_{z}(u\boldsymbol{c})$ , combined with a standard centered discretization of the parabolic terms, in a finite volume (FV) framework. This numerical technique relies on the understanding that there is a one to one correspondence between the variables $\boldsymbol{w}$ and the concentrations $\boldsymbol{c}$ , so that (4) becomes a system of conservation laws when $D_{a}=0$ . Then, a conservative discretization of the convective terms guarantees mass conservation for the conserved variables, $\boldsymbol{w}$ , and, as a consequence, the shock-capturing property, i.e. shocks (for $D_{a}=0$ ) or steep profiles (for $D_{a}>0$ ) in the numerical solution have the correct speed of propagation (the reader is referred to e.g. [5] for a complete description of conservative schemes for systems of conservation laws).

Numerical schemes that combine a conservative discretization of the convective terms with a standard discretization of the parabolic terms have been successfully used to compute numerical approximations to the solution of convection-dominated second order PDEs and systems (see e.g. [6, 7, 8]). However, since the function $\boldsymbol{C}(\boldsymbol{w})$ cannot be explicitly determined when $N>1$ , Javeed et al. propose to update the values of the vector ${\boldsymbol{c}}$ by solving instead the following linearized version of (4)

[TABLE]

using an upwind, flux-limited, high resolution, conservative discretization of the derivative of the convective flux $u\boldsymbol{c}$ . The approach in [4] is attractive because it incorporates modern shock-capturing numerical techniques in the computation of the convective fluxes in (5), leading to numerical solutions which are free of numerical oscillations. However, the need to update directly the vector $\boldsymbol{c}$ forces the authors to abandon the *conservative formulation * (4) of the ED model and, as a consequence, we shall see that the resulting scheme fails to be conservative, leading to wrong speeds in the propagating fronts.

One of the objectives of this paper is to show that there is a globally well-defined, one-to-one correspondence between the vector of concentrations $\boldsymbol{c}$ and the conserved variables $\boldsymbol{w}$ , so that (4) can be rewritten as follows:

[TABLE]

with $\boldsymbol{C}(\boldsymbol{w})$ a continuously differentiable function, satisfying $\boldsymbol{C}=\boldsymbol{W}^{-1}$ . We show that, although there is no explicit expression for the function $\boldsymbol{C}(\boldsymbol{w})$ for $N>1$ , the value of $\boldsymbol{C}(\boldsymbol{w})$ for any $\boldsymbol{w}\neq 0$ , $w_{i}\geq 0$ can be determined by computing the only positive root of a well defined rational function. Hence, the necessary transfer of information required by a conservative scheme can be carried out.

In addition, we show that the structure of the Jacobian matrix $\boldsymbol{W}^{\prime}(\boldsymbol{c})$ can be computed in a rather straightforward manner by using the secular equation, as in [9]. As a consequence, we show that all the eigenvalues of the Jacobian matrix $\boldsymbol{C}^{\prime}(\boldsymbol{w})$ are strictly positive and pair-wise different, which allows us to prove the strict hyperbolicity of the model when $D_{a}=0$ , and the well-posed character of (4)-(3) for $D_{a}>0$ .

These results provide the theoretical background to implement state-of-the-art, fully conservative numerical schemes for the ED model, with Langmuir-type adsorption isotherms. In this paper we propose to use a second order Implicit-Explicit Runge-Kutta (IMEX-RK) scheme, that incorporates an off-the shelf Weighted-Non-Oscillatory (WENO) discretization of the convective flux terms. IMEX-RK schemes for convection-dominated parabolic PDEs combine the efficiency inherent to an implicit treatment of the second order derivatives (since the stability restrictions on the time step are the same as the CFL restriction that holds for $D_{a}=0$ ), with the robustness associated to a non-linear, high order conservative discretization of the convective derivative, which is treated in an explicit fashion.

The paper is organized as follows: in section 2 we analyze the mathematical structure of the ED model. In particular, we prove that the inverse function $\boldsymbol{C}(\boldsymbol{w})$ is globally well defined and smooth in $\mathbb{R}_{+}^{N}:=\{\boldsymbol{y}\in\mathbb{R}^{N}\colon y_{i}\geq 0\}$ . Our analysis relies on the eigen-structure of the Jacobian matrix $\boldsymbol{W}^{\prime}(\boldsymbol{c})$ , which can be determined by rewriting it as a rank-one perturbation of a diagonal matrix. In section 3, we discuss the application of conservative schemes to the ED model (4)-(3). We show the effects of considering a non-conservative numerical scheme, versus a fully conservative one, and discuss the convenience of using implicit techniques when $D_{a}>0$ . Section 4 describes a simple second-order IMEX-WENO scheme and discusses the various issues required for its implementation in numerical simulations of chromatographic processes that fit the ED model (4)-(3). In section 5 we show some numerical experiments to test the performance of our WENO-IMEX-RK2 scheme. We close with some conclusions and perspectives for future work.

2 The mathematical structure of the Equilibrium Dispersive Model

The ED model (1)-(3) can be rewritten as (4), where $\boldsymbol{w}=\boldsymbol{W}(\boldsymbol{c})=(W_{i}(\boldsymbol{c}))_{i=1}^{N}$ with

[TABLE]

Since $1+\sum_{j}b_{j}c_{j}\geq 1$ , $\forall\boldsymbol{c}\in\mathbb{R}_{+}^{N}=[0,\infty)^{N}$ , it follows that the set $\mathbb{R}_{+}^{N}$ is in the interior of the domain of the function $\boldsymbol{W}\colon\mathbb{R}^{N}\to\mathbb{R}^{N}$ . Moreover, $\boldsymbol{W}(\mathbb{R}_{+}^{N})\subseteq\mathbb{R}_{+}^{N}$ .

We shall prove that $\boldsymbol{W}\colon\mathbb{R}_{+}^{N}\to\mathbb{R}_{+}^{N}$ , is a continuously differentiable bijection. The *local * invertibility of $\boldsymbol{W}$ will follow from the inverse function theorem, hence we first analyze the Jacobian Matrix $\boldsymbol{W}^{\prime}(\boldsymbol{c})$ . Since this matrix can be written as a rank-one perturbation of a diagonal matrix, the analysis of its structure can be carried out via the secular equation (see also [9] and references therein)

In what follows we shall assume that the components of the mixture are ordered so that $0<a_{1}<a_{2}<\ldots<a_{N}$ .

Theorem 2.1.

For any ${\boldsymbol{c}}\in(0,\infty)^{N}$ , the Jacobian matrix ${\boldsymbol{W}}^{\prime}({\boldsymbol{c}})$ is diagonalizable, with real, strictly positive, pairwise distinct, eigenvalues $\lambda_{1},\dots,\lambda_{N}$ satisfying

[TABLE]

Proof 1.

Consider the function $p:\mathbb{R}_{+}^{N}\rightarrow[1,\infty)$ , $p({\boldsymbol{c}}):=1+\sum_{j}b_{j}c_{j}=1+\boldsymbol{b}^{T}\boldsymbol{c}$ . For any $\boldsymbol{c}\in\mathbb{R}_{+}^{N}$ we can write, with the aid of Kronecker’s delta $\delta_{i,j}$

[TABLE]

Hence, if we define (dropping the explicit $\boldsymbol{c}$ dependence for simplicity)

[TABLE]

we can write

[TABLE]

For any fixed $\boldsymbol{c}\in(0,\infty)^{N}_{+}$ , the eigen-structure of such matrices can be easily determined (see [9]) by computing the roots of the rational function

[TABLE]

Notice that its poles, $d_{i}=1+\eta_{i}/p(\boldsymbol{c})>1$ satisfy $d_{1}<d_{2}<\cdots<d_{N}$ (since $0<a_{1}<a_{2}<\ldots<a_{N}$ ). We can easily check that

[TABLE]

Hence, there must be at least $N$ roots of $S_{\boldsymbol{c}}(\lambda)$ , $\lambda_{i}$ , $1\leq i\leq N$ , satisfying

[TABLE]

It is easy to see that these are the only roots of $S_{\boldsymbol{c}}(\lambda)$ , since $S_{\boldsymbol{c}}(\mu)=0$ implies $Q(\mu)=0$ for

[TABLE]

which is a polynomial of degree $N$ , with $N$ roots at most. Hence, the roots $\lambda_{i}$ , $i=1,\ldots,N$ in (8) must be all the roots of $Q(\lambda)$ and, as a consequence, all the roots of $S_{\boldsymbol{c}}(\lambda)$ .

The above argument shows that $\forall\boldsymbol{c}\in(0,\infty)^{N}$ , $S_{\boldsymbol{c}}(\lambda)$ has $N$ different, strictly positive roots. We claim that these roots are, precisely, the eigenvalues of the Jacobian matrix $J=\boldsymbol{W}^{\prime}(\boldsymbol{c})$ .

To prove the claim above, we observe that $D-\lambda_{i}I$ is invertible for any of the roots of $S_{\boldsymbol{c}}(\lambda)$ , hence, for $\tau$ in (9) we may define ( $\boldsymbol{c}\neq 0$ implies $\boldsymbol{\tau}\neq 0$ )

[TABLE]

and check that $J\boldsymbol{x}^{i}=(D+\boldsymbol{\tau}\boldsymbol{b}^{T})\boldsymbol{x}^{i}=\lambda_{i}\boldsymbol{x}^{i}$ . Observe that

[TABLE]

which implies $\boldsymbol{b}^{T}\boldsymbol{x}^{i}=1$ , hence

[TABLE]

Thus, $J(\boldsymbol{c})={\boldsymbol{W}}^{\prime}({\boldsymbol{c}})$ has $N$ strictly positive, pairwise distinct, eigenvalues and it is, therefore, diagonalizable.

Remark 2.1.

Theorem 2.1 implies that $\boldsymbol{W}^{\prime}(\boldsymbol{c})$ is non singular $\forall\boldsymbol{c}\in\mathbb{R}^{N}_{+}$ , hence the inverse function theorem guarantees the existence of a local inverse in a neighborhood of any $\boldsymbol{w}=\boldsymbol{W}(\boldsymbol{c})$ , $\boldsymbol{c}\in\mathbb{R}^{N}_{+}$ . We shall prove that, in fact, there is a globally defined inverse function $\boldsymbol{C}:\mathbb{R}^{N}_{+}\rightarrow\mathbb{R}^{N}_{+}$ . For this, we prove first the following result.

Lemma 2.1.

For any fixed $\boldsymbol{w}\in\mathbb{R}_{+}^{N}$ , the rational function $R_{\boldsymbol{w}}:\mathbb{R}\rightarrow\mathbb{R}$ ,

[TABLE]

has only one positive root, $\rho_{0}(\boldsymbol{w})$ . In addition, $1\leq\rho_{0}(\boldsymbol{w})\leq p(\boldsymbol{w})=1+\boldsymbol{b}^{T}\boldsymbol{w}$ .

Proof 2.

It is easy to see that $w_{i}>0$ , $i=1,\dots,N$ , may be assumed without loss of generality. The rational function $R_{\boldsymbol{w}}(y)$ satisfies

[TABLE]

hence it has at least $N+1$ real roots, $\rho_{0},\ldots,\rho_{N}$ , satisfying

[TABLE]

Note that for $p({\boldsymbol{w}})=1+\boldsymbol{b}^{T}{\boldsymbol{w}}$ we have

[TABLE]

while $R_{\boldsymbol{w}}(1)=\sum_{i}b_{i}w_{i}/(1+\eta_{i})>0$ , hence, we must have $1\leq\rho_{0}\leq 1+\boldsymbol{b}^{T}\boldsymbol{w}$ .

Since any root of $R_{\boldsymbol{w}}(y)$ is also a root of

[TABLE]

which is a polynomial of degree $N+1$ , $\rho_{0},\ldots,\rho_{N}$ must be all the roots of $Q(y)$ and, hence, of $R_{\boldsymbol{w}}(y)$ .

Theorem 2.2.

The function ${\boldsymbol{W}}:\mathbb{R}_{+}^{N}\rightarrow\mathbb{R}_{+}^{N}$ given by (7) is invertible. The inverse function $\boldsymbol{C}:=\boldsymbol{W}^{-1}\colon\mathbb{R}_{+}^{N}\rightarrow\mathbb{R}_{+}^{N}$ is continuously differentiable in $\mathbb{R}_{+}^{N}$ and is defined by

[TABLE]

where $\rho_{0}({\boldsymbol{w}})$ is the only positive root of the rational function $R_{\boldsymbol{w}}(y)$ in (10).

Proof 3.

Let $\boldsymbol{c}=(c_{1},\ldots,c_{N})^{T}\in\mathbb{R}_{+}^{N}$ and $p({\boldsymbol{c}}):=1+\boldsymbol{b}^{T}{\boldsymbol{c}}\geq 1$ . According to (7)

[TABLE]

hence we can write

[TABLE]

or, equivalently,

[TABLE]

since $p(\boldsymbol{c})\geq 1$ , $p(\boldsymbol{c})$ must be the only positive root of $R_{\boldsymbol{W}(\boldsymbol{c})}(y)$ in (10), i.e. $p(\boldsymbol{c})=\rho_{0}(\boldsymbol{W}(\boldsymbol{c}))$ . Then, for any $\boldsymbol{c}\in\mathbb{R}_{+}^{N}$ ,

[TABLE]

that is, $\boldsymbol{C}\circ\boldsymbol{W}=\text{id}_{\mathbb{R}_{+}^{N}}$ .

Consider now $\boldsymbol{w}\in\mathbb{R}_{+}^{N}$ . Taking into account Lemma 2.1:

[TABLE]

With this result, it easily follows that $\boldsymbol{W}\circ\boldsymbol{C}=\text{id}_{\mathbb{R}_{+}^{N}}$ since

[TABLE]

It follows from the inverse function theorem that the function $\boldsymbol{C}$ is continuously differentiable in $\mathbb{R}_{+}^{N}$ .

Corollary 2.1.

The ED model (6) is well posed. For $D_{a}=0$ the system of conservation laws is strictly hyperbolic, and for any $\boldsymbol{w}$ such that $w_{i}>0$ , all the eigenvalues $\mu_{j}$ of the Jacobian matrix $uC^{\prime}(\boldsymbol{w})$ are positive, pairwise distinct, and bounded above by $u$ . With the notation in Theorem 2.1, the $\mu_{j}=u/\lambda_{j}$ satisfy

[TABLE]

For $D_{a}>0$ , the system is parabolic in the sense of Petrovskii (cf. [10]), i.e., the eigenvalues of the matrix $D_{a}\boldsymbol{C}^{\prime}(w)$ are bounded below by some positive constant for any $w\in\mathbb{R}_{+}^{N}$ .

Proof 4.

Theorem 2.2 allows us to write system (4) as

[TABLE]

According to Theorem 2.1, all the eigenvalues of the Jacobian matrix $\boldsymbol{W}^{\prime}(\boldsymbol{c})$ are positive and greater than 1, i.e. $\lambda_{i}>1$ for $i=1,2,...,N$ . This implies that the eigenvalues of the Jacobian matrix $\boldsymbol{C}^{\prime}(\boldsymbol{w})$ for the inverse function are:

[TABLE]

Since in the ED model the flux is given by:

[TABLE]

all the characteristic speeds of the Jacobian matrix of the physical flux satisfy the bounds in the statement from (8).

For $D_{a}>0$ , the matrix $D_{a}\boldsymbol{C}^{\prime}(\boldsymbol{w})$ appearing in the diffusion term, has eigenvalues $D_{a}\lambda_{i}(\boldsymbol{C}^{\prime}(\boldsymbol{w}))^{-1}$ from Theorem 2.1 and Theorem 2.2. From (8) these eigenvalues are bounded below by

[TABLE]

and therefore (12) is a parabolic system in the sense of Petrovskii.

Remark 2.2.

The function ${\boldsymbol{C}}({\boldsymbol{w}})$ cannot be written explicitly for $N>1$ . However, for all practical purposes, the computation of $\boldsymbol{C}({\boldsymbol{w}})$ in (11) only requires to find $\rho_{0}(\boldsymbol{w})$ , the only positive root of the rational function (10). Since $1\leq\rho_{0}(\boldsymbol{w})\leq 1+\boldsymbol{b}^{T}\boldsymbol{w}$ , this root can easily be found using a convenient root finder.

3 Conservative Numerical Schemes

As mentioned in previous sections, discontinuous fronts (for $D_{a}=0$ ) or sharp profiles (for $0<D_{a}<<1$ ) are to be expected when computing numerical approximations to (4). It is well known that such solutions are notoriously hard to compute accurately with classical numerical schemes, which tend to produce moving fronts that display oscillatory, Gibbs-like, behavior. Conservative, shock-capturing, schemes respect a fundamental property: the conservation in time of the ’total mass’ of the conserved variables ( $\int\boldsymbol{w}(x,t)dx$ in the ED model). This fact (via the Lax-Wendroff theorem, see [5] for further details), leads to numerical solutions which have the correct behavior in terms of the velocity of shocks. Combined with appropriate discretization of the diffusive terms, a high resolution conservative discretization of the convective terms can be used to ’capture’ the sharp fronts that appear in convection dominated systems, avoiding the numerical oscillations observed in more conventional schemes.

In [4], the first attempt at using modern shock capturing techniques for the model equation (1) is carried out. The starting point in [4] is the following semi-discrete scheme for the evolution of the cell-averages of $\boldsymbol{w}$ in a uniform mesh with grid points $z_{j}=(j-\frac{1}{2})\Delta z$ , $j=1,\dots,m$ , where $\Delta z=L/m$ , $L$ being the length of the column ( $z=0$ and $z=L$ represent the beginning and the end of the column, respectively)

[TABLE]

Here $\boldsymbol{w}_{j}(t)\approx\frac{1}{\Delta z}\int_{z_{j-1/2}}^{z_{j+1/2}}\boldsymbol{w}(z,t)dz$ and $\hat{\boldsymbol{f}}_{j+1/2}=\hat{\boldsymbol{f}}({\boldsymbol{w}}_{j-p},\dots,{\boldsymbol{w}}_{j+q})$ , for some function $\hat{\boldsymbol{f}}$ of $p+q+1$ arguments, is the numerical flux at the cell-boundary $z_{j+1/2}$ . As observed in [4], the simplest choice is the first order upwind numerical flux, which for the ED model becomes $\boldsymbol{\hat{f}}_{j+1/2}=u\boldsymbol{c}_{j}$ (notice that our analysis confirms that all propagation speeds are positive, hence this is indeed the upwind choice).

Considering the first order upwind numerical fluxes, centered differences for the parabolic term and the Forward Euler method for the time evolution, we get the following conservative, fully discrete scheme

[TABLE]

which is first order in space and time. Higher order conservative schemes may be obtained from (13) by using appropriate ODE solvers for the time evolution and high-resolution numerical flux functions in the discretization of the convective derivative.

The lack of an explicit expression for the relation ${\boldsymbol{c}}=\boldsymbol{C}({\boldsymbol{w}})$ prevented the authors in [4] from using directly the conservative formulation of the numerical schemes derived from (13). As an alternative, they consider the following (non-conservative) linearization of (4)

[TABLE]

which allows the authors to construct a numerical scheme which updates directly the values of $\boldsymbol{c}$ . A second-order (in space and time) scheme is implemented in [4] by using flux-limiting techniques to construct a high-resolution numerical flux function and an appropriate ode solver.

The second order scheme proposed in [4] produces the sharp, non-oscillatory numerical profiles that characterize a state-of-the-art high resolution scheme, however, it fails to be conservative and the fundamental property of conservation of the total mass no longer holds. As a consequence, the shock fronts (or the nearly discontinuous profiles, when $D_{a}<<1$ ) do not move at the correct speed.

On the other hand, according to the results of the previous section, finding $\boldsymbol{c}=\boldsymbol{C}(\boldsymbol{w})$ , only requires the computation of the positive root of a given rational function, so that the fully conservative version of these schemes can be easily implemented. To illustrate the different behavior of the conservative versus the non-conservative numerical schemes, we consider the ED model for single-component chromatographic elution:

3.1 Single-component chromatographic elution.

The model becomes (we do not need the vector notation)

[TABLE]

For the Langmuir’s isotherm:

[TABLE]

we can readily find the analytical expression for $c(w)$ :

[TABLE]

which, as observed in Theorem 2.2, corresponds to

[TABLE]

with $\rho_{0}(w)$ the unique positive solution of the rational function

[TABLE]

For $D_{a}=0$ , it is easy to see that the quantity

[TABLE]

must be conserved at the continuous level during the time evolution (after the moment when all the solute has been injected and until it begins to leave the column). At the discrete level, we measure the conservation of total mass in the numerical solution by computing the corresponding discretized magnitude:

[TABLE]

Conservative schemes for hyperbolic conservation laws and systems maintain this quantity for all time. The following test problem shows that this property no longer holds for the numerical schemes proposed in [4].

We consider the following set of parameters $D_{a}=0$ , $a=1,b=1$ , $u=1$ , $\epsilon=0.5$ and assume that the component is injected between $t=0$ and $t=0.2$ with $c=1$ at $z=0$ .

Table 1 displays the value of $W^{n}$ in (15), corresponding to three different times, obtained by following the non-conservative strategy (NCS) proposed in [4] and a conservative scheme (CS) of the same order. Figures 1 and 2 show snapshots of the numerical solution at the times shown in the table, where the side effect of the lack of conservation can be clearly observed: when using the non-conservative scheme, the speed of the shock front is slower than expected.

The experiments carried out for the single elution experiment point out that the use of small mesh parameters in the non-conservative formulation helps to mask the effect of the loss of total mass. However, reducing the spatial step-size is not an efficient alternative, due to the stability restrictions imposed by the explicit treatment of the parabolic terms. A linearized Von-Neumann stability analysis for the first order scheme (14) leads to the following stability constraint (see [6] for details)

[TABLE]

with $\varrho(\boldsymbol{C^{\prime}(w)})$ the spectral radius of $\boldsymbol{C}^{\prime}(\boldsymbol{w})$ . Since

[TABLE]

numerical stability is obtained provided that

[TABLE]

which implies a limitation on the time step $\Delta t\leq C_{1}\Delta z^{2}$ when $u\Delta z<2D_{a}$ . It is shown in [6] that similar stability restrictions apply to higher order explicit schemes for convection-diffusion equations and systems.

In Figure 3 we display the numerical solution at $T=0.5$ obtained with the second order CS applied to the single-elution ED model with $D_{a}=0.0005$ . The simulation corresponds to $m=500$ , for which the stability bound in (16) becomes (taking $C_{0}=1$ ) 0.6666. As observed in the figure, values of the ratio $\Delta t/\Delta z$ above the stability bound produce numerical oscillations that grow in time. The plot in Figure 3 should be compared with the plots shown in Figure 2, that shows the simulation corresponding to $D_{a}=0$ , for which a ratio of $\Delta t/\Delta x=0.9$ is appropriate.

The stability restrictions that result from the explicit treatment of the parabolic terms impose time steps that can be much smaller than the those required in the purely hyperbolic case. In order to avoid these strict stability restrictions we may turn to Implicit-Explicit strategies, in which the parabolic term is handled implicitly, while the convective term is discretized using any convenient high-resolution shock capturing scheme. In [6] we proved that the stability restriction for the first order upwind IMEX scheme, applied to a convection-diffusion system similar to (12) is the same as for the purely hyperbolic case.

Applying IMEX techniques to a convection-diffusion system requires to solve a system of equations at each time step, but very often these systems have a sparse structure and the cost of its solution is offset by the increased time step allowed by their stability constraints (see e.g. [6]). In the following section, we examine the specific issues that arise in the application of IMEX-WENO schemes for the ED model.

4 IMEX-WENO schemes

As observed in the previous section, Implicit-Explicit strategies represent an interesting option for the numerical simulation of chromatographic processes. IMEX-WENO schemes have been used in fairly similar scenarios [6, 11] and extensive numerical testing shows that they provide an efficient and robust tool for the numerical simulations of convection-dominated parabolic PDEs. In this section we describe a simple second order IMEX-WENO strategy in order to illustrate the numerical issues involved in IMEX numerical schemes for convection-diffusion systems, and in particular for the ED model.

For the application of this technique to the ED model, it is appropriate to rewrite (6) as follows

[TABLE]

Considering the same computational setting as before, i.e. $z_{j}=(j-1/2)\Delta z$ , the starting point in a high order WENO scheme applied to the ED model is the following semi-discrete scheme

[TABLE]

where $\boldsymbol{w}(t)$ is an $m\times N$ matrix whose $j$ -th column, $\boldsymbol{w}_{j}(t)$ , is an approximation of $w(z_{j},t)\in\mathbb{R}^{N}$ , $j=1,\dots,m$ , $\boldsymbol{{\cal{L}}}$ represents the spatial discretization of the convective term $-\frac{\partial}{\partial z}\boldsymbol{f}(\boldsymbol{w})$ and $\boldsymbol{{\boldsymbol{{\cal D}}}}$ the spatial discretization of the diffusion term $\frac{\partial}{\partial z}\boldsymbol{g}\left(\boldsymbol{w},\frac{\partial\boldsymbol{w}}{\partial z}\right)$ in (17).

The schemes we propose are Finite Difference schemes, hence, they compute numerical approximations to the point-values of the conserved variables, $\boldsymbol{w}_{j}(t_{n})\approx\boldsymbol{w}(x_{j},t_{n})$ , $t_{n}=n\Delta t$ and are characterized by a conservative discretization of the convective and diffusive terms of the form, (dropping the dependencies for simplicity)

[TABLE]

using convective and diffusive numerical fluxes, $\boldsymbol{\hat{{f}}}_{j+1/2}$ , $\boldsymbol{\hat{{g}}}_{j+1/2}$ , respectively, that approximate the respective fluxes at the corresponding cell interface $z_{j+1/2}=z_{j}+\Delta z/2$ .

The convective numerical flux

[TABLE]

is computed by finite difference WENO schemes of order $2p+1$ (nowadays an almost black-box routine) which entail applying WENO reconstructions to split convective fluxes

[TABLE]

We refer the interested reader to [12] for more details, and simply mention here that its blind application requires to specify a numerical viscosity $\alpha_{j+\frac{1}{2}}$ at each interface, which must be a local upper-bound of the size of all the characteristic speeds of the Jacobian matrix of the physical flux. For the ED model, we know that all characteristic speeds are bounded by $u$ (see Corollary 2.1), hence we may take $u$ as the numerical viscosity at each cell interface. In our numerical experiments we shall use the WENO5 numerical flux (in a component-wise fashion [12]). We recall Remark 2.2 for the computation of the convective fluxes.

The diffusive numerical fluxes are computed by second order centered finite differences

[TABLE]

The discretization of the boundary conditions (2) at $z=0=z_{1/2}$ is given by prescribing the sum of convective and diffusive numerical fluxes as follows:

[TABLE]

The term $\mathcal{L}_{1}+\mathcal{D}_{1}$ is modified accordingly:

[TABLE]

The discretization of the boundary conditions (2) at $z=1=z_{m+1/2}$ consists in taking $\boldsymbol{\hat{g}}_{m+1/2}=0=\left.g\left(\frac{\partial\boldsymbol{w}}{\partial z}\right)\right|_{z=1}$ , so the term $\mathcal{L}_{m}+\mathcal{D}_{m}$ is modified accordingly:

[TABLE]

With all the previous comments,

[TABLE]

where $\boldsymbol{w}_{j}$ is the $j$ -th column of the $N\times m$ matrix $\boldsymbol{w}$ and $\mathcal{A}$ is the tridiagonal $m\times m$ matrix given by

[TABLE]

Furthermore, the needed convective fluxes as in (19), for $j=1,\dots,m$ such that $k=j-p<1$ or $k=j+p+1>m$ , require values at the corresponding ghost cells $x_{k}$ , whose indices are, for the first case, $k=-p+1,\cdots,0$ and for the second one, $k=m+1,\dots,m+p$ . We obtain the values at those ghost cells by using extrapolation with a linear polynomial that satisfies the boundary condition and that interpolates the data for the internal point which is symmetric with respect to the boundary. For $z=0$ and $k=1-j$ , $j=1,\dots,p$ , taking into account (2), this extrapolation yields the value

[TABLE]

Notice that for $D_{a}=0$ this reduces to

[TABLE]

For $z=1$ and $k=m+j$ , $j=1,\dots,p$ , taking into account (2), this extrapolation yields the value

[TABLE]

Fully discrete, high order, schemes are obtained by using an appropriate Runge-Kutta ODE solver on (18). To obtain a fully explicit second order scheme for approximations $w^{n}_{i,j}\approx w_{i,j}(t_{n})\approx w_{i}(x_{j},t_{n})$ , $i=1,\dots,N$ , $j=1,\dots,m$ , we may use

[TABLE]

where we use the notation $\boldsymbol{{\cal D}}^{*}=\boldsymbol{{\cal D}}(\boldsymbol{w}^{*}),\boldsymbol{{\cal L}}^{*}=\boldsymbol{{\cal L}}(\boldsymbol{w}^{*})$ . As observed in the previous section, the stability requirements of this scheme would impose the use of time steps which are proportional to the spatial step-size.

On the other hand, as shown in [6], the stability restrictions for the following implicit-explicit Runge-Kutta 2 (IMEX RK2) scheme

[TABLE]

are the same as those of the purely hyperbolic case, i.e.

[TABLE]

Hence, for the ED model considered in this paper

[TABLE]

is sufficient for stability.

4.1 IMEX-WENO schemes for the ED model

The interplay between the variables $\boldsymbol{w}$ and $\boldsymbol{c}$ has to be taken into account in order to solve the nonlinear system of equations involved in an IMEX-WENO scheme such as (22)-(23). Specifically, for a single component chromatography,(22) explicitly reads as

[TABLE]

This nonlinear equation should be solved by, e.g., Newton’s method. The difficulty that the nonlinearity in $\boldsymbol{C}^{*}(\boldsymbol{w}^{n+1/2})$ is affected by the matrix $\mathcal{A}$ can be overcome by performing a change of variables $\boldsymbol{c}^{n+1/2}:=\boldsymbol{C}^{*}(\boldsymbol{w}^{n+1/2})$ , with which (25) reads now as:

[TABLE]

where $\boldsymbol{W}^{*}(\boldsymbol{c})_{i,j}=W_{i}(c_{j})$ , $j=1,\dots,m,i=1,\dots,N$ .

For simplicity in the description, we shall consider first the case of a single component ( $N=1$ ) and drop the index $n+1/2$ . In this case, (26) becomes

[TABLE]

with (cf. (20))

[TABLE]

Thus, the vector of unknowns $\boldsymbol{c}=(c_{1},\ldots,c_{m})^{T}$ satisfies the nonlinear system

[TABLE]

where the matrix $M(\boldsymbol{c})$ is defined as

[TABLE]

with $E(\boldsymbol{c})$ diagonal

[TABLE]

We may solve this nonlinear system by the standard Newton’s method

[TABLE]

Since $\mathcal{F}_{i}(c)=e(c_{i})c_{i}-\sum_{k=1}^{m}\mathcal{A}_{i,k}c_{k}$ and $(e(c)c)^{\prime}=1+\eta(1+bc)^{-2}$

[TABLE]

Hence, solving (27) for each new iterate only involves the solution of a tridiagonal system.

For the $N$ -component ED model, with a Langmuir type adsorption isotherm, we have

[TABLE]

where $i=1,\dots,N$ refers to the component of the mixture and $j=1,\dots,m$ refers to the grid point under consideration. For an $N\times m$ matrix $A$ , we denote $\mathcal{V}(A)\in\mathbb{R}^{Nm}$ given by juxtaposition of columns, i.e., $\mathcal{V}(A)_{N(j-1)+i}=A_{i,j}$ . By applying $\mathcal{V}$ to (26) and using the identity $\mathcal{V}(BXA^{T})=(A\otimes B)\mathcal{V}(X)$ , we get the vectorial equation

[TABLE]

for the unknown ${\boldsymbol{{\cal C}}}:=\mathcal{V}(\boldsymbol{c})$ , where ${\cal M}({\boldsymbol{{\cal C}}})$ is the $Nm\times Nm$ block tridiagonal matrix:

[TABLE]

with the $N\times N$ identity matrix $I$ and ${\cal E}({\boldsymbol{{\cal C}}})$ being a block-diagonal matrix

[TABLE]

with $E^{k}$ ( $k=1,\cdots,m$ ) $N\times N$ diagonal matrices with diagonal elements:

[TABLE]

As before, Newton’s method requires the computation of $\mathcal{F}^{\prime}(\boldsymbol{{\cal C}})$ in (30). From (29)

[TABLE]

Therefore $\mathcal{F}^{\prime}(\boldsymbol{{\cal C}})=\widehat{\mathcal{E}}-\mathcal{A}\otimes I$ , where $\widehat{\mathcal{E}}$ is the block diagonal matrix, whose diagonal blocs are the $N\times N$ (full) matrices

[TABLE]

Thus, in the multi-component case, finding ${\boldsymbol{{\cal C}}}$ in (30) by Newton’s method (28) involves solving a block-tridiagonal system with small blocks of size $N\times N$ at each iteration step. This can be efficiently carried out by using a standard block tridiagonal LU factorization algorithm (see [13] for details).

The second step in the IMEX-RK2 scheme (22) is explicit, so that we directly obtain ${\boldsymbol{w}}_{j}^{n+1}$ , $j=1,\ldots,m$ . Through the computation of the only positive root of $R_{\boldsymbol{w}}(p)=0$ with $\boldsymbol{w}=\boldsymbol{w}_{j}$ (e.g. by Newton’s method) we obtain $p_{j}$ from which we easily get the vector $\boldsymbol{c}_{j}^{n+1}$ at the next time step via (11).

Applying IMEX-WENO schemes to the ED model can, thus, be carried out with a moderate computational effort. For multi-component mixtures, the larger time steps allowed by the stability restrictions compensate the additional effort with respect to a fully explicit alternative.

5 Numerical experiments

In this section we perform several numerical experiments to illustrate the behavior of the proposed IMEX-WENO scheme.

First, we consider a simple experiment that shows that the IMEX alternative is able to obtain robust and reliable results when $D_{a}>0$ with a stability restriction of the type stated in (24). Second, we consider the simulation of a three-component mixture proposed in [4].

5.1 IMEX versus Explicit Stability Restrictions

We consider the same test for the single-component ED model considered in section 3.1 ( $u=a=b=\eta=1,\epsilon=0.5$ and the same initial conditions).

In Figure 4 (left) we show the numerical solution of the purely hyperbolic problem at $T=0.5$ for $\Delta t/\Delta x=0.9$ . In Figure 4 (right) we show the solution of the ED model for $D_{a}=0.005$ obtained with the second order IMEX-WENO scheme described in the previous section, also for $\Delta t/\Delta z=0.9$ and $T=0.5$ . As expected, the simulation is free of numerical artifacts and describes correctly the effect of the parabolic term with respect to the solution obtained for the purely hyperbolic model.

As observed in section 3.1, the use of fully explicit schemes when $D_{a}>0$ lead to a stability restriction that depends on $\Delta z$ . In Figure 5 we show the numerical results obtained with the explicit two-step Runge-Kutta scheme (21) for $m=500$ ( $\Delta t/\Delta x=0.9$ ) and $m=1000$ ( $\Delta t/\Delta x=0.7$ ). In both cases, the stability restriction (16) does not hold, and numerical oscillations develop. For $m=1000$ and $\Delta t/\Delta x=0.9$ the numerical oscillations are so large that the numerical solution (not shown) is not representative of any meaningful behavior.

In Table 2 we display the maximum values for the ratio $k/h$ for which no oscillatory behavior is observed in the numerical simulation obtained with the fully explicit scheme (21) (less than 3% with respect to value of a reference solution, obtained with the IMEX scheme for $m=2000$ ), and the value of the denominator in (16). The table shows that condition (16), which was developed in [6] for the upwind first order scheme, must be enforced in order to get stable numerical solutions when using any explicit scheme for the ED model.

5.2 Three component displacement chromatography

Displacement chromatography relies on the idea that one component (the displacer) has a stronger affinity to the solid phase than any of the other components in the sample mixture, hence it has the capability to displace the other components of the mixture from the stationary phase. For a sufficiently long column and appropriate adsorption isotherms, the concentrations of the components form rectangular regions of high concentration of one component in the mixture. The series of such zones are the so-called isotachic train [14].

We consider the case of a mixture of two components and one displacer proposed in [4] (section 4.3). The values of the parameters are: $a_{1}=4,a_{2}=5,a_{3}=6,b_{1}=4,b_{2}=5,b_{3}=1$ . In addition, $N_{t}=10000$ , $\epsilon=0.5$ and $u=0.2$ .

To compute the temporal evolution of the concentrations of the two components, $c_{1}$ and $c_{2}$ , and the displacer, $c_{3}$ as they move along the column, we use the proposed WENO-IMEX-RK2 scheme described in the previous section.

Experiment 1. Components 1 and 2 are injected between $t=0$ and $t=0.1$ with $c_{1}=c_{2}=1$ at $z=0$ . Component 3, the displacer, is injected from $t=0.1$ with $c_{3}=1$ .

Figure 6 shows the results for the first experiment for a CFL= $u\Delta t/\Delta z=0.8$ and $m=1000$ (left plot) or $m=200$ (right plot). The formation of the displacement train can be clearly appreciated in both simulations. The results of the simulation for $m=1000$ can be compared with those reported in [4].

Experiment 2: The concentration injected for the displacer is $c_{3}=0.5$ and the values of the rest of the parameters are the same as in experiment 1. This fact (see [4]) prevents the formation of a rectangular pulse for component 1. As we can see in figure 7, the line from the origin intersects the isotherm of the second component, but not the corresponding to the first component. This makes that only the second component can form a rectangular pulse. Figure 8 shows the time evolution obtained with the IMEX-RK2 scheme. We see that the numerical solution behaves as expected.

Experiment 3: The concentration injected for the displacer is further reduced, to the value $c_{3}=0.1$ . The values of the rest of the parameters are the same as in experiment 1.

In this case, as can be seen in Figure 7, none of the isotherms is intersected by the operating line. According to this, both components fail to form equilibrated rectangular pulses. The results are shown in Figure 8. Again, our IMEX-RK2 scheme reproduces correctly the expected behavior.

6 Conclusions

In this paper we have examined the ED model with Langmuir-type adsorption isotherms. We have proven that for $D_{a}=0$ , the model can be written as a system of conservation laws, since the relation between the conserved variables $\boldsymbol{w}=\boldsymbol{c}+\frac{1-\epsilon}{\epsilon}\boldsymbol{q}(\boldsymbol{c})$ and the physical concentrations, $\boldsymbol{c}$ , admits a smooth, globally well defined inverse, $\boldsymbol{C}(\boldsymbol{w})$ . The properties of the function $\boldsymbol{W}(\boldsymbol{c})=\boldsymbol{c}+\frac{1-\epsilon}{\epsilon}\boldsymbol{q}(\boldsymbol{c})$ , relating the conserved variables and the physical concentrations for Langmuir-type isotherms, are exploited in order to show that the system of conservation laws is strictly hyperbolic.

The inverse function $\boldsymbol{C}(\boldsymbol{w})$ does not admit an explicit expression for $N>1$ , however we show that $\boldsymbol{C}(\boldsymbol{w})$ can be efficiently computed by finding the only positive root of a rational function. The capability to compute $\boldsymbol{C}(\boldsymbol{w})$ at a reasonable cost allows us to design fully conservative explicit schemes, which ensure the correct propagation of discontinuous, or nearly discontinuous, fronts. We discuss the advantages of the fully conservative schemes compared to the non-conservative high-resolution schemes proposed in [4].

Implicit-Explicit strategies are often considered for convection-dominated parabolic systems of PDEs, such as the ED model for $D_{a}<<1$ , since the stability restrictions of these schemes are less severe that those obtained for fully discrete schemes. We propose a second order WENO-IMEX scheme in order to illustrate the numerical difficulties arising from the non-explicit relation between the conserved and physical variables in the ED model. We show that the structure of the Langmuir adsorption isotherms can be used to set up a change of variables that allows us to compute the physical concentrations (instead of the conserved variables) in implicit steps. The solution of the resulting nonlinear system of equations that arises at each time step can be found by Newton’s method and it involves solving a block-tridiagonal linear system, with small blocks that are of the size of the number of components in the mixture, at each iteration. These systems may be efficiently solved by standard block-tridiagonal routines, making the proposed WENO-IMEX-RK technique a very efficient and robust scheme for this model. Several numerical experiments show the performance and capabilities of the proposed numerical scheme.

In this paper we have applied the scheme to the equilibrium dispersive model for a fixed bed. However, it might be adapted to more general models of chromatographic column and also to simulating moving beds (SMB) models.

Acknowledgements

This research was partially supported by Ministerio de Economía y Competitividad under grant MTM2014-54388-P with the participation of FEDER.

References

[1]

M. Mazzotti and A. Rajendran.

Equilibrium theory-based analysis of nonlinear waves in separation processes.

The Annual Review of Chemical and Biomolecular Engineering, 4:119–141, 2013.

[2]

G. Guiochon, G. Shirazi, and M. Katti.

Fundamentals of preparative and nonlinear chromatography (2nd ed.).

Elsevier, 2006.

[3]

A. Seidel-Morgenstern.

Experimental determination of single solute and competitive adsoprtion isotherms.

Journal of Chromatography A, 1037:255–272, 2004.

[4]

S. Javeed, A. Qamar, A. Seidel-Morgenstern, and G. Warnecke.

Efficient and accurate numerical simulation of nonlinear chromatographic processes.

Computers and Chemical Engineering, 35:2294–2305, 2011.

[5]

R. LeVeque.

Numerical Methods for Conservation laws.

Birkhauser-Verlag, 1990.

[6]

R. Donat, F. Guerrero, and P. Mulet.

IMEX WENO schemes for two-phase flow vertical equilibrium processes in a homogeneous porous medium.

Applied Mathematics and Information Sciences, 7(5):1865–1878, 2013.

[7]

Raimund Bürger, Pep Mulet, and Luis M. Villada.

Regularized nonlinear solvers for IMEX methods applied to diffusively corrected multispecies kinematic flow models.

SIAM J. Sci. Comput., 35(3):B751–B777, 2013.

[8]

Sebastiano Boscarino, Raimund Bürger, Pep Mulet, Giovanni Russo, and Luis M. Villada.

Linearly implicit IMEX Runge-Kutta methods for a class of degenerate convection-diffusion problems.

SIAM J. Sci. Comput., 37(2):B305–B331, 2015.

[9]

R. Donat and P. Mulet.

A secular equation for the Jacobian matrix of certain multispecies kinematic flow models.

Numerical Methods for Partial Differential Equations, 26:159–175, 2010.

[10]

S. D. Eidelman and N. V. Zhitarashu.

Parabolic boundary value problems, volume 101 of Operator Theory: Advances and Applications.

Birkhäuser Verlag, Basel, 1998.

[11]

F. Guerrero, R. Donat, and P. Mulet.

Solving a model for 1-D, three-phase flow vertical equilibrium processes in a homogeneous porous medium by means of a Weighted Essentially Non Oscillatory numerical scheme.

Computers and Mathematics with Applications, 66:1284–1298, 2013.

[12]

C.-W. Shu.

High order weighted essentially nonoscillatory schemes for convection dominated problems.

SIAM Review, 51(1):82–126, 2009.

[13]

G.H. Golub and C.F. van Loan.

Matrix Computations.

The Johns Hopkins University Press, 1996.

[14]

J. Cazes.

Encyclopedia of Chromatography.

Den New Dekker Encyclopedias. Taylor & Francis, 2001.

Bibliography14

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] M. Mazzotti and A. Rajendran. Equilibrium theory-based analysis of nonlinear waves in separation processes. The Annual Review of Chemical and Biomolecular Engineering , 4:119–141, 2013.
2[2] G. Guiochon, G. Shirazi, and M. Katti. Fundamentals of preparative and nonlinear chromatography (2nd ed.) . Elsevier, 2006.
3[3] A. Seidel-Morgenstern. Experimental determination of single solute and competitive adsoprtion isotherms. Journal of Chromatography A , 1037:255–272, 2004.
4[4] S. Javeed, A. Qamar, A. Seidel-Morgenstern, and G. Warnecke. Efficient and accurate numerical simulation of nonlinear chromatographic processes. Computers and Chemical Engineering , 35:2294–2305, 2011.
5[5] R. Le Veque. Numerical Methods for Conservation laws . Birkhauser-Verlag, 1990.
6[6] R. Donat, F. Guerrero, and P. Mulet. IMEX WENO schemes for two-phase flow vertical equilibrium processes in a homogeneous porous medium. Applied Mathematics and Information Sciences , 7(5):1865–1878, 2013.
7[7] Raimund Bürger, Pep Mulet, and Luis M. Villada. Regularized nonlinear solvers for IMEX methods applied to diffusively corrected multispecies kinematic flow models. SIAM J. Sci. Comput. , 35(3):B 751–B 777, 2013.
8[8] Sebastiano Boscarino, Raimund Bürger, Pep Mulet, Giovanni Russo, and Luis M. Villada. Linearly implicit IMEX Runge-Kutta methods for a class of degenerate convection-diffusion problems. SIAM J. Sci. Comput. , 37(2):B 305–B 331, 2015.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Taxonomy

Implicit-Explicit WENO schemes for the equilibrium dispersive model of chromatography.

Abstract

keywords:

1 Introduction

2 The mathematical structure of the Equilibrium Dispersive Model

Theorem 2.1**.**

Proof 1**.**

Remark 2.1**.**

Lemma 2.1**.**

Proof 2**.**

Theorem 2.2**.**

Proof 3**.**

Corollary 2.1**.**

Proof 4**.**

Remark 2.2**.**

3 Conservative Numerical Schemes

3.1 Single-component chromatographic elution.

4 IMEX-WENO schemes

4.1 IMEX-WENO schemes for the ED model

5 Numerical experiments

5.1 IMEX versus Explicit Stability Restrictions

5.2 Three component displacement chromatography

6 Conclusions

Acknowledgements

References

Theorem 2.1.

Proof 1.

Remark 2.1.

Lemma 2.1.

Proof 2.

Theorem 2.2.

Proof 3.

Corollary 2.1.

Proof 4.

Remark 2.2.