Iterative solvers for the Maxwell–Stefan diffusion equations: Methods and applications in plasma and particle transport

Abstract

In this paper, we are motivated to discuss a model based on a local thermodynamic equilibrium, weakly ionized plasma-mixture model used for medical and technical applications in etching processes. For studying the model, we consider a simplified model based on the Maxwell–Stefan model, which describes multicomponent diffusive fluxes in the gas mixture. The MS model is more adequate to describe complex mixtures without dominating background species. Based on additional conditions to the fluxes, we obtain an irreducible and quasi-positive diffusion matrix. Such problems result in nonlinear diffusion equations, which are more delicate to solve as simpler standard diffusion equations with Fickian’s approach. Here, we propose an efficient explicit time-discretization method, which is embedded to a fast iterative solver for the nonlinearities. Such a combination of coupling discretization and solver methods allows to simulate the delicate nonlinear differential equations more effectively. We present the efficiency and accuracy of the iterative solvers for some first ternary component gaseous mixtures and discuss the details of the numerical methods.

Keywords:

• Maxwell–Stefan approach
• plasma model
• multicomponent mixture
• explicit discretization schemes
• iterative schemes

AMS subject classifications:

• 35K25
• 35K20
• 74S10
• 70G65

Public Interest Statement

A multicomponent model based on the Stefan–Maxwell approach is presented. Iterative solver approaches are used to solve the nonlinear modelling problem.

1. Introduction

We are motivated to understand the gaseous mixtures of a normal pressure and room temperature plasma. The understanding of normal pressure and room temperature plasma applications is important for applications in medical and technical processes. Since many years, the increasing importance of plasma chemistry based on the multicomponent plasma is a key factor in understanding the gaseous mixture processes (see for low pressure plasma Senega & Brinkmann, 2006 and for atmospheric pressure regimes Tanaka, 2004).

We consider a simplified Maxwell–Stefan diffusion (MSD) equation to model the gaseous mixture of multicomponent plasma. Here, we consider of a macroscopic model, while the limits to apply to a kinetic (microscopic) model are discussed in Boudin, Grec, and Salvarani (2015). While the most classical description of the diffusion goes back to the Fickian’s approach (see Fick, 1995), we apply the modern description of the multicomponent diffusion based on the Maxwell–Stefan’s approach (see Maxwell, 1867). The novel approach considers a more detailed description of the flux and concentration, which are indeed not only proportionally coupled as in the simplified Fickian’s approach. Here, we deal with an inter-species force balance, which allows to model cross-effects, e.g. the so-called reverse diffusion (uphill diffusion in the direction of the gradients).

Such a more detailed modeling results in irreducible and quasi-positive diffusion matrices, which can be reduced by transforming with reductions or with Perron–Frobenius theorems to the solvable partial differential equations (see Bothe, 2011). The obtained system of nonlinear partial differential equations is delicate to solve with standard discretization and solver methods. Therefore, we have taken into account effective linearization methods, e.g. iterative fix-point schemes, to overcome the nonlinearities. Alternative methods exist and are explained in Böttcher (2010) and Spille-Kohoff, Preuß, and Böttcher (2012). Here, they reduce the MSD equation and solve it explicitly, but such methods are restricted to ternary or quaternary systems. Further multicomponent approaches with MSD equations are discussed for only stationary problems, while the solver methods embed the nonlinear MSD approaches into the finite volume discretization schemes (see Peerenboom, van Boxtel, Janssen, & van Dijk, 2014). Such approximations lack with respect to solve nonstationary problems.

The paper is outlined as follows.

In Section 2, we present our mathematical model. A possible reduced model for the further approximations is derived in Section 3. In Section 4, we discuss the underlying numerical schemes. The first numerical results are presented in Section 5. In the contents, that are given in Section 6, we summarize our results.

2. Mathematical model

For the full plasma model, we assume that the neutral particles can be described as the fluid dynamical model, where the elastic collision defines the dynamics and few inelastic collisions are, among other reasons, responsible for the chemical reactions.

To describe the individual mass densities, as well as the global momentum and the global energy as the dynamical conservation quantities of the system, corresponding conservation equations are derived from Boltzmann equations.

The individual character of each species is considered by mass conservation equations and the so-called difference equations.

The extension of the nonmixtured multicomponent transport model (Senega & Brinkmann, 2006) is done with respect to the collision integrals related to the right-hind side sources of the conservation laws.

The conservation laws of the neutral elements are given as$$\begin{array}{cc}& \frac{\mathit{\partial }}{\mathit{\partial }t}{\mathit{\rho }}_s+\frac{\mathit{\partial }}{\mathit{\partial }\mathit{r}}·{\mathit{\rho }}_s{\mathit{u}}_s=m_s{Q}_n^{(s)},\hfill \\ \hfill & \frac{\mathit{\partial }}{\mathit{\partial }t}\mathit{\rho }\mathit{u}+\frac{\mathit{\partial }}{\mathit{\partial }\mathit{r}}·\left({\underline{\underline{P}}}^{\ast }+\mathit{\rho }\mathit{u}\mathit{u}\right)=-Q_m^{(e)},\hfill \\ \hfill & \frac{\mathit{\partial }}{\mathit{\partial }t}{\mathcal{E}}_{tot}^{\ast }+\frac{\mathit{\partial }}{\mathit{\partial }\mathit{r}}·\left({\mathcal{E}}_{tot}^{\ast }\mathit{u}+{\mathit{q}}^{\ast }+{\underline{\underline{P}}}^{\ast }·\mathit{u}\right)=-Q_{\mathcal{E}}^{(e)},\hfill \end{array}$$

where ${\mathit{\rho }}_s$ : density of species i, $\mathit{\rho }={\sum }_{i=1}^N{\mathit{\rho }}_i$, $\mathit{u}$ : velocity, and ${\mathcal{E}}_{tot}^{\ast }$ : total energy of the neutral particles.

Further, the variable $Q_n^{(s)}$ is the collision term of the mass conservation equation, $Q_m^{(e)}$ is the collision term of the momentum conservation equation, and $Q_{\mathcal{E}}^{(e)}$ is the collision term of the energy conservation equation.

We derive the collision term with respect to the Chapmen–Enskog method (see Chapman & Cowling, 1990) and achieve for the first derivatives the following results:(1) $$\begin{array}{cc}\hfill m_s{Q}_n^{(s)}& =-\mathrm{\nabla }·({\mathit{\rho }}_i\sum _{j=0}{\mathbf{V}}_i^j),\hfill \end{array}$$(1) (2) $$\begin{array}{cc}\hfill Q_m^{(e)}& =-\sum _{i=1}^{n_s}{\mathit{\rho }}_i{F}_i,\hfill \end{array}$$(2) (3) $$\begin{array}{cc}\hfill Q_{\mathcal{E}}^{(e)}& =-\sum _{i=1}^{n_s}{\mathit{\rho }}_i\mathit{\rho }{F}_i(\mathbf{u}+\sum _{j=0}{V}_i^{(j)}),\hfill \end{array}$$(3)

where $i=1,\dots ,n_s$, $F_i$ is an external force per unit mass (see Boltzmann equation); further, the diffusion velocity is given as:(4) $$\begin{array}{cc}\hfill {\mathbf{V}}_i^0& =0,\hfill \end{array}$$(4) (5) $$\begin{array}{cc}\hfill {\mathbf{V}}_i^1& =-\sum _{j=1}^N{D}_{ij}(d_j+k_{T_j}\frac{\mathrm{\Delta }T}{T}),\hfill \end{array}$$(5)

where ${\sum }_{i=1}^N{d}_i=0$,(6) $$\begin{array}{cc}\hfill d^i& =\mathrm{\nabla }{x}_i+x_i\frac{\mathrm{\nabla }p}{p}-\frac{{\mathit{\rho }}_i}{\mathit{\rho }}{F}_i,\hfill \end{array}$$(6) (7) $$\begin{array}{cc}\hfill d_i& =d^i-y_i\sum _j{d}_j^{\ast },\hfill \end{array}$$(7)

where $x_i=\frac{n_s}{n}$ is the molar fraction of species i.

We have an additional constraint based on the mass fraction of each species:(8) $$\begin{array}{c}\hfill \frac{\mathit{\partial }}{\mathit{\partial }t}{y}_i+\mathrm{\nabla }{y}_i=R_i(y_1,\dots ,y_N),\end{array}$$(8)

where $y_i$ is the mass fraction of species i and $R_i$ is the net production rate of species i due to the reactions.

Remark 1

The full model problem considers a fully coupled system of conservation laws and Maxwell–Stefan equations. Each equation is coupled such that the gaseous mixture influences the transport equations and vice versa. In the following, we decouple the equations system and consider only the delicate Maxwell–Stefan equations.

3. Simplified model with Maxwell–Stefan diffusion equations

We discuss in the following a multicomponent gaseous mixture with three species (ternary mixture). The model problem is discussed in the experiments of Duncan and Toor (1962).

Here, they studied an ideal gaseous mixture of the following components:

(1)	Hydrogen (${\text{H}}_2$, first species),
(2)	Nitrogen (${\text{N}}_2$, second species), and
(3)	Carbon dioxide (${\text{CO}}_2$, third species).

The Maxwell–Stefan equations are given for the three species as (see also Boudin, Grec, & Salvarani, 2012):(9) $$\begin{array}{cc}& {\mathit{\partial }}_t{\mathit{\xi }}_i+\mathrm{\nabla }·N_i=0,1\le i\le 3,\hfill \end{array}$$(9) (10) $$\begin{array}{cc}& \sum _{j=1}^3{N}_j=0,\hfill \end{array}$$(10) (11) $$\begin{array}{cc}& \frac{{\mathit{\xi }}_2{N}_1-{\mathit{\xi }}_1{N}_2}{D_{12}}+\frac{{\mathit{\xi }}_3{N}_1-{\mathit{\xi }}_1{N}_3}{D_{13}}=-\mathrm{\nabla }{\mathit{\xi }}_1,\hfill \end{array}$$(11) (12) $$\begin{array}{cc}& \frac{{\mathit{\xi }}_1{N}_2-{\mathit{\xi }}_2{N}_1}{D_{12}}+\frac{{\mathit{\xi }}_3{N}_2-{\mathit{\xi }}_2{N}_3}{D_{23}}=-\mathrm{\nabla }{\mathit{\xi }}_2,\hfill \end{array}$$(12)

where the domain is given as $\mathrm{\Omega }\in \mathrm{I}{\mathrm{R}}^d,d\in {IN}^{+}$ with ${\mathit{\xi }}_i\in C^2$.

For such ternary mixture, we can rewrite the three differential Equations (9) and (11 and 12) with the help of the zero-condition (10) into two differential equations, given as:(13) $$\begin{array}{cc}& {\mathit{\partial }}_t{\mathit{\xi }}_i+\mathrm{\nabla }·N_i=0,1\le i\le 2,\hfill \end{array}$$(13) (14) $$\begin{array}{cc}& \frac{1}{D_{13}}{N}_1+\mathit{\alpha }{N}_1{\mathit{\xi }}_2-\mathit{\alpha }{N}_2{\mathit{\xi }}_1=-\mathrm{\nabla }{\mathit{\xi }}_1,\hfill \end{array}$$(14) (15) $$\begin{array}{cc}& \frac{1}{D_{23}}{N}_2-\mathit{\beta }{N}_1{\mathit{\xi }}_2+\mathit{\beta }{N}_2{\mathit{\xi }}_1=-\mathrm{\nabla }{\mathit{\xi }}_2,\hfill \end{array}$$(15)

where $\mathit{\alpha }=\left(\frac{1}{D_{12}}-\frac{1}{D_{13}}\right)$ and $\mathit{\beta }=\left(\frac{1}{D_{12}}-\frac{1}{D_{23}}\right)$.

Further, we have the relations:

(1)	Third mole fraction: ${\mathit{\xi }}_3=1-{\mathit{\xi }}_1-{\mathit{\xi }}_2$,
(2)	Third molar flux: $N_3=-N_1-N_2$.

4. Numerical methods

In the following, we discuss the numerical methods which are based on iterative schemes with embedded explicit discretization schemes (see also Geiser, 2015,in press). We apply the following methods:

(1)	Iterative scheme in time (global linearization with matrix method),
(2)	Iterative scheme in time (local linearization with Richardson’s method).

For spatial discretization, we apply finite volume or finite difference methods. The underlying time discretization is based on a first-order explicit Euler method.

4.1. Iterative scheme in time (global linearization with matrix method)

We solve the iterative scheme:(16) $$\begin{array}{cc}& {\mathit{\xi }}_1^{n+1}={\mathit{\xi }}_1^n-\mathrm{\Delta }t{D}_{+}{N}_1^n,\hfill \end{array}$$(16) (17) $$\begin{array}{cc}& {\mathit{\xi }}_2^{n+1}={\mathit{\xi }}_2^n-\mathrm{\Delta }t{D}_{+}{N}_2^n,\hfill \end{array}$$(17) (18) $$\begin{array}{cc}& \left(\begin{array}{cc}A& B\\ C& D\end{array}\right)\left(\begin{array}{c}{N}_1^{n+1}\hfill \\ N_2^{n+1}\hfill \end{array}\right)=\left(\begin{array}{c}-D_{-}{\mathit{\xi }}_1^{n+1}\hfill \\ -D_{-}{\mathit{\xi }}_2^{n+1}\hfill \end{array}\right),\hfill \end{array}$$(18)

for $j=0,\dots ,J$ , where ${\mathit{\xi }}_1^n={({\mathit{\xi }}_{1,0}^n,\dots ,{\mathit{\xi }}_{1,J}^n)}^T$, ${\mathit{\xi }}_2^n={({\mathit{\xi }}_{2,0}^n,\dots ,{\mathit{\xi }}_{2,J}^n)}^T$ and $I_J\in \mathrm{I}{\mathrm{R}}^{J+1}\times \mathrm{I}{\mathrm{R}}^{J+1}$, $N_1^n={(N_{1,0}^n,\dots ,N_{1,J}^n)}^T$, $N_2^n={(N_{2,0}^n,\dots ,N_{2,J}^n)}^T$ and $I_J\in \mathrm{I}{\mathrm{R}}^{J+1}\times \mathrm{I}{\mathrm{R}}^{J+1}$, where $n=0,1,2,\dots ,N_{end}$ and $N_{end}$ are the number of time steps, i.d. $N_{end}=T/\mathrm{\Delta }t$.

The matrices are given as:(19) $$\begin{array}{cc}& A,B,C,D\in \mathrm{I}{\mathrm{R}}^{J+1}\times \mathrm{I}{\mathrm{R}}^{J+1},\hfill \end{array}$$(19) (20) $$\begin{array}{cc}& A_{j,j}=\frac{1}{D_{13}}+\mathit{\alpha }{\mathit{\xi }}_{2,j},j=0\dots ,J,\hfill \end{array}$$(20) (21) $$\begin{array}{cc}& B_{j,j}=-\mathit{\alpha }{\mathit{\xi }}_{1,j},j=0\dots ,J,\hfill \end{array}$$(21) (22) $$\begin{array}{cc}& C_{j,j}=-\mathit{\beta }{\mathit{\xi }}_{2,j},j=0\dots ,J,\hfill \end{array}$$(22) (23) $$\begin{array}{cc}& D_{j,j}=\frac{1}{D_{23}}+\mathit{\beta }{\mathit{\xi }}_{1,j},j=0\dots ,J,\hfill \end{array}$$(23) (24) $$\begin{array}{cc}& A_{i,j}=B_{i,j}=C_{i,j}=D_{i,j}=0,i,j=0\dots ,J,i\ne J,\hfill \end{array}$$(24)

meaning that the diagonal entries given for the scale case in Equation (13) and the outer diagonal entries are zero.

The explicit form with time discretization is given as:

4.2. Iterative scheme in time (local linearization with Richardson’s method

We solve the iterative scheme given in the Richardson iterative scheme:(47) $$\begin{array}{cc}& {\mathit{\xi }}_1^{n+1,k}={\mathit{\xi }}_1^n-\mathrm{\Delta }t{D}_{+}{N}_1^{n+1},\hfill \end{array}$$(47) (48) $$\begin{array}{cc}& {\mathit{\xi }}_2^{n+1,k}={\mathit{\xi }}_2^n-\mathrm{\Delta }t{D}_{+}{N}_2^{n+1},\hfill \end{array}$$(48) (49) $$\begin{array}{cc}& \left(\begin{array}{cc}{A}^{n+1,k-1}& B^{n+1,k-1}\\ C^{n+1,k-1}& D^{n+1,k-1}\end{array}\right)\left(\begin{array}{c}{N}_1^{n+1}\hfill \\ N_2^{n+1}\hfill \end{array}\right)=\left(\begin{array}{c}-D_{-}{\mathit{\xi }}_1^{n+1,k-1}\hfill \\ -D_{-}{\mathit{\xi }}_2^{n+1,k-1}\hfill \end{array}\right),\hfill \end{array}$$(49)

for $j=0,\dots ,J$ , where ${\mathit{\xi }}_1^n={({\mathit{\xi }}_{1,0}^n,\dots ,{\mathit{\xi }}_{1,J}^n)}^T$, ${\mathit{\xi }}_2^n={({\mathit{\xi }}_{2,0}^n,\dots ,{\mathit{\xi }}_{2,J}^n)}^T$ and $I_J\in \mathrm{I}{\mathrm{R}}^{J+1}\times \mathrm{I}{\mathrm{R}}^{J+1}$, $N_1^n={(N_{1,0}^n,\dots ,N_{1,J}^n)}^T$, $N_2^n={(N_{2,0}^n,\dots ,N_{2,J}^n)}^T$ and $I_J\in \mathrm{I}{\mathrm{R}}^{J+1}\times \mathrm{I}{\mathrm{R}}^{J+1}$, where $n=0,1,2,\dots ,N_{end}$ and $N_{end}$ are the number of time steps, i.d. $N_{end}=T/\mathrm{\Delta }t$.

Further, $k=1,2,\dots ,K$ is the iteration index where ${\mathit{\xi }}_1^{n+1,0}={({\mathit{\xi }}_{1,0}^n,\dots ,{\mathit{\xi }}_{1,J}^n)}^T$, ${\mathit{\xi }}_2^{n+1,0}={({\mathit{\xi }}_{2,0}^n,\dots ,{\mathit{\xi }}_{2,J}^n)}^T$, and $I_J\in \mathrm{I}{\mathrm{R}}^{J+1}\times \mathrm{I}{\mathrm{R}}^{J+1}$ is the start solution given with the solution at $t=t^n$.

The matrices are given as:(50) $$\begin{array}{cc}& A^{n+1,k-1},B^{n+1,k-1},C^{n+1,k-1},D^{n+1,k-1}\in \mathrm{I}{\mathrm{R}}^{J+1}\times \mathrm{I}{\mathrm{R}}^{J+1},\hfill \end{array}$$(50) (51) $$\begin{array}{cc}& A_{j,j}^{n+1,k-1}=\frac{1}{D_{13}}+\mathit{\alpha }{\mathit{\xi }}_{2,j}^{n+1,k-1},j=0\dots ,J,\hfill \end{array}$$(51) (52) $$\begin{array}{cc}& B_{j,j}^{n+1,k-1}=-\mathit{\alpha }{\mathit{\xi }}_{1,j}^{n+1,k-1},j=0\dots ,J,\hfill \end{array}$$(52) (53) $$\begin{array}{cc}& C_{j,j}^{n+1,k-1}=-\mathit{\beta }{\mathit{\xi }}_{2,j}^{n+1,k-1},j=0\dots ,J,\hfill \end{array}$$(53) (54) $$\begin{array}{cc}& D_{j,j}^{n+1,k-1}=\frac{1}{D_{23}}+\mathit{\beta }{\mathit{\xi }}_{1,j}^{n+1,k-1},j=0\dots ,J,\hfill \end{array}$$(54) (55) $$\begin{array}{cc}& A_{i,j}^{n+1,i-1}=B_{i,j}^{n+1,i-1}=C_{i,j}^{n+1,i-1}=D_{i,j}^{n+1,i-1}=0,quadi,j=0\dots ,J,i\ne J,\hfill \end{array}$$(55)

meaning that the diagonal entries given for the scale case in Equation (95) and the outer diagonal entries are zero.

The explicit form with time discretization is given as:

5. Numerical experiments

In the following, we concentrate on the following three-component system, which is given as:(79) $$\begin{array}{cc}& {\mathit{\partial }}_t{\mathit{\xi }}_i+{\mathit{\partial }}_x{N}_i=0,1\le i\le 3,\hfill \end{array}$$(79) (80) $$\begin{array}{cc}& \sum _{j=1}^3{N}_j=0,\hfill \end{array}$$(80) (81) $$\begin{array}{cc}& \frac{{\mathit{\xi }}_2{N}_1-{\mathit{\xi }}_1{N}_2}{D_{12}}+\frac{{\mathit{\xi }}_3{N}_1-{\mathit{\xi }}_1{N}_3}{D_{13}}=-{\mathit{\partial }}_x{\mathit{\xi }}_1,\hfill \end{array}$$(81) (82) $$\begin{array}{cc}& \frac{{\mathit{\xi }}_1{N}_2-{\mathit{\xi }}_2{N}_1}{D_{12}}+\frac{{\mathit{\xi }}_3{N}_2-{\mathit{\xi }}_2{N}_3}{D_{23}}=-{\mathit{\partial }}_x{\mathit{\xi }}_2,\hfill \end{array}$$(82)

where the domain is given as $\mathrm{\Omega }\in \mathrm{I}{\mathrm{R}}^d,d\in {IN}^{+}$ with ${\mathit{\xi }}_i\in C^2$.

The parameters and the initial and boundary conditions are given as:

(1)	$D_{12}=D_{13}=0.833$ (means $\mathit{\alpha }=0$) and $D_{23}=0.168$ (Uphill diffusion, semi-degenerated Duncan and Toor experiment),
(2)	$D_{12}=0.0833,D_{13}=0.680$ and $D_{23}=0.168$ (asymptotic behavior, Duncan and Toor experiment (see Duncan & Toor, 1962)),
(3)	$J=140$ (spatial grid points),
(4)	The time step restriction for the explicit method is given as: $\mathrm{\Delta }t\le \frac{{(\mathrm{\Delta }x)}^2}{2max\{D_{12},D_{13},D_{23}\}}$,
(5)	The spatial domain is $\mathrm{\Omega }=[0,1]$, the time-domain $[0,T]=[0,1]$,
(6)	The initial conditions are: (i) Uphill example(83) $$\begin{array}{cc}\hfill {\mathit{\xi }}_1^{in}(x)& =\left\{\begin{array}{cc}0.8\hfill & \text{if}0\le x<0.25,\hfill \\ 1.6(0.75-x)\hfill & \text{if}0.25\le x<0.75,\hfill \\ 0.0\hfill & \text{if}0.75\le x\le 1.0,\hfill \end{array}\right.\hfill \end{array}$$(83) (84) $$\begin{array}{cc}\hfill {\mathit{\xi }}_2^{in}(x)& =0.2,\text{for all}x\in \mathrm{\Omega }=[0,1].\hfill \end{array}$$(84) (ii) Diffusion example (Asymptotic behavior)(85) $$\begin{array}{cc}\hfill {\mathit{\xi }}_1^{in}(x)& =\left\{\begin{array}{cc}0.8\hfill & \text{if}0\le x\in 0.5,\hfill \\ 0.0\hfill & \text{else},\hfill \end{array}\right.\hfill \end{array}$$(85) (86) $$\begin{array}{cc}\hfill {\mathit{\xi }}_2^{in}(x)& =0.2,\text{for all}x\in \mathrm{\Omega }=[0,1].\hfill \end{array}$$(86)
(7)	The boundary conditions are of no-flux type:(87) $$\begin{array}{cc}& N_1=N_2=N_3=0,\text{on}\mathit{\partial }\mathrm{\Omega }\times [0,1].\hfill \end{array}$$(87)

We could reduce to a simpler model problem as:(88) $$\begin{array}{cc}& {\mathit{\partial }}_t{\mathit{\xi }}_i+{\mathit{\partial }}_x·N_i=0,1\le i\le 2,\hfill \end{array}$$(88) (89) $$\begin{array}{cc}& \frac{1}{D_{13}}{N}_1+\mathit{\alpha }{N}_1{\mathit{\xi }}_2-\mathit{\alpha }{N}_2{\mathit{\xi }}_1=-{\mathit{\partial }}_x{\mathit{\xi }}_1,\hfill \end{array}$$(89) (90) $$\begin{array}{cc}& \frac{1}{D_{23}}{N}_2-\mathit{\beta }{N}_1{\mathit{\xi }}_2+\mathit{\beta }{N}_2{\mathit{\xi }}_1=-{\mathit{\partial }}_x{\mathit{\xi }}_2,\hfill \end{array}$$(90)

where $\mathit{\alpha }=\left(\frac{1}{D_{12}}-\frac{1}{D_{13}}\right)$, $\mathit{\beta }=\left(\frac{1}{D_{12}}-\frac{1}{D_{23}}\right)$.

We rewrite into:(91) $$\begin{array}{cc}& {\mathit{\partial }}_t{\mathit{\xi }}_1+{\mathit{\partial }}_x·N_1=0,\hfill \end{array}$$(91) (92) $$\begin{array}{cc}& {\mathit{\partial }}_t{\mathit{\xi }}_2+{\mathit{\partial }}_x·N_2=0,\hfill \end{array}$$(92) (93) $$\begin{array}{cc}& \left(\begin{array}{cc}\frac{1}{D_{13}}+\mathit{\alpha }{\mathit{\xi }}_2& -\mathit{\alpha }{\mathit{\xi }}_1\\ -\mathit{\beta }{\mathit{\xi }}_2& \frac{1}{D_{23}}+\mathit{\beta }{\mathit{\xi }}_1\end{array}\right)\left(\begin{array}{c}{N}_1\hfill \\ N_2\hfill \end{array}\right)=\left(\begin{array}{c}-{\mathit{\partial }}_x{\mathit{\xi }}_1\hfill \\ -{\mathit{\partial }}_x{\mathit{\xi }}_2\hfill \end{array}\right),\hfill \end{array}$$(93)

and we have(94) $$\begin{array}{cc}& {\mathit{\partial }}_t{\mathit{\xi }}_1+{\mathit{\partial }}_x·N_1=0,\hfill \end{array}$$(94) (95) $$\begin{array}{cc}& {\mathit{\partial }}_t{\mathit{\xi }}_2+{\mathit{\partial }}_x·N_2=0,\hfill \end{array}$$(95) (96) $$\begin{array}{cc}& \left(\begin{array}{c}{N}_1\hfill \\ N_2\hfill \end{array}\right)=\frac{D_{13}{D}_{23}}{1+\mathit{\alpha }{D}_{13}{\mathit{\xi }}_2+\mathit{\beta }{D}_{23}{\mathit{\xi }}_1}\left(\begin{array}{cc}\frac{1}{D_{23}}+\mathit{\beta }{\mathit{\xi }}_1& \mathit{\alpha }{\mathit{\xi }}_1\\ \mathit{\beta }{\mathit{\xi }}_2& \frac{1}{D_{13}}+\mathit{\alpha }{\mathit{\xi }}_2\end{array}\right)\left(\begin{array}{c}-{\mathit{\partial }}_x{\mathit{\xi }}_1\hfill \\ -{\mathit{\partial }}_x{\mathit{\xi }}_2\hfill \end{array}\right).\hfill \end{array}$$(96)

The next step is to apply the semi-discretization of the partial differential operator $\frac{\mathit{\partial }}{\mathit{\partial }x}$.

We apply the first differential operator in Equations (93) and (94) as an forward upwind scheme given as(97) $$\begin{array}{c}\hfill \frac{\mathit{\partial }}{\mathit{\partial }x}=D_{+}=\frac{1}{\mathrm{\Delta }x}·\left(\begin{array}{ccccc}-1& 0& \dots & & 0\\ 1& -1& 0& \dots & 0\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& & 1& -1& 0\\ 0& \dots & 0& 1& -1\end{array}\right)\in \mathrm{I}{\mathrm{R}}^{(J+1)\times (J+1)},\end{array}$$(97)

and the second differential operator in Equation (13) as an backward upwind scheme given as(98) $$\begin{array}{c}\hfill \frac{\mathit{\partial }}{\mathit{\partial }x}=D_{-}=\frac{1}{\mathrm{\Delta }x}·\left(\begin{array}{ccccc}-1& 1& 0& \dots & 0\\ 0& -1& 1& 0& \dots \\ ⋮& \ddots & \ddots & \ddots & \ddots \\ 0& \dots & 0& -1& 1\\ 0& & \dots & 0& -1\end{array}\right)\in \mathrm{I}{\mathrm{R}}^{(J+1)\times (J+1)}.\end{array}$$(98)

5.1. Experiments with the iterative scheme in time (global linearization)

In the first experiments, we test the first iterative scheme (iterative scheme in time (global linearization)).

We test the different experiments (uphill and diffusion examples) and obtain the results as shown in Figure 1 for the uphill example.

Figure 1. Results of the mole fraction (concentration) ${\mathit{\xi }}_1$, ${\mathit{\xi }}_2$ and ${\mathit{\xi }}_3$.

Remark 2

In Figure 1, we obtain a typical complex mixture, without dominating background gas. While the concentration ${\mathit{\xi }}_1$ increases and ${\mathit{\xi }}_3$ decreases, the important concentration ${\mathit{\xi }}_2$ decreases and increases around the value 0.2. Such a behavior cannot be produced with simple standard Fickian’s approach and here it is important to deal with the MS approach.

The concentration and their fluxes are given in Figure 2.

Figure 2. Concentration and flux of iterative scheme in time (global linearization). Notes: The upper figures present the results of the flux $N_1$ and concentration gradient $-{\mathit{\partial }}_x{\mathit{\xi }}_1$. The lower figures present the results of flux $N_2$ and concentration gradient $-{\mathit{\partial }}_x{\mathit{\xi }}_2$.

Remark 3

In Figure 2, we obtain a fine resolved behavior of the important second species ${\mathit{\xi }}_2$. Here, we see detailed the oscillatory behavior of the flux $N_2$ and the concentration gradient $-{\mathit{\partial }}_x{\mathit{\xi }}_2$ in the time interval $t=[0,1]$. For the next time interval $t=[1,2]$, we obtained a stabilized behavior of the complex mixture. Here, it is important to resolve the nonlinearity very accurately for unstable initialization of the complex mixture.

The full plots in time and space of the concentrations and their fluxes are given in Figure 3.

Figure 3. Results of the 3D plots in time and space. Notes: The upper figures present the results of the flux $N_1$ and concentration gradient $-{\mathit{\partial }}_x{\mathit{\xi }}_1$. The lower figures present the results of the flux $N_2$ and concentration gradient $-{\mathit{\partial }}_x{\mathit{\xi }}_2$.

The space–time regions where $-N2{\mathit{\partial }}_x{\mathit{\xi }}_2\ge 0$ for the uphill diffusion and asymptotic diffusion, given in Figure 4.

Figure 4. Asymptotic diffusion (left-hand side) and uphill diffusion (right-hand side) in the space–time region.

Remark 4

In Figure 3, we present the delicate unstable behavior of the mixture around time scale $t=[0,1]$ and the spatial scale $x=[0,1]$ in a 3D plot with the first and second concentration. The steep gradients of the concentrations have to be resolved with a very fine time step. We also obtain the same delicate results of the initialization process in Figure 4. Here, we present the space–time regions for the uphill and asymptotic diffusions. In both experiments, we see that we have a backflow of the mixture, meaning the second concentration ${\mathit{\xi }}_2$ has both decreasing and increasing gradients. For both experiments, we resolve the initial process with fine time steps.

Remark 5

In the first numerical method, we apply global linearization based on the time steps. Meaning, we deal with an explicit time discretization that solves the linearized equation forward. All effects are resolved by taking into account the CFL condition. Therefore, we achieve better results with finer time steps, e.g. $\mathrm{\Delta }{t}_{CFl}/8$, such that the global linearization, via the time step, is important for our first numerical method.

5.2. Iterative schemes in time (local linearization)

In the next series of experiments, we apply the more refined linearization scheme, meaning the iterative approximation in a single time step. We deal with a backward idea and resolve the nonlinearity with a fix-point method (see Geiser, 2011) in each local time step. Such a modified method allows to enlarge the local time steps while we have an additional fix-point scheme for the nonlinearity (see Kelley, 1995).

We apply the numerical convergence of the schemes with the reference solution of the explicit method by ${\mathrm{\Xi }}_{ref}=({\mathit{\xi }}_1,{\mathit{\xi }}_2)$, where the time step is $\mathrm{\Delta }{t}_{CFL}/8$ for this refined solution, with the error only marginal.

Based on the reference solution, we deal with the following errors:(99) $$\begin{array}{cc}\hfill E_{L_{1,\mathrm{\Delta }x}}(t)& ={\int }_{\mathrm{\Omega }}|{\mathrm{\Xi }}_{method,J,\mathrm{\Delta }x}(x,t)-{\mathrm{\Xi }}_{ref}(x,t)|dx\hfill \\ \hfill & =\mathrm{\Delta }x\sum _{i=1}^N|{\mathrm{\Xi }}_{method,J,\mathrm{\Delta }x}(x_i,t)-{\mathrm{\Xi }}_{ref}(x_i,t)|,\hfill \end{array}$$(99)

where $method,J$ is the Richardson method (second method, see Section 4.1) with J iterative steps and $\mathrm{\Delta }t=\mathrm{\Delta }{t}_{CFL},\mathrm{\Delta }{t}_{CFL}/2,\mathrm{\Delta }{t}_{CFL}/4$.

Further, $method,expl$ is the explicit method (first method, see Section 4.2) with $\mathrm{\Delta }t=\mathrm{\Delta }{t}_{CFL},\mathrm{\Delta }{t}_{CFL}/2,\mathrm{\Delta }{t}_{CFL}/4$.

We apply different versions of time steps and iterative steps and a reference solution is obtained with a fine time step $\mathrm{\Delta }t=\mathrm{\Delta }{t}_{CFL}/4$. We see improvements in Figure 5 and the errors in Figure 6.

Figure 5. Solutions of different time steps and iterative steps of the Richardson method (left-hand side: concentration ${\mathit{\xi }}_1$ and right-hand side: concentration ${\mathit{\xi }}_2$).

Remark 6

In Figure 5, we compare a reference solution (first method) with fine time steps with flexible time step and iterative step solutions (second method). We obtain some accurate solutions with the second method with more larger time steps and more iterative steps. Based on the large time steps, the computational time decreases and additional number of more iterative cycles did not influence the amount of computational work.

Figure 6. Errors of the different time step and iterative step solutions of the second method (left-hand side: error with the reference solution at the full time interval $t\in [0,1.4]$ and right-hand side: error with the reference solutions at the initial time interval $t\in [0,0.1]$).

Remark 7

In Figure 6, we present the $L_1$ errors of the second method (local linearization) compared with the reference solution (first method with very fine time steps). Here, we see the benefit of large time steps with $N=100$ and $K=800$, meaning we have only 100 time steps and $800$ iterative steps, which are not expensive. Therefore, we could gain the same results as with many small time steps $N=80000$ and only one iterative step $K=1$, such that the relaxation method benefits with the iterative cycles and we could enlarge the time steps. We obtain the same accurate results for larger time steps, instead of a computational intensive fine resolution with small time steps. Such a novel treatment with a local linearization and additional iterative cycles reduces the computational amount of the novel second scheme.

Remark 8

The second method applies a linear linearization based on the iterative approaches in each single time step. The spatial and time discretizations are embedded into the iterative solver method. We have the benefit of relaxation in each local time step with the high resolution of the spatial and timescales. Therefore, we see a more accurate solution also with larger time steps than in the global linearization method.

6. Conclusions and discussion

We present a fluid model based on a delicate mixture of components. Such a model can be resolved by a MSD equation, while we can embed the complex mixture processes. The underlying problems for such a more delicate diffusion matrix are discussed. Based on a delicate nonlinear partial differential equation, which has to be reduced to a solvable linearized partial differential equation, we have to discuss two linearization approaches. The first approach deals with a global linearization, which can be controlled by the time steps. The second approach deals with a local linearization and applies an iterative scheme. Such a novel approach is more flexible and can be controlled by the time steps and additionally with the iterative steps. Therefore, we obtain the same accuracy with much more larger time steps and moderate iterative steps and reduce the computational amount. For the first test examples, we achieve more accurate results for a new local linearized scheme and optimize their computational amount with larger time steps. In future, we concentrate on numerical convergence analysis of the local linearization schemes and generalize our results to real-life applications.

Additional information

Funding

Funding. The author has received no direct funding for this research.

Notes on contributors

Jürgen Geiser

Jüergen Geiser (http://homepage.ruhr-uni-bochum.de/Juergen.Geiser/), researcher and lecturer at the Ruhr-University of Bochum, Germany, has been involved in teaching and research projects and has collaborated with engineering and physicist groups on numerical modeling of technical and physical models. The research activity refers to the mathematical modeling, numerics, and analysis of transport and flow problems in engineering applications, e.g. groundwater modeling and plasma modelling. He is a specialist in multiscale solvers and iterative solvers and most of the topics of the special issue. Moreover, Juergen Geiser is the author of scientific books and editor of various scientific journals, and thus able to manage the editorial activity.