Mathematical modeling of erosion and deposition in porous media

Erosion and deposition are represented as the evolution of solid bodies due to the forces exerted by the fluid or air on the contact surfaces, which both often lead to reconfiguration and change of the topology and structure of the porous media. These processes are notably very complicated and challenging to study. In this work, we formulate simplified and idealized mathematical models to examine the internal evolution of flow networks in the setting of cylindrical channels, undergoing a unidirectional flow, by using asymptotic and numerical techniques. Starting from the Stokes equations combined with the advection-diffusion equation for solid transport, we propose a model to construct a complete analysis of both the erosion and deposition. The considered approach is of the form of threshold laws: the fluid-solid interface erosion and deposition occur when the total shear stress is, respectively, greater or lower than some specified critical values, depending on the solid material. As a consequence of the erosion and deposition, the radii of the channels in the structure expand and shrink, respectively, due to several key parameters, which we find and investigate in this paper. We also perform a parametric study to quantify the correlation between these threshold values and the particle concentration in the flow. A comprehensive parametric study of the constructed model reveals that the final configuration of the structure can be predicted from the system parameters.

Figure 1

Schematic showing a structure and channels with the channel radii $A(X,Y,Z,T)$ and particle concentration $C(X,Y,Z,T)$.

Figure 2

Erosion and deposition system diagram.

Figure 3

Channel erosion: comparison of the channel radius $a\left(t\right)$ vs time $t$, based on our model and the dimensionless version of model proposed by Bonnelli et al. [17] given in (36) and (38), respectively, for three dimensionless erosion coefficients ${\kappa }_{\mathrm{er}}$. The following parameters are used: $W=0.5 \mu \text{m}$, $\rho =1000$ kg ${\mathrm{m}}^{-3}, {\rho }_s=1600$ kg ${\mathrm{m}}^{-3}, {\mathrm{\Sigma }}_e=0$ kg ${\mathrm{m}}^{-1} {\mathrm{s}}^{-2}, {\mathrm{\Sigma }}_d=0$ kg ${\mathrm{m}}^{-1} {\mathrm{s}}^{-2}, P_{\mathrm{in}}-P_{\mathrm{out}}=10$ kg ${\mathrm{m}}^{-1} {\mathrm{s}}^{-2}$, and $B_d=0 {\mathrm{mol}}^{-1} {\mathrm{kg}}^{-1} {\mathrm{m}}^5$ s. For all simulations, the corresponding dimensionless parameters are ${\tau }_e=0, {\tau }_d=0, {\beta }_e=1,1.5,2 [{\beta }_e={\kappa }_{\mathrm{er}}\mu /\left(2\rho W\right)], {\beta }_d=0$, and $a_{\mathrm{in}}=0.5$.

Figure 4

Channel erosion: comparison of the channel radius $a\left(t\right)$ vs time $t$, based on our model and the dimensionless version of model proposed by Bonnelli et al. [17] given in (36) and (38), respectively, for two dimensionless erosion coefficients ${\kappa }_{\mathrm{er}}$. The following parameters are used: $W=0.5 \mu \text{m}$, $\rho =1000$ kg ${\mathrm{m}}^{-3}, {\rho }_s=1600$ kg ${\mathrm{m}}^{-3}, {\mathrm{\Sigma }}_e=0$ kg ${\mathrm{m}}^{-1} {\mathrm{s}}^{-2}, {\mathrm{\Sigma }}_d=0$ kg ${\mathrm{m}}^{-1} {\mathrm{s}}^{-2}, P_{\mathrm{in}}-P_{\mathrm{out}}=10$ kg ${\mathrm{m}}^{-1} {\mathrm{s}}^{-2}$, and $B_d=0 {\mathrm{mol}}^{-1} {\mathrm{kg}}^{-1} {\mathrm{m}}^5$ s. For all simulations, the corresponding dimensionless parameters are ${\tau }_e=0, {\tau }_d=0, {\beta }_e=0.001,0.01 [{\beta }_e={\kappa }_{\mathrm{er}}\mu /\left(2\rho W\right)], {\beta }_d=0$, and $a_{\mathrm{in}}=0.5$.

Figure 5

Channel erosion: comparison of the channel radius $a\left(t\right)$ vs time $t$, based on our model and the dimensionless version of model proposed by Bonnelli et al. [17] given in (36) and (38), respectively, for three dimensionless erosion thresholds ${\mathrm{\Sigma }}_e$. The following parameters are used: $W=0.5 \mu \text{m}$, $\rho =1000$ kg ${\mathrm{m}}^{-3}, {\rho }_s=1600$ kg ${\mathrm{m}}^{-3}, {\kappa }_{\mathrm{er}}=4 {\mathrm{m}}^{-1}s, {\mathrm{\Sigma }}_d=0$ kg ${\mathrm{m}}^{-1} {\mathrm{s}}^{-2}, P_{\mathrm{in}}-P_{\mathrm{out}}=10$ kg ${\mathrm{m}}^{-1} {\mathrm{s}}^{-2}$, and $B_d=0 {\mathrm{mol}}^{-1} {\mathrm{kg}}^{-1} {\mathrm{m}}^5$ s. For all simulations, the corresponding dimensionless parameters are ${\tau }_d=0, {\tau }_e=0.0017,0.1667,0.5$, (${\beta }_e=8$), ${\beta }_d=0$, and $a_{\mathrm{in}}=0.5$.

Figure 6

Channel radius evolution $a(x,t)$ [black curves in (a)–(c)] and the particle concentration profile $c(x,t)$ [red curves in (d)–(f)], at several different times, for three different erosion shear stress thresholds ${\tau }_e$ with their corresponding erosion coefficients ${\lambda }_e$: ${\tau }_e=0.8$ and ${\lambda }_e=0.6$; ${\tau }_e=0.4$ and ${\lambda }_e=0.3$; and ${\tau }_e=0.1$ and ${\lambda }_e=0.075$. For all simulations, the following parameters are used: $\mathrm{Pe}=2, {\alpha }_e=0.75, {\alpha }_d=0.75 {\lambda }_d=0, {\tau }_d=0, {\beta }_e=1, {\beta }_d=0, c_{\mathrm{in}}=0$, and $a_{\mathrm{in}}=\frac{0.65}{0.55}\left(\frac{6.5+cos\left(20\pi x\right)}{10}\right)e^{-x}$.

Figure 7

Channel radius evolution $a(x,t)$ vs $x$, for several different erosion coefficients ${\beta }_e$. For all simulations, the following parameters are considered: $\mathrm{Pe}=2, {\alpha }_e={\alpha }_d=0.75, {\lambda }_e=0.0375, {\lambda }_d=0, {\tau }_e=0.05, {\tau }_d=0, {\beta }_d=0, c_{\mathrm{in}}=0$, and $a_{\mathrm{in}}=\frac{0.65}{0.55}\left(\frac{6.5+cos\left(20\pi x\right)}{10}\right)e^{-x}$.

Figure 8

Channel deposition: comparison of the channel radius $a\left(t\right)$ vs time $t$, based on our model and the dimensionless version of model proposed by Jäger et al. [13] given in (36) and (48), respectively, for three dimensionless deposition coefficients ${\kappa }_{\mathrm{dep}}$. The following parameters are used: $W=0.5 \mu \text{m}$, $\rho =1000$ kg ${\mathrm{m}}^{-3}, {\rho }_s=1600$ kg ${\mathrm{m}}^{-3}, {\mathrm{\Sigma }}_e=0$ kg ${\mathrm{m}}^{-1} {\mathrm{s}}^{-2}, {\mathrm{\Sigma }}_d=0.4$ kg ${\mathrm{m}}^{-1} {\mathrm{s}}^{-2}, P_{\mathrm{in}}-P_{\mathrm{out}}=10$ kg ${\mathrm{m}}^{-1} {\mathrm{s}}^{-2}$, and $B_e=0 {\mathrm{kg}}^{-1} {\mathrm{m}}^2$ s. For all simulations, the corresponding dimensionless parameters are ${\tau }_e=0, {\tau }_d=24, {\beta }_e=0, {\beta }_d=1,1.5,2 [{\beta }_d={\kappa }_{\mathrm{dep}}\mu C_{\mathrm{typ}}/\left(2\rho W\right)]$, and $a_{\mathrm{in}}=0.5$. The deposition coefficient in Jäger et al. [13], ${\kappa }_{\mathrm{dep}}$, is related to deposition coefficient in our model, $B_d$, via ${\kappa }_{\mathrm{dep}}={\rho }_s{B}_d$.

Figure 9

Channel deposition: comparison of the channel radius $a\left(t\right)$ vs time $t$, based on our model and the dimensionless version of model proposed by Jäger et al. [13] given in (36) and (48), respectively, for three deposition threshold coefficients ${\mathrm{\Sigma }}_d$. The following parameters are used: $W=0.5 \mu \text{m}$, $\rho =1000$ kg ${\mathrm{m}}^{-3}, {\rho }_s=1600$ kg ${\mathrm{m}}^{-3}, {\mathrm{\Sigma }}_e=0$ kg ${\mathrm{m}}^{-1} {\mathrm{s}}^{-2}, {\kappa }_{\mathrm{dep}}=10 {\mathrm{m}}^{-1}$ s, $P_{\mathrm{in}}-P_{\mathrm{out}}=10$ kg ${\mathrm{m}}^{-1} {\mathrm{s}}^{-2}$, and $B_e=0 {\mathrm{kg}}^{-1} {\mathrm{m}}^2$ s. For all simulations, the corresponding dimensionless parameters are ${\tau }_e=0, {\tau }_d=0.1667$, 0.3333, 0.6667, ${\beta }_e=0, {\beta }_d=2$, and $a_{\mathrm{in}}=0.5$. The deposition coefficient in Jäger et al. [13], ${\kappa }_{\mathrm{dep}}$, is related to deposition coefficient in our model, $B_d$, via ${\kappa }_{\mathrm{dep}}={\rho }_s{B}_d$.

Figure 10

Illustration of the pore radius evolution $a(x,t)$ (black curves) and its respective particles concentration profile $c(x,t)$ (red curves), at several different times, and for three different deposition shear stress thresholds ${\tau }_d$ and their corresponding stickiness coefficients ${\lambda }_d$: (a) and (d) ${\tau }_d=0.4$ and ${\lambda }_d=0.4$, (b) and (e) ${\tau }_d=0.6$ and ${\lambda }_d=0.45$, (c) and (f) ${\tau }_d=0.9$ and ${\lambda }_d=0.675$. For all simulations, the following parameters and initial channel profile are considered: ${\tau }_e=2$ (high enough to cancel erosion), $\mathrm{Pe}=2, {\alpha }_e=0.75, {\alpha }_d=0.75, {\lambda }_e=1.5, {\beta }_e=0, {\beta }_d=0.5, c_{\mathrm{in}}=0.5, a_{\mathrm{in}}=\frac{0.65}{0.55}\left(\frac{6.5+cos\left(20\pi x\right)}{10}\right)e^{-x}$.

Figure 11

Illustration of the channel volume change, $\mathrm{\Delta }V={\int }_0^1(a^2(x,t_f)-a^2(x,0))\mathrm{d}x$, in terms of a parametric study of the shear stress thresholds ${\tau }_e$ and ${\tau }_d$. The following parameters and initial channel profile are considered: $\mathrm{Pe}=2, {\alpha }_e={\alpha }_d=0.75, {\lambda }_e=0.4, {\lambda }_d=0.4, {\beta }_e=0.3, {\beta }_d=0.3, c_{\mathrm{in}}=1$, and $a_{\mathrm{in}}=0.3$.

Figure 12

(a): Illustration of the volume change, $\mathrm{\Delta }V={\int }_0^1[a^2(x,t_f)-a^2(x,0)]\mathrm{d}x$, of the channel in terms of a parametric study of the shear threshold values ${\tau }_e$ and ${\tau }_d$, with an initial concentration at the channel inlet $c_{\mathrm{in}}=0.3$. The channel pore radius evolution $a(x,t)$ [black curves in (b)–(d)] and its respective shear stress profile evolution $\tau (x,t)$ [red curves in (e)–(g)] are simulated for three pairs of shear threshold values: in (b) and (e) ${\tau }_e=0.27, {\tau }_d=0.11$, in (c) and (f) ${\tau }_e=0.59, {\tau }_d=0.21$, in (d) and (g) ${\tau }_e=1.23, {\tau }_d=0.84$. For all simulations, the following parameters are considered: $\mathrm{Pe}=2, {\alpha }_e={\alpha }_d=0.75, {\lambda }_e=0.4, {\lambda }_d=0.4, {\beta }_e=0.3, {\beta }_d=0.3$, and $a_{\mathrm{in}}=\frac{0.65}{0.55}\left(\frac{6.5+cos\left(20\pi x\right)}{10}\right)e^{-x}$.