Khuzhayorov Bakhtiyor Khuzhayorovich
| Dzhiyanov Tursunpulot Ortikovich* | Kholliev Fakhriddin Bokhodirovich | Mamatov Shakhobiddin Sayfitdinovich
© 2024 The authors. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).
OPEN ACCESS
The processes of transport of various pollutants in porous media are of great practical importance. Such contaminants can be solid colloidal particles suspended in the carrier fluid. In the process of transfer, particles can be deposited in pores, which significantly change the permeability and porosity characteristics of the medium. The heterogeneity of porous media considerably affects the transfer of these particles. One of the macroscopically inhomogeneous media is zonal inhomogeneous media, consisting of several zones with different characteristics. In such media, generalized mathematical models have not yet been developed that take into account the zonal inhomogeneity of the medium, various linear and nonlinear, reversible and irreversible kinetics of deposition of solid particles from the liquid into the pore space, etc. In this work, a model is generalized for the transfer of solid particles in a two-zone porous medium. In this work, a mathematical model is considered for colloidal particles transport process in a two-zone porous medium and both the zones having the reversible retentions of particles with different characteristics (parameters). It is shown that the nonlinear kinetics of particle deposition, other parameters being equal, leads to an intensification of particle deposition in pores. As the index n decreases from unity, the rate of particle deposition increases in both zones of the medium. As a consequence, the concentration of suspended particles in the mobile fluid in both zones decreases.
approximation, colloidal particles, difference schemes, mathematical model, particles transport, porous medium
A fractured porous media (FPM) is a typical example of structured media [1, 2]. The main paths of movement of a fluid with suspended solid particles or dissolved substances of fracture is assumed. Porous blocks are treated as impermeable for fluid in a simplified model and particles or solute can diffuse into them. The zones are divided in to two in which one is with a mobile fluid (fractures) and the second is with stationary one (porous blocks) having the mass transport process between these two zones.
Spreading of substances in the porous medium in advance manner may raise the results in many factors. Using this phenomenon there may be certain difficulties of making mathematical model and few models in this field were presented by researchers [3-6]. Of course, each approach has its own advantages and disadvantages. If the first-order kinetics is used as a rule without specifying the geometric characteristics of the medium, then in the diffusion approach it is necessary to determine the geometry of the zones and consider the mass transfer between the zones as a diffusion process. Тwo-zone models, where the liquid is mobile in both zones, but with different scales, are more preferable. In addition, in each of the zones there may be different sections where the deposition of particles can occur in completely different ways: equilibrium, non-equilibrium, linear, nonlinear, etc.
Leij and Bradford [7] presented a transport model of colloidal particles in a medium with double porosity is presented, which adds the reversible and irreversible retention of particles, as well as the first-order mass transfer between fractures and porous blocks. The obtained analytical solution was implemented to analyze the results of experiments [8]. Good agreement between theoretical and experimental results was obtained.
Under the assumption of the reversible retention with different characteristics (parameters) in both the zones of the medium it is considered the transfer process of colloidal particles. This method proves to be more fitting because of the irreversible capture of particles, a phenomenon observed solely during the initial phase of the procedure when a monolayer of colloidal particles establishes itself on the medium's rock surface. Importantly, there exists no energy obstacle to hinder this retention [9-11]. Thus, the irreversible retention of particles becomes reversible beyond the certain time limit. Therefore, the double reversible retention kinetics is preferred in this study. It possible to consider the irreversible retention kinetics as a particular case of reversible kinetics. First, we formulate a mathematical model for solid particle transport in two-zone, two-ite fractured porous media. Two kinetics of reversible deposition of solid particles are considered (case I and II): linear and non-linear. Then we pose a transport problem for this model and numerically solve it with using two different finite difference schemes. Relative advantages of these schemes are estimated. On the basis of numerical experiments with discretized model we present results and make some conclusions.
Examine a medium featuring dual porosity as illustrated in Figure 1, revealing the mobility of fluid within one zone (fractures) and the stationary nature of the other zone (porous matrix blocks). In this context, the fluid's mobility is acknowledged in both zones. In such a mobile zone, it has double porosity or double permeability, the transfer of matter occurs with different intensities sometimes contrasting, like the movement of fluid. Similarly, it finds application in macroscopically non-uniform media, where convective matter transfer is possible in both zones.
The first zone, denoted by index 1, exhibits elevated permeability, while the second zone, denoted by index 2, demonstrates reduced permeability. Each zone is further divided into subsections, denoted by indices a and s in the notation.
Suppose two zones form coexisting continua. This allows us to consider one medium, at each point of which two characteristics are determined: Concentration of particles suspended in a fluid, concentration of deposited particles, porosity, diffusion coefficients, fluid filtration rates, etc. Mass transfer of particles occurs between the continua, which is considered proportional to the concentration difference in the continuum. This approach is known as interpenetrating or coexisting continua. Models of two coexisting continua have been effectively used in studying the processes of filtration of homogeneous liquids in fractured-porous media [12]. Now this methodology has become widely used in problems of transport of various substances in inhomogeneous porous media.
Figure 1. A scheme of particle transfer in a two-zone medium
The equations of transfer of particles in the one-dimensional case are written in the form [13]:
$\begin{aligned} & \rho \frac{\partial S_{b i}}{\partial t}+\rho \frac{\partial S_{d i}}{\partial t}+\theta_i \frac{\partial C_i}{\partial t} =\theta_i D_i \frac{\partial^2 C_i}{\partial x^2}-\theta_i v_i \frac{\partial C_i}{\partial x}+\alpha\left(C_n-C_i\right),(i=1,2 ; n=3-i),\end{aligned}$ (1)
where, t is the time, s; x-coordinate, m; Di is longitudinal dispersion coefficients, m2⁄s; vi is the fluid velocity, m⁄s; v1>v2; C1, C2 is the volumetric concentration of particles in a fluid; Sbi and Sdi-concentration of deposited particles, m3⁄kg; θi is porosity of zones; ρ is the bulk density of the medium, kg/m3; α is coefficient of mass transfer between zones, s-1.
Two cases considered for the kinetics of particle deposition.
Case 1. The deposition of particles occurs in each of the sections of the zones reversibly with considering the linear kinetic equations are written as:
$\rho \frac{\partial S_{b i}}{\partial t}=\theta_i k_{a i} C_i-\rho k_{a d i} S_{a i},(i=1,2)$ (2)
$\rho \frac{\partial S_{d i}}{\partial t}=\theta_i k_{s i} C_i-\rho k_{s d i} S_{s i},(i=1,2)$ (3)
Here, kai and ksi represent the deposition coefficients of particles from the fluid phase i onto the solid phase, specifically, onto the solid framework of the environment, measured in s-1. Additionally, kadj and ksdj denote the coefficients for the reattachment of particles from the solid phase and their transition back to the fluid phase, measured in s-1 .
Case 2. The deposition of particles occurs reversibly in each section of the zones with considering the nonlinear kinetic equations:
$\rho \frac{\partial S_{b i}}{\partial t}=\theta_i k_{a i} C_i^n-\rho k_{a d i} S_{b i},(i=1,2)$ (4)
$\rho \frac{\partial S_{d i}}{\partial t}=\theta_i k_{s i} C_i^n-\rho k_{s d i} S_{d i},(i=1,2)$ (5)
where, n(<1) is a real indicator.
The initial and boundary conditions are in the form:
$C_i(0, \mathrm{x})=0, S_{b i}(0, x)=0, S_{d i}(0, x)=0$, (6)
$C_i(t, 0)=c_0$, (7)
$\frac{\partial C_i}{\partial x}(t, \infty)=0, i=1,2.$ (8)
The concentration fields are Ci(t, x), Sbi(t, x), Sdi(t, x) (i=1, 2.) from the Eqs. (1)-(5) with additional conditions Eqs. (6)-(8).
To solve problem Eqs. (1)-(8), we use the finite difference method [14].
Scheme I. In this scheme, in the balance Eq. (1), the mass transfer term $\alpha\left(C_n-C_i\right)$, as well as $\left(C_i\right)$ in Eqs. (2) and (3), are approximated at the lower time layer j. Difference schemes have the form:
$\begin{aligned} \rho \frac{\left(S_{b i}\right)_l^{j+1}-\left(S_{b i}\right)_l^j}{\tau} & +\rho \frac{\left(S_{d i}\right)_l^{j+1}-\left(S_{d i}\right)_l^j}{\tau}+\theta_i \frac{\left.\left(C_i\right)\right)_l^{j+1}-\left(C_i\right)_l^j}{\tau} =\theta_i D_i \frac{\left(C_i\right)_{l-1}^{j+1}-2\left(C_i\right)_l^{j+1}+\left(C_i\right)_{l+1}^{j+1}}{h^2}-\\ &\theta_i v_i \frac{\left(C_i\right)_l^{j+1}-\left(C_i\right)_{l-1}^{j+1}}{h} +\alpha\left(\left(C_n\right)_l^j-\left(C_i\right)_l^j\right), (i=1,2 ; n=3-i),\end{aligned}$ (9)
$\begin{aligned} \rho \frac{\left(S_{b i}\right)_l^{j+1}-\left(S_{b i}\right)_l^j}{\tau}= & \theta_i k_{a i}\left(C_i\right)_l^j-\rho k_{a d i}\left(S_{b i}\right)_l^{j+1} (i=1,2),\end{aligned}$ (10)
$\begin{aligned} & \rho \frac{\left(S_{d i}\right)_l^{j+1}-\left(S_{d i}\right)_l^j}{\tau}=\theta_i k_{s i}\left(C_i\right)_l^j-\rho k_{s d i}\left(S_{d i}\right)_l^{j+1}(i=1,2),\end{aligned}$ (11)
where, $\left(C_i\right)_l^j,\left(S_{b i}\right)_l^j,\left(S_{d i}\right)_l^j$ are the grid values of the functions $C_i(t, x), S_{b i}(t, x), S_{d i}(t, x),(i=1,2)$ at the point $\left(t_j, x_l\right)$.
From the grid Eqs. (10) and (11) we determine $\left(S_{b i}\right)_l^{j+1},\left(S_{d i}\right)_l^{j+1}$,
$\left(S_{b i}\right)_l^{j+1}=p_{i 1}\left(S_{b i}\right)_l^j+p_{i 2},(i=1,2)$ (12)
$\left(S_{d i}\right)_l^{j+1}=q_{i 1}\left(S_{d i}\right)_l^j+q_{i 2},(i=1,2)$ (13)
where,
$\begin{gathered}p_{i 1}=\frac{1}{1+\tau k_{a d i}}, p_{i 2}=\frac{\tau \theta_i k_{a i}}{\rho\left(1+k_{a d i}\right)}\left(C_i\right)_l^j, \\ q_{i 1}=\frac{1}{1+\tau k_{s d i}}, q_{i 2}=\frac{\tau \theta_i k_{s i}}{\rho\left(1+\tau k_{s d i}\right)}\left(C_i\right)_l^j, (i=1,2) .\end{gathered}$
Eq. (9) is reduced in the form:
$\begin{gathered}A_i\left(C_i\right)_{l-1}^{j+1}-B_i\left(C_i\right)_l^{j+1}+E_i\left(C_i\right)_{l+1}^{j+1}=-\left(F_i\right)_l^j\,\,(i=1,2),\end{gathered}$ (14)
where,
$\begin{gathered}A_i=\frac{\theta_i D_i}{h^2}+\frac{\theta_i v_i}{h}, B_i=\frac{\theta_i}{\tau}+\frac{2 \theta_i D_i}{h^2}+\frac{\theta_i v_i}{h}, E_i=\frac{\theta_i D_i}{h^2}, \\ \left(F_i\right)_l^j=\frac{\theta_i}{\tau}\left(C_i\right)_l^j-\frac{\rho}{\tau}\left(\left(S_{b i}\right)_l^{j+1}+\left(S_{d i}\right)_l^{j+1}-\left(S_{b i}\right)_l^j-\left(S_{d i}\right)_l^j\right)+\alpha\left(\left(c_n\right)_l^j-\left(C_i\right)_l^j\right),\,\,(i=1,2 ; n=3-i) .\end{gathered}$
The initial and boundary conditions in a discrete manner take the following form:
$\begin{gathered}\left(C_i\right)_l^0=0,\left(S_{b i}\right)_l^0=0,\left(S_{d i}\right)_l^0=0,\left(S_{d i}\right)_l^0=0, l=\overline{1, L}, \\ \left(C_i\right)_0^j=c_0,\left(C_i\right)_{L-1}^j=\left(C_i\right)_L^j=0, j=\overline{0, J}, i=1,2 .\end{gathered}$
The values of $\left(S_{b i}\right)_l^{j+1},\left(S_{d i}\right)_l^{j+1}$ are obtained using Eqs. (12) and (13). Subsequently, the Thomas’ algorithm is applied to solve the systems of linear algebraic equations given in Eq. (14), yielding the values of $\left(C_i\right)_l^{j+1}$, where $(i=1,2)$. The stability of schemes (12) and (13) is ensured by the fact that $p_{i 1}, q_{i 1}<1$. Moreover, the stability conditions of the Thomas’ algorithm are satisfied for solving (14) for both $(i=1,2)$ [14, 15].
When examining Eq. (1) in conjunction with Case 2, involving Eqs. (4) and (5), only the approximations of the latter equations undergo changes. These approximations follow a similar approach to Eqs. (10) and (11).
$\begin{aligned} & \rho \frac{\left(S_{b i}\right)_l^{j+1}-\left(S_{b i}\right)_l^j}{\tau}= \theta_i k_{a i}\left(C_i^n\right)_l^j-\rho k_{a d i}\left(S_{b i}\right)_l^{j+1},\,\,(i=1,2),\end{aligned}$ (15)
$\begin{gathered}\rho \frac{\left(S_{d i}\right)_l^{j+1}-\left(S_{d i}\right)_l^j}{\tau}=\theta_i k_{s i}\left(C_i^n\right)_l^j-\rho k_{s d i}\left(S_{d i}\right)_l^{j+1}, \,\, (i=1,2),\end{gathered}$ (16)
that can be written as Eqs. (12) and (13) with to determine $\left(S_{b i}\right)_l^{j+1},\left(S_{d i}\right)_l^{j+1}$ the same difference Eqs. (12) and (13) with:
$\begin{gathered}p_{i 1}=\frac{1}{1+\tau k_{a d}}, p_{i 2}=\frac{\tau \theta_i k_{a i}}{\rho\left(1+\tau k_{a d i}\right.}\left(C_i^n\right)_l^j, \\ q_{i 1}=\frac{1}{1+\tau k_{s d}}, q_{i 2}=\frac{\tau \theta_i k_{s i}}{\rho\left(1+\tau k_{s d i}\right)}\left(C_i^n\right)_l^j,(i=1,2) .\end{gathered}$
Eq. (14) retains its original form. The methodology for computing solutions mirrors that employed for Eqs. (12)-(14).
Scheme II. Now we approximate the mass transfer term $\alpha\left(C_n-C_i\right)$ in Eq. (1) and $C_i, C_i^n$ in kinetic Eqs. (2)-(5) on the time layer $j+1$. This leads to SLAE in the Case 1 and -systems of nonlinear algebraic equations in the Case 2.
Eqs. (1)-(3) were approximated as follows:
$\begin{aligned} & \rho \frac{\left(S_{b i}\right)_l^{j+1}-\left(S_{b i}\right)_l^j}{\tau}+\rho \frac{\left(S_{d i}\right)_l^{j+1}-\left(S_{d i}\right)_l^j}{\tau} +\theta_i \frac{\left(C_i\right)_l^{j+1}-\left(C_i\right)_l^j}{\tau} =\theta_i D_i \frac{\left(C_i\right)_{l-1}^{j+1}-2\left(C_i\right)_i^{j+1}+\left(C_i\right)_{l+1}^{j+1}}{h^2}\\ &- \theta_i v_i \frac{\left(C_i\right)_l^{j+1}-\left(C_i\right)_{l-1}^{j+1}}{h}+\alpha\left(\left(C_n\right)_l^{j+1}-\left(C_i\right)_l^{j+1}\right), (i=1,2 ; n=3-i),\end{aligned}$ (17)
$\begin{aligned} \rho \frac{\left(S_{b i}\right)_l^{j+1}-\left(S_{b i}\right)_l^j}{\tau}= & \theta_i k_{a i}\left(C_i\right)_l^{j+1}-\rho k_{a d i}\left(S_{b i}\right)_l^{j+1}\,\, (i=1,2),\end{aligned}$ (18)
$\begin{aligned} \rho \frac{\left(S_{d i}\right)_l^{j+1}-\left(S_{d i}\right)_l^j}{\tau}= & \theta_i k_{s i}\left(C_i\right)_l^{j+1}-\rho k_{s d i}\left(S_{d i}\right)_l^{j+1} \,\,(i=1,2) .\end{aligned}$ (19)
Extracting $\left(S_{b i}\right)_l^{j+1},\left(S_{d i}\right)_l^{j+1}$ are accomplished through the solutions provided by grid Eqs. (18) and (19):
$\left(S_{b i}\right)_l^{j+1}=p_{i 1}\left(S_{b i}\right)_l^j+p_{i 2},(i=1,2)$, (20)
$\left(S_{d i}\right)_l^{j+1}=q_{i 1}\left(S_{d i}\right)_l^j+q_{i 2},(i=1,2)$ (21)
where,
$\begin{aligned} & p_{i 1}=\frac{1}{1+\tau k_{a d i}}, p_{l 2}=\frac{\tau \theta_i k_{a i}}{\rho\left(1+\tau k_{a d i}\right)}\left(C_i\right)_l^{j+1}, \\ & q_{i 1}=\frac{1}{1+\tau k_{s d i}}, q_{i 2}=\frac{\tau \theta_i k_{s i}}{\rho\left(1+\tau k_{s d i}\right)}\left(C_i\right)_l^{j+1},(i=1,2) .\end{aligned}$
Grid Eq. (17) are reduced to the form:
$\begin{gathered}A_i\left(C_i\right)_{l-1}^{j+1}-B_i\left(C_i\right)_l^{j+1}+E_i\left(C_i\right)_{l+1}^{j+1}+G_i\left(C_n\right)_l^{j+1} =-\left(F_i\right)_l^j, \,\, (l=1,2 ; m=3-l),\end{gathered}$ (22)
where,
$\begin{aligned} & A_i=\frac{\theta_i D_i}{h^2}+\frac{\theta_i v_i}{h}, B_i=\frac{\theta_i}{\tau}+\frac{2 \theta_i D_i}{h^2}+\frac{\theta_i v_i}{h}+\frac{\theta_i k_{a i}}{1+\tau k_{a d i}}+\frac{\theta_i k_{s i}}{1+\tau k_{s d i}}+\alpha, \\ & E_i=\frac{\theta_i D_i}{h^2}, G_i=\alpha, \\ & \left(F_i\right)_l^j=\frac{\rho}{\tau}\left(1-\frac{1}{1+\tau k_{a d i}}\right)\left(S_{b i}\right)_l^j+\frac{\rho}{\tau}\left(1-\frac{1}{1+\tau k_{s d i}}\right)\left(S_{d i}\right)_l^j+\frac{\theta_l}{\tau}\left(C_i\right)_l^j \,\, (i=1,2 ; n=3-i) .\end{aligned}$
In contrast to Eq. (14), SLAEs (22) on the $j+1$ layer are connected for i=1 and i=2. To solve these systems, we use the iterative method [14-16]. We construct the iterative scheme in such a way that at $G_1\left(C_2\right)_l^{j+1}$, the term is taken from the lower iteration layer i=1. In the iterative scheme for i=2 the term $G_2\left(C_1\right)_l^{j+1}$ is already taken from the upper iteration layer. Thus, iterative schemes are composed like the Seidel procedure and have the form:
$\begin{gathered}A_1\left(C_1^{(k+1)}\right)_{l-1}^{j+1}-B_1\left(C_1^{(k+1)}\right)_l^{j+1}+E_1\left(C_1^{(k+1)}\right)_{l+1}^{j+1} =-G_1\left(C_2\right)_l^{j+1}-\left(F_1\right)_l^j,\end{gathered}$ (23)
$\begin{array}{r}A_2\left(C_2^{(k+1)}\right)_{l-1}^{j+1}-B_2\left(C_2^{(k+1)}\right)_l^{j+1}+E_2\left(C_2^{(k+1)}\right)_{l+1}^{j+1} =-G_2\left(C_1\right)_l^{j+1}-\left(F_2\right)_l^j,\end{array}$ (24)
where, k is the iteration number, k=0, 1, ....
For Eqs. (20), (21), an iterative scheme can be omitted. The calculation scheme, in contrast to the previous one (Eqs. (12)-(14)), will change. Here, at each time level, Eqs. (23) and (24) are first solved by the Thomas’ algorithm, after reaching the required accuracy of the iterative scheme according to Eqs. (20) and (21) are determined $\left(S_{b i}\right)_l^{j+1}$ and $\left(S_{d i}\right)_l^{j+1}$ using the coefficients $\left(S_{b i}\right)_l^{j+1}$ and $\left(S_{d i}\right)_l^{j+1}$
$p_{i 2}=\frac{\tau \theta_i k_{a i}}{\rho\left(1+\tau k_{a d i}\right)}\left(\stackrel{(k+1)}{C_i}\right)_l^j$, (25)
$q_{i 2}=\frac{\tau \theta_i k_{s i}}{\rho\left(1+\tau k_{s d i}\right)}\left(C_i^{(k+1)}\right)_l^j,(i=1,2)$ (26)
The iterative process ends when the following conditions are satisfied:
$\begin{gathered}\max _l\left|\binom{(k+1)}{C_1}_l^{j+1}-\binom{(k)}{C_1}_l^{j+1}\right|<\varepsilon \text { and } \max _l\left|\binom{(k+1)}{C_2}_l^{j+1}-\binom{(k)}{C_2}_l^{j+1}\right|<\varepsilon,\end{gathered}$ (27)
where, $\varepsilon$ is the specified accuracy of calculations.
In the Case 2, the approximations of Eqs. (4) and (5) are similar to Eqs. (18) and (19) and we have the same Eqs. (20) and (21) with $p_{i 1}, q_{i 1}$ and
$p_{i 2}=\frac{\tau \theta_i k_{a i}}{\rho\left(1+\tau k_{a d i}\right)}\left(C_i^n\right)_l^{j+1}, q_{i 2}=\frac{\tau \theta_i k_{s i}}{\rho\left(1+\tau k_{s d i}\right)}\left(C_i^n\right)_l^{j+1}$, (28)
and in Eq. (22) the coefficients take the form:
$\begin{aligned} & A_i=\frac{\theta_i D_i}{h^2}+\frac{\theta_i v_i}{h}, B_i=\frac{\theta_i}{\tau}+\frac{2 \theta_i D_i}{h^2}+\frac{\theta_i v_i}{h}+\alpha, \\ & E_i=\frac{\theta_i D_i}{h^2}, G_i=\alpha, \\ & \left(F_i\right)_l^j=\frac{\rho}{\tau}\left(\left(S_{b i}\right)_l^j-\left(S_{b i}\right)_l^{j+1}\right)+\frac{\rho}{\tau}\left(\left(S_{d i}\right)_l^j-\left(S_{d i}\right)_l^{j+1}\right)+\frac{\theta}{\tau}\left(C_i\right)_l^j, \,\, (i=1,2) .\end{aligned}$
For $\left(C_i\right)_l^{j+1}, i=1,2$ a system of nonlinear algebraic equations is obtained:
$\begin{gathered}A_i\left(C_i\right)_{l-1}^{j+1}-B_i\left(C_i\right)_l^{j+1}+E_i\left(C_i\right)_{l+1}^{j+1}+K_i\left(C_i^n\right)_l^{j+1} \,\, l=\overline{1, L-1}, j=\overline{0, J-1},\end{gathered}$ (29)
where, $A_i, B_i, E_i, G_i$ remain unchanged, i.e. as in Eq. (22), and
$\begin{gathered}K_i=\frac{\tau \theta_i}{\rho}\left(\frac{k_{a i}}{1+\tau k_{a d i}}+\frac{k_{s i}}{1+\tau k_{s d i}}\right), \\ \left(F_i\right)_l^j=\frac{\rho}{\tau}\left(\left(1-p_{i 1}\right)\left(S_{b i}\right)_l^j+\left(1-q_{i 1}\right)\left(S_{d i}\right)_l^j\right)+\frac{\theta}{\tau}\left(C_i\right)_l^j,(l=1,2) .\end{gathered}$
The iterative scheme for (29) is constructed as follows:
$A_1(\stackrel{(k+1)}{C_1})_{l-1}^{j+1}-B_1(\stackrel{(k+1)}{C_1})_l^{j+1}+E_1(\stackrel{(k+1)}{C_1})_{l+1}^{j+1}=-K_1(\stackrel{(k)}{C^n_1})_l^{j+1}-G_1(\stackrel{(k)}{C_2})_l^{j+1}-\left(F_1\right)_l^j$, (30)
$A_2\left(\stackrel{(k+1)}{C_2}\right)_{l-1}^{j+1}-B_2\left(\stackrel{(k+1)}{C_2}\right)_l^{j+1}+E_2\left(\stackrel{(k+1)}{C_2}\right)_{l+1}^{j+1}=-K_2\left(\stackrel{(k)}{C_2^n}\right)_l^{j+1}-G_2\left(\stackrel{(k+1)}{C_1}\right)_l^{j+1}-\left(F_2\right)_l^j$, (31)
In the k+1 iteration layer, the systems represented by Eqs. (30) and (31) are linear and are effectively addressed using the Thomas’ algorithm. The stability criteria for the Thomas’ algorithm are duly met.
As an initial iterative approximation, we can take $\left(C_i\right)_l^j, i=1,2$. Note that in the Scheme II, the presence of the mass transfer term $\alpha\left(C_n-C_i\right)$ plays a positive role in ensuring the stability conditions for the Thomas’ algorithm when solving systems (23), (24), (30), (31). The scheme of calculating solutions is similar to that for Eqs. (20)-(24). First, from (30), (31), we determine $\left(C_i\right)_l^{j+1}, i=1,2$. Then $\left(S_{b i}\right)_l^{j+1}$, $\left(S_{d i}\right)_l^{j+1}, i=1,2$ from equations similar to Eqs. (20), (21) with coefficients $p_{i 1}, q_{i 1}$ and Eq. (28).
Computations were conducted employing the specified initial parameter values in accordance with the numerical solution outlined above.
$v_1=10^{-4} \mathrm{~m} / \mathrm{s}, v_2=10^{-5} \mathrm{~m} / \mathrm{s}, D_1=v_1 \cdot \beta, D_2=v_2 \cdot \beta, \beta=0,005 \mathrm{~m}, \theta_1=0,1, \theta_2=0,4, k_{a 1}=3 \cdot 10^{-4} \mathrm{~s}^{-1}, k_{a 2}=4 \cdot 10^{-4} \mathrm{~s}^{-1}, \quad k_{s 1}=4 \cdot 10^{-4} \mathrm{~s}^{-1}, k_{s 2}=5 \cdot 10^{-4} \mathrm{~s}^{-1}, \quad \rho=1800 \mathrm{~kg} / \mathrm{m}^3, c_0$$=0,1 \text { and several values of } k_{\text {adl }}, k_{s d l} \text { and } \alpha \text {. }$ The following grid steps are taken $h=0,01 \mathrm{~m}, \tau=0,5 \mathrm{~s}$.
Figures 2-4 depict the outcomes associated with Case 2. In Figure 2, it illustrates the reduction in the indicator n from unity, while holding the remaining parameters constant, resulting in a deceleration of the development of profiles for suspended particle concentrations. In this scenario, the concentration of deposited particles exhibits an accelerated evolution (Figures 3 and 4). To put it differently, a decrease in the n indicator, with the other model parameters remaining constant, leads to an enhancement in the particle deposition in both regions. Consequently, there is a delay in the distribution of particle concentrations in the mobile fluid of both zones.
Comparable outcomes were achieved through the utilization of Scheme II. Without presenting the comprehensive set of results, we will provide comparative visual representations for a single instance (Figure 5), denoting the concentrations calculated by Scheme II as $\bar{C}_i$ Examining the figures reveals the primary distinction in solutions is particularly noticeable for i=2, corresponding to $C_2$. Notably, Scheme II yields underestimated values for $C_2$. Meanwhile, the values for $C_1$ and $\bar{C}_1$ closely align across the entire range of concentration variations along the x-coordinate at a given time.
The devised iterative schemes exhibit favorable convergence. Achieving a solution with the precision of $\varepsilon=10^{-4}$ proved to be satisfactory with a completion of 8-10 iterations for Eqs. (23) and (24), and 28-30 iterations for Eqs. (30) and (31).
Figure 2. Concentration profiles Ci at t=1800s, $\begin{gathered}\alpha=10^{-5} \mathrm{~s}^{-1}, k_{a d 1}=2,5 \cdot 10^{-4} s^{-1}, k_{a d 2}=2 \cdot 10^{-4} s^{-1}, k_{s d 1}=2 \cdot 10^{-4} s^{-1}, k_{s d 2}=10^{-4} s^{-1}\end{gathered}$
Note: …… n=0,8, ----- n=0,9, —— n=1,0
To attain a more elaborate presentation of outcomes for both strategies, namely Scheme I and Scheme II, the subsequent values were computed:
$\begin{aligned} \delta_{i s} & =\left|C_i-\bar{C}_i\right|, \delta_{i S a}=\left|S_{b i}-\bar{S}_{b i}\right|, \\ & \delta_{i S s}=\left|S_{d i}-\bar{S}_{d i}\right|, i=1,2,\end{aligned}$
where, $\bar{S}_{b i}, \bar{S}_{d i}$ represent the values of $S_{b i}$ and $S_{d i}$ obtained in accordance with Scheme II, respectively.
Figure 3. Concentration profiles $\begin{gathered}\alpha=10^{-5} \mathrm{~s}^{-1}, k_{a d 1}=2,5 \cdot 10^{-4} s^{-1}, k_{a d 2}=2 \cdot 10^{-4} s^{-1} k_{s d 1}=2 \cdot 10^{-4} s^{-1}, k_{s d 2}=10^{-4} s^{-1}\end{gathered}$
Note: …… n=0,8, ----- n=0,9, —— n=1,0
Figure 4. Concentration profiles $\begin{gathered}S_{d i}, i=1 ; 2, \text { at } t=1800 \mathrm{~s} \alpha=10^{-5} \mathrm{~s}^{-1}, k_{a d 1}=2,5 \cdot 10^{-4} \mathrm{~s}^{-1}, k_{a d 2}=2 \cdot 10^{-4} \mathrm{~s}^{-1}, k_{s d 1}=2 \cdot 10^{-4} \mathrm{~s}^{-1}, k_{s d 2}=10^{-4} \mathrm{~s}^{-1}\end{gathered}$
Note: …… n=0,8, ----- n=0,9, —— n=1,0
Figure 5. Concentration profiles $\begin{gathered}C_i \text { (Scheme I) and } \bar{C}_i\left(\text { Scheme II) at } \alpha=10^{-5} \mathrm{~s}^{-1}, k_{a d 1}=2,5 \cdot 10^{-4} \mathrm{~s}^{-1}, k_{a d 2}=2 .\right. 10^{-4} \mathrm{~s}^{-1}, k_{s d 1}=2 \cdot 10^{-4} \mathrm{~s}^{-1}, k_{s d 2}=10^{-4} \mathrm{~s}^{-1} t=1800 \mathrm{~s}, n=0,8, l=1,2\end{gathered}$
Figure 6. Surfaces $\begin{gathered}\delta_{1 c}(\mathrm{a}), \delta_{2 c} \text { (b), } \delta_{1 S a} \text { (c), } \delta_{2 S a}(\mathrm{~d}), \delta_{1 S s} \text { (e), } \delta_{2 S s} \text { (f) при } n= 1, \alpha=10^{-5} \mathrm{~s}^{-1} k_{a d 1}=2,5 \cdot 10^{-4} s^{-1}, k_{a d 2}=2 \cdot 10^{-4} s^{-1}, k_{s d 1}=2 . 10^{-4} s^{-1}, k_{s d 2}=10^{-4} s^{-1}\end{gathered}$
Figures 6 and 7 illustrate some representative surfaces, namely $\delta_{i c}, \delta_{i S a}, \delta_{i S s}, i=1,2$. These surfaces reveal discernible maxima in their dependencies on both x and t. Notably, the absolute disparities between the solutions for $S_{b i}$ and $S_{d i}$ are negligible. Substantial distinctions in outcomes are solely observed for $C_i(t, x), i=1,2$. The computations indicate that $\delta_{i c}$ can attain maximum values up to 0.0015 (Figure 7), which is notably significant. In essence, both schemes yield closely aligned results, and the observed discrepancies in solutions can be deemed acceptable.
Figure 7. Surfaces $\begin{gathered}\delta_{1 c} \text { (a), } \delta_{2 c} \text { (b), } \delta_{1 S a} \text { (c), } \delta_{2 S a} \text { (d), } \delta_{1 S s} \text { (e), } \delta_{2 S s} \text { (f) at } n= \alpha=10^{-5} \mathrm{~s}^{-1}, k_{a d 1}=2,5 \cdot 10^{-4} \mathrm{~s}^{-1}, k_{a d 2}=2 \cdot 10^{-4} s^{-1} k_{s d 1}=2 \cdot 10^{-4} s^{-1}, k_{s d 2}=10^{-4} s^{-1}\end{gathered}$
To form a judgment regarding the comparative merits of Scheme I and Scheme II, conducting test computations for a problem featuring a known exact solution is imperative. The numerical solutions should then be juxtaposed with the exact solution. In this context, a preliminary evaluation of Scheme II can be provided, guided by the overarching inclination favoring implicit schemes over explicit ones, even though, computationally speaking, Scheme II is relatively more resource-intensive.
As noted above, the deposition of particles in the pores changes the structure of the pore space and, as a consequence, the porosity and permeability of the medium. The model needs to be improved to account for these phenomena. In this case, you can use the approaches [16-19].
In this paper, a model of solid particle transport in a two-zone porous medium is generalized. In this case, there are two sections with significantly different particle settling characteristics in each section. Here the sedimentation of particles in both parts of the zones is assumed to be reversed and the process is expressed by kinetic equations. Two different schemes were used to solve the generalized model numerically. The solutions obtained using this model were analyzed.
In nonlinear kinetics, it is shown that the particle settling increases when the index decreases to one, without changing the other parameters. As a result, there is a delay in the development of the particle concentration distribution in the mobile fluid. As further step in this work, we determined some parameters of the model used here by solving the inverse problems. In addition, it was necessary to study the possibility of replacing the bicontinuum considered here with a monocontinuum and such work was carried out [13, 20].
[1] Sahimi, M. (2011). Flow and Transport in Porous Media and Fractured Rock: From Classical Methods to Modern Approaches. John Wiley & Sons.
[2] van Golf-Racht, T.D. (1982). Fundamentals of Fractured Reservoir Engineering. Elsevier, 12: 732.
[3] Gerke, H.H., van Genuchten, M.T. (1996). Macroscopic representation of structural geometry for simulating water and solute movement in dual-porosity media. Advances in Water Resources, 19(6): 343-357. https://doi.org/10.1016/0309-1708(96)00012-7
[4] Selim, H.M., Ma, L. (1998). Physical Nonequilibrium in Soils: Modeling and Application. CRC Press.
[5] Šimůnek, J., van Genuchten, M.T. (2008). Modeling nonequilibrium flow and transport processes using HYDRUS. Vadose Zone Journal, 7(2): 782-797. https://doi.org/10.2136/vzj2007.0074
[6] Toride, N., Leij, F.J., Van Genuchten, M.T. (1995). The CXTFIT code for estimating transport parameters from laboratory or filed tracer experiments. Riverside, CA: US Salinity Laboratory.
[7] Leij, F.J., Bradford, S.A. (2013). Colloid transport in dual-permeability media. Journal of Contaminant Hydrology, 150: 65-76. https://doi.org/10.1016/j.jconhyd.2013.03.010
[8] Bradford, S.A., Bettahar, M., Simunek, J., Van Genuchten, M.T. (2004). Straining and attachment of colloids in physically heterogeneous porous media. Vadose Zone Journal, 3(2): 384-394. https://doi.org/10.2113/3.2.384
[9] Bradford, S.A., Torkzaban, S. (2008). Colloid transport and retention in unsaturated porous media: A review of interface‐, collector‐, and pore‐scale processes and models. Vadose Zone Journal, 7(2): 667-681. https://doi.org/10.2136/vzj2007.0092
[10] Elimelech, M., Gregory, J., Jia, X. (2013). Particle Deposition and Aggregation: Measurement, Modelling and Simulation. Butterworth-Heinemann.
[11] Gitis, V., Rubinstein, I., Livshits, M., Ziskind, G. (2010). Deep-bed filtration model with multistage deposition kinetics. Chemical Engineering Journal, 163(1-2): 78-85. https://doi.org/10.1016/j.cej.2010.07.044
[12] Barenblatt, G.I., Entov, V.M., Ryzhik, V.M. (1990). Theory of Fluid Flows Through Natural Rocks. Springer Dordrecht. Boston, London.
[13] Khuzhayorov, B.K., Dzhiyanov, T.O., Eshdavlatov, Z.Z. (2022). Numerical investigation of solute transport in a non-homogeneous porous medium using nonlinear kinetics. International Journal of Mechanical Engineering and Robotics Research, 11(2): 79-85. https://doi.org/10.18178/ijmerr.11.2.79-85
[14] Samarskii, A.A. (2001). The Theory of Difference Schemes. CRC Press. https://doi.org/10.1201/9780203908518
[15] Samarskii, A.A., Nikolaev, E.S. (1989). Methods for Grid Equations. Volume II Iterative Methods Translated from the Russian by Stephen G. Nash. Basel, Boston, Berlin, Birkhauser, p. 506.
[16] Samarskiy A.A., Vabishchevich P.N. (1995) Computational Heat Transfer. Volume 2: Finite Difference Methodology. John Wiley & Sons, p. 421.
[17] Khuzhaerov, B. (1990). Effects of blockage and erosion on the filtration of suspensions. Journal of Engineering Physics, 58: 185-190.
[18] Khuzhayorov, B.K. (1999). A model of colmatation suffosion filtration. Journal of Porous Media, 2(2). https://doi.org/10.1615/JPorMedia.v2.i2.40
[19] Khuzhayorov, B.K., Djiyanov, T.O., Yuldashev, T.R. (2019). Anomalous nonisothermal transfer of a substance in an inhomogeneous porous medium. Journal of Engineering Physics and Thermophysics, 92: 104-113. https://doi.org/10.1007/s10891-019-01912-y
[20] Khuzhayorov, B.K., Mustafokulov, J.A., Dzhiyanov, T.O., Zokirov, M.S. (2022). Solute transport with non-equilibrium adsorption in a non-homogeneous porous medium. WSEAS Transactions on Fluid Mechanics, 17: 181-188. https://doi.org/10.37394/232013.2022.17.18
