Competitive Adsorption and Diffusion of Gases in a Microporous Solid

Abstract

The experimental and theoretical study of the co-adsorption and co-diffusion of several gases through a microporous solid and the instantaneous (out of equilibrium) distribution of the adsorbed phases is particularly important in many fields, such as gas separation, heterogeneous catalysis, purification of confined atmospheres, reduction of exhaust emissions contributing to global warming, etc. The original NMR imaging technique used gives a signal characteristic of each adsorbed gas at each instant and at each level of the solid and therefore the distribution of several gases in competitive diffusion and adsorption. But it does not allow to determine separately the inter- and intra-crystallite quantities. A new fast and accurate analytical method for the calculation of the coefficients of co-diffusing gases in the intra- and inter-crystallite spaces of microporous solid (here ZSM 5 zeolite) is developed, using high-performance methods (iterative gradient methods of residual functional minimization and analytical methods of influence functions) and mathematical co-adsorption models, as well as the NMR spectra of each adsorbed gas in the bed. These diffusion coefficients and the gas concentrations in the inter- and intra-crystallite spaces are obtained for each position in the bed and for different adsorption times.

1. Introduction

Knowledge of the co-diffusion and co-adsorption coefficients of reactants and products is essential when a heterogeneous catalytic reaction is performed by flowing gases through a microporous catalyst bed. But generally the distribution of the various reactants adsorbed on the catalyst is very inhomogeneous and, moreover, very variable from one reactant to another. It is therefore necessary to determine at every moment the diffusion coefficient of each reactant in the presence of the others and its instantaneous distribution along the length of the catalyst bed.

Classical H-MRI should be a good technique for monitoring the co-diffusion and co-adsorption of several gases flowing through a microporous bed. However, since the signal obtained is not specific for each gas, this requires that each experiment be performed several times under identical conditions and each time with only one not deuterated gas. To remedy the drawbacks of classical imaging, we have used the NMR imaging technique, named slice selection procedure, to follow the diffusion and adsorption of a gas in a microporous bed [1, 2, 3]. The sample is displaced vertically, step-by-step, relative to a very thin coil detector during the adsorption of the gas. The bed is assumed to consist of N very thin layers of solid, and the region probed is limited to each layer, so that the variation of the concentration of gas absorbed at the level of each layer is obtained as a function of time. An interesting feature of this technique is its ability to visualize directly the co-diffusion of several gases. Indeed, the NMR signals are quantitatively characteristic of the adsorbed gases. They can therefore provide directly, at every moment and at every level of the bed, the distribution of several gases competing in diffusion and adsorption. We have presented in a previous paper the experimental results of the co-diffusion of benzene and hexane through a silicalite bed [4]. In [5, 6] we have developed a mathematical methodology for efficient linearization of similar models. Using Heaviside’s operational method and Laplace’s integral transformation method, we have built solutions allowing fast calculations for two-component co-adsorption in a heterogeneous zeolite bed and for the dehydration of natural gas [7]. In this chapter we have improved the methods previously used to compute the diffusion coefficients against time, increasing the accuracy and speed of calculations by significantly reducing the iteration number. This made it possible to use them for the co-adsorption of several gases diffusing along such a column.

2. Experimental

The NMR imaging technique, the sample-holder bulb containing the liquid phase in equilibrium with the gas phase, and the narrow zone monitored by the detector have been described in [1, 2, 3], respectively.

The upper face of the cylindrical bed of zeolite crystallites is exposed to a constant pressure of each gas (Figure 1). The diffusion of the two gases is axial in the macropores of the inter-crystallite space (z direction along the height, l, of the bed) and radial in the micropores of the zeolite. According to the experimental conditions, the zeolite bed consists of a large number, N, of very thin layers of solid, of thickness $\mathrm{\Delta }{l}_k=l_k-l_{k-1}$ , perpendicular to the propagation of the gas in the z direction. The corresponding coefficients of inter- and intra-crystallite space are D_inter,k and D_intra,k , respectively.

Figure 1.

Distribution of the layers (left) and corresponding parameters (right).

3. Experimental results: gaseous benzene and hexane adsorption curves

The experimental results have been summarized in [4]: the spectrum of each gas at every instant and every level of the solid and the benzene and hexane concentrations along the sample, for each diffusion time. Here we shall only use the evolution, as a function of time, of the benzene and hexane concentrations at different levels of the sample, on which the calculations of the diffusion coefficients and the instantaneous inter- and intra-crystallite concentrations are based [8]. Figure 2 shows clearly that, under the chosen experimental conditions, benzene hinders the diffusion of hexane, and this at every moment. Moreover, it can be noticed that, at equilibrium, the amount of benzene within the zeolite is twice that of hexane, indicating quantitatively the relative affinity to the two adsorbates.

Figure 2.

Evolution vs. time of the benzene and hexane concentrations (arbitrary units) at different levels of the sample (continuous, experimental curves; dotted, their approximations used for simulation) [from Ref. [8], reprinted with permission from ACS].

These curves display modulations as a function of time, which must be averaged for all subsequent mathematical representations. These modulations are weak at the lower layers of the tube and can be due to errors in the measurement of small amounts. Those closer to the arrival of the gas are greater and are similar for the two gases. We suggested that these fluctuations may be due to the fact that inter-crystallite adsorption at levels close to the gas phase is fast compared to the liquid-gas equilibrium, which is not as instantaneous for a mixture as for a single component [8]. Each slight decrease of the gas pressure could correspond a slight fast desorption.

4. A mathematical model of competitive co-adsorption and co-diffusion in microporous solids

4.1 Co-adsorption model in general formulation

The model presented is similar to the bipolar model [2, 3, 8, 9]. By developing the approach described by Ruthven and Kärger [10, 11] and Petryk et al. [5] concerning the elaboration of a complex process of co-adsorption and co-diffusion, it is necessary to specify the most important hypotheses limiting the process.

The general hypothesis adopted to develop the model presented in the most general formulation is that the interaction between the co-adsorbed molecules of several gases and the adsorption centers on the surface in the nanoporous crystallites is determined by the nonlinear competitive equilibrium function of the Langmuir type, taking into account physical assumptions [10]:

1. Co-adsorption is caused by the dispersion forces whose interaction is established by Lennard-Jones and the electrostatic forces of gravity and repulsion described by van der Waals [11].

2. The co-diffusion process involves two types of mass transfer: diffusion in the macropores (inter-crystallite space) and diffusion in the micropores of crystallites (intra-crystallite space).

3. During the evolution of the system toward equilibrium, there is a concentration gradient in the macropores and/or in the micropores.

4. Co-adsorption occurs on active centers distributed over the entire inner surface of the nanopores (intra-crystallite space) [10, 11]. All crystallites are spherical and have the same radius R; the crystallite bed is uniformly packed.

5. Active adsorption centers adsorb molecules of the ith adsorbate, forming molecular layers of adsorbate on their surfaces.

6. Adsorbed molecules are held by active centers for a certain time, depending on the temperature of the process.

Taking into account these hypotheses, we have developed a nonlinear co-adsorption model. The meaning of the symbols is given in the nomenclature:

$$\frac{{\mathit{дС}}_s\left(t,Z\right)}{\mathit{дt}}=\frac{D_{{\text{inter}}_s}}{l^2}\frac{{\partial }^2{С}_s}{\partial Z^2}-e_{\text{inter}}{\tilde{K}}_s\frac{D_{{\text{intra}}_s}}{R^2}{\left(\frac{\partial Q_s}{\partial X}\right)}_{X=1}$$ E1

$$-H\frac{\partial T\left(t,z\right)}{\partial t}-{\mathit{uh}}_g\frac{\partial T}{\partial z}-\sum _{s=1}^m\mathrm{\Delta }{\overline{H}}_s\frac{\partial {\overline{Q}}_s}{\partial t}-2\frac{{\alpha }_h}{R_{\mathit{\text{column}}}}T+\mathrm{\Lambda }\frac{{\partial }^{2}T}{\partial z^2}=0$$ E2

$$\frac{{\mathit{дQ}}_s\left(t,X,Z\right)}{\mathit{дt}}=\frac{D_{{\text{intra}}_s}}{R^2}\left(\frac{{\partial }^2{Q}_s}{\partial X^2}+\frac{2}{X}\frac{\partial Q_s}{\partial X}\right)$$ E3

with initial conditions

$$C_s\left(t=0,Z\right)=0;Q_s\left(t=0,X,Z\right)=0;Z\in \left(0,1\right),X\in \left(0,1\right),s=\overline{1,\mathrm{m}}$$ E4

boundary conditions for coordinate X of the crystallite

$$\frac{\partial }{\partial X}{Q}_s\left(t,X=0,Z\right)=0\left(\text{symmetry conditions}\right),$$ E5

$$Q_s\left(t,X=1,Z\right)=\frac{K_s\left(T\right)C_s\left(t,Z\right)}{1+\sum _{s_1=1}^m{K}_{s_1}\left(T\right)C_{s_1}\left(t,Z\right)},s=\overline{1,\mathrm{m}}\left(\text{Langmuir equilibrium}\right),$$ E6

boundary and interface conditions for coordinate Z

$$C_s\left(t,1\right)=1,\frac{{\mathit{дС}}_s}{\mathit{дZ}}\left(t,Z=0\right)=0,t>0$$ E7

$$T\left(t,\mathrm{Z}\right){\left|{}_{Z=1}=T_{\mathit{\text{initial}}},\frac{\partial }{\partial z}T(t,\mathrm{Z})\right|}_{Z=0}=0.$$ E8

with $K_s\left(T\right)=k_{0s}exp\left(-\frac{\mathrm{\Delta }{H}_s}{R_{g}T}\right)$ .

Here the activation energy is the heat of adsorption defined as $\mathrm{\Delta }{H}_s=\overline{\varphi }-\left(U_{g_s}-U_{{\mathit{ads}}_s}\right)-R_{g}T$ , where $U_{g_s}-U_{{\mathit{ads}}_s}$ —the difference between the kinetic energies of the molecule of the ith component of the adsorbate in the gaseous and adsorbed states is the magnitude of the Lennard-Jones potential, averaged over the pore volume of the adsorbent [11].

The non-isothermal model (1)–(8) can easily be transformed into isothermal model, removing the temperature Eq. (2) and condition (8) and replacing the functions $K_s\left(T\right)$ with the corresponding equilibrium constants $K_s$ . The competitive diffusion coefficients $D_{{\text{intra}}_s}$ and $D_{{\text{inter}}_s}$ can be considered as functions of the time and the position of the particle in the zeolite bed.

4.2 The inverse model of co-diffusion coefficient identification: application to the benzene-hexane mixture

On the basis of a developed nonlinear co-adsorption model (1)–(8), we construct an inverse model for the identification of the competitive diffusion coefficients $D_{{\text{intra}}_s}$ and $D_{{\text{inter}}_s}$ as a function of time and coordinate in the zeolite bed.

The mathematical model of gas diffusion kinetics in the zeolite bed is defined in domains: ${\mathrm{\Omega }}_{k_t}=\left(0,t^{\mathit{\text{total}}}\right)\times {\mathrm{\Omega }}_k,\left({\mathrm{\Omega }}_k=\left(L_{k-1},L_k\right),k=\overline{1,N+1},L_0=0<L_1<\dots <L_{N+1}=1\right)$ by the solutions of the system of differential equations

$$\frac{{\mathit{дС}}_{s_k}\left(t,Z\right)}{\mathit{дt}}=\frac{D_{{\text{inter}}_{s_k}}}{l^2}\frac{{\partial }^2{С}_{s_k}}{\partial Z^2}-e_{{\text{inter}}_k}{\tilde{K}}_{s_k}\frac{D_{{\text{intra}}_{s_k}}}{R^2}{\left(\frac{\partial Q_{s_k}}{\partial X}\right)}_{X=1}$$ E9

$$\frac{{\mathit{дQ}}_{s_k}\left(t,X,Z\right)}{\mathit{дt}}=\frac{D_{{\text{intra}}_{s_k}}}{R^2}\left(\frac{{\partial }^2{Q}_{s_k}}{\partial X^2}+\frac{2}{X}\frac{\partial Q_{s_k}}{\partial X}\right)$$ E10

with initial conditions

$$C_{s_k}\left(t=0,Z\right)=0;Q_{s_k}\left(t=0,X,Z\right)=0;X\in \left(0,1\right),Z\in {\mathrm{\Omega }}_k,k=\overline{1,N+1},$$ E11

boundary and interface conditions for coordinate Z

$$C_{s_1}\left(t,L_1\right)=1,\frac{{\mathit{дС}}_{s_1}}{\mathit{дZ}}\left(t,Z=0\right)=0,t\in \left(0,t^{\mathit{\text{total}}}\right);$$ E12

$$\begin{array}{l}{\left[C_{s_k}\left(t,Z\right)-C_{s_k}\left(t,Z\right)\right]}_{Z=L_k}=0,\frac{\partial }{\partial Z}{\left[D_{{\text{inter}}_{s_{k-1}}}{C}_{s_{k-1}}\left(t,Z\right)-D_{{\text{inter}}_{s_k}}{C}_{s_k}\left(t,Z\right)\right]}_{Z=L_k}=0,k\\ =\overline{1,N},t\in \left(0,t^{\mathit{\text{total}}}\right);\end{array}$$ E13

boundary conditions for coordinate X in the particle

$$\frac{\partial }{\partial X}{Q}_{s_k}\left(t,X=0,Z\right)=0,Q_{s_k}\left(t,X=1,Z\right)=K_s{C}_{s_k}\left(t,Z\right)\left(\text{equilibrium conditions}\right),Z\in {\mathrm{\Omega }}_k,k=\overline{1,N+1}.$$ E14

Additional condition (NMR-experimental data)

$${\left.\left[C_{s_k}\left(t,Z\right)+{\overline{Q}}_{s_k}\left(t,Z\right)\right]\right|}_{h_k}={\left.M_{s_k}\left(t,Z\right)\right|}_{h_k},s=\overline{1,2};h_k\in {\mathrm{\Omega }}_k.$$ E15

The problem of the calculation (9)–(15) is to find unknown functions $D_{{\text{intra}}_{\mathrm{s}}}\in {\mathrm{\Omega }}_t,D_{{\text{inter}}_{\mathrm{s}}}\in {\mathrm{\Omega }}_t$ ( $D_{{\text{intra}}_{\mathrm{s}}}>0,D_{{\text{inter}}_{\mathrm{s}}}>0$ , $s=\overline{1,2}$ ), when absorbed masses $C_{s_k}\left(t,Z\right)+{\overline{Q}}_{s_k}\left(t,Z\right)$ satisfy the condition (15) for every point $h_k\subset {\mathrm{\Omega }}_k$ of the kth layer [8, 12].

Here

$${\mathrm{e}}_{{\text{inter}}_k}=\frac{{\mathrm{\epsilon }}_{{\text{inter}}_k}{c}_{\mathit{sk}}}{{\mathrm{\epsilon }}_{{\text{inter}}_k}{c}_{\mathit{sk}}+\left(1-{\mathrm{\epsilon }}_{{\text{inter}}_k}\right)q_{\mathit{sk}}}\approx \frac{{\mathrm{\epsilon }}_{{\text{inter}}_k}}{\left(1-{\mathrm{\epsilon }}_{{\text{inter}}_k}\right){\tilde{\mathrm{K}}}_{\mathit{sk}}},{\tilde{\mathrm{K}}}_{s_k}=\frac{{\mathrm{q}}_{s_k\infty }}{{\mathrm{c}}_{s_k\infty }},$$

where ${\overline{Q}}_s\left(t,Z\right)=\underset{0}{\overset{1}{\int }}{Q}_s\left(t,X,Z\right)\mathit{XdX}$ is the average concentration of adsorbed component s in micropores and ${\left.M_s\left(t,Z\right)\right|}_{h_k}$ is the experimental distribution of the mass of the sth component absorbed in macro- and micropores at $h_k\subset {\mathrm{\Omega }}_k$ (results of NMR data, Figure 2).

4.3 Iterative gradient method of co-diffusion coefficient identification

The calculation of $D_{{\text{intra}}_{s_k}}$ and $D_{{\text{inter}}_{s_k}}$ is a complex mathematical problem. In general, it is not possible to obtain a correct formulation of the problem (9)–(15) and to construct a unique analytical solution, because of the complexity of taking into account all the physical parameters (variation of temperature and pressure, crystallite structures, nonlinearity of Langmuir isotherms, etc.), as well as the insufficient number of reliable experimental data, measurement errors, and other factors.

Therefore, according to the principle of Tikhonov and Arsenin [13], later developed by Lions [14] and Sergienko et al. [15], the calculation of diffusion coefficients requires the use of the model for each iteration, by minimizing the difference between the calculated values and the experimental data.

The calculation of the diffusion coefficients (9)–(15) is reduced to the problem of minimizing the functional of error (16) between the model solution and the experimental data, the solution being refined incrementally by means of a special calculation procedure which uses fast high-performance gradient methods [6, 8, 12, 15].

According to [12, 15], and using the error minimization gradient method for the calculation of $D_{{\text{intra}}_{s_k}}$ and $D_{{\text{inter}}_{s_k}}$ of the sth diffusing component, we obtain the iteration expression for the n + 1th calculation step:

$$\begin{array}{c}{D}_{{\text{intra}}_{s_k}}^{n+1}\\ \left(t\right)=D_{int{\mathit{ra}}_{s_k}}^n\left(t\right)-\nabla J_{D_{{\text{intra}}_{s_k}}}^n\left(t\right)\frac{{\left[C_{s_k}\left(D_{{\text{inter}}_{s_k}}^n,D_{{\text{intra}}_{s_k}}^n;t,h_k\right)+{\overline{Q}}_{s_k}\left(D_{{\text{inter}}_{s_k}}^n,D_{{\text{intra}}_{s_k}}^n;t,h_k\right)-M_{s_k}\left(t\right)\right]}^2}{{‖\nabla J_{D_{{\text{intra}}_{s_k}}}^n\left(t\right)‖}^2+{‖\nabla J_{D_{{\text{inter}}_{s_k}}}^n\left(t\right)‖}^2},\\ \begin{array}{l}{D}_{{\text{inter}}_{s_k}}^{n+1}\left(t\right)=D_{{\text{inter}}_{s_m}}^n\left(t\right)-\nabla J_{D_{{\text{inter}}_{s_k}}}^n\left(t\right)\frac{{\left[C_{s_k}\left(D_{{\text{inter}}_{s_k}}^n,D_{{\text{intra}}_{s_k}}^n;t,h_k\right)+{\overline{Q}}_{s_k}\left(D_{{\text{inter}}_{s_k}}^n,D_{{\text{intra}}_{s_k}}^n;t,h_k\right)-M_{s_k}\left(t\right)\right]}^2}{{‖\nabla J_{D_{{\text{intra}}_{s_k}}}^n\left(t\right)‖}^2+{‖\nabla J_{D_{{\text{inter}}_{\mathit{sk}}}}^n\left(t\right)‖}^2},\\ t\in \left(0,t^{\mathit{\text{total}}}\right)\end{array}\end{array}$$ E16

where $J\left(D_{{\text{inter}}_{{\mathrm{s}}_k}},D_{{\text{intra}}_{{\mathrm{s}}_k}}\right)$ is the error functional, which describes the deviation of the model solution from the experimental data on $h_k\in {\mathrm{\Omega }}_k$ , which is written as

$$\begin{array}{l}J\left(D_{{\text{inter}}_{{\mathrm{s}}_k}},D_{{\text{intra}}_{{\mathrm{s}}_k}}\right)\\ =\frac{1}{2}\underset{0}{\overset{T}{\int }}{\left[C_s\left(\tau ,Z,D_{{\text{inter}}_{{\mathrm{s}}_k}},D_{{\text{intra}}_{{\mathrm{s}}_k}}\right)+{\overline{Q}}_s\left(t,Z,D_{{\text{inter}}_{{\mathrm{s}}_k}},D_{{\text{intra}}_{{\mathrm{s}}_k}}\right)-M_{s_k}\left(t\right)\right]}_{h_k}^2\mathit{d\tau },h_k\in {\mathrm{\Omega }}_k,k=\overline{1,N+1},\end{array}$$ E17

$\nabla J_{D_{{\text{intter}}_{s_k}}}^n\left(t\right)$ , $\nabla J_{D_{{\text{intra}}_{s_k}}}^n\left(t\right)$ (the gradients of the error functional), $J\left(D_{{\text{inter}}_{{\mathrm{s}}_k}},D_{{\text{intra}}_{{\mathrm{s}}_k}}\right)$ .

${‖\nabla J_{D_{{\text{inter}}_{\mathit{sk}}}}^n\left(t\right)‖}^2=\underset{0}{\overset{T}{\int }}{\left[\nabla J_{D_{{\text{inter}}_{\mathit{sk}}}}^n\left(t\right)\right]}^2\mathit{dt}$ , ${‖\nabla J_{D_{{\text{intra}}_{\mathit{sk}}}}^n\left(t\right)‖}^2=\underset{0}{\overset{T}{\int }}{\left[\nabla J_{D_{{\text{intra}}_{\mathit{sk}}}}^n\left(t\right)\right]}^2\mathit{dt}$ .

4.4 Analytical method of co-diffusion coefficient identification

With the help of iterative gradient methods on the basis of the minimization of the residual functional, very precise and fast analytical methods have been developed making it possible to express the diffusion coefficients in the form of time-dependent analytic functions (16). For their efficient use, it is necessary to have an extensive experimental database, with at least two experimental observation conditions for the simultaneous calculation of $D_{{\text{intra}}_{s_k}}$ and $D_{{\text{inter}}_{s_k}}$ coefficients. Our experimental studies were carried out for five Z positions of the swept zeolite layer for each of the adsorbed components. The data were not sufficient to fully apply this simultaneous identification method to these five sections. We therefore used a combination of the analytical method and the iterative gradient method for determining the co-diffusion coefficients.

Using Eqs. (9)–(15), it is possible to calculate $D_{{\text{intra}}_{s_k}},D_{{\text{inter}}_{\mathit{sk}}}$ as a function of time using the experimental data obtained by NMR scanning. In particular, in Eqs. (9) and (10), the co-diffusion coefficients can be set directly as functions of the time t: $D_{{\text{intra}}_{s_k}}\left(t\right),D_{{\text{inter}}_{\mathit{sk}}}\left(t\right)$ . In this case, the boundary condition (11) can be given in a more general form—also as a function of time:

$$C_{s_1}\left(t,1\right)=C_s^{\mathit{\text{initial}}}\left(t\right).$$ E18

Experimental NMR scanning conditions are defined simultaneously for all P observation surfaces:

$${\left[C_{\mathit{sk}}\left(t,Z\right)+\underset{0}{\overset{1}{\int }}{Q}_{\mathit{sk}}\left(t,X,Z\right)\mathit{dX}\right]}_{Z=h_i}={\left.M_{{\mathit{sk}}_i}\left(t,Z\right)\right|}_{h_i},i=\overline{1,P},s=\overline{1,2};h_i\in \underset{k=1}{\overset{N+1}{\cup }}{\mathrm{\Omega }}_k$$ E19

For simplicity we design $u\left(t,Z\right)=C_{\mathit{sk}}\left(t,Z\right),v\left(t,\mathrm{X},Z\right)=Q_{\mathit{sk}}\left(t,\mathrm{X},Z\right)$ ,

$$b\left(t\right)=D_{{\text{intra}}_{\mathit{sk}}}\left(t\right)/{\mathrm{R}}^2,{\chi }_i\left(t\right)=M_{{\mathit{sk}}_i}\left(t\right),i=\overline{1,\mathrm{P}},$$

and considering Eq. (10) in flat form, its solution can be written as [16]:

$$v\left(t,\mathrm{X},\mathrm{Z}\right)=-\underset{0}{\overset{t}{\int }}{\mathcal{H}}_{4\xi }^{\left(2\right)}\left(t,\tau ,\mathrm{X},1)\right)b\left(\tau \right)u\left(\tau ,\mathrm{Z}\right)\mathit{d\tau }$$ E20

where ${\mathcal{H}}_{4\xi }^{\left(2\right)}\left(t,\tau ,\mathrm{X},\xi \right)=-2\sum _{m=0}^{\infty }{e}^{-{\eta }_m^2\left({\theta }_2\left(\mathrm{t}\right)-{\theta }_2\left(\tau \right)\right)}{\eta }_{m}cos{\eta }_{m}X\cdot {\left(-1\right)}^m$ .

Here the Green influence function of the particle ${\mathcal{H}}_k^{\left(2\right)},\overline{k=1,4}$ is used; it has the form [17].

${\mathcal{H}}_4^{\left(2\right)}\left(t,\tau ,\mathrm{X},\xi \right)=2\sum _{m=0}^{\infty }{e}^{-{\eta }_m^2\left({\theta }_2\left(\mathrm{t}\right)-{\theta }_2\left(\tau \right)\right)}cos{\eta }_{m}Xcos{\eta }_m\xi ,{\eta }_m=\frac{2m+1}{2}\pi ,$

${\mathcal{H}}_3^{\left(2\right)}\left(t,\tau ,\mathrm{X},\xi \right)=2\sum _{m=0}^{\infty }{e}^{-{\eta }_m^2\left({\theta }_2\left(\mathrm{t}\right)-{\theta }_2\left(\tau \right)\right)}sin{\eta }_{m}Xsin{\eta }_m\xi ,\begin{array}{cc}\begin{array}{cc}& \end{array}{\eta }_m=\frac{2m+1}{2}\pi ,& \end{array}$

${\mathcal{H}}_2^{\left(2\right)}\left(t,\tau ,\mathrm{X},\xi \right)=1+2\sum _{m=1}^{\infty }{e}^{-{\eta }_m^2\left({\theta }_2\left(\mathrm{t}\right)-{\theta }_2\left(\tau \right)\right)}cos{\eta }_{m}X\cdot cos{\eta }_m\xi ,\begin{array}{cc}\begin{array}{cc}& {\eta }_m=\mathit{m\pi },\end{array}& \end{array}$ where ${\theta }_2\left(t\right)=\underset{0}{\overset{t}{\int }}b\left(s\right)\mathit{ds}$ .

The notation ${\mathcal{H}}_{4\mathit{\tau \xi }}^{\left(2\right)},{\mathcal{H}}_{4\mathit{\xi \xi }}^{\left(2\right)}$ means partial derivatives of the influence function ${\mathcal{H}}_4^{\left(2\right)}$ relative to the definite variables $\tau$ and $\xi$ , respectively.

Based on formula (20), we calculate

$$v_X\left(t,\mathrm{X},\mathrm{Z}\right)=-\underset{0}{\overset{t}{\int }}{\mathcal{H}}_{4\mathit{\xi X}}^{\left(2\right)}\left(t,\tau ,\mathrm{X},1\right)b\left(\tau \right)u\left(\tau ,\mathrm{Z}\right)\mathit{d\tau }$$ E21

Integrating parts (21), taking into account the relations.

$${\mathcal{H}}_{4X}^{\left(2\right)}\left(t,\tau ,\mathrm{X},\xi \right)=-{\mathcal{H}}_{3\xi }^{\left(2\right)}\left(t,\tau ,X,\xi \right),{\mathcal{H}}_{3\tau }^{\left(2\right)}\left(t,\tau ,X,\xi \right)=-b\left(\tau \right){\mathcal{H}}_{3\mathit{\xi \xi }}^{\left(2\right)}\left(t,\tau ,X,\xi \right),$$

and the initial condition $u_{\left|t=0\right.}=0$ , we find

$$v_X\left(t,\mathrm{X},\mathrm{Z}\right)=\underset{0}{\overset{t}{\int }}{\mathcal{H}}_3^{\left(2\right)}\left(t,\tau ,\mathrm{X},1)\right)u_{\tau }\left(\tau ,\mathrm{Z}\right)\mathit{d\tau }$$ E22

We substitute the expression $v\left(t,\mathrm{X},\mathrm{Z}\right)$ (20) in the observation conditions (19):

$$u\left(t,h_i\right)-\underset{0}{\overset{1}{\int }}\mathit{XdX}\underset{0}{\overset{t}{\int }}{\mathcal{H}}_{4\xi }^{\left(2\right)}\left(t,\tau ,\mathrm{X},1\right)b\left(\tau \right)u\left(\tau ,h_i\right)\mathit{d\tau }={\chi }_i\left(\mathrm{t}\right),i=\overline{1,\mathrm{P}}$$ E23

Integrating parts (23) and taking into account equality

$${\mathcal{H}}_{4\xi }^{\left(2\right)}\left(t,\tau ,X,1\right)=-{\mathcal{H}}_{3X}^{\left(2\right)}\left(t,\tau ,X,1\right),$$

we obtain [16]

$$u\left(t,h_i\right)={\chi }_i\left(\mathrm{t}\right)-\underset{0}{\overset{t}{\int }}{\mathcal{H}}_3^{\left(2\right)}\left(t,\tau ,1,1\right)b\left(\tau \right)u\left(\tau ,h_i\right)\mathit{d\tau }+\underset{0}{\overset{t}{\int }}\underset{0}{\overset{1}{\int }}{\mathcal{H}}_3^{\left(2\right)}\left(t,\tau ,\mathrm{X},1\right)b\left(\tau \right)u\left(\tau ,h_i\right)\mathit{\text{dXdτ}},i=\overline{1,\mathrm{P}}$$ E24

Let us first put $u\left(t,h_P\right)={\mu }_{\mathit{sP}}\left(t\right)=C_s^{\mathit{\text{initial}}}\left(t\right)$ , where $Z=h_P$ is the observation surface, approaching the point of entry into the work area Z = 1.

Then Eq. (24) for i = P will be

$$\underset{0}{\overset{t}{\int }}{\mathcal{H}}_3^{\left(2\right)}\left(t,\tau ,1,1\right)b\left(\tau \right){\mu }_{\mathit{sP}}\left(t\right)\mathit{d\tau }={\chi }_{\mathit{sP}}\left(t\right)-{\mu }_{\mathit{sP}}\left(t\right)+\underset{0}{\overset{t}{\int }}\underset{0}{\overset{1}{\int }}{\mathcal{H}}_3^{\left(2\right)}\left(t,\tau ,\mathrm{x},1\right)b\left(\tau \right){\mu }_{\mathit{sP}}\left(t\right)\mathit{\text{dxdτ}}$$ E25

Applying to Eq. (25) the formula $\underset{\tau }{\overset{t}{\int }}b\left(\sigma \right){\mathcal{H}}_2^{\left(2\right)}\left(t,\sigma ,0,0)\right){\mathcal{H}}_4^{\left(2\right)}\left(t,\sigma ,0,0)\right)\mathit{d\sigma }=1$ ,

obtained by Ivanchov [16], and taking into account.

${\mathcal{H}}_3^{\left(2\right)}\left(t,\sigma ;1,1\right)={\mathcal{H}}_4^{\left(2\right)}\left(t,\sigma ;0,0\right)$ , we obtain

$$\begin{array}{l}\underset{0}{\overset{t}{\int }}b\left(\tau \right){\mu }_{\mathit{sP}}\left(\tau \right)\mathit{d\tau }=\underset{0}{\overset{t}{\int }}{\mathcal{H}}_2^{\left(2\right)}\left(t,\sigma ,0,0\right)b\left(\sigma \right)\left({\chi }_{\mathit{sP}}\left(\sigma \right)-{\mu }_{\mathit{sP}}\left(\sigma \right)\right)\mathit{d\sigma }\\ +\underset{0}{\overset{t}{\int }}{\mathcal{H}}_2^{\left(2\right)}\left(t,\sigma ,0,0\right)b\left(\sigma \right)\mathit{d\sigma }\underset{0}{\overset{\sigma }{\int }}\underset{0}{\overset{1}{\int }}{\mathcal{H}}_3^{\left(2\right)}\left(\sigma ,\tau ,\mathrm{x},1\right)b\left(\tau \right){\mu }_{\mathit{sP}}\left(t\right)\mathit{\text{dXdτ}},t\in \left[0,t^{\mathit{\text{total}}}\right]\end{array}$$ E26

Differentiating Eq. (26) by t, after the transformations series, we obtain

$${\mu }_{\mathit{sP}}\left(t\right)=\underset{0}{\overset{t}{\int }}{\mathcal{H}}_2^{\left(2\right)}\left(t,\sigma ,0,0\right)\left(b\left(\sigma \right){\mu }_{\mathit{sP}}\left(\sigma \right)+{\chi }_{\mathit{sP}}\left(\sigma \right)-{\mu }_{\mathit{sP}}\left(\sigma \right)\right)\mathit{d\sigma }$$ E27

After multiplying Eq. (27) on the expression ${\mathcal{H}}_4^{\left(2\right)}\left(t,\sigma ,0,0)\right)b\left(\sigma \right)$ , the integration by τ and the differentiation by t

$$b\left(t\right){\mu }_{\mathit{sP}}\left(t\right)+{\chi }_{\mathit{sP}}\left(\sigma \right)-{\mu }_{\mathit{sP}}\left(\sigma \right)=b\left(t\right)\underset{0}{\overset{t}{\int }}{\mathcal{H}}_2^{\left(2\right)}\left(t,\tau ,0,0)\right){\mu }_{\mathit{sP}}^{\prime }\left(\tau \right))\mathit{d\tau }.$$

So we obtain the expression for calculating the co-diffusion coefficient in the intra-crystallite space:

$$D_{int{\mathit{ra}}_{\mathit{sP}}}\left(t\right)\equiv R^{2}b\left(t\right)=R^2\frac{{\chi }_{\mathit{sP}}^{\prime }\left(t\right)-{\mu }_{\mathit{sP}}^{\prime }\left(t\right)}{\underset{0}{\overset{t}{\int }}{\mathcal{H}}_2^{\left(2\right)}\left(t,\tau ,0,0)\right){\mu }_{\mathit{sP}}^{\prime }\left(\tau \right)\mathit{d\tau }-{\mu }_{\mathit{sP}}\left(t\right)},t\in \left[0,t^{\mathit{\text{total}}}\right]$$ E28

Using calculated $D_{{\text{intra}}_{\mathit{sP}}}\left(t\right)$ with the formula (28) on the observation limit $h_P$ , we define the gradient method $D_{{\text{inter}}_{\mathit{sP}}}\left(t\right)$ in the same way. With $D_{{\text{intra}}_{\mathit{sP}}}\left(t\right)$ and $D_{{\text{inter}}_{\mathit{sP}}}\left(t\right)$ in $h_P$ , we calculate ${С}_{\mathit{sk}}\left(t,h_P\right)$ , substituting it in ${\mu }_{\mathit{sP}-1}\left(t\right)={С}_{\mathit{sk}}\left(t,h_P\right)$ for the next coefficient $D_{{\text{inter}}_{\mathit{si}}}\left(t\right),i=\overline{P-1,1}$ calculations. All subsequent coefficients $D_{int{\mathit{ra}}_{\mathit{si}}}\left(t\right)$ will be calculated by the formula

$$D_{int{\mathit{ra}}_{\mathit{si}}}\left(t\right)\equiv R^2{b}_{\mathit{si}}\left(t\right)=R^2\frac{{\chi }_{\mathit{si}}^{\prime }\left(t\right)-{\mu }_{\mathit{si}}^{\prime }\left(t\right)}{\underset{0}{\overset{t}{\int }}{\mathcal{H}}_2^{\left(2\right)}\left(t,\tau ,0,0)\right){\mu }_{\mathit{si}}^{\prime }\left(\tau \right)\mathit{d\tau }-{\mu }_{\mathit{si}}\left(t\right)},i=\overline{P-1,1}$$ E29

with parallel computing $D_{{\text{inter}}_{\mathit{si}}}\left(t\right),i=\overline{P-1,1}$ .

5. Numerical simulation and analysis: co-diffusion coefficients: concentration profiles in inter- and intra-crystallite spaces

The variation against time of the benzene and hexane intra-crystallite diffusion coefficients, $D_{{\text{intra}}_{1,\mathrm{k}}}$ аnd $D_{{\text{intra}}_{2,\mathrm{k}}}$ , respectively, is presented in Figure 3 for the five coordinate positions: 6, 8, 10, 12, and 14 mm, defined now from the top of the bed. The curves for benzene $D_{{\text{intra}}_{1,\mathrm{k}}}$ are pseudo exponentials. $D_{{\text{intra}}_{1,\mathrm{k}}}$ decreases from 9.0 E−13 to about 1.0 E−14 a.u. (equilibrium) depending on the position of the crystallite and the time, as well as on the amount of adsorbed gas. The shapes of the variations of $D_{{\text{intra}}_{2,\mathrm{k}}}$ for hexane are roughly the same, but the diffusion coefficients are higher, from about 9.0 E−12 to 3.0 E−13 a.u.

Figure 3.

Variation of intra-crystallite diffusion coefficients (arbitrary units) for benzene ${\mathit{D}}_{{\mathit{\text{intra}}}_{\mathit{1},\mathit{k}}}$ (left) and hexane ${\mathit{D}}_{{\mathit{\text{intra}}}_{\mathit{2},\mathit{k}}}$ (right) against time, at different positions in the bed. (Top) time range 6–240 mn, (bottom) time range 100–240 mn.

Figure 4 presents the variation against time of the benzene and hexane diffusion coefficients in inter-crystallite space, $D_{{\text{inter}}_{1,\mathrm{k}}}$ аnd $D_{{\text{inter}}_{2,\mathrm{k}}}$ , for the same positions. These coefficients decrease with time from 2.0 E−6 to 1.0 E−7 a.u. (equilibrium) for benzene and from 3 E−5 to 1.0 E−6 a.u. for hexane, depending on the bed position, and increase adsorbed concentrations.

Figure 4.

Variation of inter-crystallite diffusion coefficients (a.u.) for benzene (left) and hexane (right) against time at different positions in the bed.

Figure 5 shows the variation against time of the calculated concentrations C for benzene and hexane in the inter-crystallite space. As can be seen, these concentrations approach the equilibrium values for a diffusion time around 250 min. But the variations of the concentrations with time are rather different for the two gases.

Figure 5.

Variation of the inter-crystallite concentration (a.u.) calculated for benzene (left) and hexane (right) against time and at different positions in the bed.

Figure 6 shows the variation of the concentrations Q(t,X,z) of adsorbed benzene (left) and hexane (right) in the micropores of the intra-crystallite space from the surface (abscissa-1) to the center (abscissa-0) of the crystallites located between 6 and 14 mm from the top of the bed and after 25–200 min. of diffusion (a, b, c, and d, respectively). The gradients increase, and the mean concentrations decrease with the increasing distance of the particles from the arrival of the gases. The particles at 6 and 8 mm are saturated with benzene after 100 min, but not yet with hexane.

Figure 6.

Distribution of the benzene (left) and hexane (right) concentrations in the intra-crystallite space from the surface (abscissa 1) to the center (abscissa 0) of the crystallites, at different times: (1) dark blue, t = 25 min; (2) green, t = 50 min; (3) brown, t = 100 min; and (4) red, t = 200 min.

6. Conclusion

The main result of this work is the possibility, from a single experiment, of simultaneously distributing several co-diffusing gases in a porous solid and of using the methods of mathematical modeling to analyze for each of them the distribution of their concentrations in the intra- and inter-crystallite spaces.

Using the experimental NMR data and proposed co-adsorption models, the identification procedures for calculating the co-diffusion coefficients for two or more components in intra- and inter-crystallite spaces are developed. These procedures use the iterative gradual identification methods on minimizing of the Lagrange error function and rapid analytic methods based on the influence function. The co-diffusion coefficients were obtained as a function of time for different positions along the catalyst bed. In particular, those in the intra-crystallite space were computed by the analytical method which allowed a calculation with a relatively high degree of discretization over time and to reduce practically twice the volume of iterative calculations. Using these results, the concentrations of co-diffusing benzene and hexane in the inter- and intra-crystallite spaces were calculated for each time and each position in the bed.

Nomenclature

$k=\overline{1,N+1}$

layer number, subscript k will be added to all the following symbols to specify that they are characteristic of the kth layer

c

adsorbate concentration in macropores

c_∞

equilibrium adsorbate concentration in macropores

C = c/c_∞

dimensionless adsorbate concentration in macropores

D_inter

diffusion coefficient in macropores, m²/s

D_intra

diffusion coefficient in micropores, m²/s

$\tilde{K}$

adsorption equilibrium constant

l

bed length, m

$\mathrm{\Delta }l=l_k-l_{k-1};k=\overline{1,N+1}$

layer thickness (all layers have the same thickness)

L

dimensionless bed length (L = 1)

q

adsorbate concentration in micropores

q_∞

equilibrium adsorbate concentration in micropores

Q = q/q_∞

dimensionless concentration of adsorbate in micropores

$T$

temperature of gas phase flow, K, and time total, s

M

mass total

$u$

velocity of gas phase flow, m/s²

$\mathrm{\Lambda }$

coefficient of thermal diffusion along the columns

$h_g$

gas heat capacity, kJ/(kg K)

μ

molecular mass of adsorbate, kg/mol

$H$

total heat capacity of the adsorbent and gas, kJ/(kg K)

${\alpha }_h$

heat transfer coefficient

R_column

column radius, m

R_g

gas constant, kJ mol/(m³ K)

$\mathrm{\Delta }{H}_i$

activation energy ( $\mathrm{\Delta }{H}_i=\mathrm{\Delta }{\overline{H}}_i/\mu$ ), kJ/mol

$\mathrm{\Delta }{\overline{H}}_i$

adsorption heat, kJ/kg

$k_{0i}$

empirical equilibrium coefficient for the i adsorbate, depending on the adsorbent properties and the diffusing adsorbate component ( $k_{0i}$ equal to the ratio of the desorption and adsorption rate constants)

x

distance from crystallite center, m

R

mean crystallite radius, m (we assume that the crystallites are spherical)

X = x/R

dimensionless distance from crystallite center

z

distance from the bottom of the bed for mathematical simulation, m

Z = z/l

dimensionless distance from the bottom of the bed

t

time

$\tau ,\xi$

variables of integration

$t^{\mathit{\text{total}}}$

total duration of co-adsorption, mn

L_k

dimensionless position of the kth layer

h_k

(L_k –L _k-1)/2

ε_inter

inter-crystallite bed porosity

e_inter

value utilized in Eq. (9)

n

iteration number of identification

m

number of adsorbed components

P

number of NMR observation surfaces

s

index of adsorbate component

i

index of NMR observation surface

initial

index of initial value (concentrations, temperature)

macro

index of extended Lagrange functional component for inter-crystallite space

micro

index of extended Lagrange functional component for intra-crystallite space

Appendix

A. Iterative gradient method of the identification of co-diffusion coefficients

The methodology for solving the direct boundary problem (9)–(15), which describes the diffusion process in a heterogeneous nanoporous bed, is developed in [9, 12, 15]. According to [12] the procedure for determining the diffusion coefficients (16) requires a special technique for calculating the gradients $\nabla J_{D_{{\text{intra}}_{s_k}}}^n\left(t\right)$ and $\nabla J_{D_{{\text{inter}}_{s_k}}}^n\left(t\right)$ of the residual functional (17). This leads to the problem of optimizing the extended Lagrange function [12, 15]:

$$\mathrm{\Phi }\left(D_{{\text{inter}}_{\mathrm{sk}}},D_{{\text{intra}}_{\mathrm{sk}}}\right)=J_s+I_{s_{\mathit{\text{macro}}}}+I_{s_{\mathit{\text{micro}}}},$$ EA.1

where $I_{s_{\mathit{\text{macro}}}},I_{s_{\mathit{\text{micro}}}}$ are the components given by Eqs. (A.2) and (A.3), corresponding to the macro- and microporosity, respectively

$$I_{s_{\mathit{\text{macro}}}}=\underset{0}{\overset{T}{\int }}\underset{L_{k-1}}{\overset{L_k}{\int }}{\varphi }_{s_k}\left(t,Z\right)\left(\frac{\partial C_{s_k}}{\partial t}-\frac{D_{{\text{inter}}_{\mathit{sk}}}}{l^2}\frac{{\partial }^2{C}_{s_k}}{\partial Z^2}+e_{{\text{inter}}_k}{K}_{s_k}\frac{D_{{\text{intra}}_{s_k}}}{R^2}{\left(\frac{\partial Q{\left(t,X,Z\right)}_{s_k}}{\partial X}\right)}_{X=1}\right)\mathit{\text{dZdt}}$$ EA.2

$$I_{s_{\mathit{\text{micro}}}}=\underset{0}{\overset{T}{\int }}\underset{0}{\overset{1}{\int }}\underset{L_{k-1}}{\overset{L_k}{\int }}{\psi }_{s_k}\left(t,X,Z\right)\left(\frac{\partial Q_{s_k}\left(t,X,Z\right)}{\partial t}-\frac{D_{int{\mathit{ra}}_{\mathit{sk}}}}{R^2}\left(\frac{{\partial }^2{Q}_{s_k}}{\partial X^2}+\frac{2}{X}\frac{\partial Q_{s_k}}{\partial X}\right)\right)\mathit{\text{XdXdZdt}}$$ EA.3

$J_s$ is the residual functional (17), ${\varphi }_{s_k},{\psi }_{s_k}{,}_{s,=\overline{1,2}}$ —unknown factors of Lagrange, to be determined from the stationary condition of the functional $\mathrm{\Phi }\left(D_{{\text{inter}}_{\mathrm{sk}}},D_{{\text{intra}}_{\mathrm{sk}}}\right)$ [9, 15]:

$$\mathrm{\Delta \Phi }\left(D_{{\text{inter}}_{\mathrm{sk}}},D_{{\text{intra}}_{\mathrm{sk}}}\right)\equiv \mathrm{\Delta }{J}_s+\mathrm{\Delta }{I}_{s_{\mathit{\text{macro}}}}+\mathrm{\Delta }{I}_{s_{\mathit{\text{micro}}}}=0$$ EA.4

The calculation of the components in Eq. (A.4) is carried out by assuming that the values $D_{{\text{inter}}_{\mathrm{sk}}},D_{{\text{intra}}_{\mathrm{sk}}}$ are incremented by $\mathrm{\Delta }{D}_{{\text{inter}}_{\mathrm{sk}}},\mathrm{\Delta }{D}_{{\text{intra}}_{\mathrm{sk}}}$ . As a result, concentration $C_{s_k}\left(t,Z\right)$ changes by increment $\mathrm{\Delta }{C}_{\mathit{sk}}\left(t,Z\right)$ and concentration $Q_{s_k}\left(t,X,Z\right)$ by increment $\mathrm{\Delta }{Q}_{s_k}\left(t,X,Z\right),s=\overline{1,2}$ .

Conjugate problem. The calculation of the increments $\mathrm{\Delta }{J}_s$ , $\mathrm{\Delta }{J}_{s_{\mathit{\text{macro}}}}$ , and $\mathrm{\Delta }{J}_{s_{\mathit{\text{micro}}}}$ in Eq. (A.4) (using integration by parts and the initial and boundary conditions of the direct problem (9)–(15)) leads to solving the additional conjugate problem to determine the Lagrange factors ${\varphi }_{s_k}$ and ${\psi }_{s_k}$ of the functional (A.1) [15]:

$$\frac{\partial {\varphi }_{s_k}\left(t,Z\right)}{\partial t}+\frac{D_{{\text{inter}}_{\mathit{sk}}}}{l^2}\frac{{\partial }^2{\varphi }_{s_k}}{\partial Z^2}+e_{{\text{inter}}_k}{K}_{s_k}\frac{D_{{\text{intra}}_{s_k}}}{R^2}{\left.\frac{\partial {\psi }_{s_k}\left(t,X,Z\right)}{\partial X}\right|}_{X=1}=E_{s_k}^n\left(t\right)\delta \left(Z-h_k\right)$$ EA.5

where $E_{s_k}^n\left(t\right)=C_{s_k}\left({D_{int{\mathit{ra}}_{\mathit{sk}}}}_n^{},{D_{int{\mathit{er}}_{\mathit{sk}}}}_n^{};t,h_k\right)+{\overline{Q}}_{s_k}\left({D_{int{\mathit{ra}}_{\mathit{sk}}}}_n^{},{D_{int{\mathit{er}}_{\mathit{sk}}}}_n^{};t,h_k\right)-M_{\mathit{sk}}\left(t\right)$ ,

$\delta \left(Z-h_k\right)$ (function of Dirac) [15].

$$\frac{\partial {\psi }_{s_k}\left(t,X,Z\right)}{\partial t}+\frac{D_{int{\mathit{ra}}_{\mathit{sk}}}}{R^2}\left(\frac{{\partial }^2{\psi }_{s_k}}{\partial X^2}+\frac{2}{X}\frac{\partial {\psi }_{s_k}}{\partial X}\right)=E_{s_k}^n\left(t\right)\delta \left(Z-h_k\right)$$ EA.6

$${\varphi }_{s_k}{\left(t,Z\right)}_{\left|t=t^{\mathit{\text{total}}}\right.}=0;{\psi }_{s_k}{\left(t,X,Z\right)}_{\left|t=t^{\mathit{\text{total}}}\right.}=0\left(\text{conditions}\mathrm{at}t=t^{\mathit{\text{total}}}\right);$$ EA.7

$$\frac{\partial }{\partial X}{\psi }_{s_k}{\left(t,X,Z\right)}_{\left|X=0\right.}=0;{\psi }_{s_k}{\left(t,X,Z\right)}_{\left|X=1\right.}={\phi }_{s_k}\left(t,Z\right)$$ EA.8

$${\varphi }_{s_k}\left(t,Z=L_k\right)=0,{\varphi }_{s_{k-1}}\left(t,Z=L_{k-1}\right)=0;s=\overline{1,2},k=\overline{N,2},$$ EA.9

$${\varphi }_{s_1}\left(t,L_1\right)=0,\frac{д{\varphi }_{s_1}}{\mathit{дZ}}\left(t,Z=0\right)=0$$ EA.10

We have obtained the solutions ${\varphi }_{\mathit{sk}}$ and ${\psi }_{\mathit{sk}}$ of problem (A.5)–(A.10) using Heaviside operational method in [15].

Substituting in the direct problem (9)–(15) $D_{{\text{inter}}_{\mathrm{sk}}},D_{{\text{intra}}_{\mathrm{sk}}}$ , ${С}_{s_k}\left(t,Z\right)$ , and $Q_{s_k}\left(t,X,Z\right)$ by the corresponding values with increments $D_{{\text{inter}}_{\mathrm{sk}}}+\mathrm{\Delta }{D}_{{\text{inter}}_{\mathrm{sk}}},D_{{\text{intra}}_{\mathrm{sk}}}+\mathrm{\Delta }{D}_{{\text{intra}}_{\mathrm{sk}}}$ , $C_{s_k}\left(t,Z\right)+\mathrm{\Delta }{C}_{s_k}\left(t,z\right)$ , and $Q_{s_k}\left(t,X,Z\right)+\mathrm{\Delta }{Q}_{s_k}\left(t,X,Z\right)$ , subtracting the first equations from the transformed ones and neglecting second-order terms of smallness, we obtain the basic equations of the problem (9)–(15) in terms of increments $\mathrm{\Delta }{C}_{s_k}\left(t,Z\right)$ and $\mathrm{\Delta }{Q}_{s_k}\left(t,X,Z\right),s=\overline{1,2}$ in the operator form

$$\mathcal{L}{w}_{s_k}\left(t,X,Z\right)={\mathrm{X}}_{s_k},w_{s_k}\in \left(0,1\right)\cup {\mathrm{\Omega }}_{\mathit{kt}},k=\overline{1,N+1}$$ EA.11

Similarly, we write the system of the basic equations of conjugate boundary problem (A.5)–(A.10) in the operator:

$${\mathcal{L}}^{\ast }{\mathrm{\Psi }}_{s_k}\left(t,X,Z\right)=E_{s_k}\left(t\right)\delta \left(Z-h_k\right),{\mathrm{\Psi }}_{s_k}\in \left(0,1\right)\cup {\mathrm{\Omega }}_{\mathit{kt}},k=\overline{1,N+1}$$ EA.12

where $\mathcal{L}=\left[\begin{array}{cc}\frac{\partial }{\partial t}-\frac{\partial }{\partial Z}\left(D_{{\text{inter}}_{\mathrm{sk}}}\frac{\partial }{\partial Z}\right)& e_{int{\mathit{er}}_k}\frac{D_{{\text{intra}}_{\mathrm{sk}}}}{R}\frac{\partial }{\partial X}\left|{}_{X=1}\right.\\ 0& \frac{\partial }{\partial t}-\frac{D_{{\text{intra}}_{\mathrm{sk}}}}{R^2}\left(\frac{{\partial }^2}{\partial X^2}+\frac{2}{X}\frac{\partial }{\partial X}\right)\end{array}\right]$ , ${\mathcal{L}}^{\ast }=\left[\begin{array}{cc}\frac{\partial }{\partial t}+\frac{\partial }{\partial Z}\left(D_{{\text{inter}}_{\mathrm{sk}}}\frac{\partial }{\partial Z}\right)& e_{int{\mathit{er}}_k}\frac{D_{{\text{intra}}_{\mathrm{sk}}}}{R^2}\frac{\partial }{\partial X}\left|{}_{X=1}\right.\\ 0& \frac{\partial }{\partial t}+\frac{D_{{\text{intra}}_{\mathrm{sk}}}}{R^2}\left(\frac{{\partial }^2}{\partial X^2}+\frac{2}{X}\frac{\partial }{\partial X}\right)\end{array}\right]$ , $w_{s_k}\left(t,X,Z\right)=\left[\begin{array}{c}\mathrm{\Delta }{C}_{s_k}\left(t,Z\right)\\ \mathrm{\Delta }{Q}_{s_k}\left(t,X,Z\right)\end{array}\right]$ , ${\mathrm{\Psi }}_{s_k}\left(t,X,Z\right)=\left[\begin{array}{c}{\varphi }_{s_k}\left(t,Z\right)\\ {\psi }_{s_k}\left(t,X,Z\right)\end{array}\right]$ .

$${\mathrm{X}}_{s_k}\left(t,X,Z\right)=\left[\begin{array}{c}\frac{\partial }{\partial Z}\left(\mathrm{\Delta }{D}_{{\text{inter}}_{\mathrm{sk}}}\frac{\partial }{\partial Z}{C}_{s_k}\right)-e_{int{\mathit{er}}_k}\frac{\mathrm{\Delta }{D}_{{\text{intra}}_{\mathrm{sk}}}}{R^2}\frac{\partial }{\partial X}{Q}_{s_k}{\left(t,X,Z\right)}_{X=1}\\ \frac{\mathrm{\Delta }{D}_{{\text{intra}}_{\mathrm{sk}}}}{R^2}\left(\frac{{\partial }^2}{\partial X^2}+\frac{2}{X}\frac{\partial }{\partial X}\right)Q_{s_k}\left(t,X,Z\right)\end{array}\right]$$ EA.13

where ${\mathcal{L}}^{\ast }$ is the conjugate Lagrange operator of operator $\mathcal{L}$ .

The calculated increment of the residual functional (17), neglecting second-order terms, has the form

$$\begin{array}{ll}\mathrm{\Delta }{J}_s\left(D_{{\text{intra}}_{\mathrm{sk}}},D_{{\text{inter}}_{\mathrm{sk}}}\right)& =\underset{0}{\overset{T}{\int }}\underset{L_{k-1}}{\overset{L_k}{\int }}{\mathcal{L}}^{-1}{\mathrm{X}}_{s_{k,1}}\left(t,Z\right)\cdot E_{s_k}\left(t\right)\delta \left(Z-h_k\right)\mathit{\text{dZdt}}\\ & +\underset{0}{\overset{T}{\int }}\underset{L_{k-1}}{\overset{L_k}{\int }}\underset{0}{\overset{1}{\int }}{\mathcal{L}}^{-1}{\mathrm{X}}_{s_{k2}}\left(t,X,Z\right)\cdot E_{s_k}\left(t\right)\delta \left(Z-h_k\right)\mathit{\text{XdXdZdt}}\end{array}$$ EA.14

wherе $w_{s_k}={\mathcal{L}}^{-1}{\mathrm{X}}_{s_k}$ and ${\mathcal{L}}^{-1}$ is the inverse operator of operator $\mathcal{L}$ .

Defining the scalar product

$$\left(\mathcal{L}{w}_{s_k}\left(t,X,Z\right),{\mathrm{\Psi }}_{s_k}\left(t,X,Z\right)\right)=\left[\begin{array}{c}\underset{{\mathrm{\Omega }}_{\mathit{kT}}}{\iint }\mathcal{L\Delta }{C}_{s_k}\left(t,Z\right){\varphi }_{s_k}\left(t,Z\right)\mathit{\text{dZdt}}\\ \underset{\left(0,R\right)\cup {\mathrm{\Omega }}_{\mathit{kT}}}{\iiint }\mathcal{L\Delta }{Q}_{s_k}\left(t,X,Z\right){\psi }_{s_k}\left(t,X,Z\right)\mathit{\text{XdXdZdt}}\end{array}\right]$$ EA.15

and taking into account (A.19) Lagrange’s identity [12, 15]

$$\left(\mathcal{L}{w}_{s_k}\left(t,X,Z\right),{\mathrm{\Psi }}_{s_k}\left(t,X,Z\right)\right)=\left(w_{s_k}\left(t,X,Z\right),{\mathcal{L}}^{\ast }{\mathrm{\Psi }}_{s_k}\left(t,X,Z\right)\right)$$ EA.16

and the equality ${\mathcal{L}}^{-1\ast }\left[E_{s_k}\left(t\right)\delta \left(Z-h_k\right)\right]={\mathrm{\Psi }}_{s_k}$ , we obtain the increment of the residual functional expressed by the solution of conjugate problem (A.6)–(A.10) and the vector of the right-hand parts of Eq. (A.13):

$$\mathrm{\Delta }{J}_s\left(D_{{\text{inter}}_{\mathrm{sk}}},D_{{\text{intra}}_{\mathrm{sk}}}\right)=\left({\mathrm{\Psi }}_{s_k}\left(t,X,Z\right),{\mathrm{X}}_{s_k}(t,X,Z)\right)$$ EA.17

where ${\varphi }_{s_k}\left(t,Z\right)$ and ${\psi }_{s_k}\left(t,X,Z\right)$ belong to ${\overline{\mathrm{\Omega }}}_{\mathit{kt}}$ and $\left[0,1\right]\cup {\overline{\mathrm{\Omega }}}_{\mathit{kt}}$ , respectively, ${\mathcal{L}}^{-1\ast }$ is the conjugate operator to inverse operator ${\mathcal{L}}^{-1}$ , and ${\mathrm{\Psi }}_{s_k}$ is the solution of conjugate problem (A.5)–(A.10).

Reporting in Eq. (A.17) the components ${\mathrm{X}}_{s_k}\left(t,X,Z\right)$ taking into account Eq. (A.18), we obtain the formula which establishes the relationship between the direct problem (9)–(15) and the conjugate problem (A.6)–(A.10) which makes it possible to obtain the analytical expressions of components of the residual functional gradient:

$$\mathrm{\Delta }{J}_s\left(D_{{\text{intra}}_{\mathrm{sk}}},D_{{\text{inter}}_{\mathrm{sk}}}\right)=\begin{array}{c}\left({\varphi }_{s_k}\left(t,Z\right),\frac{\partial }{\partial Z}\left(\mathrm{\Delta }{D}_{{\text{inter}}_{\mathrm{sk}}}\frac{\partial }{\partial Z}{C}_{s_k}\right)-e_{int{\mathit{er}}_k}\frac{\mathrm{\Delta }{D}_{{\text{intra}}_{\mathrm{sk}}}}{R^2}\frac{\partial }{\partial X}{Q}_{s_k}{\left(t,X,Z\right)}_{X=1}\right)\\ +\left({\psi }_{s_k}\left(t,X,Z\right),\frac{\mathrm{\Delta }{D}_{{\text{intra}}_{\mathrm{sk}}}}{R^2}\left(\frac{{\partial }^2}{\partial X^2}+\frac{2}{X}\frac{\partial }{\partial X}\right)Q_{s_k}\left(t,X,Z\right)\right)\hfill \end{array}$$ EA.18

Differentiating expression (A.18), by $\mathrm{\Delta }{D}_{{\text{intra}}_{\mathrm{sk}}}$ and $\mathrm{\Delta }{D}_{{\text{inter}}_{\mathrm{sk}}}$ , respectively, and calculating the scalar products according to Eq. (A.15), we obtain the required analytical expressions for the gradient of the residual functional in the intra- and inter-crystallite spaces, respectively:

$$\begin{array}{l}\nabla J_{D_{{\text{intra}}_{\mathrm{sk}}}}\left(t\right)=-\frac{e_{{\text{inter}}_k}}{R^2}\underset{L_{k-1}}{\overset{L_k}{\int }}\frac{\partial }{\partial X}{Q}_{s_k}\left(t,1,Z\right){\varphi }_{s_k}\left(t,Z\right)\mathit{dZ}\\ +\frac{1}{R^2}\underset{L_{k-1}}{\overset{L_k}{\int }}\underset{0}{\overset{1}{\int }}\left(\frac{{\partial }^2}{\partial X^2}+\frac{2}{X}\frac{\partial }{\partial X}\right)Q_{s_k}\left(t,X,Z\right){\psi }_{s_k}\left(t,X,Z\right)\mathit{\text{XdXdZ}}\end{array}$$ EA.19

$$\nabla J_{D_{{\text{inter}}_{\mathrm{sk}}}}\left(t\right)=\underset{L_{k-1}}{\overset{L_k}{\int }}\frac{{\partial }^2{C}_{s_k}\left(t,Z\right)}{\partial Z^2}{\varphi }_{s_k}\left(t,Z\right)\mathit{dZ}$$ EA.20

The formulas of gradients $\nabla J_{D_{{\text{intra}}_{s_k}}}^n\left(t\right)$ and $\nabla J_{D_{{\text{inter}}_{s_k}}}^n\left(t\right)$ include analytical expressions of the solutions of the direct problem (9)–(14) and inverse problem (A.5)–(A.10). They provide high performance of computing process, avoiding a large number of inner loop iterations by using exact analytical methods [2, 15].

B. The linearization schema of the nonlinear co-adsorption model: system of linearized problems and construction of solutions

The linearization schema of nonlinear co-adsorption (1)–(8) is shown in order to demonstrate the simplicity of implementation for the case of two diffusing components ( $m=2$ ) and isothermal adsorption. The simplified model (1)–(8) for the case of $m=2$ is converted into the form

$$\frac{{\mathit{дС}}_s\left(t,Z\right)}{\mathit{дt}}=\frac{D_{{\text{inter}}_s}}{l^2}\frac{{\partial }^2{С}_s}{\partial Z^2}-e_{\text{inter}}{\tilde{K}}_s\frac{D_{{\text{intra}}_s}}{R^2}{\left(\frac{\partial Q_s}{\partial X}\right)}_{X=1}$$ EA.21

$$\frac{{\mathit{дQ}}_s\left(t,X,Z\right)}{\mathit{дt}}=\frac{D_{{\text{intra}}_s}}{R^2}\left(\frac{{\partial }^2{Q}_s}{\partial X^2}+\frac{2}{X}\frac{\partial Q_s}{\partial X}\right),s=\overline{1,2}$$ EA.22

with initial conditions

$$C_s\left(t=0,Z\right)=0;Q_s\left(t=0,X,Z\right)=0;X\in \left(0,1\right),s=\overline{1,2}$$ EA.23

boundary conditions for coordinate X of the crystallite

$$\frac{\partial }{\partial X}{Q}_s\left(t,X=0,Z\right)=0$$ EA.24

$$\begin{array}{l}{Q}_1\left(t,X=1,Z\right)=\frac{K_1{C}_2\left(t,Z\right)}{1+K_1{C}_1\left(t,Z\right)+K_2{C}_2\left(t,Z\right)}\left(\text{Langmuir equilibrium}\right),\\ Q_2\left(t,X=1,Z\right)=\frac{K_2{C}_2\left(t,Z\right)}{1+K_1{C}_1\left(t,Z\right)+K_2{C}_2\left(t,Z\right)}\end{array}$$ EA.25

boundary and interface conditions for coordinate Z

$$C_s\left(t,1\right)=1,\frac{{\mathit{дС}}_s}{\mathit{дZ}}\left(t,Z=0\right)=0,t\in \left(0,t^{\mathit{\text{total}}}\right)$$ EA.26

$K_1=\frac{{\overline{\theta }}_1}{p_1\left(1-{\overline{\theta }}_1-{\overline{\theta }}_2\right)},K_2=\frac{{\overline{\theta }}_2}{p_2\left(1-{\overline{\theta }}_1-{\overline{\theta }}_2\right)}$ , where $p_1,p_2$ are the co-adsorption equilibrium constants and partial pressure of the gas phase for 1-th and 2-th component and ${\overline{\theta }}_1,{\overline{\theta }}_2$ , are the intra-crystallite spaces occupied by the corresponding adsorbed molecules. The expression ${\phi }_s\left(C_1,C_2\right)=\frac{C_s\left(t,Z\right)}{1+K_1{C}_1\left(t,Z\right)+K_2{C}_2\left(t,Z\right)}$ is represented by the series of Tailor [5]:

$${\phi }_s\left(C_1,C_2\right)={\phi }_s\left(0,0\right)+\left({\left.\frac{\partial {\phi }_s}{\partial C_1}\right|}_{\left(0,0\right)}{\mathrm{C}}_1+{\left.\frac{\partial {\phi }_s}{\partial C_2}\right|}_{\left(0,0\right)}{\mathrm{C}}_2\right)+\frac{1}{2!}\left(\begin{array}{l}{\left.\frac{{\partial }^2{\phi }_s}{\partial C_1^2}\right|}_{\left(0,0\right)}{C}_1^2+2{\left.\frac{{\partial }^2{\phi }_s}{\partial C_1{C}_2}\right|}_{\left(0,0\right)}{C}_1{C}_1+\\ \begin{array}{cccc}\begin{array}{ccc}& & \end{array}& & & \end{array}+{\left.\frac{{\partial }^2{\phi }_s}{\partial {\mathrm{C}}_2^2}\right|}_{\left(0,0\right)}{\mathrm{C}}_2^2\end{array}\right)+\dots$$ EA.27

As a result of transformations, limiting to the series not higher than the second order, we obtain

$$\begin{array}{c}\frac{K_1{C}_1\left(t,Z\right)}{1+K_1{C}_1\left(t,Z\right)+K_2{C}_2\left(t,Z\right)}=K_1{C}_1\left(t,Z\right)-\left(K_1^2{C}_1^2\left(t,Z\right)+\frac{1}{2}{K}_1{K}_2{C}_1\left(t,Z\right)C_2\left(t,Z\right)\right),\\ \frac{K_2{C}_2\left(t,Z\right)}{1+K_1{C}_1\left(t,Z\right)+K_2{C}_2\left(t,Z\right)}=K_2{C}_2\left(t,Z\right)-\left(\frac{1}{2}{K}_1{K}_2{C}_1^2\left(t,Z\right)C_2^1\left(t,Z\right)+K_2^2{C}_2^2\left(t,Z\right)\right)\end{array}$$ EA.28

Substituting the expanded expression (A.28) in Eq. (A.25) of nonlinear systems (A.20)–(A.26), we obtain

$$\begin{array}{l}{Q}_1{\left(t,\mathrm{X},Z\right)}_{X=1}=K_1{C}_1\left(t,Z\right)-\epsilon \left(C_1^2\left(t,Z\right)+\frac{1}{2}\frac{K_2}{K_1}{C}_1\left(t,Z\right)C_2\left(t,Z\right)\right),\\ Q_2{\left(t,\mathrm{X},Z\right)}_{X=1}=K_2{C}_2\left(t,Z\right)-\epsilon \left(\frac{1}{2}\frac{K_2}{K_1}{C}_1\left(t,Z\right)C_2\left(t,Z\right)+{\left(\frac{K_2}{K_1}\right)}^2{C}_2^2\left(t,Z\right)\right)\end{array}$$ EA.29

where $\epsilon =K_1^2<<1$ is the small parameter.

Taking into account the approximate equations of the kinetics of co-adsorption (A.29) containing the small parameter $\epsilon$ , we search for the solution of the problem (A.21)–(A.26) by using asymptotic series with a parameter $\epsilon$ in the form [6, 7]

$$C_s\left(t,\mathrm{Z}\right)=C_{s_0}\left(t,\mathrm{Z}\right)+\epsilon C_{s_1}\left(t,\mathrm{Z}\right)+{\epsilon }^2{C}_{s_2}\left(t,\mathrm{Z}\right)+\dots ,$$ EA.30

$$Q_s\left(t,X,Z\right)=Q_{s_0}\left(t,X,Z\right)+\epsilon Q_{s_1}\left(t,X,Z\right)+{\epsilon }^2{Q}_{s_2}\left(t,X,Z\right)+\dots ,s=\overline{1,2}$$ EA.31

As the result of substituting the asymptotic series (A.30)–(A.31) into the equations of the nonlinear boundary problem (A.21)–(A.26) considering Eq. (A.28), the problem (A.21)–(A.26) will be parallelized into two types of linearized boundary problems [6]:

The problem $A_{s_0},s=\overline{1,2}$ (zero approximation with initial and boundary conditions of the initial problem): to find a solution in the domain $D=\left\{\left(t,\mathrm{X},\mathrm{Z}\right):t>0,X\in \left(0,1\right),Z\in \left(0,1\right)\right\}$ of a system of partial differential equations

$$\frac{{\mathit{дС}}_{s_0}\left(t,Z\right)}{\mathit{дt}}=\frac{D_{{\text{inter}}_s}}{l^2}\frac{{\partial }^2{С}_{s_0}}{\partial Z^2}-e_{\text{inter}}{\tilde{K}}_s\frac{D_{{\text{intra}}_s}}{R^2}{\left(\frac{\partial Q_{s_0}}{\partial X}\right)}_{X=1}$$ EA.32

$$\frac{{\mathit{дQ}}_{s_0}\left(t,X,Z\right)}{\mathit{дt}}=\frac{D_{{\text{intra}}_s}}{R^2}\left(\frac{{\partial }^2{Q}_{s_0}}{\partial X^2}+\frac{2}{X}\frac{\partial Q_{s_0}}{\partial X}\right)$$ EA.33

with initial conditions.

$$C_{s_0}\left(t=0,Z\right)=0;Q_{s_0}\left(t=0,X,Z\right)=0;X\in \left(0,1\right),s=\overline{1,2}$$ EA.34

boundary conditions for coordinate X of the crystallite

$$\frac{\partial }{\partial X}{Q}_{s_0}\left(t,X=0,Z\right)=0$$ EA.35

$$Q_{s_0}\left(t,X=1,Z\right)=K_s{C}_{s_0}\left(t,Z\right),s=\overline{1,2}$$ EA.36

boundary and interface conditions for coordinate Z.

$$C_{s_0}\left(t,1\right)=1,\frac{{\mathit{дС}}_{s_0}}{\mathit{дZ}}\left(t,Z=0\right)=0,t\in \left(0,T\right)$$ EA.37

The problem $A_n;n=\overline{1,\infty }$ (nth approximation with zero initial and boundary conditions): to construct in the domain D a solution of a system of equations

$$\frac{{\mathit{дС}}_{s_n}\left(t,Z\right)}{\mathit{дt}}=\frac{D_{{\text{inter}}_s}}{l^2}\frac{{\partial }^2{С}_{s_n}}{\partial Z^2}-e_{\text{inter}}{\tilde{K}}_s\frac{D_{{\text{intra}}_s}}{R^2}{\left(\frac{\partial Q_{s_n}}{\partial X}\right)}_{X=1}$$ EA.38

$$\frac{{\mathit{дQ}}_{s_n}\left(t,X,Z\right)}{\mathit{дt}}=\frac{D_{{\text{intra}}_s}}{R^2}\left(\frac{{\partial }^2{Q}_{s_n}}{\partial X^2}+\frac{2}{X}\frac{\partial Q_{s_n}}{\partial X}\right)$$ EA.39

with initial conditions.

$$C_{s_n}\left(t=0,Z\right)=0;Q_{s_n}\left(t=0,X,Z\right)=0;s=\overline{1,2}$$ EA.40

boundary conditions for coordinate X of the crystallite.

$$\frac{\partial }{\partial X}{Q}_{s_n}\left(t,X=0,Z\right)=0$$ EA.41

$$\begin{array}{l}{Q}_{1_n}{\left(t,X,Z\right)}_{X=1}=K_1{C}_{1_n}\left(t,Z\right)-\sum _{\nu =0}^{n-1}{C}_{1_v}\left(t,\mathrm{Z}\right)\left(C_{1,n-1-\nu }\left(t,\mathrm{Z}\right)+\frac{1}{2}\frac{K_2}{K_1}{C}_{2,n-1-\nu }\left(t,\mathrm{Z}\right)\right),\\ Q_{2_n}{\left(t,X,Z\right)}_{X=1}=K_2{C}_{2_n}\left(t,Z\right)-\sum _{\nu =0}^{n-1}{C}_{2_s}\left(t,\mathrm{Z}\right)\left(\frac{1}{2}\frac{K_2}{K_1}{C}_{1,n-1-\nu }\left(t,\mathrm{Z}\right)+{\left(\frac{K_2}{K_1}\right)}^2{C}_{2,n-1-\nu }\left(t,\mathrm{Z}\right)\right)\end{array}$$ EA.42

boundary and interface conditions for coordinate Z.

$$C_{s_n}\left(t,1\right)=0,\frac{{\mathit{дС}}_{s_0}}{\mathit{дZ}}\left(t,Z=0\right)=0,t\in \left(0,t^{\mathit{\text{total}}}\right)$$ EA.43

The problems $A_{s_0},s=\overline{1,2}$ are linear with respect to zero approximation $C_{s_0},Q_{s_0}$ ; the problems $A_{s_n};\mathrm{n}=\overline{1,\infty }$ are linear with respect to the n-th approximation $C_{s_n},Q_{s_n}$ and nonlinear with respect to all previous n-1 approximations $C_{s_0},\dots ,C_{s_{n-1}}$ .

As demonstrated for the two-component adsorption model (A.21)–(A.26), our proposed methodology can easily be developed and applied to the co-adsorption of any number of gases.