Numerical Study of Suspension Filtration Model in Porous Medium with Modified Deposition Kinetics

Abstract

Filtration is one of the most used technologies in chemical engineering. Development of computer technology and computational mathematics made it possible to explore such processes by mathematical modeling and computational methods. The article deals with the study of suspension filtration in a porous medium with modified deposition kinetics. It is suggested that deposition is formed in two types, reversible and irreversible. The model of suspension filtration in porous media consists of the mass balance equation and kinetic equations for each type of deposition. The model includes dynamic factors and multi-stage deposition kinetics. By using the symmetricity of porous media, the higher dimensional cases are reduced to the one-dimensional case. To solve the problem, a stable, effective and simple numerical algorithm is proposed based on the finite difference method. Sufficient conditions for stability of schemes are found. Based on numerical results, influences of dynamic factors on solid particle transport and deposition characteristics are analyzed. It is shown that the dynamic factors mainly affect the profiles of changes in the concentration of deposition of the active zone.

1. Introduction

The problem of filtering disperse systems in a porous medium is of big practical importance. Many technological processes, natural phenomena, production processes, etc. are associated with the flow of dispersed systems in porous [1,2] and fractured-porous media [3,4]. In contrast to the filtration of homogeneous liquids, when filtering disperse systems, a number of new phenomena arise, the study of which is very important for understanding the mechanisms of the filtration process. If we consider suspensions as dispersed, then the most important question is to study the processes of transport of the dispersed phase of solid particles in a porous medium, the interaction of these particles with the surface of the skeleton of a porous medium, deposition in the pore space, the release of pore channels from deposited particles, etc. Note that such questions were studied very poorly from both the theoretical and experimental sides. In industrial and natural conditions, there are only a few studies. Perhaps one of the reasons for the poor knowledge of the problem is the extreme complexity of the physical and mathematical modeling of particle transport processes in pore space [5]. Fast development of computer technology and computational mathematics in recent years allow us to use more complex models to describe the process.

It is obvious that the study of the filtration of disperse systems taking into account the deposition and release of particles in addition to the experimental direction should be based on mathematical models [6]. However, the specifics of filtering dispersed systems in comparison with filtering conventional homogeneous systems are not always possible to take into account in mathematical models. Moreover, the mechanism of the deposition of particles in pores and their release by the filtration flow has not yet been sufficiently studied, although experimental data convincingly show clogging of the pores and it is possible to fix the distribution of the impurity concentration along the length of the filtration zone [7,8]. The theoretical basis of filtration of inhomogeneous fluids can be found in the books of Ives [9] and Shekhtman [10].

From the above, we see that one of the most important tasks in the study of dispersed system filtration is the developing of mathematical models of the process that adequately reflect the basic characteristic properties of the process under study. Various approaches have been developed to simulate the suspended particle transport and deposition process in porous media or in literature, the so-called deep bed filtration, which can be arbitrarily divided into the so-called empirical models [11,12,13], network models [14,15,16], models for analyzing particle trajectories [17,18,19] and stochastic models [20,21,22].

In the empirical (phenomenological) approach, the mathematical model for suspension filtration in porous media consists of mass balance equation and capture kinetics [11,12,23,24,25,26,27,28,29,30]. Balance equation for three dimensional case is in the following form

$$\frac{\partial (mc+\rho )}{\partial t}+v\left(\frac{\partial c}{\partial x}+\frac{\partial c}{\partial y}+\frac{\partial c}{\partial z}\right)=0,$$

(1)

where c is the concentration of the suspension, v is the filtration velocity, m is the porosity of the medium, $\rho$ is the concentration of deposition, t is the time, $x, y, z$ are coordinates of space. This equation is symmetric with respect to $x, y, z$.

As mentioned in [31] the filters used in water purification have two main categories: symmetric (homogeneous) and asymmetric (anisotropic). The pores in symmetric filters have close to a uniform diameter throughout the depth of the filter. So, in homogeneous porous media all characteristics are the same in all directions of coordinates $(x, y, z)$. One can assume that particle concentration in the stream c and the velocity v is uniform across the filter cross-section. Accordingly, the expected depositions are uniform across the cross-section, too. These common assumptions mean that the model can be developed as one-dimensional in space. Balance equation in one-dimensional case takes the form [11,12,23,24,25,26,27,28,29,30]

$$\frac{\partial (mc+\rho )}{\partial t}+v\frac{\partial c}{\partial x}=0,$$

(2)

Kinetic equation is in following form [11,23,24,25,26,27,28]

$$\frac{\partial \rho }{\partial t}=\lambda \left(\rho \right)c,$$

(3)

where $\lambda \left(\rho \right)$ is the coefficient characterizing the kinetics.

The first kinetic equation for deep bed filtration is given by Iwasaki [32]. Kinetic equations given in the literature can be divided into two approaches [6,12]: irreversible deposition, which means deposited particles never get back to flow again and reversible deposition, meaning some of the deposited particles may be detached from the media grain by the flow. Irreversible kinetic equations are usually given in the form (3), where $\lambda \left(\rho \right)$ can be changed during the filtration. $\lambda \left(\rho \right)$ can be increasing [32] or decreasing function [33], but in some cases a combination of both kinds of behaviour mentioned above [34] are used. First reversible kinetic equations given by Mints [35] mentioned and experimentally showed that during the filtration due to the increasing of pressure gradient, detachment of less strongly linked particles from grain occurs. This kind of process can be described in the form $\partial \rho /\partial t={\beta }_{att}c-{\beta }_{det}\rho ,$ where ${\beta }_{att}$ and ${\beta }_{det}$ phenomenological attachment and detachment coefficients respectively. In [12], given multistage deposition kinetics, in the first stage an irreversible layer is formed, then in the operable stage reversible deposition is formed.

In [12,28,29] mass balance equation is given with diffusion equation

$$\frac{\partial (mc+\rho )}{\partial t}+v\frac{\partial c}{\partial x}=D\frac{{\partial }^{2}c}{\partial x^2},$$

(4)

where D is diffusion (dispersion) coefficient.

In the last decade fractional diffusion equation is used for fractural porous media [3,4,36].

In the present paper, we consider the problem of filtering a suspension in a porous medium, taking into account the deposition of solid particles of the suspension in the pore volume and their release. For this we use the well-known model of Venitsianov [23,24] and suggest its modification. The problem is considered taking into account dynamic factors in the kinetics of deposition. To solve the problem, a numerical method is developed based on finite difference scheme. The stability of the method is proved. By using the symmetry properties of porous media, the model was reduced to a one-dimensional case. The numerical algorithm developed in the present paper can be easily used for 2D and 3D cases by using “alternating direction method”, which means the values of functions in new time layer will be found separately for each coordinate. Based on the numerical solution of the problem, the role of dynamic factors on the filtration and transport characteristics is estimated.

2. Formulation of the Problem

The deposition in the porous space of the deep bed filters has two forms—washable and nonwashable. Accordingly, the filter zones are called active and passive [23,24]. The active zones washed by the flow form a reversible deposition with a concentration of ${\rho }_a$, passive zones that are stagnant form an irreversible deposition do the same with a concentration of ${\rho }_p$. Denote the total filter capacity by ${\rho }_0$. From the foregoing

$${\rho }_0={\rho }_{a0}+{\rho }_{p0}$$

(5)

where ${\rho }_{a0}$ and ${\rho }_{p0}$ are the capacities of the active and passive zones, respectively. The indicated capacities have dynamic characteristics. They depend not only on the “quality” of the dispersed phase but also on the speed and structure of the flow, as well as the geometry of the layer [23].

We consider a semi-infinite homogeneous media with an initial porosity $m_0$, filled with a homogeneous liquid. At the point $x=0$, starting from $t=0$ to when the reservoir enters the suspension with a concentration $c_0$ and filtration velocity $v\left(t\right)=v_0=\mathrm{const}$.

The system of equations for suspension filtration with two deposition zones with a constant speed regime consists of the equation of balance and kinetic equations for each zone, which in the one-dimensional case is represented as

$$m_0\frac{\partial c}{\partial t}+v\frac{\partial c}{\partial x}+\frac{\partial {\rho }_a}{\partial t}+\frac{\partial {\rho }_p}{\partial t}=0,$$

(6)

$$\frac{\partial {\rho }_a}{\partial t}={\beta }_a\left(c-\frac{{\rho }_a}{{\rho }_{a0}}{c}_0\right),$$

(7)

$$\frac{\partial {\rho }_p}{\partial t}={\beta }_p\left({\rho }_p\right)c,$$

(8)

where c is the concentration of the suspension (m${}^3$/m${}^3$), v is the filtration velocity (m/s), $m_0$ is the initial porosity of the medium, ${\rho }_a$ is the concentration of deposition in the active zone (m${}^3$/m${}^3$), ${\rho }_p$ is the concentration of deposition in the passive zone (m${}^3$/m${}^3$), ${\beta }_a$ is the coefficient characterizing the kinetics in the active zone (1/s), ${\beta }_p$ is the coefficient (1/s) associated with the effect of compaction (aging) of the deposition, which was proposed as ${\beta }_p=\alpha \left({\rho }_p\right){\beta }_{p0}$ [23],

$$\alpha \left({\rho }_p\right)=\left\{\begin{array}{ccc}1\hfill & \mathrm{for}\hfill & {\rho }_p\le {\rho }_{p1},\hfill \\ {\rho }_{p1}/{\rho }_p\hfill & \mathrm{for}\hfill & {\rho }_{p1}<{\rho }_p<{\rho }_{p0} ,\hfill \\ 0\hfill & \mathrm{for}\hfill & {\rho }_p={\rho }_{p0} ,\hfill \end{array}\right.$$

(9)

where ${\rho }_{p1}$ is the value of ${\rho }_p$ at which “aging” begins.

Therefore, at the beginning of deposition formation, $\alpha =1$. Starting with a certain concentration of ${\rho }_{p1}$, the value of $\alpha$ becomes less than 1, and a further decrease in $\alpha$ goes inversely with the amount of deposition ${\rho }_p$. Finally, when the deposition concentration is close to saturation, $\alpha$ decreases more rapidly. This section is approximated by a step.

The system of Equations (6)–(8) is solved taking into account the above dependence (9) and with the following initial and boundary conditions

$$c(x,0)=0,{\rho }_a(x,0)={\rho }_p(x,0)=0,c(0,t)=c_0=\mathrm{const}.$$

(10)

3. Finite Difference Schemes for Model

To solve problem (6)–(10), we apply the finite difference method [37,38,39,40,41]. In the area $D=\left\{0\le x<\infty ,0\le t\le T\right\}$ we introduce a net, where T is the maximum time during which the process is studied. For this, we divide the interval $\left[0, \infty \right)$ in steps h and $\left[0, T\right]$ into J parts in steps of $\tau$. As a result we have a net

$${\omega }_{h\tau }=\left\{\left(x_i, t_j\right), x_i=ih,i=0,1,\dots ,t_j=j\tau , j=0,1,\dots ,J,\tau =T/\phantom{TJ}J\right\}.$$

(11)

Instead of functions $c\left(t, x\right)$, ${\rho }_a\left(t, x\right)$, ${\rho }_p\left(t, x\right)$ we consider the net functions and denote their values at the nodes $\left(x_i,t_j\right)$ by $c_i^j$, ${\rho }_{a,i}^j$, ${\rho }_{p,i}^j$, respectively.

Equation (6) is approximated on a grid ${\omega }_{h\tau }$ by the equation

$$m_0\frac{c_i^{j+1}-c_i^j}{\tau }+v\frac{c_i^{j+1}-c_{i-1}^{j+1}}{h}+\frac{{\rho }_{a,i}^{j+1}-{\rho }_{a,i}^j}{\tau }+\frac{{\rho }_{p,i}^{j+1}-{\rho }_{p,i}^j}{\tau }=0.$$

(12)

The difference scheme for Equation (7) will have the form

$$\frac{{\rho }_{a,i}^{j+1}-{\rho }_{a,i}^j}{\tau }={\beta }_a\left(c_i^j-\frac{{\rho }_{a,i}^j}{{\rho }_{a0}}{c}_0\right).$$

(13)

The difference scheme for Equation (8) is

$$\frac{{\rho }_{p,i}^{j+1}-{\rho }_{p,i}^j}{\tau }=\alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}{c}_i^j.$$

(14)

The initial and boundary conditions (10) are also presented in a net form as follows

$$\begin{array}{c}{\rho }_{a,i}^j=0,i=\overline{0,I},j=0,\hfill \\ {\rho }_{p,i}^j=0,i=\overline{0,I},j=0,\hfill \\ c_i^j=0,i=\overline{0,I},j=0,\hfill \\ c_i^j=c_0,i=0,j=\overline{0,J},\hfill \end{array}$$

(15)

where I is a sufficiently large number for which the equation $c_I^j=0$ is approximately satisfied.

4. Stability of the Finite Difference Schemes

Theorem 1.

The sufficient stability conditions for the initial data of the schemes (12) and (13) are

$$\left\{\begin{array}{c}\tau \le min\left\{\frac{2{\rho }_{a0}}{{\beta }_a{c}_0}; \frac{2m_{0}h}{{\beta }_{a}h-v+h\alpha \left({\rho }_{p,i}^j\right){\beta }_0}\right\},\hfill \\ v<h({\beta }_a+\alpha \left({\rho }_{p,i}^j\right){\beta }_0).\hfill \end{array}\right.$$

(16)

and the scheme (14) is unconditionally stable.

Proof of Theorem 1.

Schemes (12)–(14) are two-layer schemes. When studying the stability of two-layer schemes, we will use their canonical form [37]

$$B{y}_t+Ay=\phi \left(t\right),t=t_n=n\tau \in {\omega }_{\tau },y\left(0\right)=y_0.$$

(17)

where $y=y_n=y\left(t_n\right),\widehat{y}=y_{n+1}=y\left(t_{n+1}\right)=y(t_n+\tau ), y_t=\frac{\widehat{y}-y}{\tau }$, $y_n$ is vector-solution defined on the layer $t_n$,

$$y_n=\left\{y_0^n, y_1^n, \dots , y_I^n\right\}.$$

(18)

Let $H_h$ be a finite-dimensional real space, ( , ) be the scalar product, $∥x∥=\sqrt{(x,x)}$ is the norm in $H_h$.

It is necessary to find sufficient conditions for the stability of the scheme (17) and obtain a priori estimates for its solution, expressing the stability of the scheme from the initial data.

The solution to the problem (17) can be represented in the form $y=\overline{y}+\tilde{y},$ where $\overline{y}$ is a solution to the homogeneous equation with the initial condition $\overline{y} \left(0\right)=y\left(0\right)=y_0:$

$$B{y}_t+Ay=0,t\in {\omega }_{\tau },y\left(0\right)=y_0,$$

(19)

and $\tilde{y}$ is a solution of a nonhomogeneous equation with zero initial condition

$$B{y}_t+Ay=\phi \left(t\right),t={\omega }_{\tau },y\left(0\right)=0.$$

(20)

Estimation of the solution of problem (19)

$${∥y(t+\tau )∥}_{\left(1\right)}\le M_1{∥y\left(0\right)∥}_{\left(1\right)}$$

(21)

means that the scheme (17) is stable according to the initial data, and the estimate of the solution to the problem (20)

$${∥y(t+\tau )∥}_{\left(1\right)}\le M_2\underset{0\le t^{\prime }\le t}{max}{∥\phi \left(t^{\prime }\right)∥}_{\left(2\right)}$$

(22)

expresses the stability of the scheme (17) on the right hand side.

It is known [37] that the uniform stability of the scheme (19) from the initial data implies the stability of the scheme (20) on the right-hand side. Therefore, it is enough to show the uniform stability of the scheme (19). For this, it is necessary to obtain an a priori estimate such as

$$∥y_{n+1}∥\le \delta ∥y_n∥,$$

for all n and for all $y_n$, where ${\delta }^n\le M_1$, $\delta$ and $M_1$ are independent of the net steps and n.

To study stability, a method based on an estimate of the norm of the transition operator from layer to layer can be applied.

We write the difference scheme (19) in the form

$$\widehat{y}=Sy,S=E-\tau B^{-1}A,$$

(23)

where S is transition operator.

If $A=A^{*}>0$, we can obtain the stability condition in the energy norms ${∥·∥}_A$. From the inequality

$${∥\widehat{y}∥}_A={∥Sy∥}_A\le {∥S∥}_A{∥y∥}_A,$$

it can be seen that the scheme (23) is stable in $H_A$, if the norm of the transition operator does not exceed unity:

$${∥S∥}_A\le 1.$$

(24)

We will check the inequality (24) for transition operators for (12)–(14).

The scheme (12) can be represented as

$$m_0\frac{c_i^{j+1}-c_i^j}{\tau }+v\frac{c_i^{j+1}-c_{i-1}^{j+1}}{h}+{\beta }_a\left(c_i^j-\frac{{\rho }_{a,i}^j}{{\rho }_{a0}}{c}_0\right)+\alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}{c}_i^j=0.$$

We transform this equation to obtain

$$m_{0}h(c_i^{j+1}-c_i^j)+v\tau (c_i^{j+1}-c_{i-1}^{j+1})+\tau h{\beta }_a\left(c_i^j-\frac{{\rho }_{a,i}^j}{{\rho }_{a0}}{c}_0\right)+\tau h\alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}{c}_i^j=0,$$

$$(m_{0}h+v\tau )c_i^{j+1}=v\tau c_{i-1}^{j+1}+\left(m_{0}h-\tau h{\beta }_a-\tau h\alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}\right)c_i^j-\tau h{\beta }_a\frac{{\rho }_{a,i}^j}{{\rho }_{a0}}{c}_0,$$

$$c_i^{j+1}=\frac{\left(m_{0}h-\tau h{\beta }_a-\tau h\alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}\right)}{m_{0}h+v\tau }{c}_i^j+f,$$

(25)

where $f=\frac{1}{m_{0}h+v\tau }\left(v\tau c_{i-1}^{j+1}-\tau h{\beta }_a\frac{{\rho }_{a,i}^j}{{\rho }_{a0}}{c}_0\right).$ Introducing the solution vector $c_j=\left\{c_0^j, c_1^j, \dots , c_I^j\right\}$, from (25) discarding f (since stability is checked against the initial data) we get

$$∥c_{j+1}^{}∥\le ∥S∥∥c_j∥,$$

where

$$∥S∥=\left|\frac{m_{0}h-\tau h{\beta }_a-\tau h\alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}}{m_{0}h+v\tau }\right|\le 1.$$

This inequality is split down into two inequalities. The first one is

$$\frac{m_{0}h-\tau h{\beta }_a-\tau h\alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}}{m_{0}h+v\tau }\le 1,$$

$$m_{0}h-\tau h{\beta }_a-\tau h\alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}\le m_{0}h+v\tau ,$$

$$-\tau h{\beta }_a-\tau h\alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}\le v\tau ,$$

$$-h\left({\beta }_a+\alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}\right)\le v.$$

The validity of this inequality is obvious, since the left hand side is negative and the right hand side is positive. So let us look at the second inequality

$$-1\le \frac{m_{0}h-\tau h{\beta }_a-\tau h\alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}}{m_{0}h+v\tau },$$

$$-m_{0}h-v\tau \le m_{0}h-\tau h{\beta }_a-\tau h\alpha \left({\rho }_{p,i}^j\right){\beta }_{p0},$$

$$\tau \left(h{\beta }_a+h\alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}-v\right)\le 2m_{0}h.$$

By the second inequality in (16) $v<h({\beta }_a+\alpha \left({\rho }_{p,i}^j\right){\beta }_0)$, and hence

$$\tau \le \frac{2m_{0}h}{{\beta }_{a}h-v+h\alpha \left({\rho }_{p,i}^j\right){\beta }_0},$$

We transform (13) to obtain

$${\rho }_{a,i}^{j+1}={\rho }_{a,i}^j+\tau {\beta }_a\left(c_i^j-\frac{{\rho }_{a,i}^j}{{\rho }_{a0}}{c}_0\right),$$

whence

$${\rho }_{a,i}^{j+1}={\rho }_{a,i}^j\left(1-\frac{\tau {\beta }_a{c}_0}{{\rho }_{a0}}\right)+\tau {\beta }_a{c}_i^j.$$

Ignoring the term $\tau {\beta }_a{c}_i^j$ and introducing the solution vector ${\rho }_{a,j}=\left\{{\rho }_{a,0}^j, {\rho }_{a,1}^j, \dots , {\rho }_{a,I}^j\right\}$ the scheme is written as

$${\rho }_{a,j+1}^{}=S·{\rho }_{a,j}^{},$$

where $S=\left(1-\frac{\tau {\beta }_a{c}_0}{{\rho }_{a0}}\right)E$, E is unit operator.

Condition $∥S∥\le 1$ gives $\left|1-\frac{\tau {\beta }_a{c}_0}{{\rho }_{a0}}\right|\le 1,$ from which we can find $\tau \le \frac{2{\rho }_{a0}}{{\beta }_a{c}_0}$.

From (14) we obtain

$${\rho }_{p,i}^{j+1}={\rho }_{p,i}^j+\tau \alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}{c}_i^j$$

If ${\rho }_{p,i}^j\le {\rho }_{p1}$, then $\alpha \left({\rho }_{p,i}^j\right)=1$ and hence ${\rho }_{p,i}^{j+1}={\rho }_{p,i}^j+\tau {\beta }_{p0}{c}_i^j$. Ignoring the term $\tau {\beta }_{p0}{c}_i^j$ and by introducing the solution vector ${\rho }_{p,j}=\left\{{\rho }_{p,0}^j, {\rho }_{p,1}^j, \dots , {\rho }_{p,I}^j\right\}$ last formula can be written in the form

$${\rho }_{p,j+1}^{}=S·{\rho }_{p,j}^{},$$

where $S=E$, $∥S∥=1$.

If ${\rho }_{p1}<{\rho }_{p,i}^j<{\rho }_{p0}$, then $\alpha \left({\rho }_{p,i}^j\right)=\frac{{\rho }_{p1}}{{\rho }_{p,i}^j}$,

$${\rho }_{p,i}^{j+1}={\rho }_{p,i}^j+\tau \frac{{\rho }_{p1}}{{\rho }_{p,i}^j}{\beta }_{p0}{c}_i^j.$$

Since ${\rho }_{p1}<{\rho }_{p,i}^j<{\rho }_{p0}$, then $0<\frac{{\rho }_{p1}}{{\rho }_{p,i}^j}<1$, we will take the maximum value,

$${\rho }_{p,i}^{j+1}={\rho }_{p,i}^j+\tau {\beta }_{p0}{c}_i^j.$$

$${\rho }_{p,j+1}^{}=S·{\rho }_{p,j}^{},S=E,∥S∥=1.$$

If ${\rho }_{p,i}^j={\rho }_{p0}$, then $\alpha \left({\rho }_{p,i}^j\right)=0$ and ${\rho }_{p,i}^{j+1}={\rho }_{p,i}^j$, ${\rho }_{p,j+1}^{}=S·{\rho }_{p,j}^{}$,

$$S=E,∥S∥=1.$$

□

Let us look to another difference scheme. Let us introduce implicit schemes for (7) and (8) in the form

$$\frac{{\rho }_{a,i}^{j+1}-{\rho }_{a,i}^j}{\tau }={\beta }_a\left(c_i^{j+1}-\frac{{\rho }_{a,i}^{j+1}}{{\rho }_{a0}}{c}_0\right),$$

(26)

$$\frac{{\rho }_{p,i}^{j+1}-{\rho }_{p,i}^j}{\tau }=\alpha \left({\rho }_{p,i}^{j+1}\right){\beta }_{p0}{c}_i^j$$

(27)

instead of explicit schemes (13) and (14).

Theorem 2.

Scheme (26) is unconditionally stable, but scheme (27) is unstable.

Proof Theorem 2.

Since

$${\rho }_{a,i}^{j+1}\left(1+\frac{\tau c_0}{{\rho }_{a0}}\right)={\rho }_{a,i}^j+\tau {\beta }_a{c}_i^{j+1},$$

we obtain

$${\rho }_{a,i}^{j+1}={\left(1+\frac{\tau c_0}{{\rho }_{a0}}\right)}^{-1}{\rho }_{a,i}^j+\tau {\beta }_a{c}_i^{j+1}.$$

Clearly,

$${\left(1+\frac{\tau c_0}{{\rho }_{a0}}\right)}^{-1}<1.$$

Therefore scheme (26) is unconditionally stable.

Next, for (27), if

$${\rho }_{p1}<{\rho }_{p,i}^{j+1}<{\rho }_{p0},then\alpha \left({\rho }_{p,i}^{j+1}\right)=\frac{{\rho }_{p1}}{{\rho }_{p,i}^{j+1}}.$$

Using the Taylor series, we obtain

$$\alpha \left({\rho }_{p,i}^{j+1}\right)=\frac{{\rho }_{p1}}{{\rho }_{p,i}^j}+\frac{{\rho }_{p,i}^{j+1}-{\rho }_{p,i}^j}{1}\left(-\frac{{\rho }_{p1}^{}}{{\left({\rho }_{p,i}^j\right)}^2}\right),$$

and so

$${\rho }_{p,i}^{j+1}={\rho }_{p,i}^j+\left(\frac{{\rho }_{p1}}{{\rho }_{p,i}^j}+\frac{{\rho }_{p,i}^{j+1}-{\rho }_{p,i}^j}{1}\left(-\frac{{\rho }_{p1}^{}}{{\left({\rho }_{p,i}^j\right)}^2}\right)\right)\tau {\beta }_{p0}{c}_i^j.$$

Finally,

$${\rho }_{p,i}^{j+1}=\left(\frac{{\left({\rho }_{p,i}^j\right)}^2+2{\rho }_{p1}^{}\tau {\beta }_{p0}{c}_i^j}{{\left({\rho }_{p,i}^j\right)}^2+{\rho }_{p1}^{}\tau {\beta }_{p0}{c}_i^j}\right)·{\rho }_{p,i}^j.$$

It is obvious that

$$\frac{{\left({\rho }_{p,i}^j\right)}^2+2{\rho }_{p1}^{}\tau {\beta }_{p0}{c}_i^j}{{\left({\rho }_{p,i}^j\right)}^2+{\rho }_{p1}^{}\tau {\beta }_{p0}{c}_i^j}>1.$$

meaning that scheme (27) is unstable. □

It is unexpected that, in some cases explicit difference schemes are “better” than implicit ones.

5. Numerical Algorithm

Transforming the difference schemes (12)–(14) we get

$$c_i^{j+1}=\frac{h}{v\tau +h{m}_0}\left(\frac{v\tau }{h}{c}_{i-1}^{j+1}+m_0{c}_i^j-\left({\rho }_{a,i}^{j+1}-{\rho }_{a,i}^j+{\rho }_{p,i}^{j+1}-{\rho }_{p,i}^j\right)\right),$$

(28)

$${\rho }_{a,i}^{j+1}={\rho }_{a,i}^j+\tau {\beta }_a\left(c_i^j-\frac{{\rho }_{a,i}^j}{{\rho }_{a0}}{c}_0\right),$$

(29)

$${\rho }_{p,i}^{j+1}={\rho }_{p,i}^j+\tau \alpha \left({\rho }_{p,i}^j\right){\beta }_{p0}{c}_i^j,$$

(30)

where

$$\alpha \left({\rho }_{p,i}^j\right)=\left\{\begin{array}{ccc}1\hfill & \mathrm{for}\hfill & {\rho }_{p,i}^j\le {\rho }_{p1},\hfill \\ {\rho }_{p1}/\phantom{{\rho }_{p1}{\rho }_{p,i}^j}{\rho }_{p,i}^j\hfill & \mathrm{for}\hfill & {\rho }_{p1}<{\rho }_{p,i}^j<{\rho }_{p0} ,\hfill \\ 0\hfill & \mathrm{for}\hfill & {\rho }_{p,i}^j={\rho }_{p0} .\hfill \end{array}\right.$$

The system (28)–(30) is solved under the known initial and boundary conditions (15). The calculations are carried out in the following order. According to (29) and (30), the values ${\rho }_{a,i}^{j+1}$ and ${\rho }_{p,i}^{j+1}$ are determined through the known quantities ${\rho }_{a,i}^j$, ${\rho }_{p,i}^j$ and $c_i^j$ of the lower layer at the corresponding points, from (28) we determine $c_i^{j+1}$.

6. Numerical Experiments and Their Analyses

To carry out numerical experiments we developed a program in Python.

As the initial parameters, we take the following numerical values as [23,24]: $c_0=0.05$, $m_0=0.3$, $v_0={10}^{-4}$ m/s, ${\rho }_0=0.1$, ${\rho }_{a0}=0.01$, ${\rho }_{p0}=0.09$, ${\beta }_{a0}=0.005$ s${}^{-1}$, ${\beta }_{p0}=0.005$ s${}^{-1}$.

Let us analyse the numerical results. Over time, the values of c, ${\rho }_a$ and ${\rho }_p$ at fixed points in the reservoir increase (Figure 1). With an increase in the parameter ${\rho }_{p1}$ to $0.05$ at $t=450$ s (Figure 1), near the point $x=0$, the capacitance of the passive zone is completely saturated with deposition and this leads to an increase in the concentration of suspended solids in the liquid c and deposition in the active zone ${\rho }_a$. Increasing the parameter ${\rho }_{p1}$ leads to a delay in the front of the advance c, ${\rho }_a$ and ${\rho }_p$.

Figure 2 shows graphs of the characteristics of the solute transport for various ${\rho }_{p1}$, which allow us to estimate the effect of this parameter. As can be seen from the figures, nonmonotonic changes in all characteristics occur. An increase in ${\rho }_{p1}$ leads to a progressive dynamics of $c/\phantom{c{c}_0}{c}_0$, ${\rho }_a$, ${\rho }_p$ for small x and for relatively large x the opposite is observed.

It can be seen from Figure 3 that increasing the value of the parameter ${\beta }_{p0}$ leads to increased “ageing”, that is the concentration of irreversible deposition increases near the point $x=0$. This leads to an increase in the concentration of reversible deposition and concentrations of suspended particles. However, at the points away from $x=0$ we can see the opposite scene. Most of the particles are deposited near the point $x=0$ and this does not lead to distribution of suspended particles further. Figure 4 shows the time of the beginning of ageing according to the coordinate for various values of ${\rho }_{p1}$. For ${\rho }_{p1}=$ 0.01, ageing is started at the point $x=0.02$ approximately at $t=150$ s, and for ${\rho }_{p1}=0.03$ at $t=350$ s, for ${\rho }_{p1}=0.05$ at $t=550$ s (Figure 4).

7. Improving the Model

Over the time when deposition is increasing it leads to change in porosity and permeability of porous media. Changes in porosity can be written as $m=m_0-\left({\rho }_a+{\rho }_p\right)$. It also leads to change in filtration coefficient $K\left(m\right)$. Several formulas have been obtained for filtration coefficient in terms of porosity. Of these formulas the Carman–Kozeny formula [42,43,44] is widely used, which was developed for laminar flow in a packed bed of spherical grains. In our case porous media is homogenous and filtration velocity is small enough. So, for $K\left(m\right)$, we use Carman–Kozeny equation [42,43,44]

$$K\left(m\right)={k_0{m}^3/\phantom{k_0{m}^3\left(1-m\right)}\left(1-m\right)}^2,$$

(31)

where $k_0=const$ is initial permeability coefficient of media before deposition occurs.

In constant flow regime, increasing of porosity and permeability leads to change in pressure gradient also. Darcy’s law basically used an equation describing the flow of a fluid through a porous medium. To find pressure gradient, we use Darcy’s law [26,27], which was first determined experimentally by Darcy, but it can be derived from the Navier–Stokes equations for small Reynolds numbers when inertial effects are ignored. In our case filtration velocity is small enough and we can ignore inertial effects, and so we use Darcy’s law in the form

$$v=K\left(m\right)\left|\nabla p\right|,$$

(32)

where $\left|\nabla p\right|$ is module of pressure gradient.

In [42], kinetic equations of the process of deposition (capture) of solid particles of a suspension as well as their release were generalised, taking into account dynamic factors. In accordance with [42], the process of particle deposition and release depends on the pressure gradient, and the greater module of pressure gradient the smaller probability of the deposition of the particles and greater the probability of the release of the particles from pores. In accordance with [42,45], the kinetic Equation (7) characterising the deposition and release of solid particles in the active zone of the porous medium is written as

$$\frac{\partial {\rho }_a}{\partial t}=\frac{{\beta }_a}{1+\gamma \left|\nabla p\right|}c-{\beta }_a\frac{{\rho }_a\left(1+\omega \left|\nabla p\right|\right)}{{\rho }_{a0}}{c}_0.$$

(33)

where, constant parameters $\gamma$ and $\omega$ characterise the intensity of the influence of $\left|\nabla p\right|$ on attachment and detachment of particles, respectively.

Thus, model of suspension filtration in porous media taking into account changes in porous medium and dynamical factors consists of (6), (8), (31)–(33). The system of partial differential equations is solved under initial and boundary conditions (10) by using the finite difference schemes.

Some results are presented in Figure 5, Figure 6 and Figure 7. Dynamic factors mainly affect the profiles of changes in the concentration of the deposition in active zone. In the case when dynamic factors are not taken into account over time, the capacity of active zone ${\rho }_{a0}$ is completely saturated with deposition, here near the point $x=0$ over time first ${\rho }_a$ increases, then decreases, and finally $\partial {\rho }_a/\phantom{\partial {\rho }_a\partial t}\partial t=0$, even ${\rho }_{a0}$ is not completely saturated with deposition (Figure 5b). Why the capacity of active zone is not completely saturating with deposition can be explained by the fact that the capacity of active zone is dynamical. Which means by increasing the pressure gradient the structure of porous media can be partially destroyed and some previously deposited particles can get back to flow. This conclusion is quite consistent with conclusions of Mints [35]. The dynamical equilibrium, meaning that the number of deposited particles within the unit time is equal to the detached particles in this time in the active zone, occurs before the deposited particles volume get the capacity of the active zone. Dynamical equilibrium in active zone leads to an increase in suspended particles concentration which in turn leads to an increase in concentration of deposited particles in the passive zone. For ${\rho }_{p1}=$ 0.01 from $t\approx 1700 s$ and for ${\rho }_{p1}=0.03$ from $t\approx 600 s$ the dynamical equilibrium in active zone occurs, that is $\partial {\rho }_a/\phantom{\partial {\rho }_a\partial t}\partial t=0$ (Figure 6a). We can say that an increase in the parameters $\gamma$ and $\omega$ leads to a decrease in ${\rho }_a$ at the corresponding points in the reservoir. Figure 7 shows the graphics of $\left|\nabla p\right|$ for various ${\rho }_{p1}$. For large values of ${\rho }_{p1}$ we can see the presence of large and almost constant values of $\left|\nabla p\right|$ for small x, as in the previous cases, the nature of the $\left|\nabla p\right|$ changes in greater values of x.

8. Discussion

On the basis of numerical analyses it has been established that for small values of ageing parameter, the deposition concentration in the passive zone still does not reach the capacitance of the passive zone in the studied times. With an increase in the ageing parameter near the point of inlet, the capacity of the passive zone is completely saturated with deposition and this leads to an increase in the concentration of suspended solids in fluid and the deposition in the active zone. It is shown that the consideration of dynamic factors mainly affects the profiles of changes in the concentration of the deposition in active zone. In the case when dynamic factors are not taken into account over time, the capacity of active zone is completely saturated with deposition. Here, near the inlet point over the time, deposition in active zone first increases, then decreases, and finally stops changing; even the capacity of active zone is not completely saturated with deposition. This can be explained as follows: over the time deposition is growing, it leads to a decrease of porosity and permeability. In constant, the filtration velocity regime decrease of porosity leads to an increase of pressure gradient. So, the probability of deposited particles detachment in the active zone increases and dynamical equilibrium (number of deposited and released particles in unity of time) occurs before the deposition reaches the capacity.

9. Conclusions

In this paper, the problem of filtration of a single-component suspension in a porous medium is posed. A mathematical model of suspension filtration in porous media is considered as a system of partial differential equations, which consists of mass balance equation, kinetic equations, Darcy’s law and Carman–Kozeny equation. The model takes into account the effects of ageing and dynamical factors. The influence of module of pressure gradient have been taken into account in kinetic equations directly. This allows us to estimate how the change in porous media characteristics affects particle deposition and release processes. To solve the problem, a numerical method has been developed on the basis of finite difference schemes. Though there are many finite difference schemes examined for partial differential equations, there are only a few examined for the system of partial differential equations. In addition, stability of difference schemes has been proved. As we know, implicit schemes are “better” than explicit schemes in stability, but in our case it has been shown that in some cases explicit schemes are more stable than implicit ones. To carry out numerical experiments, a program has been developed in Python.