2.1. Description of the Dynamics of Complex Multicomponent Diffusion
Let us assume that in a microporous medium with a geometry defined by the shape factor
m, diffusion involving
K components occurs under isothermal conditions. The mechanisms contributing to mass transport in a complex diffusion process are: molecular diffusion and Knudsen diffusion in pores, viscous flow, and surface diffusion. The mathematical model of such diffusion for unsteady conditions is given by the system of equations [
24,
25]:
where
The system of Equation (1) represents the relationship between the changes in the amount of mass in the gas phase and those in the solid, and the molar flux densities in the pores
Np and by surface diffusion
Ns. In the first term of the right-hand side of Equation (1), the porosity of the solid ε does not appear, as this quantity is incorporated in the definitions of the effective diffusion coefficients. The mathematical model (1) was developed assuming local thermodynamic equilibrium between the adsorbed phase and the bulk phase in the pores at each value of the spatial coordinate
z. If the fugacity coefficients of the gas phase components
fi are close to unity, then the aforementioned equilibrium is given by the following relationships:
Surface diffusion occurs when the adsorption of components on the solid surface takes place. The sorption isotherms in the
K-component system are generally defined as:
where a fractional surface occupancy of component
i is defined as:
In light of the assumption of the thermodynamic equilibrium mentioned above, that is, when the concentrations in the gas phase and on the surface of the solid are unambiguously dependent on each other, the derivatives of surface concentrations with respect to time appearing in Equation (1) can be written as:
The term
J(
p) is the derivative of vector
q with respect to vector
p. It is therefore a Jacobi matrix, which is calculated using the information on sorption isotherms (5). Finally, the mathematical model of multicomponent complex diffusion can be expressed as follows:
where Jacobi matrix
J is defined as:
According to the dusty gas model [
16,
27], the vector of molar flux densities in the pores
Np, appearing in Equation (1), is given by the equation:
where the elements of matrices
Bp,
E and
p are defined according to Equations (11) and (12):
and
The effective molecular and Knudsen diffusion coefficients are calculated respectively as:
It is believed [
28] that the coefficients of molecular diffusion in pores
Dij may be calculated according to the correlation of Fuller and co-authors [
29,
30].
It is worth mentioning here that the method of estimating the values of Knudsen diffusion coefficients was the subject of detailed discussion in the literature [
31,
32,
33,
34,
35,
36,
37,
38,
39]. Krishna pointed out that the Knudsen diffusion coefficients can be correlated with Henry constant He for adsorption [
28]. They increase with increases in this constant. In the case of many industrial processes such as adsorption or chemical reactions in the porous catalysts, the dominant mechanism of mass transport in the pores may be Knudsen diffusion [
2,
3,
4,
5,
27]. Therefore, the contribution of this transport mechanism should be present in the general model (8). Knudsen diffusion is said to occur when collisions of molecules with pore walls are more frequent than with each other. The occurrence of the Knudsen mechanism is favored by low pressures. This is because then the mean free path of molecules increases.
The parameter
e is the so-called effective coefficient of permeability of the porous material:
It turns out that in multicomponent mixtures containing compounds with different molar masses, a total pressure gradient develops in the porous material [
26,
40]. This arises as a consequence of the difference in the Knudsen diffusion rates of the individual components. For this reason, the definition of molar flux density (10) needs to include an expression for viscous flow. This term introduces an additional mass transport mechanism which, unlike Knudsen diffusion, is not a separative mechanism with respect to individual substances.
An extremely significant step in the development of the Maxwell–Stefan formalism was its application to the modelling of the diffusive mass transport of species adsorbed on the internal surfaces of porous materials. A great contribution to the extension of the Maxwell–Stefan model to the description of surface diffusion was made by the works of Krishna [
19,
20,
21]. It was proved that the framework of the Maxwell–Stefan model holds well also for process simulations and calculations of surface diffusion.
According to the mentioned suggestions of Krishna, the mass transport due to surface diffusion can be described by a system of equations:
After assuming thermodynamic equilibrium (4), the expression for the vector of molar flux densities,
Ns is obtained:
where matrix
q* is given by:
and
The expression for
Ns similarly to
Np appears in Equation (1). The values of the thermodynamic factors are calculated from the following relationship:
It should be remarked that the information about the thermodynamic equilibrium (4) is necessary only when surface diffusion occurs. Moreover, the adoption of model (18) implies the need to determine inverse functions with respect to isotherms (5), i.e.,
This factor makes computations more difficult, especially when functions (5) are not given explicitly, e.g., when
qi values are derived from numerical calculations e.g., according to ideal adsorption solution theory (IAST), real adsorption solution theory (RAST) [
17] or are obtained on the basis of vacancy solution theory (VST) or from molecular simulations. For this reason, an alternative model, which does not contain a matrix of thermodynamic factors
Γ, is presented below, and then a method for the determination of steady states for the discussed complex multicomponent diffusion is proposed.
2.2. An Alternative Model of Surface Diffusion
It can be easily shown that for the surface diffusion of a single species (
K = 1), its molar flux density is:
After applying the Maxwell–Stefan formalism, one can derive the formula for surface diffusion for
K components:
In matrix notation we have:
After utilizing expression (9) in Equation (25) one gets:
where
Bs matrix determines surface diffusion:
The matrix
Δ present in Equations (25) and (26) has a diagonal form
The determination of the inverse of the matrix
Δ is therefore an elementary problem. In surface diffusion models, there are coefficients of two types. They are the diffusion coefficients of the given
i-th component
and the counter-sorption diffusivities
between component A
i and component A
j. In this paper, the effect of the fractional surface occupancy on the values of the surface diffusion coefficients is addressed based on the following relationship:
where
is the diffusion coefficient corresponding to zero coverage, i.e., θ
t = 0 and α is the degree of confinement. Furthermore, the following formula holds:
The Maxwell–Stefan counter-sorption diffusivities were calculated according to the recommendations available in the literature [
22]:
2.3. Method for Determining Steady States
Mathematical formalism implies that the determination of the steady state should consist in taking the partial pressure derivatives with respect to time in Equation (8) to be equal to zero and solving the resulting system of equations. Following in this manner, from Equations (1), (10) and (25) the following formula will be obtained:
It is a system of K second order ordinary differential equations due to K partial pressures pi (i = 1, 2, …, K). With given values of partial pressures for z = 0 and for z = L these equations can be solved obtaining K functions pi(z). Equation (32) was obtained assuming thermodynamic equilibrium between the gas and the adsorbed phase. Therefore knowing isotherms qi(p) one can calculate functions qi (z) for i = 1, 2, …, K.
Since the transport of mass in the pores and via surface diffusion are parallel mechanisms, both of these diffusion routes add up algebraically and give the total molar flux density per unit cross section. Due to the assumption of interphase thermodynamic equilibrium (5) and the resulting Equation (7), another approach for the determination of steady states can be proposed. To the best of the authors’ knowledge, such a method has not yet been published. The proposed method takes advantage of the transformed Equations (10) and (26). We illustrate it here using examples of diffusion through flat or cylindrical membranes of small thickness
L in comparison with the membrane radius. In these cases, the cylindrical membrane can be treated as flat and one can assume
m = 0. Then, the following system of 2 ×
K nonlinear equations is obtained that needs to be solved:
where ξ is a dimensionless spatial coordinate. The pore volume fraction has been explained previously (Equations (14)–(16)). The expression
Bs [
q(
p)] indicates that the values of
q in the matrix
Bs (Equation (27)) are determined using the partial pressures calculated from Equation (33). In Equation (34) the matrix of thermodynamic factors
Γ does not appear, and the only matrix to be inverted is the diagonal matrix
Δ. Information about the adsorption equilibrium is incorporated in the matrices
J and
Δ. Using the mathematical formalism, we state that a system of
K second order differential Equation (32) correspond to 2 ×
K first order differential Equations (33) and (34).
There are boundary conditions associated with the systems of differential Equations (33) and (34). Their form depends on the mode in which the process is carried out. Let us assume that the boundary conditions are of the first type, i.e.,
The contribution of the external resistance to mass transport, i.e., from the gas to the membrane surface, decreases with increasing membrane thickness L. For a sufficiently large thickness L almost all the resistance to mass transport is localized in the membrane.
The method for determining the steady state of a specific diffusion system can be summarized by the following algorithm:
- (a)
adopt tentative values of Np and Ns;
- (b)
assume the boundary conditions (35), (37) for ξ = 0;
- (c)
integrate systems of ordinary differential Equations (33) and (34) from ξ = 0 to ξ = 1;
- (d)
verify the fulfilment of the boundary conditions (36), (38) for ξ = 1;
- (e)
if the boundary conditions are not met, improve the values of Np and Ns e.g., with the help of Newton’s overriding algorithm and return to point (b);
- (f)
if the boundary conditions are met, it means that the profiles of all 2 × K variables, i.e., p(ξ) and q(ξ) = q [p(ξ)] and the molar flux densities corresponding to pore diffusion Np and surface diffusion Ns describe the determined steady state. Thence it is possible to calculate Ntot = Np + Ns.
From the above algorithm, it follows that the determination of the steady state of the diffusion process under consideration reduces to solving a system of 2 ×
K algebraic equations, which can be written as:
where
λ is a vector of model parameters.
The steady-state characteristics of considered diffusion processes include concentration profiles in the gas y(ξ) and in the solid θ(ξ), as well as molar flux densities Np, Ns, Ntot. The above algorithm illustrates how to determine all these quantities.
The systems of differential Equations (33) and (34) with boundary conditions (35)–(38) constitute a mathematical model of steady-state complex multicomponent diffusion. On the basis of this model, process analysis can be performed and the influence of the individual constituents, i.e., molecular and Knudsen diffusion, viscous flow and surface diffusion, can be assessed depending on the process conditions. Nevertheless, in order to assess the validity of the model (33)–(38) its reliability should be verified. The method for determining steady states was verified by comparing the results obtained from Equations (33) and (34) with the results obtained according to the dynamic model (8) for sufficiently long time t.
Let λ
k be the
kth element of the vector of parameters
λ. Applying the continuation of solutions of the system of Equations (39) with respect to this chosen parameter, one gets a branch of steady states. The use of a continuation algorithm, e.g., the local parameterization method [
41], gives the possibility to determine the parametric dependence of the diffusion model (33)–(38).