Splitting Algorithm for Numerical Solving of Parabolic Nonlinear Multi - Component Convection – Diffusion Mass Transfer Equations

1. Introduction

Any process modelled by time dependent PDEs can be naturally split into simpler additive subprocesses. The splitting (or fractional step) methods, [1,2,3], decouple a complex, high-dimensional problem into simpler sub-problems (the classical divide et impera strategy). The main advantage offered by this approach consists of: each sub-problem can be time integrated using a solver designed for its particular characteristics. The error generated by operator splitting is the main drawback of this method.

Today, splitting methods are widely employed to construct geometric integrators for many problems, [4]. They are commonly used in environmental sciences, combustion and life science processes. Due to the present framework, only the articles that deal with splitting methods applied to non-linear convection (advection) – diffusion – reaction equations will be mentioned in this section.

For an advection-diffusion-reaction equation from atmospheric pollution, Lanser and Verwer [5] studied the Strang splitting error using Lie operators. The convergence analysis for a high order splitting scheme (Richardson extrapolation and Strang splitting) applied to a non-linear diffusion-reaction equation can be viewed in Descombes [6]. The splitting error for both linear and non-linear splitting methods was also analysed by Hundsdorfer and Verwer, [7]. Fully implicit, semi-implicit and operator-splitting algorithms are the time-integration methods investigated in [8,9] for numerical solving of non-linear reaction-diffusion, [8], and radiation-diffusion, [9], equations. For reaction-diffusion equations, the stability of the splitting methods is analysed in [10]. The results obtained in [10] were extended to advection-diffusion-reaction equations in [11]. The case with indefinite reaction terms was also analysed in [11]. A brief review of the splitting methods used in air pollution modeling can be viewed in [12]. Numerical strategies developed to improve the splitting techniques applied to non-linear diffusion- reaction and convection-diffusion-reaction equations are presented in [13,14]. These strategies use estimates of the splitting error. Recent splitting techniques developed for non-linear convection-diffusion-reaction equations can be viewed in [15,16,17]. In references [5,6,7,8,9,10,11,12,13,14,15,16,17] the reaction / radiation term is non-linear.

The mathematical model for mass transfer in multi-component mixtures, [18,19], is a strongly coupled system of partial differential equations. A review of the algorithms employed for numerical solving of the parabolic multi-component mass transfer equations can be seen in [20,21]. Until [20,21], the splitting method was not applied to multi-component mass transfer equations. In the present work, the splitting strategy developed [20,21] for linear multi-component convection – diffusion equations is extended to non-linear, multi-component (the diffusion coefficients depend on chemical species concentrations), convection – diffusion equations. It must be mentioned that the present nonlinearity is not similar to those usually encountered in PDEs. Generally accepted explicit relations for the dependence diffusion coefficients versus concentrations are not available for any number of chemical species, n. The splitting error is analysed theoretically. Numerical experiments that investigate the discrete error approximation are presented.

Section 2 shows the analysis of the splitting error and the present algorithm. Section 3 shows the test problem. The numerical simulations and results can be viewed in Sect. 4. The conclusions of this work are mentioned in Sect. 5.

2. Splitting Method

2.1. Numerical algorithm description

For a constant mass density multi-component mixture with n chemical species, the balance equations for the chemical species concentrations, defined on Ω × (0, ∞), Ω ⊂ R^N, are as follows:

$$\frac{\partial w_1}{\partial t}+v ·{\nabla w}_1=\sum _{k=1}^{n-1}\nabla ·\left(D_{1,k}\left(w\right)\nabla w_k\right)$$

$$\frac{\partial w_{n-1}}{\partial t}+v·\nabla w_{n-1}={\sum }_{k=1}^{n-1}\nabla ·\left(D_{n-1,k}\left(w\right)\nabla w_k\right)$$

(1)

where w_k is the mass concentration of component k, $w=\left(w_1,w_2,....\right)$, v is the velocity vector, $D_{i,k} \left(w\right)$ is the Fick diffusion coefficient and ∇ the del (nabla) operator. In equations (1) the chemical reactions are not taken into consideration. In the next paragraphs, for the brevity of writing, the Fick diffusion coefficient will be symbolized by D_i,k.

Discretizing the spatial derivatives from (1), one obtains,

$$\frac{\partial w_1}{\partial t}+{\Lambda }_h{w}_1=\sum _{k=1}^{n-1}{{\nabla }_h D}_{1,k}{\nabla }_h{w}_k$$

$$\frac{\partial w_{n-1}}{\partial t}+{\Lambda }_h{w}_{n-1}={\sum }_{k=1}^{n-1}{{\nabla }_{h}D}_{n-1,k}{\nabla }_h{w}_k$$

(2)

where ∇_h and Λ_h are the discrete approximations of the continuous operators. The present splitting method is defined by:

$$\begin{gathered} \left[\begin{array}{cccc}{{\nabla }_h D}_{11}{\nabla }_h-{\Lambda }_h& {{\nabla }_h D}_{12}{\nabla }_h& \cdots & {{\nabla }_h D}_{1,n-1}{\nabla }_h\\ {{\nabla }_{ h}D}_{21}{\nabla }_h& {{\nabla }_{ h}D}_{22}{\nabla }_h-{\Lambda }_h& \cdots & {{\nabla }_{h}D}_{2,n-1}{\nabla }_h\\ ⋮& ⋮& \ddots & ⋮\\ {{\nabla }_{h}D}_{n-\mathrm{1,1}}{\nabla }_h& \cdots & \cdots & {\nabla }_h{D}_{n-1,n-1}{\nabla }_h-{\Lambda }_h\end{array}\right]= \\ \left[\begin{array}{cccc}\frac{1}{2}\left({\nabla }_h{D}_{11}{\nabla }_h-{\Lambda }_h\right)& {\nabla }_h{D}_{12}{\nabla }_h& \cdots & {{\nabla }_{h}D}_{1,n-1}{\nabla }_h\\ 0& \frac{1}{2}\left({\nabla }_h{D}_{22}{\nabla }_h-{\Lambda }_h\right)& \cdots & {{\nabla }_{h}D}_{2,n-1}{\nabla }_h\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \frac{1}{2}\left({{\nabla }_{h}D}_{n-1,n-1}{\nabla }_h-{\Lambda }_h\right)\end{array}\right]+ \\ \left[\begin{array}{cccc}\frac{1}{2}\left({\nabla }_h{D}_{11}{\nabla }_h-{\Lambda }_h\right)& 0& \cdots & 0\\ {{\nabla }_{h}D}_{21}{\nabla }_h& \frac{1}{2}\left({{\nabla }_{h}D}_{22}{\nabla }_h-{\Lambda }_h\right)& 0& 0\\ ⋮& ⋮& \ddots & ⋮\\ {{\nabla }_{h}D}_{n-\mathrm{1,1}}{\nabla }_h& \cdots & \cdots & \frac{1}{2}\left({\nabla }_h{D}_{n-1,n-1}{\nabla }_h-{\Lambda }_h\right)\end{array}\right] \end{gathered}$$

(3)

Relation (3) shows that in the present splitting algorithm, the numerical solving of a non-linear block algebraic system is reduced to the numerical solving of two non-linear block-triangular algebraic systems. Obviously, the previous statement is valid if implicit or semi-implicit methods are used for time integration.

The test cases analyzed in this work have 2-D spatial operators, i.e. R^N = R². Consider the geometric splitting of ∇_h and Λ_h as:

$${\nabla }_h={\nabla }_{h,x}+{\nabla }_{h,y},{\Lambda }_h={\Lambda }_{h,x}+{\Lambda }_{h,y}$$

(4)

and denote by $L_{h,x}^s,\hspace{0.17em}{L}_{h,y}^s,\hspace{0.17em}{L}_{h,x}^i,\hspace{0.17em}{L}_{h,y}^i$, the matrices

$$\begin{gathered} L_{h,x}^s=\left[\begin{array}{cccc}\frac{1}{2}\left({{\nabla }_{h,x}D}_{11}{\nabla }_{h,x}-{\Lambda }_{h,x}\right)& {\nabla }_{h,x}{D}_{12}{\nabla }_{h,x}& \cdots & {{\nabla }_{h,x}D}_{1,n-1}{\nabla }_{h,x}\\ 0& \frac{1}{2}\left({\nabla }_{h,x}{D}_{22}{\nabla }_{h,x}-{\Lambda }_{h,x}\right)& \cdots & {\nabla }_{h,x}{D}_{2,n-1}{\nabla }_{h,x}\\ ⋮& ⋮& \ddots & ⋮\\ 0& \cdots & 0& \frac{1}{2}\left({{\nabla }_{h,x}D}_{n-1,n-1}{\nabla }_{h,x}-{\Lambda }_{h,x}\right)\end{array}\right], \\ {L}_{h,y}^s=\left[\begin{array}{cccc}\frac{1}{2}\left({{\nabla }_{h,y}D}_{11}{\nabla }_{h,y}-{\Lambda }_{h,y}\right)& {{\nabla }_{h,y}D}_{12}{\nabla }_{h,y}& \cdots & {{\nabla }_{h,y}D}_{1,n-1}{\nabla }_{h,y}\\ 0& \frac{1}{2}\left({{\nabla }_{h,y}D}_{22}{\nabla }_{h,y}-{\Lambda }_{h,y}\right)& \cdots & {{\nabla }_{h,y}D}_{2,n-1}{\nabla }_{h,y}\\ ⋮& ⋮& \ddots & ⋮\\ 0& \cdots & 0& \frac{1}{2}\left({\nabla }_{h,y}{D}_{n-1,n-1}{\nabla }_{h,y}-{\Lambda }_{h,y}\right)\end{array}\right] \\ {L}_{h,x}^i=\left[\begin{array}{cccc}\frac{1}{2}\left({{\nabla }_{h,x}D}_{11}{\nabla }_{h,x}-{\Lambda }_{h,x}\right)& 0& \cdots & 0\\ {{\nabla }_{h,x}D}_{21}{\nabla }_{h,x}& \frac{1}{2}\left({\nabla }_{h,x}{D}_{22}{\nabla }_{h,x}-{\Lambda }_{h,x}\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ {\nabla }_{h,x}{D}_{n-\mathrm{1,1}}{\nabla }_{h,x}& ....& \cdots & \frac{1}{2}\left({{\nabla }_{h,x}D}_{n-1,n-1}{\nabla }_{h,x}-{\Lambda }_{h,x}\right)\end{array}\right] \\ {L}_{h,y}^i=\left[\begin{array}{cccc}\frac{1}{2}\left({{\nabla }_{h,y}D}_{11}{\nabla }_{h,y}-{\Lambda }_{h,y}\right)& 0& \cdots & 0\\ {{\nabla }_{h,y}D}_{21}{\nabla }_{h,y}& \frac{1}{2}\left({\nabla }_{h,y}{D}_{22}{\nabla }_{h,y}-{\Lambda }_{h,y}\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ {\nabla }_{h,y}{D}_{n-\mathrm{1,1}}{\nabla }_{h,y}& \cdots & \cdots & \frac{1}{2}\left({{\nabla }_{h,y}D}_{n-1,n-1}{\nabla }_{h,y}-{\Lambda }_{h,y}\right)\end{array}\right] \end{gathered}$$

(5)

One step for time integration is the Strang integrator with four pieces, [2]:

$$\frac{\partial w^1}{\partial t}+L_{h,x}^s{w}^1=0,t_{i-1}\le t\le t_{i-1/2}$$

(6a)

$$\frac{\partial w^2}{\partial t}+L_{h,y}^s{w}^2=0,t_{i-1}\le t\le t_{i-1/2}$$

(6b)

$$\frac{\partial w^3}{\partial t}+L_{h,x}^i{w}^3=0,t_{i-1}\le t\le t_{i-1/2}$$

(6c)

$$\frac{\partial w^4}{\partial t}+L_{h,y}^i{w}^4=0,t_{i-1}\le t\le t_i$$

(6d)

$$\frac{\partial w^5}{\partial t}+L_{h,x}^i{w}^5=0,t_{i-1/2}\le t\le t_i$$

(6e)

$$\frac{\partial w^6}{\partial t}+L_{h,y}^s{w}^6=0,t_{i-1/2}\le t\le t_i$$

(6f)

$$\frac{\partial w^7}{\partial t}+L_{h,x}^s{w}^7=0,t_{i-1/2}\le t\le t_i$$

(6g)

where $w=\left(w_1,w_2,....\right)$. The initial conditions are: $w^1$ ( t_i-1 ) = w ( t_i-1 ), $w^k$ ( t_i-1 ) = $w^{k-1}$ ( t_i-1/2 ) ( for k = 2, 3, 4 ) and $w^k$ ( t_i-1/2 ) = $w^{k-1}$ ( t_i-1/2 ) ( for k = 5, 6, 7 ). In relations (6), the superscript refers to each fractional step. The reasons for the selection of the Strang splitting integrator are presented in the sub-section 2.2.

2.2. Splitting Error Analysis

The truncation error for any time-split algorithm is given by the sum of the following two terms: (1) the splitting error and (2) the error of the discrete approximation (finite differences, finite elements, finite volumes, etc.). The splitting error will be analysed in this section.

Consider the general, non-linear, autonomous system of differential equations,

$$\frac{\partial w}{\partial t}=A \left(w\right),$$

(7)

with the initial condition,

w (x, t₀) = w₀, x ∈ Rⁿ.

(8)

where A is a smooth function of w. For well-posed problems, the solution of (7) at any time t₁ > t₀, (t₁ = t₀ + δ t) is uniquely determined by the initial condition. This solution can be written as,

w ( t₁ ) = S ( t₁, t₀ ) w₀ = S ( δ t ) w₀.

(9)

where S ( δ t ) is the solution operator.

The operator A can be additively split into two or more pieces. Consider the case of two pieces,

A ( w ) = A₁ (w) + A₂ ( w ).

(10)

where A₁ and A₂ are smooth functions of w. The resulted subproblems

$$\frac{\partial w}{\partial t}=A_1 \left(w\right),$$

(11a)

and

$$\frac{\partial w}{\partial t}=A_2 \left(w\right),$$

(11b)

should be easier to solve than the full problem. Denote by S₁ ( δ t ) and S₂ ( δ t ) the solution operators for the subproblems (11a) and (11b). The usual or Lie-Trotter splitting algorithm corresponding to (11) is,

w ( t₁ ) = S₂ ( δ t ) S₁ ( δ t ) w₀

(12)

For the algorithm (12), the splitting error operator E_split ( δ t ) is defined by,

E_split ( δ t ) = S₂ ( δ t ) S₁ ( δ t ) – S ( δ t)

Using Taylor expansions, one can prove (see LeVeque, [22]) that the splitting error for the algorithm (12) verifies the relation:

$$E_{split} \left(\delta t\right)= S_2 \left(\delta t\right) S_1 \left(\delta t\right)-S \left(\delta t\right)= \frac{{\delta t}^2}{2} \left[A_1 , A_2\right]+ \mathcal{O} \left({\delta t}^3\right)$$

(13)

where [A₁, A₂] is the Lie bracket (or Lie-Jacobi bracket or commutator) defined by

$$\left[A_1, A_2\right] \left(w\right)=A_{2w} \left(w\right) A_1\left(w\right)-A_{1w}\left(w\right) A_2\left(w\right)$$

(14)

and A_iw is

$$A_{iw}=\frac{\partial A_i}{\partial w}, i=1, 2.$$

Relation (13) shows that if [A₁,A₂] = 0, the Lie-Trotter splitting is second order accurate. If [A₁,A₂] does not vanish, the Lie-Trotter splitting is only first order accurate. To compute the Lie bracket for the splitting scheme used in this work, we will consider, for brevity of writing, the case of a ternary system. For a ternary system, the present splitting can be written as,

$$\begin{gathered} \left[\begin{array}{cc}\nabla D_{11}\nabla w_1-v \nabla w_1& \nabla D_{12}\nabla w_2\\ \nabla D_{21}\nabla w_1& \nabla D_{22}\nabla w_2-v \nabla w_2\end{array}\right]= \\ (A) \\ \left[\begin{array}{cc}\frac{1}{2}\left(\nabla D_{11}\nabla w_1-v \nabla w_1\right)& \nabla D_{12}\nabla w_2\\ 0& \frac{1}{2}\left(\nabla D_{22}\nabla w_2-v \nabla w_2\right)\end{array}\right]+\left[\begin{array}{cc}\frac{1}{2}\left(\nabla D_{11}\nabla w_1-v \nabla w_1\right)& 0\\ \nabla D_{21}\nabla w_1& \frac{1}{2}\left(\nabla D_{22}\nabla w_2-v \nabla w_2\right)\end{array}\right] \\ (A_1)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}(A_2) \end{gathered}$$

(15)

In this case, A_iw is a three – order tensor (2×2×2) which can be written synthetically as,

$$A_{iw}=\left[\frac{\partial A_i}{\partial w_1}⌈ \frac{\partial A_i}{\partial w_2}\right], i=1, 2,$$

(16)

where ${\partial A}_i/{\partial w}_1, {\partial A}_i/{\partial w}_2,$ are 2 × 2 matrices. The products A_2w A₁ and A_1w A₂ are also three – order tensors. To evaluate the Lie bracket [A₁ , A₂], the matrices

$$\frac{\partial A_2}{\partial w_1} A_1, \frac{\partial A_2}{\partial w_2} A_1, \frac{\partial A_1}{\partial w_1} A_2, \frac{\partial A_1}{\partial w_2} A_2.$$

were computed. Denote by a₁₁, a₁₂, a₂₁ and a₂₂ the elements of matrix A. The matrices

$$\frac{\partial A_2}{\partial w_1} A_1=\left[\begin{array}{cc}\frac{1}{2}\left(\nabla \frac{\partial D_{11}}{\partial w_1} \nabla w_1+\nabla D_{11} \nabla -v \nabla \right) \frac{1}{2} a_{11}& \frac{1}{2}\left(\nabla \frac{\partial D_{11}}{\partial w_1} \nabla w_1+\nabla D_{11} \nabla -v \nabla \right) a_{12}\\ \left(\nabla \frac{\partial D_{21}}{\partial w_1} \nabla w_1+\nabla D_{21} \nabla \right) \frac{1}{2} a_{11}& \left(\nabla \frac{\partial D_{21}}{\partial w_1} \nabla w_1+\nabla D_{21} \nabla \right) a_{12}+\frac{1}{2} \left(\nabla \frac{\partial D_{22}}{\partial w_1} \nabla w_2\right) \frac{1}{2} a_{22}\end{array}\right]$$

and

$$\frac{\partial A_1}{\partial w_1} A_2= \left[\begin{array}{cc}\frac{1}{2}\left(\nabla \frac{\partial D_{11}}{\partial w_1} \nabla w_1+\nabla D_{11} \nabla -v \nabla \right) \frac{1}{2} a_{11}+\left(\nabla \frac{\partial D_{12}}{\partial w_1} \nabla w_2\right) a_{21}& \left(\nabla \frac{\partial D_{12}}{\partial w_1} \nabla w_2\right) \frac{1}{2}{a}_{22}\\ \frac{1}{2}\left(\nabla \frac{\partial D_{22}}{\partial w_1} \nabla w_2\right) a_{21}& \left(\nabla \frac{\partial D_{22}}{\partial w_1} \nabla w_2\right)\frac{1}{2}{a}_{22}\end{array}\right]$$

show that [A₁,A₂] ≠ 0. Thus, the splitting defined by (12) and (15) is only first-order accurate.

For the cases when the number of chemical species n is, n > 3, the computation of the Lie bracket leads to similar results. Only the elements from the main diagonal of the matrices have additional terms. The other elements of the matrices are similar to those previously presented. It must be also mentioned that the following relation,

D_i,_j ≠ D_j,_i, i ≠ j

is generally valid for the Fick diffusion coefficients. Under these conditions, it is obvious that the Lie bracket does not vanish when n > 3.

A second – order accurate splitting, regardless the values of the Lie bracket [A₁,A₂], is the Strang splitting [2]. The Strang splitting is defined by,

w ( t₁ ) = S₁ ( δ t /2 ) S₂ ( δ t ) S₁ ( δ t /2 ) w₀

(17)

The relation that describes the Strang splitting is a palindromic structure (Bou-Rabee and Sanz-Serna, [23]). For the Strang splitting, the solution operators S, S₁ and S₂ satisfy the relation [22,23],

$${S_1\left(\delta t/2\right) S}_2\left(\delta t\right) S_1\left(\delta t/2\right)-S\left(\delta t\right)=\frac{{\delta t}^3}{12}\left[A_1 , \left[A_1 ,A_2\right]\right]+\frac{{\delta t}^3}{24}\left[A_2 , \left[A_1 ,A_2\right]\right]+\mathcal{O} \left({\delta t}^4\right)$$

(18)

Relation (18) shows that Strang splitting is a second-order integrator. If the operator A is split in many pieces (m pieces), it can be proved by induction [23,24], that the Strang splitting

w (t₁) = S₁(δ t /2) S₂(δ t /2)…..S_m-1 (δ t /2) S_m (δ t) S_m-1 (δ t /2) …. S₂ (δ t /2) S₁ (δ t /2) w₀

(19)

is second – order accurate. However, it must be mentioned that the increase in the number of pieces increases the splitting error. The present splitting algorithm (6) belongs to the classes of algorithms defined by (19).

3. Test Problem

The present test problem is the same as in [20]. It models the mass transfer inside a sphere of constant radius R when the resistance to the mass transfer outside the sphere is negligible (the sphere surface concentrations are equal to the free – stream concentrations), [25,26,27]. In the absence of chemical reactions and neglecting the Marangoni and buoyancy effects, the dimensionless convection - diffusion equations read as follows:

$$\frac{\partial W_i}{\partial \tau }+\frac{Pe}{2}\left(V_R\frac{\partial W_i}{\partial r}+\frac{V_{\theta }}{r}\frac{\partial W_i}{\partial \theta }\right)={\sum }_{k=1}^{n-1}{\nabla D}_{ik}\nabla W_k, i=1, 2, \dots , n\text{-}1$$

(20)

where

$$\begin{gathered} \nabla D_{ik}\nabla =\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{D}_{ik}\frac{\partial }{\partial r}\right)+\frac{1}{r^2\mathit{sin}\theta }\frac{\partial }{\partial \theta }\left(\mathit{sin}\theta D_{ik}\frac{\partial }{\partial \theta }\right), \\ Pe=\frac{d U_0}{D_{ref}} , \tau =\frac{t D_{ref}}{R^2} , \end{gathered}$$

d is the diameter of the sphere, D_ref the reference diffusion coefficient, r the dimensionless radial coordinate, t the time, U₀ the free stream velocity, W the dimensionless concentration of the chemical species and θ the polar angle. For the computation of the dimensionless velocities field (V_R, V_θ) we assume incompressible, steady, axisymmetric flow inside and outside the sphere and constant physical properties.

The boundary conditions are:

- θ = 0, π

$${\sum }_{j=1}^{n-1}{D}_{ij}\frac{1}{r}\frac{\partial W_j}{\partial \theta }=0 ⟺ \frac{\partial W_i}{\partial \theta }=0, i=1,2,\dots , n\text{-}1$$

(21a)

- r = 1

$$W_i=W_i^R, i=1,2,\dots , n\text{-}1$$

(21b)

- r = 0

W_i = finite, i = 1,2,…, n-1

(21c)

The dimensionless initial conditions are:

$$\mathsf{\tau } = 0, W_{\mathrm{i}} = W_i^0, i=1,2,\dots , n\text{-}1$$

(22)

Applying the (well-known) transformation,

$$W_i\hspace{0.17em}=\hspace{0.17em}\frac{Z_i}{r}, i=1,2,\dots , n\text{-}1$$

the mathematical model equations become:

$$\frac{\partial Z_i}{\partial \tau }+\frac{Pe}{2}\left[V_R\left(\frac{\partial Z_i}{\partial r}-\frac{Z_i}{r}\right)+\frac{V_{\theta }}{r}\frac{\partial Z_i}{\partial \theta }\right]={\sum }_{k=1}^{n-1}{\widehat{\nabla } D}_{ik}\widehat{\nabla }{Z}_k, i=1,2,\dots , n\text{-}1$$

(23)

where

$$\widehat{\nabla } D_{ik}\widehat{\nabla }=\frac{\partial }{\partial r}\left( D_{ik}\frac{\partial }{\partial r}\right)-\frac{1}{r} \frac{\partial D_{ik}}{\partial r}+\frac{1}{r^2} \left[ \frac{\partial }{\partial \theta } \left( D_{ik} \frac{\partial }{\partial \theta }\right)+\mathit{cot}\theta D_{ik}\frac{\partial }{\partial \theta }\right],$$

The boundary and initial conditions are:

- θ = 0, π

$$\frac{\partial \hspace{0.17em}{Z}_i}{\partial \hspace{0.17em}\theta }\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}0, i=1,2,\dots , n\text{-}1$$

(24a)

- r = 1

$$Z_i=Z_i^R, i=1,2,\dots , n\text{-}1$$

(24b)

- r = 0

Z_i = 0, i = 1,2,…, n-1

(24c)

-τ = 0

$${\mathrm{Z}}_{\mathrm{i}}=Z_i^0, i=1,2,\dots , \mathrm{n}\text{-}1$$

(24d)

4. Results and Discussions

Uniform meshes with N × N points,

$$\begin{gathered} 0=r_1<r_2<\mathrm{...}<r_{N-1}<r_N=1,r_k=\hspace{0.17em}(k\hspace{0.17em}-\hspace{0.17em}1\hspace{0.17em})\hspace{0.17em}{h}_r,0={\theta }_1<{\theta }_2<\mathrm{...}<{\theta }_{N-1}<{\theta }_N=\pi , \\ {\theta }_l=\hspace{0.17em}(\hspace{0.17em}l\hspace{0.17em}-\hspace{0.17em}1\hspace{0.17em})\hspace{0.17em}{h}_{\theta },\hspace{0.33em}k,l=\overline{1,N}, \end{gathered}$$

h_r = 1 / ( N – 1 ) and h_θ = π / ( N – 1 ) being the mesh steps size, were used for the discretization of the spatial derivatives of equations (23). The diffusion terms were discretized with the second-order finite difference scheme as follows,

$$\begin{array}{c}\frac{\partial }{\partial r}\left( D_{ik}\frac{\partial Z_k}{\partial r}\right)+\frac{1}{r^2} \frac{\partial }{\partial \theta } \left( D_{ik} \frac{\partial Z_k}{\partial \theta }\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\cong \frac{D_{i,j}^{k+\frac{1}{2},l}\left(Z_j^{k+1,l}-Z_j^{k,l}\right)-D_{i,j}^{k-\frac{1}{2},l}\left(Z_j^{k,l}-Z_j^{k-1,l}\right)}{h_r^2}+\frac{1}{r^2} \frac{D_{i,j}^{k,l+\frac{1}{2}}\left(Z_j^{k,l+1}-Z_j^{k,l}\right)-D_{i,j}^{k,l-\frac{1}{2}}\left(Z_j^{k,l}-Z_j^{k,l-1}\right)}{h_{\theta }^2}\end{array}$$

(25)

The diffusion coefficients from (25) are the arithmetic averages of the mesh point values. Appendix A shows the methodology used to calculate the multi-component Fick diffusion coefficients. Appendix B shows the Stefan – Maxwell diffusion coefficients and the choice of the reference diffusion coefficients. The convection term and the boundary conditions (24a) were also discretized with the central, second-order accurate scheme.

The time integration of each fractional step from (6) was done by an implicit method. In each fractional step, a non-linear, tri-diagonal matrix system must be solved. The successive substitution method solves the non-linear system. The following error criterion was used for the non-linear iteration: the difference between the numerical solutions of two consecutive non-linear iterations, in the discrete L_∞ norm, must be smaller than 10^-8.

The values considered for the dimensionless parameter Pe are, 0 ≤ Pe ≤ 10³. Only the creeping flow solutions [28,29] for the dimensionless velocity field were used. Two mass transfer cases were considered in this work. The first is the usual mass transfer for which $W_k^0$ ≠ $W_k^R$ for all k. The second is the so-called „Multi-component effects without an imposed concentration difference of a particular component”, [27]. In this case, for at least one of the chemical species, the condition $W_k^0$ = $W_k^R$ holds. This case is discussed in detail in [20].

The error for the discrete approximation is the first issue analysed in this section. The following methodology is used to investigate the accuracy of the spatial discretization:

• compute the numerical solutions on the following spatial meshes: 33 × 33, 65 × 65, 129 × 129 and 257 × 257, numbered by m = 1, 2, 3, 4;
• denote by Z^(m-1) and Z^(m) the solutions computed on the two consecutive meshes m – 1 and m, respectively; compute the difference between Z^(m-1) and Z^(m) by,

  $$Er{r}^{(m)}\hspace{0.17em}=\hspace{0.17em}{‖\hspace{0.17em}{Z}^{(m)}\hspace{0.17em}-\hspace{0.17em}{Z}^{(m-1)}\hspace{0.17em}‖}_p$$

  (26)

  where ${‖,‖}_p$is the discrete p norm; the difference between Z^(m-1) and Z^(m) is computed for the common grid points;
• the spatial convergence rate is computed as,

  $$R^{(m)}\hspace{0.17em}\approx \hspace{0.17em}{\mathit{log}}_2\hspace{0.17em}\left(\hspace{0.17em}\frac{Er{r}^{(m-1)}}{Er{r}^{(m)}}\right)\hspace{0.17em}.$$

  (27)

From the numerical simulations made in order to compute the spatial convergence rate, a selection is presented in Table 1 and Table 2. The discrete Euclidean norm (p = 2) was used in relation (26). The time step used to compute the results presented in Table 1 and Table 2 was constant and equal to ∆ τ = 10^-4. Similar results were obtained for values of the dimensionless time τ smaller or greater than 0.085 (Table 1) and 0.075 (Table 2). The numerical simulations made and the data presented in tables 1 and 2 show that:

• R^(m)≈ 2 for all Pe values used in this work; this means second – order spatial accuracy for the present numerical solutions;
• the number of chemical species, n, does not affect the convergence rate;
• the values of $W_k^0$, $W_k^R$and the Stefan - Maxwell diffusion coefficients do not influence the convergence rate.

The error for the time discretization was investigated following the approach presented in [10], (because the present test problem does not have an analytical solution). To apply the approach presented in [10], a reference solution must be first calculated. The reference solution was computed numerically as “the Richardson extrapolation of two solutions obtained with the smallest ∆ τ values”,[10]. The error used to approximate the convergence rate is the discrete Euclidean norm of the difference between the current solution and the reference solution.

Figure 1 shows one case selected from the present computations. This selection can be viewed as a typical one for the present results. The spatial mesh has 129 × 129 points. Figure 1 shows that, as expected, the increase in Pe increases the numerical error. However, regardless the Pe numbers values, the numerical solutions converge with the decrease in the time step and are second – order accurate in time.

A stable solution was obtained in all numerical simulations for all the combinations of parameters, initial and boundary values for W_k. The numerical simulations also shown that the number of chemical species n and the Stefan - Maxwell diffusion coefficients do not affect the stability of the present algorithm.

To verify the splitting error, the present test problem was also solved numerically using the fully implicit method from [30,31] combined with the defect – correction iteration [32] (in order to achieve second-order accuracy). The values of the relative differences between the present results (average concentrations) and those provided by the fully implicit method are smaller than 1%.

Figure 2 and Figure 3 show the time variation of the present numerical solutions. The quantity depicted in these figures is the sphere dimensionless average (geometric average) concentration. The cases presented in figures 2 and 3 refer only to usual mass transfer, i.e. mass transfers from the surrounding fluid into the sphere when $W_k^R$> $W_k^0$ and mass transfers from the sphere to the surrounding fluid when $W_k^R$< $W_k^0$. Time variation of the dimensionless concentrations for “mass transfer without an imposed concentration difference of a particular component”, [27], is not presented here. This case was fully discussed in [20].

Figure 2 and Figure 3 show that the dimensionless sphere average concentrations vary monotonously, from the initial value to the surface value. As expected, the mass transfer rate increases with the increase in Pe. The effect of Pe is stronger in the range 100 ≤ Pe ≤ 1000. The results presented in figures 2 and 3 show that the numerical solutions converge to the correct solution.

5. Conclusions

In this work, the splitting algorithm developed in [20,21] for linear convection – diffusion equations was extended to non-linear convection – diffusion equations. For the present algorithm, only block-triangular algebraic systems must be solved. Each block contains a tridiagonal system. The stability and accuracy of this algorithm were analysed both theoretically and numerically.

Theoretically, we showed that a second-order splitting error can be provided by Strang splitting. Also, the numerical simulations made have shown that the discrete approximation is stable and second – order accurate in space and time.

The splitting techniques applied to multi-component mass transfer equations can be further developed. First, the chemical reaction term must be added. An operator (physical) splitting can be applied to the convection-diffusion-reaction equations. The present splitting algorithm can be applied to the diffusion step.