Adsorption on Fractal Surfaces: A Non Linear Modeling Approach of a Fractional Behavior

Abstract

This article deals with the random sequential adsorption (RSA) of 2D disks of the same size on fractal surfaces with a Hausdorff dimension $1<d<2$. According to the literature and confirmed by numerical simulations in the paper, the high coverage regime exhibits fractional dynamics, i.e., dynamics in $t^{-1/d}$ where d is the fractal dimension of the surface. The main contribution this paper is that it proposes to capture this behavior with a particular class of nonlinear model: a driftless control affine model.

1. Introduction

Adsorption is a widely encountered phenomenon, especially in physics, chemistry or biology, e.g., to capture pollutants [1,2], for water purification [3], or in gas sensors [4,5]. From a mathematical point of view, the idealized adsorption process is known as random sequential adsorption (RSA) which produces fractional behaviors.

RSA has been widely studied in the literature especially in 2 dimensions. There are results on the final jamming value [6,7,8], and on the kinetic of the density of free spaces which exhibits a $t^{-1/2}$ behavior [6,9,10]. Numerous models have been proposed to capture this kinetic. An interesting analysis of these models is proposed in [11]. If $q\left(t\right)$ denotes the amount of adsorbed particles, most of the models discussed express the adsorption rate $\frac{\mathrm{d} q\left(t\right)}{\mathrm{d} t}$ as a function of the difference $q_e-q\left(t\right)$, where $q_e$ denotes the amount of particles adsorbed at equilibrium. The following models have been proposed: by Lagergren in 1898 [12]: $\frac{\mathrm{d} q\left(t\right)}{\mathrm{d} t}=k_1(q_e-q\left(t\right))$; by Kopelman in 1988 [13]: $\frac{\mathrm{d} q\left(t\right)}{\mathrm{d} t}=k_2{t}^{-h}(q_e-q\left(t\right))\mathrm{with}0\le h\le 1$; by Ho and Mckay in 2000 [14]: $\frac{\mathrm{d} q\left(t\right)}{\mathrm{d} t}=k_3{(q_e-q\left(t\right))}^2$; and by Brouers and Sotolongo-Costain in 2006 [15]: $\frac{{\mathrm{d}}^{\alpha }q\left(t\right)}{\mathrm{d} t^{\alpha }}=k_4{(q_e-q\left(t\right))}^n.$

If $\theta \left(t\right)=\frac{q\left(t\right)}{q_m}$ denotes the coverage density, where $q_m$ is the maximum amount of adsorbed particles and if c denotes the concentration of particles close to the surface (bulk), the following model was also proposed by Haerifar and Azizian in 2012 [16] and Bashiri and Shajari in 2014 [11]: $\frac{\mathrm{d} \theta \left(t\right)}{\mathrm{d} t}=k_5{t}^{-h}c(1-\theta \left(t\right))\mathrm{with}0\le h\le 1.$

However, the RSA process is non-linear. If the flow of particles is doubled, the filling dynamic of the surface is almost doubled, but the final value is the same as for a single flow. Moreover, if the flow is stopped then the filling is obviously stopped as well. Furthermore, if the flow restarts, the filling will restart at the same point. Thus, none of the time varying and fractional models previously cited permit a consistent modeling of the RSA phenomenon. That is why a nonlinear model is proposed by the authors of [17].

In practice, the surface on which the adsorption occurs is not absolutely flat, which leads some authors to consider RSA on fractal surfaces [18]. Analysis of RSA on fractal surfaces has produced a major result: the high coverage regime exhibits a fractional behavior depending on the surface dimension d, i.e., dynamics in $t^{-1/d}$ [6,9,18]. Such a property is in accordance with the analysis proposed in [19,20], which leads the author of [20] to claim that a kinetic in $t^{\nu }$ on space of dimension 1 produces a kinetic in $t^{\nu /d}$ on space of dimension d. This is precisely what we observe with the RSA process whose kinetic of free places density for the high coverage regime exhibits a $t^{-1}$ behavior in 1D [21]. Physical examples of adsorption on fractal surfaces are also considered in numerous studies [22,23,24,25,26].

In line with reference [17], the main contribution of the present paper is that it proposes to capture the RSA kinetics on fractal surfaces with driftless control affine nonlinear models. An example of a driftless control-affine system is the nonholonomic integrator introduced in [27]. This class of system some times referred to as Brockett’s system, or the Heisenberg system, as it also arises in quantum mechanics [28]. Analysis and control results for this class of systems can be found in [29,30]. Two models are proposed: a simple one which is proved analytically to produce a fractional dynamic behavior and a more complex and accurate one. The accuracy of the latter model was assessed via numerical simulations on three types of fractal surfaces.

2. Evidence of the Fractional Asymptotic Behavior of Some Fractal Surfaces

This article deals with the random sequential adsorption (RSA) of 2D disks on a fractal surface. A definition of RSA is now given.

Definition 1 (RSA).

Random sequential adsorption is a process in which particles are constantly trying to attach themselves to a random location on a surface. If the incoming particle does not overlap any previously attached particles, then it binds irreversibly.

In order to model the phenomenon of physical adsorption, the random sequential adsorption model was used. For any surface, it is suggested in the literature [6,18] that the long time coverage density of the RSA follows a power law:

$${\theta }_{\infty }-\theta \left(t\right)\sim t^{-1/d},$$

(1)

where d is the Hausdorff dimension of the surface on which the 2D disks arise. The fractal surfaces considered are a set of points contained in the square $[0,1]\times [0,1]$.

In order to simulate the RSA process of 2D disks on a given fractal surface, the RSA algorithm used is the same as in [17] and is recalled in Algorithm 1:

Algorithm 1 Random Sequential Adsorption

A random point $c=(p,q)$ of the fractal is selected at each iteration of the process. A disc of radius R and center c will fix on the surface if: A part of the disc does not lie outside the surface $[0,1]\times [0,1]$ containing the fractal, which is true if the following two conditions are satisfied: $p,q\ge R$ and $p,q\le 1-R$; There is no overlap between the current disc and a previously fixed disc, that is if $\mathrm{d} (c,c_k)\ge 2R$, where for all $k\in \mathbb{N}$, $c_k=(p_k,q_k)$ is the center of a disc previously fixed on S and where $\mathrm{d}(c,c_k)=\sqrt{{(p-p_k)}^2+{(q-q_k)}^2}$ is the distance between discs.

The fractal surfaces considered are illustrated in Figure 1.

For each of the fractal surfaces previously cited, a RSA process was performed. As shown in [18], the asymptotic value of the density depends on the number of iterations in the fractal building process. The number of iterations needed to have relevant results seems to be greater than 6. The final values given in

Table 1. Final values of the density of adsorbed disks.

Fractal	Final Value
Vicsek	∼0.68
Sierpinski triangle	∼0.62
Sierpinski carpet	∼0.58

are those proposed in [18] and are close to those plotted in Figure 2. The differences can be explained by differences in the simulation conditions (the number of iterations in the fractal building process, the number of trials, the ratio between the size of the disks and the size of the fractal).

The densities of adsorbed disks as a function of time (number of trials) for each of the considered fractals are shown in Figure 2.

Figure 3 shows ${\theta }_{\infty }-\theta$ as a function of $t^{-1/d}$, where d is the dimension of the fractal. This figure provides evidence that the density of adsorbed disks has a power-law asymptotic behavior of the form

$${\theta }_{\infty }-\theta \left(t\right)\sim K{t}^{-1/d},$$

(2)

where d is the fractal surface dimension (Hausdorff dimension). The Hausdorff dimension of the Vicsek fractal is $d_V=\frac{log\left(5\right)}{log\left(3\right)}\approx 1.4649$, $d_{ST}=\frac{log\left(3\right)}{log\left(2\right)}\approx 1.585$ for the Sierpinski triangle, $d_{SC}=\frac{log\left(8\right)}{log\left(3\right)}\approx 1.8928$ for the Sierpinski carpet. The straight line on each subfigure of Figure 3 is the plot of $K{t}^{-1/d}$ to highlight the fractional behavior.

Figure 3 confirms that the results of the simulations are consistent with the expected power-law behaviors.

3. Power-Law Non Linear Dynamical Modeling

3.1. Detailed Modeling Approach on the Vicsek Fractal

In this section, the class of model considered is called a driftless control affine nonlinear model [27,28] and the first one is defined by the relation

$$\begin{array}{c}\hfill \dot{x}\left(t\right)=A{(x\left(t\right)+C)}^{1-\frac{1}{\alpha }}·u\left(t\right).\end{array}$$

(3)

$u\left(t\right)$ can be viewed as an input and x(t) is the model state and output. $A\in \mathbb{R}$, $C\in \mathbb{R}$ and $-1<\alpha <1$ ($\alpha \in {\mathbb{R}}^{*}$). This model produces trajectories of the form

$$x\left(t\right)=\frac{1}{\alpha }{\left({\int }_1^{t}Au\left(t\right)\mathrm{d} t+K\right)}^{\alpha }-C,$$

(4)

where $K\in \mathbb{R}$ is an integration constant. For $u\left(t\right)=1$, $A=1$, $K=0$ and $C=0$, Equation (4) becomes

$$x\left(t\right)=\frac{1}{\alpha }{t}^{\alpha },$$

(5)

and the corresponding trajectories are given in Figure 4 for various values of $\alpha$.

The trajectories in Figure 4 have shapes very similar to that of the densities variations of fixed particles considered above in Figure 2, i.e., with very fast growths for short times, followed by very slow progressions towards a steady state. Figure 5 shows the response of model (3) for various values of $\alpha$ with $A=-1$, $u\left(t\right)=1$, $K=0$ and $C=0$. These responses look like the one of a fractional integrator impulse response.

This model was thus used in a first approach to model the RSA density on the fractals considered in the previous section. An optimization was performed on parameters A and C. Since the power-law behavior only holds for the high coverage regime (long times) of the RSA process, two optimization criteria were used:

$$\begin{array}{cc}\hfill \mathrm{criterion}1:& {\displaystyle \epsilon =\sum _{t=0}^{t_{max}}{(\theta \left(t\right)-{\theta }_{model}\left(t\right))}^2}\end{array}$$

(6)

$$\begin{array}{cc}\hfill \mathrm{criterion}2:& {\displaystyle {\epsilon }_{TL}=\sum _{t=t_1}^{t_{max}}{(\theta \left(t\right)-{\theta }_{model}\left(t\right))}^2.}\end{array}$$

(7)

In (6) and (7), $t_{max}$ denotes the last time point in the considered RSA data and $t_1$ denotes the first time value considered for the computation of the criterion.

For the Vicsek fractal, parameter $\alpha$ was imposed equal to $-\frac{1}{d_V}$, and the result of this optimization with criterion $\epsilon$ is given in Figure 6. The following parameters were obtained: $A=4.5743\times {10}^{-4}$ and $C=-0.7137$.

The same optimization was also done using criterion ${\epsilon }_{TL}$ and $t_1=1\times {10}^5$. A comparison of the RSA density data and of the models obtained with criterion $\epsilon$ and ${\epsilon }_{LT}$ is shown in Figure 7.

Even if the model parameters were computed using criterion ${\epsilon }_{LT}$, criterion $\epsilon$ was computed for the resulting model, and conversely, criterion ${\epsilon }_{LT}$ was computed with the model obtained through the minimization of criterion $\epsilon$. The results are reported in

Table 2. Error criterion for the three modeling approach.

	Model (3)	Model (3)	Model (8)
	with $\mathit{\epsilon }$	with ${\mathit{\epsilon }}_{\mathit{LT}}$	with $\mathit{\epsilon }$
${\epsilon }_{LT}$	$2.0669$	$0.3012$	$1.6765$
$\epsilon$	$5.5433$	$11.9854$	$2.2896$

. This table shows that ${\epsilon }_{LT}$ is lower when the parameters are computed with ${\epsilon }_{LT}$. Model (3) is thus a very good candidate to model fractional dynamics (high coverage times) and achieve a compromise when computed with criterion $\epsilon$.

In order to reduce this error, the more general model

$$\dot{x}\left(t\right)=f\left(x\left(t\right)\right)·u\left(t\right),$$

(8)

is now considered, where f is a polynomial function:

$$f\left(x\right)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3+a_4{x}^4.$$

(9)

A theoretical justification of this class of model appears in [31]. The proposed model for RSA kinetics of anisotropic particles is based on the available surface function (ASF) concept and is defined by

$$\frac{\mathrm{d} \theta \left(t\right)}{\mathrm{d} t}=\frac{1}{2\pi }\int \Phi (\theta \left(t\right),\Omega )d\Omega$$

(10)

where $\Phi \left(\theta \right(t),\Omega )$ is the probability of adding a new particle with orientation $\Omega$ to the surface when the coverage is $\theta \left(t\right)$. The function $\Phi \left(\theta \right(t),\Omega )$ cannot be obtained exactly, but approximations can be computed under the form of series expansions for low and high coverage regimes. These two approximations are then combined to provide an approximate description of the kinetics over the entire coverage range. The following two interpolation formulas are proposed in [31]:

$$\Phi \left(\zeta \right)={(1-\zeta )}^4(1+c_1\zeta +c_2{\zeta }^2)\mathrm{with}\zeta =\frac{\theta \left(t\right)}{{\theta }_{\infty }}$$

(11)

and

$$\Phi \left(\zeta \right)=\frac{{(1-\zeta )}^4}{(1+c_1\zeta +c_2{\zeta }^2)}.$$

(12)

Parameters $c_1, c_2, d_1$ and $d_2$ are then computed to fit the series expansions of the function $\Phi \left(\theta \right(t),\Omega )$ in (10). An expansion similar to (11) is used in [32,33]:

$$\Phi \left(\zeta \right)={(1-\zeta )}^4(1+c_1\zeta +c_2{\zeta }^2+c_3{\zeta }^3).$$

(13)

In [32], parameters $c_1$, $c_2$, and $c_3$ are computed as in [31] and to fit insulin adsorption data as [33]. These works, and particularly [31], fully justify the interest of nonlinear models for RSA kinetic modeling and by extension for adsorption kinetics in physical time. However, these models have a constrained form to meet the asymptotic behaviors for low and high surface coverage, which reduces the accuracy of the model outside these coverage regimes.

In our approach a more general expansion is considered in the form of model (8) and (9). If $u\left(t\right)=1$, and according to relation (8), this polynomial evaluated in $\theta \left(t\right)$ is equal to the derivative of $\theta \left(t\right)$. An optimization was performed on the parameters $a_i$ in order to minimize the quadratic error $\epsilon$.

The function $\dot{\theta }\left(t\right)$ is computed numerically from data of Figure 2 and the optimized function f are plotted in Figure 8 as functions of $\theta \left(t\right)$.

The diagram in Figure 9 is an implementation of relation (8) where function f is given by relation (9). Such a diagram is used for simulation of the density of adsorbed disks as a function of time.

A comparison of RSA density data and the response of model (8) is proposed in Figure 10.

The parameters $a_0$, $a_1$, $a_2$, $a_3$, and $a_4$ given in

Table 3. Parameters of f in relation (9).

${\mathit{a}}_0$	${\mathit{a}}_1$	${\mathit{a}}_2$	${\mathit{a}}_3$	${\mathit{a}}_4$
$8.2105\times {10}^{-5}$	$-2.4909\times {10}^{-4}$	$9.6858\times {10}^{-5}$	$2.9337\times {10}^{-4}$	$-2.3250\times {10}^{-4}$

were obtained using the Matlab function fmincon.

Both on the long time and on all the data, the errors are lower than the first approach with $\epsilon$. Even if the error on the long time is greater than when the optimization done only on the long time (see Table 2), the compromise is better. Considering this, model (8) is better than model (3).

3.2. Result for the Other Fractals

For the Sierpinski triangle, the result of the modeling with model (8) is shown in Figure 11.

The parameters $a_0$, $a_1$, $a_2$, $a_3$, and $a_4$ given in

Table 4. Parameters of f in relation (9) for the modeling of the Sierpinski triangle.

${\mathit{a}}_0$	${\mathit{a}}_1$	${\mathit{a}}_2$	${\mathit{a}}_3$	${\mathit{a}}_4$
$4.8939\times {10}^{-4}$	$-1.7073\times {10}^{-3}$	$3.0119\times {10}^{-4}$	$4.1781\times {10}^{-3}$	$-3.6734\times {10}^{-3}$

were obtained.

The error with criterion $\epsilon$ (relation (6)) is

$$\epsilon =2.7037.$$

(14)

For the Sierpinski carpet, the result of the modeling with model (8) is shown in Figure 12.

The parameters $a_0$, $a_1$, $a_2$, $a_3$, and $a_4$ given in

Table 5. Parameters of f in relation (9) for the modeling of the Sierpinski carpet.

${\mathit{a}}_0$	${\mathit{a}}_1$	${\mathit{a}}_2$	${\mathit{a}}_3$	${\mathit{a}}_4$
$1.3388\times {10}^{-4}$	$-5.6145\times {10}^{-4}$	$4.1346\times {10}^{-4}$	$8.5644\times {10}^{-4}$	$-1.0131\times {10}^{-3}$

were obtained.

The error with criterion $\epsilon$ (relation (6)) is

$$\epsilon =1.2337.$$

(15)

4. Conclusions

Random sequential adsorption on fractal surfaces is considered in this paper. Numerical simulations confirm that the studied phenomenon exhibits a fractional dynamic linked to the dimension of the fractal surfaces considered. In the literature, fractional behaviors are almost exclusively modeled using fractional models. However, this implicit link between fractional behaviors and fractional models does not result from physical justification. Moreover, as fractional models are doubly infinite dimensional models and thus of infinite memory [34], several limitations are associated with them [35,36]. However, many other classes of model permit to capture fractional behaviors [37,38,39,40]. Here, two driftless control-affine nonlinear models are proposed to capture this dynamic. For the first one, its ability to produce a fractional behavior is proved analytically. It is accurate on the high coverage regime but not on the beginning of the process. For this reason a more complex model but with a limited number of parameters is proposed. Numerical simulations on several fractal surfaces confirmed its efficiency.

The authors now intend to apply this modeling approach to real adsorption phenomena such as those encountered in Love wave sensors [4], nitrogen adsorption [41] or oxygen adsorption [24].