Prerpints.org logo
Preprint
Article

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

Altmetrics

Downloads

125

Views

31

Comments

0

Submitted:

23 March 2023

Posted:

24 March 2023

You are already at the latest version

Alerts
Abstract
The splitting method developed in a previous study for linear, multi – component, convection-diffusion mass transfer equations is extended to non-linear, multi-component, convection-diffusion mass transfer equations. The diffusion term is non-linear (the diffusion coefficients depend on chemical species concentrations). Theoretical results about the splitting error are presented, as well as numerical simulations for mixtures with n ≤ 5. The numerical experiments show second-order accuracy and stability for the present splitting algorithm.
Keywords: 
Subject: Computer Science and Mathematics  -   Mathematics
2010 Mathematics Subject Classification: 80A20, 80A32

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, ∞), Ω ⊂ RN, are as follows:
∂ w 1 ∂ t + v   · ∇ w 1 = ∑ k = 1 n − 1 ∇ · D 1 , k w ∇   w k
∂ w n − 1 ∂ t + v · ∇ w n − 1 = ∑ k = 1 n − 1 ∇ · D n − 1 , k w ∇ w k
where wk is the mass concentration of component k,  w = w 1 , w 2 , . . . . , v is the velocity vector, D i , k   w 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 Di,k.
Discretizing the spatial derivatives from (1), one obtains,
∂ w 1 ∂ t + Λ h w 1 = ∑ k = 1 n − 1 ∇ h   D 1 , k ∇ h w k
∂ w n − 1 ∂ t + Λ h w n − 1 = ∑ k = 1 n − 1 ∇ h D n − 1 , k ∇ h w k
where ∇h and Λh are the discrete approximations of the continuous operators. The present splitting method is defined by:
∇ h   D 11 ∇ h − Λ h ∇ h   D 12 ∇ h ⋯ ∇ h   D 1 , n − 1 ∇ h ∇   h D 21 ∇ h ∇   h D 22 ∇ h − Λ h ⋯ ∇ h D 2 , n − 1 ∇ h ⋮ ⋮ ⋱ ⋮ ∇ h D n − 1,1 ∇ h ⋯ ⋯ ∇ h D n − 1 , n − 1 ∇ h − Λ h = 1 2 ∇ h D 11 ∇ h − Λ h ∇ h D 12 ∇ h ⋯ ∇ h D 1 , n − 1 ∇ h 0 1 2 ∇ h D 22 ∇ h − Λ h ⋯ ∇ h D 2 , n − 1 ∇ h ⋮ ⋮ ⋱ ⋮ 0 0 0 1 2 ∇ h D n − 1 , n − 1 ∇ h − Λ h + 1 2 ∇ h D 11 ∇ h − Λ h 0 ⋯ 0 ∇ h D 21 ∇ h 1 2 ∇ h D 22 ∇ h − Λ h 0 0 ⋮ ⋮ ⋱ ⋮ ∇ h D n − 1,1 ∇ h ⋯ ⋯ 1 2 ∇ h D n − 1 , n − 1 ∇ h − Λ h
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. RN = R2. Consider the geometric splitting of ∇h and Λh as:
∇ h = ∇ h , x + ∇ h , y , Λ h = Λ h , x + Λ h , y
and denote by L h , x s ,   L h , y s ,   L h , x i ,   L h , y i , the matrices
L h , x s = 1 2 ∇ h , x D 11 ∇ h , x − Λ h , x ∇ h , x D 12 ∇ h , x ⋯ ∇ h , x D 1 , n − 1 ∇ h , x 0 1 2 ∇ h , x D 22 ∇ h , x − Λ h , x ⋯ ∇ h , x D 2 , n − 1 ∇ h , x ⋮ ⋮ ⋱ ⋮ 0 ⋯ 0 1 2 ∇ h , x D n − 1 , n − 1 ∇ h , x − Λ h , x , L h , y s = 1 2 ∇ h , y D 11 ∇ h , y − Λ h , y ∇ h , y D 12 ∇ h , y ⋯ ∇ h , y D 1 , n − 1 ∇ h , y 0 1 2 ∇ h , y D 22 ∇ h , y − Λ h , y ⋯ ∇ h , y D 2 , n − 1 ∇ h , y ⋮ ⋮ ⋱ ⋮ 0 ⋯ 0 1 2 ∇ h , y D n − 1 , n − 1 ∇ h , y − Λ h , y L h , x i = 1 2 ∇ h , x D 11 ∇ h , x − Λ h , x 0 ⋯ 0 ∇ h , x D 21 ∇ h , x 1 2 ∇ h , x D 22 ∇ h , x − Λ h , x ⋯ 0 ⋮ ⋮ ⋱ ⋮ ∇ h , x D n − 1,1 ∇ h , x . . . . ⋯ 1 2 ∇ h , x D n − 1 , n − 1 ∇ h , x − Λ h , x L h , y i = 1 2 ∇ h , y D 11 ∇ h , y − Λ h , y 0 ⋯ 0 ∇ h , y D 21 ∇ h , y 1 2 ∇ h , y D 22 ∇ h , y − Λ h , y ⋯ 0 ⋮ ⋮ ⋱ ⋮ ∇ h , y D n − 1,1 ∇ h , y ⋯ ⋯ 1 2 ∇ h , y D n − 1 , n − 1 ∇ h , y − Λ h , y
One step for time integration is the Strang integrator with four pieces, [2]:
∂ w 1 ∂ t + L h , x s w 1 = 0 , t i − 1 ≤ t ≤ t i − 1 / 2
∂ w 2 ∂ t + L h , y s w 2 = 0 , t i − 1 ≤ t ≤ t i − 1 / 2
∂ w 3 ∂ t + L h , x i w 3 = 0 , t i − 1 ≤ t ≤ t i − 1 / 2
∂ w 4 ∂ t + L h , y i w 4 = 0 , t i − 1 ≤ t ≤ t i
∂ w 5 ∂ t + L h , x i w 5 = 0 , t i − 1 / 2 ≤ t ≤ t i
∂ w 6 ∂ t + L h , y s w 6 = 0 , t i − 1 / 2 ≤ t ≤ t i
∂ w 7 ∂ t + L h , x s w 7 = 0 , t i − 1 / 2 ≤ t ≤ t i
where w = w 1 , w 2 , . . . . . The initial conditions are: w 1 ( ti-1 ) = w ( ti-1 ), w k ( ti-1 ) = w k − 1 ( ti-1/2 ) ( for k = 2, 3, 4 ) and w k ( ti-1/2 ) = w k − 1 ( ti-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,
∂   w ∂   t = A   w ,
with the initial condition,
w (x, t0) = w0, x ∈ Rn.
where A is a smooth function of w. For well-posed problems, the solution of (7) at any time t1 > t0, (t1 = t0 + δ t) is uniquely determined by the initial condition. This solution can be written as,
w ( t1 ) = S ( t1, t0 ) w0 = S ( δ t ) w0.
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 ) = A1 (w) + A2 ( w ).
where A1 and A2 are smooth functions of w. The resulted subproblems
∂   w ∂   t = A 1   w ,
and
∂   w ∂   t = A 2   w ,
should be easier to solve than the full problem. Denote by S1 ( δ t ) and S2 ( δ t ) the solution operators for the subproblems (11a) and (11b). The usual or Lie-Trotter splitting algorithm corresponding to (11) is,
w ( t1 ) = S2 ( δ t ) S1 ( δ t ) w0
For the algorithm (12), the splitting error operator Esplit ( δ t ) is defined by,
Esplit ( δ t ) = S2 ( δ t ) S1 ( δ t ) – S ( δ t)
Using Taylor expansions, one can prove (see LeVeque, [22]) that the splitting error for the algorithm (12) verifies the relation:
E s p l i t   δ t =   S 2   δ t   S 1   δ t − S   δ t =   δ t 2 2   A 1   ,   A 2 +   O   δ t 3
where [A1, A2] is the Lie bracket (or Lie-Jacobi bracket or commutator) defined by
A 1 ,   A 2   w = A 2 w   w   A 1 w − A 1 w w   A 2 w
and Aiw is
A i w = ∂   A i ∂   w ,   i = 1 ,   2 .
Relation (13) shows that if [A1,A2] = 0, the Lie-Trotter splitting is second order accurate. If [A1,A2] 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,
∇   D 11 ∇   w 1 − v   ∇   w 1 ∇   D 12 ∇   w 2 ∇   D 21 ∇   w 1 ∇   D 22 ∇   w 2 − v   ∇   w 2 = ( A ) 1 2 ∇   D 11 ∇   w 1 − v   ∇   w 1 ∇   D 12 ∇   w 2 0 1 2 ∇   D 22 ∇   w 2 − v   ∇   w 2 + 1 2 ∇   D 11 ∇   w 1 − v   ∇   w 1 0 ∇   D 21 ∇   w 1 1 2 ∇   D 22 ∇   w 2 − v   ∇   w 2 ( A 1 )                                       ( A 2 )
In this case, Aiw is a three – order tensor (2×2×2) which can be written synthetically as,
A i w = ∂   A i ∂   w 1   ∂   A i ∂   w 2 ,   i = 1 ,   2 ,
where   ∂   A i / ∂   w 1 ,     ∂   A i / ∂   w 2 , are 2 × 2 matrices. The products A2w A1 and A1w A2 are also three – order tensors. To evaluate the Lie bracket [A1 , A2], the matrices
∂   A 2 ∂   w 1   A 1 ,     ∂   A 2 ∂   w 2   A 1 ,     ∂   A 1 ∂   w 1   A 2 ,     ∂   A 1 ∂   w 2   A 2 .  
were computed. Denote by a11, a12, a21 and a22 the elements of matrix A. The matrices
∂   A 2 ∂   w 1   A 1 = 1 2 ∇   ∂   D 11 ∂   w 1   ∇   w 1 + ∇   D 11   ∇ − v   ∇   1 2   a 11 1 2 ∇   ∂   D 11 ∂   w 1   ∇   w 1 + ∇   D 11   ∇ − v   ∇   a 12 ∇   ∂   D 21 ∂   w 1   ∇   w 1 + ∇   D 21   ∇   1 2   a 11 ∇   ∂   D 21 ∂   w 1   ∇   w 1 + ∇   D 21   ∇   a 12 + 1 2   ∇   ∂   D 22 ∂   w 1   ∇   w 2   1 2   a 22
and
∂   A 1 ∂   w 1   A 2 =   1 2 ∇   ∂   D 11 ∂   w 1   ∇   w 1 + ∇   D 11   ∇ − v   ∇   1 2   a 11 + ∇   ∂   D 12 ∂   w 1   ∇   w 2   a 21 ∇   ∂   D 12 ∂   w 1   ∇   w 2   1 2 a 22 1 2 ∇   ∂   D 22 ∂   w 1   ∇   w 2   a 21 ∇   ∂   D 22 ∂   w 1   ∇   w 2 1 2 a 22
show that [A1,A2] ≠ 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,
Di,j ≠ Dj,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 [A1,A2], is the Strang splitting [2]. The Strang splitting is defined by,
w ( t1 ) = S1 ( δ t /2 ) S2 ( δ t ) S1 ( δ t /2 ) w0
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, S1 and S2 satisfy the relation [22,23],
S 1 δ t / 2   S 2 δ t   S 1 δ t / 2 − S δ t = δ t 3 12 A 1   ,   A 1   , A 2 + δ t 3 24 A 2   ,   A 1   , A 2 + O   δ t 4
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 (t1) = S1(δ t /2) S2(δ t /2)…..Sm-1 (δ t /2) Sm (δ t) Sm-1 (δ t /2) …. S2 (δ t /2) S1 (δ t /2) w0
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:
∂ W i ∂ τ + P e 2 V R ∂ W i ∂ r + V θ r ∂ W i ∂ θ = ∑ k = 1 n − 1 ∇   D i k ∇   W k ,   i = 1 ,   2 ,   … ,   n - 1
where
∇   D i k ∇ = 1 r 2 ∂ ∂ r r 2 D i k ∂ ∂ r + 1 r 2 sin ⁡ θ ∂ ∂ θ sin ⁡ θ   D i k ∂ ∂ θ , P e = d   U 0 D r e f   ,   τ = t   D r e f R 2   ,
d is the diameter of the sphere, Dref the reference diffusion coefficient, r the dimensionless radial coordinate, t the time, U0 the free stream velocity, W the dimensionless concentration of the chemical species and θ the polar angle. For the computation of the dimensionless velocities field (VR, Vθ) we assume incompressible, steady, axisymmetric flow inside and outside the sphere and constant physical properties.
The boundary conditions are:
- θ = 0, π
∑ j = 1 n − 1 D i j 1 r ∂ W j ∂ θ = 0   ⟺   ∂   W i ∂   θ = 0 ,   i = 1 , 2 , … ,   n - 1
- r = 1
W i = W i R ,   i = 1 , 2 , … ,   n - 1
- r = 0
Wi = finite, i = 1,2,…, n-1
The dimensionless initial conditions are:
τ   =   0 ,   W i   =   W i 0 ,   i = 1 , 2 , … ,   n - 1
Applying the (well-known) transformation,
W i   =   Z i r ,   i = 1 , 2 , … ,   n - 1
the mathematical model equations become:
∂ Z i ∂ τ + P e 2 V R ∂ Z i ∂ r − Z i r + V θ r ∂ Z i ∂ θ = ∑ k = 1 n − 1 ∇ ^   D i k ∇ ^ Z k ,   i = 1 , 2 , … ,   n - 1
where
∇ ^   D i k ∇ ^ = ∂ ∂   r   D i k ∂ ∂   r − 1 r   ∂   D i k ∂   r + 1 r 2     ∂ ∂   θ     D i k   ∂ ∂   θ + cot ⁡ θ D i k ∂ ∂ θ ,
The boundary and initial conditions are:
- θ = 0, π
∂   Z i ∂   θ     =   0 ,   i = 1 , 2 , … ,   n - 1
- r = 1
Z i = Z i R ,   i = 1 , 2 , … ,   n - 1
- r = 0
Zi = 0, i = 1,2,…, n-1
-τ = 0
Z i = Z i 0 ,   i = 1 , 2 , … ,   n - 1

4. Results and Discussions

Uniform meshes with N × N points,
0 = r 1 < r 2 < ... < r N − 1 < r N = 1 , r k =   ( k   −   1   )   h r , 0 = θ 1 < θ 2 < ... < θ N − 1 < θ N = π , θ l =   (   l   −   1   )   h θ ,   k , l = 1 , N ¯ ,
hr = 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,
∂ ∂   r   D i k ∂   Z k ∂   r + 1 r 2   ∂ ∂   θ     D i k   ∂ Z k ∂   θ               ≅   D i , j k + 1 2 , l Z j k + 1 , l − Z j k , l − D i , j k − 1 2 , l Z j k , l − Z j k − 1 , l h r 2 + 1 r 2   D i , j k , l + 1 2 Z j k , l + 1 − Z j k , l − D i , j k , l − 1 2 Z j k , l − Z j k , l − 1 h θ 2
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 ≤ 103. 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,
    E r r ( m )   =     Z ( m )   −   Z ( m − 1 )   p
    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 )   ≈   log 2     E r r ( m − 1 ) E r r ( m )   .
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 Wk. 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.

Appendix A.

The following steps must be followed in order to compute D (mass-base diffusion coefficients matrix with the elements Di,j), [33]:
1. Compute the matrix B (the so-called drag effects matrix) by,
B i i   =   x i D i n   +   ∑ k   =   1 k   ≠   i n   x k D i k           i   =   1 , .... , n − 1
B i j   =   −   x i     1 D i j   −   1 D i n             i , j   =   1 , .... , n − 1 ,     i   ≠   j
In the previous relations, D i j is the Stefan – Maxwell diffusion coefficient and xi the molar fraction.
2. Compute the matrix DM by,
DM = B-1 Γ
where Γ (the matrix of the thermodynamic effects) is equal to the identity matrix for ideal mixtures. DM is the mole-based diffusion coefficients matrix.
3. Compute the matrix D by,
D = C [ w ] [ x ]-1 DM [ x ] [ w ]-1 C-1
where [ x ] and [ w ] are diagonal matrices whose the elements on the main diagonal are the mole fractions xi and the mass fractions wi, respectively. The matrix C is computed as
C i j   =   δ i j   −   w i     x j w j   −   x n w n   .
where δij is the Kronecker symbol.
In the present computations, one assumed that the molar mass of the chemical species is the same. In this situation, xi = wi and D = DM. Also, constant values for the Stefan – Maxwell diffusion coefficients were assumed.

Appendix B.

In section 4, the Stefan – Maxwell diffusion coefficients used are:
- n = 4;
set   Q 1 D 12 = 0.647 , D 13 = 0.415 , D 14 = 1.0 ,   D 23 = 0.628 , D 24 = 1.4897 , D 34 = 0.969 . set   Q 2 D 12 = 0.915 , D 13 = 0.782 , D 14 = 1.0 ,   D 23 = 972 , D 24 = 1.344 , D 34 = 1.165 .
The diffusion coefficient D 14 is equal to D 14 = 1 because D 14 was considered the reference diffusion coefficient.
- n = 5;
set   P 1 D 12   =   0.22 ,         D 13   =   0.31 ,         D 14   =   0.25 ,         D 15   =   1.0 ,         D 23   =   0.35 ,         D 24   =   0.1 ,             D 25   =   1.18 ,   D 34   =   0.43 ,         D 35   =   1.2 ,       D 45   =   1.3 .       set   P 2 D 12   =   6.326 ,         D 13   =   4.2523 ,         D 14   =   0.247 ,         D 15   =   1.0 ,       D 23   =   4.978 ,               D 24   =   0.189 ,         D 25   =   5.178 ,   D 34   =   3.325 ,         D 35   =   1.27 ,       D 45   =   0.987 .    
The diffusion coefficient D 15 is equal to D 15 = 1 because, for n = 5, D 15 is the reference diffusion coefficient. It must be mentioned that the following relation
D i j = D j i , i ≠ j ,
holds for the Stefan – Maxwell diffusion coefficients,

References

  1. N.N. Yanenko, The Method of Fractional Steps, Vol. 1, Springer, New York, 1967.
  2. G. Strang, On the construction and comparison of difference schemes, SIAM J. Num. Anal. 5 (3) (1968) 506 – 517.
  3. G.I. Marchuk, Splitting and Alternating Direction Methods, in Handbook of Numerical Analysis Vol. I, Ed. P.G. Ciarlet and J.L. Lions, pp. 199 – 462, Elsevier, North – Holland, 1990.
  4. R.I. McLachlan, G. Reinout W. Quispel, Splitting methods, Acta Numerica 11 (2002) 341 – 434.
  5. D. Lanser, J.G. Verwer, Analysis of operator splitting for advection-diffusion-reaction problems from air pollution modelling, J. Comput. Appl. Math. 111 (1999) 201 – 216. [CrossRef]
  6. S. Descombes, Convergence of a splitting method of high order for reaction-diffusion systems, Math. Comput. 70 (2001) 1481-1501. [CrossRef]
  7. W. Hundsdorfer, J. Verwer, Numerical Solution of Advection-Diffusion-Reaction Equations, pp. 325 – 417, Springer, Berlin, 2003. [CrossRef]
  8. D.L. Ropp, J.N. Shadid, C.C. Ober, Studies on the accuracy of time-integration methods for the reaction-diffusion equations, J. Comput. Phys. 194 (2004) 544 – 574. [CrossRef]
  9. C.C. Ober, J.N. Shadid, Studies on the accuracy of time-integration methods for the radiation-diffusion equations, J. Comput. Phys. 195 (2004) 743 – 772. [CrossRef]
  10. D.L. Ropp, J.N. Shadid, Stability of operator splitting methods for systems with indefinite operators: reaction-diffusion systems, J. Comput. Phys. 203 (2005) 449 – 466. [CrossRef]
  11. D.L. Ropp, J.N. Shadid, Stability of operator splitting methods for systems with indefinite operators: Advection – diffusion – reaction systems, J. Comput. Phys. 228 (2009) 3508 – 3516. [CrossRef]
  12. B. Sportisse, A review of current issues in air pollution modeling and simulation, Comput. Geosci. 11 (2007) 159 – 181. [CrossRef]
  13. S. Descombes, M. Duarte, T. Dumont, F. Laurent, V. Louvet, M. Massot, Analysis of operator splitting in the nonasymptotic regime for nonlinear reaction-diffusion equations, SIAM J. Numer. Anal., 52 (2014) 1311-1334. [CrossRef]
  14. S. Descombes, M. Duarte, M. Massot, Operator splitting methods with error estimator and adaptive time-stepping. Application to the simulation of combustion phenomena, In Splitting Methods in Communication, Imaging, Science and Engineering, Ed. R. Glowinski, S. Osher. W. Yin, pp. 1-13, Springer, Berlin, 2015. [CrossRef]
  15. G. Bertoli, G. Vilmart, Strang splitting method for semilinear parabolic problems with inhomegeneous boundary conditions: a correction based on the flow of the nonlinearity, SIAM J. Sci. Comput. 42 (3) (2020). [CrossRef]
  16. C. Beck, S. Becker, P. Cheridito, A. Jentzen, A. Neufeld, Deep splitting method for parabolic PDEs, SIAM J. Sci. Comput. 43 (5) (2021). [CrossRef]
  17. C. Liu, C. Wang, Y. Wang, A structure -preserving splitting scheme for reaction-diffusion equations with detailed balance, J. Comput. Phys. 436 (2021) 110253. [CrossRef]
  18. S.R. de Groot and P. Mazur, Non-Equilibrium Thermodynamics, chap. XI, 2nd ed., Dover, New York, 1984.
  19. R. Taylor and R. Krishna, Multicomponent Mass Transfer, Wiley, New York, 1993.
  20. Gh. Juncu, A. Nicola, C. Popa, and E. Stroila, Numerical solution of the parabolic multi-component convection- diffusion mass transfer equations by a splitting method, Numer. Heat Transfer A 71 (1) (2017) 72 – 90. [CrossRef]
  21. Gh. Juncu, A. Nicola and C. Popa, Splitting methods for the numerical solution of multi-component mass transfer problems, Mathematics and Computers in Simulations 152 (2018) 1 – 14. [CrossRef]
  22. R.J. LeVeque, Time-split methods for partial differential equations, Ph.D. Thesis, Stanford University, 1982.
  23. N. Bou-Rabee, J.M. Sanz-Serna, Geometric integrators and the Hamiltonian Monte Carlo method, Acta Numerica 27 (2018) 113 – 206. [CrossRef]
  24. D. Gottlieb, Strang-type difference schemes for multidimensional problems, SIAM J. Numer. Anal. 36 (1981) 603-626. [CrossRef]
  25. S. Ubal, C.H. Harrison, P. Grassia, W.J. Korchinsky, Numerical simulation of mass transfer in circulating drops, Chem. Eng. Sci. 65, 2934 – 2956, 2010. [CrossRef]
  26. S. Ubal, P. Grassia, C.H. Harrison, W.J. Korchinsky, Numerical simulation of multi-component mass transfer in rigid or circulating drops: multi-component effects even in the presence of weak coupling, Colloids and Surfaces A: Physicochem. Eng. Aspects vol. 380, pp. 6 – 15, 2011. [CrossRef]
  27. S. Ubal, P. Grassia, C.H. Harrison, W.J. Korchinsky, Reprint of Numerical simulation of multi-component mass transfer in rigid or circulating drops: multi-component effects even in the presence of weak coupling, Colloids and Surfaces A: Physicochem. Eng. Aspects 382, 251 – 260, 2011. [CrossRef]
  28. R. Clift, J.R. Grace, M.E. Weber, Bubbles, Drops and Particles, chap. 3, Academic Press, New York, 1978.
  29. L.G. Leal, Advanced Transport Phenomena, chap. 7, Cambridge University Press, Cambridge, 2007.
  30. Gh. Juncu, C. Popa, G. Sarbu, On Numerical Solution of Nonlinear Parabolic Multicomponent Diffusion–Reaction Problems, Anal. Univ. Ovidius Constanta, Seria Matematica 29 (3) 183 – 200, 2021.
  31. Gh. Juncu, C. Popa, G. Sarbu, On Numerical Solution of Nonlinear Parabolic Multicomponent Convection-Diffusion–Reaction Problems, submitted for publication, 2022.
  32. P. Wesseling, Principles of Computational Fluid Dynamics, Springer, Berlin, 2001. [CrossRef]
  33. A. Leahy – Dios, A. Firoozabadi, Unified model for nonideal multicomponent molecular diffusion coefficients, A.I.Ch.E. J., 53, 2932 – 2939, 2007. [CrossRef]
Preprints 70103 g001
Figure 1. Analysis of the temporal convergence on a spatial mesh with hr = 1 / 128, hθ = π / 128; n = 5, matrix P1 and Pe in the range 10 ÷ 103; the final time is τ = 1.0; W 1 R = 0.2 , W 2 R = 0.05 , W 3 R = 0.2 , W 4 R = 0.1 , W 1 0 = 0.05 , W 2 0 = 0.3 , W 3 0 = 0.01 , W 4 0 = 0.3.
Figure 1. Analysis of the temporal convergence on a spatial mesh with hr = 1 / 128, hθ = π / 128; n = 5, matrix P1 and Pe in the range 10 ÷ 103; the final time is τ = 1.0; W 1 R = 0.2 , W 2 R = 0.05 , W 3 R = 0.2 , W 4 R = 0.1 , W 1 0 = 0.05 , W 2 0 = 0.3 , W 3 0 = 0.01 , W 4 0 = 0.3.
Preprints 70103 g001
Preprints 70103 g002aPreprints 70103 g002b
Figure 2. Bulk dimensionless concentrations vs. dimensionless time for a quaternary system at various Pe numbers and set Q2 of Stefan - Maxwell diffusion coefficients; W 1 R = 0.01, W 1 0 = 0.3, W 2 R = 0.3 ,   W 2 0 = 0.05   , W 3 R =   0.1 , W 3 0 =   0.05 ; (a) W ¯ 1 ( τ ); (b) W ¯ 2 ( τ ); (c) W ¯ 3 ( τ ).
Figure 2. Bulk dimensionless concentrations vs. dimensionless time for a quaternary system at various Pe numbers and set Q2 of Stefan - Maxwell diffusion coefficients; W 1 R = 0.01, W 1 0 = 0.3, W 2 R = 0.3 ,   W 2 0 = 0.05   , W 3 R =   0.1 , W 3 0 =   0.05 ; (a) W ¯ 1 ( τ ); (b) W ¯ 2 ( τ ); (c) W ¯ 3 ( τ ).
Preprints 70103 g002aPreprints 70103 g002b
Preprints 70103 g003aPreprints 70103 g003bPreprints 70103 g003c
Figure 3. Time evolution of the bulk dimensionless concentrations for n = 5, various Pe numbers, set P1 of Stefan - Maxwell diffusion coefficients; W 1 R = 0.1 , W 1 0 = 0.01 , W 2 R = 0.01 , W 2 0 = 0.1 , W 3 R = 0.01, W 3 0 = 0.2, W 4 R = 0.2, W 4 0 = 0.01; (a) W ¯ 1 ( τ ); (b) W ¯ 2 ( τ ); (c) W ¯ 3 ( τ ); (d) W ¯ 3 ( τ ).
Figure 3. Time evolution of the bulk dimensionless concentrations for n = 5, various Pe numbers, set P1 of Stefan - Maxwell diffusion coefficients; W 1 R = 0.1 , W 1 0 = 0.01 , W 2 R = 0.01 , W 2 0 = 0.1 , W 3 R = 0.01, W 3 0 = 0.2, W 4 R = 0.2, W 4 0 = 0.01; (a) W ¯ 1 ( τ ); (b) W ¯ 2 ( τ ); (c) W ¯ 3 ( τ ); (d) W ¯ 3 ( τ ).
Preprints 70103 g003aPreprints 70103 g003bPreprints 70103 g003c
Table 1. Analysis of the spatial convergence rate; n = 4, W 1 R = 0.02 , W 2 R = 0.3 , W 3 R = 0.01 , W 1 0 = 0.3 , W 2 0 = 0.02 , W 3 0 = 0.01 , τ = 0.085, set Q1 of Stefan – Maxwell diffusion coefficients.
Table 1. Analysis of the spatial convergence rate; n = 4, W 1 R = 0.02 , W 2 R = 0.3 , W 3 R = 0.01 , W 1 0 = 0.3 , W 2 0 = 0.02 , W 3 0 = 0.01 , τ = 0.085, set Q1 of Stefan – Maxwell diffusion coefficients.
m N Err(m) R(m)
Pe=1 Pe=10 Pe=102 Pe=103 Pe=1 Pe=10 Pe=102 Pe=103
1 33 - - - - - - - -
2 65 0.6×10-3 0.67×10-3 0.29×10-2 0.22×10-2 - - - -
3 129 0.16×10-3 0.18×10-3 0.79×10-3 0.57×10-3 1.9 1.89 1.87 1.96
4 257 0.45×10-4 0.7×10-4 0.21×10-3 0.15×10-3 1.85 1.87 1.89 1.95
Table 2. Analysis of the spatial convergence rate; n = 5, W 1 R = 0.2 , W 2 R = 0.05 , W 3 R = 0.2 , W 4 R = 0.1 , W 1 0 = 0.05 , W 2 0 = 0.3 , W 3 0 = 0.01 , W 4 0 = 0.3 , τ = 0.075, set P2 of Stefan-Maxwell diffusion coefficients.
Table 2. Analysis of the spatial convergence rate; n = 5, W 1 R = 0.2 , W 2 R = 0.05 , W 3 R = 0.2 , W 4 R = 0.1 , W 1 0 = 0.05 , W 2 0 = 0.3 , W 3 0 = 0.01 , W 4 0 = 0.3 , τ = 0.075, set P2 of Stefan-Maxwell diffusion coefficients.
m N Err(m) R(m)
Pe=1 Pe=10 Pe=102 Pe=103 Pe=1 Pe=10 Pe=102 Pe=103
1 33 - - - - - - - -
2 65 0.14×10-2 0.16×10-2 0.68×10-2 0.82×10-2 - - - -
3 129 0.36×10-3 0.42×10-3 0.18×10-2 0.21×10-2 1.97 1.96 1.95 1.95
4 257 0.93×10-4 0.11×10-3 0.46×10-3 0.56×10-3 1.96 1.95 1.95 1.93
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

© 2024 MDPI (Basel, Switzerland) unless otherwise stated