**1. Introduction**

Based on the experimental and numerical simulation results, available experience indicates the importance of considering the Basset force in studying the processes of deposition and sedimentation of small particles moving close to rigid boundaries [1]. In this regard, differential formulation of the Basset force for its numerical calculation is an urgent problem in the field of computational mechanics [2]. Recently, the solution of this problem has been associated with the fractional calculus application [3].

The fractional-order differential equations allows us to generalize existing approaches in the fields of mechanical and chemical engineering, particularly for modeling hydromechanical processes, such as the sedimentation of particles in a viscous fluid, deposition of aerosols in separation channels, pneumatic classification of fine particles, nutrient release from mineral fertilizers and migration of mineral components in soil, and gas-cleaning.

However, despite the existing numerical approaches for solving fractional-order differential equations, their analytical solutions' approaches remain an incompletely studied problem. Therefore, this article aims to develop analytical techniques for solving the fractional-order differential equation of particle sedimentation considering the fractional origin of the Basset force.

To achieve this goal, the following objectives were set: substantiation of fractional origin for the Basset force, obtaining the fractional-order equation of the sedimentation of particles, solving the obtained equation analytically, validating the obtained general solution analytically, and approximating the numerically obtained case studies by obtained analytical dependencies.

**Citation:** Pavlenko, I.; Ochowiak, M.; Agarwal, P.; Olszewski, R.; Michałek, B.; Krupi ´nska, A. Improvement of Mathematical Model for Sedimentation Process. *Energies* **2021**, *14*, 4561. https://doi.org/10.3390/ en14154561

Academic Editor: Pouyan Talebizadeh Sardari

Received: 22 June 2021 Accepted: 23 July 2021 Published: 28 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The problems of particle deposition and sedimentation in a heterogeneous dispersed system have not been thoroughly investigated. The first works devoted to solving small particle motion equations in a viscous fluid and a nonuniform flow were presented by Leal [4], Maxey, and Riley [5], respectively. Mainly, weak inertia and non-Newtonian effects on the dynamics of rigid particles in an unbounded fluid were considered. However, the proposed analytical approaches were not considered the fractional origin of the Basset force. In addition, Loussaief et al. [6] investigated a spherical particle's motion in a viscous fluid and applied the Thomas algorithm to determine the characteristics of motion for a spherical particle along a slip wall. However, the proposed methodology does not allow one to obtain the general solution of the particle sedimentation equation analytically.

Coimbra and Kobayashi [7] studied the small particle motion in a viscous medium in the rotating cylinder considering the Saffman fractional-order lift force. However, this force acts on particles in vortex flows. Moreover, the authors only presented a comparison with the results obtained by considering the drag force as a dominant one. Oppenheimer et al. [8] studied the coupled thermal and hydromechanical particle motion problem in a viscous fluid at low Reynolds numbers. As a result, a general analytical expression to determine the force and its torque on a particle was derived based on the Lorentz reciprocal theorem. However, this presented analytical solution does not consider the fractional-order Basset force.

Moreno-Casas and Bombardelli [9] presented a general numerical approach for the calculation of the Basset force. A method for approximating the Basset force was proposed by van Hinsberg et al. [10]. However, such an approach does not allow one to estimate the sedimentation velocity analytically and, therefore, the sedimentation time accurately.

A number of recent research works have aimed at applying fractional calculus in the field of engineering. Particularly, Chung [11] described a general approach for solving fractional Newton's mechanics problems based on fractional-order differential equations. As a result, an approach to solve fractional-order differential equations of Newton dynamics approximately has been proposed based on infinite power series. However, such an approach cannot be applied to obtain the analytical solution of the fractional-order particle sedimentation equation rigorously.

Agila [12] proposed applying fractional Euler–Lagrange equations for solving the problem of free damped oscillations. He et al. [13] analyzed the dynamic response of viscosimeters based on fractional-order differential equations. As a result, the memory-free Yuan–Agrawal's approach for numerical integrating fractional-order differential equations was developed and proven for the particular case studies. The presented approach discovered perspectives in measuring the viscosity of fluids.

Tomovski and Sandev [14] proposed an analytical treatment of the wave equation considering fractional friction and obtained the solution based on the Mittag–Leffler-type functions. Rossikhin and Shitikova [15] studied the dynamic behavior of nonlinear oscillatory systems described by fractional-order time derivatives. The proposed approach allowed them to obtain approximations of particular oscillatory modes for the fixed equilibrium position case study.

Therefore, the following research gaps in particle deposition and sedimentation in a heterogeneous dispersed system should be stated due to the critical review mentioned above. First, the fractional-order particles sedimentation equation considering the Basset force's fractional origin should be substantiated and solved analytically. Second, the general solution of the proposed equation should be found numerically. Finally, the obtained solution should be validated analytically for the available case study and proven numerically.

#### **2. Materials and Methods**

*2.1. The Particle Motion Equation Considering the Basset Force*

According to the equation of motion for a small rigid sphere in a nonuniform flow [16], the particle motion equation moving with a time-varying velocity in projection to the positive tangential direction has the following form:

$$m\_p \frac{dv\_p(t)}{dt} = F\_\% - F\_A - F\_d - F\_f - F\_{B\prime} \tag{1}$$

where *mp*—mass of a particle (kg), *vp*—particle velocity (m/s), *t*—time (s), *Fg* = *mpg* gravity force (N), *g*—acceleration gravity (m/s), *FA* = *mfg*–Archimedes' force, (N), *mf* added mass of a flow (kg), *Fd* = 6*πμavp*(*t*)—drag force (N), *μ*—dynamic viscosity (Pa·s), *a*—radius of a particle (m), *Ff* = <sup>1</sup> <sup>2</sup>*mf dvp*(*t*) *dt* —force of the flow added mass (N), and *FB*— Basset force, determined as follows [17]:

$$F\_B = 6\pi\mu a^2 \int\_0^t \frac{\frac{dv\_p(\tau)}{d\tau}}{\sqrt{\pi\nu(t-\tau)}} d\tau,\tag{2}$$

where *ν*—kinematic viscosity (m2/s) and *τ*—time parameter (s).

The expressions mentioned above allow us to rewrite Equation (1) as follows:

$$\left(m\_p + \frac{1}{2}m\_f\right)\frac{dv\_p(t)}{dt} + 6\sqrt{\pi\rho\mu}a^2 \int\_0^t \frac{dv\_p(\tau)}{\sqrt{t-\tau}}d\tau + 6\pi\mu a v\_p(t) = \left(m\_p - m\_f\right)\mathbf{g}\_{,}\tag{3}$$

where *ρ*—fluid density (kg/m3).

The introduction of the parameters

$$n = \frac{3\pi\sqrt{\rho\mu}a^2}{m\_p + \frac{1}{2}m\_f};\ n = \frac{6\pi\mu a}{m\_p + \frac{1}{2}m\_f};\ v\_{p\infty} = \frac{\left(m\_p - m\_f\right)g}{6\pi\mu a} \tag{4}$$

allows us to rewrite integro-differential Equation (3) in the following form:

$$\frac{dv\_p(t)}{dt} + \frac{2\alpha}{\sqrt{\pi}} \int\_0^t \frac{\frac{dv\_p(\tau)}{d\tau}}{\sqrt{t - \tau}} d\tau + nv\_p(t) = nv\_{p\infty\prime} \tag{5}$$

where *<sup>n</sup>*—relaxation factor (s–1), *<sup>α</sup>*—coefficient of the Basset force (s–1/2), and *vp*<sup>∞</sup> <sup>=</sup> lim*t*→∞*v*(*t*) stationary velocity (m/s).

#### *2.2. The Fractional-Order Differential Equation of Particles Sedimentation*

For further analytically solving Equation (5), the Riemann–Liouville integral and derivative [18] and Grünwald–Letnikov fractional derivative are considered [19]:

$$\begin{split} I\_{t\_0, t}^{\beta} v\_p(t) &= \frac{1}{\Gamma(\beta)} \int\_{t\_0}^t v\_p(\tau) (t - \tau)^{\beta - 1} d\tau; \\ D\_{t\_0, t}^{\beta} v\_p(t) &= \frac{1}{\Gamma(|\beta| - \beta)} \frac{d^{|\beta|}}{dt^{|\beta|}} \int\_{t\_0}^t v\_p(\tau) (t - \tau)^{\beta - [\beta] + 1} d\tau; \\ D\_{t\_0, t}^{\beta} v\_p(t) &= \lim\_{N \to \infty} \left( \frac{t - t\_0}{N} \right)^{\beta} \sum\_{i = 0}^N \frac{\Gamma(a + i)}{i! \Gamma(a)} v\_p \left[ t - \frac{i}{N} (t - t\_0) \right] = I\_{t\_0, t}^{-\beta} v\_p(t), \end{split} \tag{6}$$

where *β*—fractional-order, *t*0—initial time (s), and Γ(*β*)—gamma function [20]:

$$\Gamma(\beta) = \int\_0^\infty e^{-x} x^{\beta - 1} dx, \text{ } \text{Re}(\beta) > 0 \tag{7}$$

Considering the connection between fractional derivatives [21]:

$$D\_{t\_0, t}^{\gamma} \left[ \frac{d^r}{dt^r} v\_p(t) \right] = \frac{d^r}{dt^r} \left[ D\_{t\_0, t}^{\gamma} v\_p(t) \right] - \frac{(t - t\_0)^{-\gamma - r + i}}{\Gamma(-\gamma - r + i + 1)} \left[ \frac{d^i}{dt^i} v\_p(t) \right]\_{t = t\_0},\tag{8}$$

in the case of *γ* = <sup>1</sup> <sup>2</sup> , *r* = 1 for initial time *t*<sup>0</sup> = 0, and considering Γ 1 2 <sup>=</sup> <sup>√</sup>*π*, the following equation can be written:

$$I\_{0,t}^{\frac{1}{2}}\left[\frac{d\upsilon\_p(t)}{dt}\right] = D\_{0,t}^{\frac{1}{2}}\left[\upsilon\_p(t)\right] - \frac{\upsilon\_p(0)}{\sqrt{\pi t}}.\tag{9}$$

Comparing Equations (2) and (9) with the Riemann–Liouville integro-differential (6) for the value of *β* = <sup>1</sup> <sup>2</sup> and initial time *t*<sup>0</sup> = 0, as well as for the case of zero initial condition (*vp*(0) = 0), allows us to determine the Basset force in terms of fractional calculus:

$$F\_B = 6\pi a^2 \sqrt{\mu \rho} I\_{0,t}^{\frac{1}{2}} \left[ \frac{dv\_p(t)}{dt} \right] = 6\pi a^2 \sqrt{\mu \rho} D\_{0,t}^{\frac{1}{2}} \left[ v\_p(t) \right]. \tag{10}$$

Therefore, Equation (5) takes the following form:

$$\frac{d\upsilon\_p(t)}{dt} + 2\alpha \frac{d^{\frac{1}{2}}\upsilon\_p(t)}{dt^{\frac{1}{2}}} + n\upsilon\_p(t) = n\upsilon\_{p\infty}.\tag{11}$$

Notably, this formula is the fractional-order differential equation of particle sedimentation that considers the Basset force. This equation can be applied for more precise modeling of the processes of fine particles sedimentation in liquids, aerosol deposition in a gas flow, and particle deposition in gas-dispersed systems.

#### **3. Results**

*3.1. The General Solution of the Fractional-Order Differential Equation of Particles Sedimentation*

For solving Equation (11) analytically, the Laplace transform [22] is used for zero initial condition (*vp*(0) = 0). In this case, the following equational operation can be obtained:

$$\left(\mathrm{s} + 2\mathrm{n}\sqrt{s} + n\right)V\_p(\mathrm{s}) = \frac{n\upsilon\_{p\infty}}{\mathrm{s}},\tag{12}$$

where *s*—complex frequency parameter (s–1) and *Vp*(*s*)—Laplace transform of the particle velocity (m):

$$V\_{\mathcal{P}}(\mathbf{s}) = \int\_0^{\infty} \upsilon\_{\mathcal{P}}(t) e^{-st} dt,\tag{13}$$

which can be determined from the algebraic Equation (12):

$$V\_{\mathbb{P}}(\mathbf{s}) = \frac{n\upsilon\_{p\infty}}{\mathrm{s}\left(\mathbf{s} + 2\alpha\sqrt{\mathbf{s}} + n\right)} = n\upsilon\_{p\infty} \left(\frac{a\_1}{\sqrt{\mathbf{s}}} + \frac{a\_2}{\mathbf{s}} + \frac{a\_3}{\sqrt{\mathbf{s}} - \theta\_1} + \frac{a\_4}{\sqrt{\mathbf{s}} - \theta\_2}\right),\tag{14}$$

where *<sup>θ</sup>*1,2 <sup>=</sup> <sup>−</sup>*<sup>α</sup>* <sup>±</sup> <sup>√</sup>*α*<sup>2</sup> <sup>−</sup> *<sup>n</sup>*—a couple of roots of the square equation *<sup>θ</sup>*<sup>2</sup> + 2*αθ* <sup>+</sup> *<sup>n</sup>* = 0.

Notably, *<sup>θ</sup>* <sup>=</sup> <sup>√</sup>*s*. The unknown parameters *<sup>a</sup>*1, *<sup>a</sup>*2, *<sup>a</sup>*3, and *<sup>a</sup>*<sup>4</sup> are obtained from the following matrix equation:

$$
\begin{bmatrix}
1 & 0 & 1 & 1 \\
\theta\_1 \theta\_2 & -(\theta\_1 + \theta\_2) & 0 & 0 \\
0 & \theta\_1 \theta\_2 & 0 & 0
\end{bmatrix}
\begin{Bmatrix}
a\_1 \\
a\_2 \\
a\_3 \\
a\_4
\end{Bmatrix} = \begin{Bmatrix}
0 \\ 0 \\ 0 \\ 1
\end{Bmatrix}.\tag{15}
$$

Based on the inverse matrix method, the following dependencies can be obtained after identical transformation:

$$a\_1 = -\frac{2\alpha}{n^2};\ a\_2 = \frac{1}{n};\ a\_{3,4} = \pm \frac{1}{2\left(\alpha \mp \sqrt{\alpha^2 - n}\right)}.\tag{16}$$

Introduction of the parameters *θ* = –*θ*<sup>1</sup> and *θ*/ = –*θ*2, after considering the equality *L*−<sup>1</sup> √1 *s* = <sup>√</sup><sup>1</sup> *<sup>π</sup><sup>t</sup>* allows us to apply inverse Laplace transform [23] to Equation (14):

$$v\_p(t) = nv\_{p\infty} \left[ \frac{a\_1}{\sqrt{\pi t}} + a\_2 H(t) + a\_3 L^{-1} \left( \frac{1}{\sqrt{s} + \theta} \right) + a\_4 L^{-1} \left( \frac{1}{\sqrt{s} + \theta'} \right) \right],\tag{17}$$

where *H*(*t*)—Heaviside step function.

Notably, the absence of singularity for this solution is proven below. For further consideration, Mittag–Leffler function [24] is used:

$$E\_{\gamma,\theta}(z) = \sum\_{i=0}^{\infty} \frac{z^i}{\Gamma(\gamma i + \beta)}. \tag{18}$$

It has the following peculiarity in terms of its Laplace transform [25]:

$$L\left[t^{\theta-1}E\_{\gamma,\theta}(\theta t^{\gamma})\right] = \frac{s^{\gamma-\theta}}{s^{\gamma}-\theta}.\tag{19}$$

Consequently, the following inverse Laplace transform for the case of *γ* = <sup>1</sup> <sup>2</sup> and *β* = 1 can be written as follows:

$$L\left[E\_{\frac{1}{2},1}\left(\theta\sqrt{t}\right)\right] = \frac{1}{\sqrt{s}\left(\sqrt{s}-\theta\right)} = \frac{1}{\theta}\left(\frac{1}{\sqrt{s}-\theta} - \frac{1}{\sqrt{s}}\right). \tag{20}$$

Therefore, the following inverse Laplace transform can be obtained:

$$L^{-1}\left(\frac{1}{\sqrt{s}-\theta}\right) = \frac{1}{\sqrt{\pi t}} + \theta E\_{\frac{1}{\pi},1}^{\frac{1}{\pi}} \left(\theta \sqrt{t}\right). \tag{21}$$

In addition, the decomposition

$$\frac{1}{s - \theta^2} = \frac{1}{2\theta} \left( \frac{1}{\sqrt{s} - \theta} - \frac{1}{\sqrt{s} + \theta} \right) \tag{22}$$

with equality *L*−<sup>1</sup> 1 *<sup>s</sup>*−*θ*<sup>2</sup> = *eθ*2*<sup>t</sup>* with Equation (21) allows us to obtain the following inverse Laplace transform:

$$L^{-1}\left(\frac{1}{\sqrt{s}+\theta}\right) = \frac{1}{\sqrt{\pi t}} + \theta E\_{\frac{1}{2},1}\left(\theta\sqrt{t}\right) - 2\theta e^{\theta^2 t}.\tag{23}$$

Therefore, Equation (16) can be rewritten as follows:

$$w\_p(t) = n\upsilon\_{p\infty} \left[ \frac{a\_1 + a\_3 + a\_4}{\sqrt{\pi t}} + a\_2 H(t) + a\_3 \theta E\_{\frac{1}{2},1} \left( \theta \sqrt{t} \right) + a\_4 \theta' E\_{\frac{1}{2},1} \left( \theta' \sqrt{t} \right) - 2 \left( a\_3 \theta e^{\theta^2 t} + a\_4 \theta' e^{\theta^2 t} \right) \right],\tag{24}$$

Notably, *a*<sup>1</sup> + *a*<sup>3</sup> + *a*<sup>4</sup> ≡ 0 due to Equation (16), and *H*(*t*) = 1 for *t* > 0. Therefore, Equation (24) can be rewritten in the following form:

$$\upsilon\_{\mathcal{P}}(t) = \upsilon\_{\mathcal{P}^{\otimes}} \left\{ 1 - n \left[ 2 \left( a\_3 \theta e^{\theta^2 t} + a\_4 \theta' e^{\theta'^2 t} \right) - a\_3 \theta E\_{\frac{1}{2},1} \left( \theta \sqrt{t} \right) - a\_4 \theta' E\_{\frac{1}{2},1} \left( \theta' \sqrt{t} \right) \right] \right\}. \tag{25}$$

Remarkably, in practical applications, the inequality *α*<sup>2</sup> < *n*, the obtained solution, seems complex because of the complexity of parameters *θ*1,2. Consequently, the belonging

of this solution to the real range of values must be proven strictly theoretically. For this purpose, the complex parameters *θ*1,2 = −*α* ± *j* <sup>√</sup>*<sup>n</sup>* <sup>−</sup> *<sup>α</sup>*<sup>2</sup> are considered (*j*–imaginary unit: *j* <sup>2</sup> = –1). The substitution of these parameters to Equation (25), as well as the use of Newton's binomial theorem [26] and properties of hyperbolic functions [27], allows us to obtain the following expressions:

$$\begin{split} a\_{3}\theta e^{\theta^{2}t} + a\_{4}\theta^{\prime}e^{\theta^{2}t} &= \frac{1}{n} \sum\_{i=0}^{\infty} \frac{i^{\frac{1}{2}}}{\Gamma\left(\frac{1}{2}+1\right)} \sum\_{r=0}^{i} \binom{i}{r} a^{i-r} \left(n-a^{2}\right)^{\frac{r}{2}} \left(\cos\frac{\pi i}{2} - \frac{a}{\sqrt{n-a^{2}}} \sin\frac{\pi i}{2}\right); \\ a\_{3}\theta E\_{\frac{1}{2},1}\left(\theta\sqrt{t}\right) + a\_{4}\theta^{\prime}E\_{\frac{1}{2},1}\left(\theta^{\prime}\sqrt{t}\right) &= \frac{1}{n} e^{-(n-2a^{2})t} \left(\cos 2a\sqrt{n-a^{2}}t - \frac{a}{\sqrt{n-a^{2}}} \sin 2a\sqrt{n-a^{2}}t\right). \end{split} \tag{26}$$

Therefore, the Velocity (25) belongs to the real range of values.

Finally, the general Equation (25) of the fractional-order differential Equation (11) of particles sedimentation that considers the Basset force takes the following form:

$$\begin{split} \upsilon\_{P}(t) &= \upsilon\_{\text{pro}} \left[ 1 - 2e^{-(n-2a^2)t} \left( \cos 2a\sqrt{n-a^2}t - \frac{a}{\sqrt{n-a^2}} \sin 2a\sqrt{n-a^2}t \right) + \\ &+ \sum\_{i=0}^{\infty} \frac{t^{\frac{i}{2}}}{\Gamma\left(\frac{i}{2} + 1\right)} \sum\_{r=0}^{i} \binom{i}{r} a^{i-r} \left(n-a^2\right)^{\frac{r}{2}} \left(\cos \frac{\pi i}{2} - \frac{a}{\sqrt{n-a^2}} \sin \frac{\pi i}{2} \right) \right]. \end{split} \tag{27}$$

#### *3.2. The Particular Case Study*

For validating the obtained general Equation (27) of the fractional-order differential Equation (11), the particular case study for *α* = 0 is considered. In this case, the solution takes the following form:

$$w\_{p0}(t) = v\_{p\infty} \left[ 1 - 2e^{-nt} + E\_{\frac{1}{2},1}' \left( \sqrt{nt} \right) \right],\tag{28}$$

where the following modified Mittag–Leffler function is introduced:

$$E\_{\gamma,\theta}^{\langle c\rangle}(z) = \sum\_{i=0}^{\infty} \frac{z^i}{\Gamma(\gamma i + \beta)} \cos \frac{\pi i}{2},\tag{29}$$

which differs from the traditional one (18) by the multiplier cos *<sup>π</sup><sup>i</sup>* 2 .

Remarkably, the modified Mittag–Leffler function has the following peculiarity:

$$E\_{\frac{1}{2},1}^{\langle c\rangle}(z) \equiv e^{-z^2}. \tag{30}$$

In this case, a particular Equation (28) takes the form

$$
v\_{p0}(t) = v\_{p\infty}(1 - e^{-nt}),\tag{31}$$

which corresponds to the traditional one [28].

Thus, the general Equation (27) discovers new areas in studying the process of particle sedimentation in a fluid flow.

#### *3.3. Analysis of Leading Orders of the General Solution*

To analyze the leading orders of the general Equation (27), the dimensionless velocity *vp*(*τ*) <sup>=</sup> *vp*(*τ*) *vp*<sup>∞</sup> , the dimensionless time *τ* = *nt*, and the dimensionless coefficient of the Basset force *ε* = <sup>√</sup>*<sup>α</sup> <sup>n</sup>* are introduced. Accordingly, Equation (27) takes the following dimensionless form:

$$\begin{split} \overline{v}\_{p}(\tau) &= 1 - 2\varepsilon^{-\tau} \varepsilon^{2\varepsilon^{2}\tau} \Big( \cos 2\varepsilon \sqrt{1 - \varepsilon^{2}} \tau - \frac{\varepsilon}{\sqrt{1 - \varepsilon^{2}}} \sin 2\varepsilon \sqrt{1 - \varepsilon^{2}} \tau \Big) + \\ &+ \sum\_{i=0}^{\infty} \frac{\tau^{\frac{i}{2}}}{\Gamma\left(\frac{i}{2} + 1\right)} \sum\_{r=0}^{i} \binom{i}{r} \varepsilon^{i-r} (1 - \varepsilon^{2})^{\frac{\tilde{r}}{2}} \Big( \cos \frac{\pi i}{2} - \frac{\varepsilon}{\sqrt{1 - \varepsilon^{2}}} \sin \frac{\pi i}{2} \Big) . \end{split} \tag{32}$$

For the first approximation, in case of relatively small values of *ε* << 1, after expansion into a Maclaurin series leaving the first-order terms *ε* only, Equation (27) considering identity Equation (30) and dependence Equation (31) can be rewritten approximately as follows:

$$
\overline{\sigma}\_p(\tau) = 1 - \varepsilon^{-\tau} - \varepsilon E\_{\frac{1}{2},1}^{\langle s \rangle}(\sqrt{\tau}) = \overline{\sigma}\_{p0}(\tau) - \delta \overline{\sigma}\_p(\tau), \tag{33}
$$

where the first-order variation has been introduced:

$$
\delta\overline{\varpi}\_{\mathcal{P}}(\tau) = \varepsilon E\_{\frac{1}{2},1}^{\langle s\rangle}(\sqrt{\tau}).\tag{34}
$$

It contains the modified Mittag–Leffler function

$$E\_{\gamma,\theta}^{\langle s\rangle}(z) = \sum\_{i=0}^{\infty} \frac{z^i}{\Gamma(\gamma i + \beta)} \sin \frac{\pi i}{2} \tag{35}$$

for values *γ* = <sup>1</sup> <sup>2</sup> , *<sup>β</sup>* = 1 and argument *<sup>z</sup>* <sup>=</sup> <sup>√</sup>*τ*.

Due to the abovementioned results, the leading terms of the Equation (27) are the first terms in the double sum for the case of *i* = *r* = 0, where *εi*–*<sup>r</sup>* = *ε*<sup>0</sup> = 1:


Notably, due to Equation (33), an increase in the Basset force (the dimensionless parameter *ε*) leads to decreased sedimentation velocity.

For practical purposes, Simplification (33) can be applied for values of the dimensionless coefficient *ε* ≤ 0.23 because the maximum relative deviation of the sedimentation velocity from the accurate analytical Solution (27) is less than 10%. Therefore, if the dimensionless coefficient *ε* is more than this critical value, Equation (27) or (32) should be given preference.

Figure 1 presents the dimensionless solution Equation (32) of the particle sedimentation Equation (11) for different dimensionless parameters *ε* in a wide range (from 0 to 1).

**Figure 1.** Comparative analysis of the dimensionless solutions for the fractional-order integrodifferential equations of particle sedimentation.

Notably, all the solutions asymptotically approach the horizontal lines of the maximum sedimentation velocities.

#### *3.4. Comparison with the Numerical Simulation Results*

The next step for proving the obtained general solution Equation (27) is its comparison with the numerical simulation results.

The direct integration of the fractional-order differential Equation (11) can be realized with approach based on the S-approximation method [29] using the block-pulse operational matrix (*S*(*t*)) [30], the elements of which are as follows:

$$S\_i(t) = \frac{1}{2} \{ \operatorname{sign}(t - i\Delta t) - \operatorname{sign}[t - (i+1)\Delta t] \},\tag{36}$$

where *i*—plot index and Δ*t* = *tmax*/*N*—timestep (*tmax*—the maximum time range; *N* number of plots for numerical integration).

In this case, the initial fractional-order differential Equation (11) is transformed to the following fractional-order integral equation:

$$\upsilon\_p(t) + 2\varkappa I\_{0,t}^{\frac{1}{2}}[\upsilon\_p(t)] + nI\_{0,t}^1[\upsilon\_p(t)] = n\upsilon\_{p\infty}t\_\prime \tag{37}$$

which satisfies the zero initial condition, as *vp*(0) = 0.

The numerical solution of this integral equation is based on the following operational analog [31]:

$$n\left(\left[E\right] + 2\mathfrak{a}\left[P^{(0.5)}\right] + n\left[P^{(1)}\right]\right)\{Y\} = \{F\},\tag{38}$$

where {*Y*}—operational column-vector of the solution, [*E*]—the identity matrix with dimensions of *<sup>N</sup>* <sup>×</sup> *<sup>N</sup>*, and [*P*<0.5>], [*P*<1>]—elements of the operational matrix [*P<<sup>β</sup>*>] (*<sup>β</sup>* = 0.5, and *β* = 1, respectively), the elements of which are determined as follows:

$$P\_{i,j}^{\langle \beta \rangle} = p(\beta, i - j)\_{\prime} \tag{39}$$

where *p*(*β*, *r*)—the following function (|*r*| = 0, 1, . . . , *N*–1):

$$P\_{i,j}^{\langle \beta \rangle} = \frac{\Delta t^{\beta}}{\Gamma(\beta + 2)} \cdot \begin{cases} 1, & i = j; \\ (r+1)^{\beta+1} - 2r^{\beta+1} + (r-1)^{\beta+1}, & i > j \\ 0, & i < j. \end{cases} \tag{40}$$

Elements *Fi* (*i* = 1, 2, ... , *N*) of the operational column-vector of external impact {*F*} are determined as follows:

$$F\_{\bar{l}} = \frac{1}{\bar{l}t} \int\_{i\Delta t}^{(i+1)\Delta t} n v\_{p\infty} t dt = \left(i + \frac{1}{2}\right) \Delta t \cdot n v\_{p\infty}.\tag{41}$$

Using the inverse matrix method [32], the operational column-vector of the solution

$$\mathbb{E}\{Y\} = \left( [E] + 2a \left[ P^{\langle 0.5 \rangle} \right] + n \left[ P^{\langle 1 \rangle} \right] \right)^{-1} \{ F \} \tag{42}$$

allows us to obtain the following numerical solution of the initial Equation (33):

$$w\_p(t) = \sum\_{i=0}^{N-1} \mathcal{Y}\_i \mathcal{S}\_i(t). \tag{43}$$

Figure 1 also presents the numerically obtained results of solving the particle sedimentation Equation (37). Approximation lines for all the data completely coincide with the analytically obtained function Equation (32). This fact additionally proves the analytical approach's adequacy for solving the fractional-order differential equation of particle sedimentation.

Remarkably, the Basset force significantly increases the dimensionless sedimentation time determined from the condition *vp* = (1 − *δ*)*vp*∞, where *δ*—deviation of the sedimentation velocity from the stationary one *vp*∞, usually chosen as equal to *δ* = 0.05.

Particularly, if the Basset force is not considered (*ε* = 0), then the dimensionless sedimentation time is equal to 3.0 (Figure 1). However, even considering the parameter *ε* = 0.05, this dimensionless time parameter increases by 43% and becomes equal to 4.3.

The fact mentioned above explains the use of significant correction factors in traditional calculations of apparatuses for the sedimentation and deposition of the particles. However, using the fractional-order differential equation of particle sedimentation avoids unjustified correction factors in modeling processes and designing related chemical technology equipment.

#### **4. Discussion**

The fractional-order differential Equation (11) and its analytical solution Equation (27) eliminate the gap in studying the particle sedimentations and deposition processes.

Notably, Wan [33] obtained the fractional-order differential equation's analytical solution for the free-falling process. However, the considered hydraulic drag as a linear combination of the velocity *vp* and its square *v*<sup>2</sup> *<sup>p</sup>* was extended to a more general case of velocity analog based on the fractional derivative of particle displacement with the order in a range from 1 to 2. However, this approach does not consider the Basset force at all. Moreover, consideration of this force makes it possible to use the fractional derivative of particle displacement with the order of 3/2 only.

Visitskii et al. [34] investigated spherical particles' sedimentation in a viscous fluid considering the Basset force. However, the solution was obtained numerically for particular case studies. In addition, numerical simulation results can be obtained using the operational matrix based on the Legendre polynomials proposed by Saadatmandi and Dehghan [35]. However, the resulting curves presented in Figure 1 remain unchanged in terms of their approximation by the obtained analytical solution Equation (27).

Sobral et al. [36] applied the "Maple" software for obtaining the particular solution of the sedimentation equation. However, this solution is quite complicated for the analysis due flaws, including complex variables that need to reduce hyperbolic functions by trigonometric ones with imaginary coefficients. Moreover, in the case of zero Basset force, the corresponding solution should be reduced to the solution of the relaxation equation.

All these disadvantages were eliminated using the Mittag–Leffler function and its particular properties. In addition, the validity of the proposed solution for the imaginary arguments was justified. Finally, the singularity of the fractional-order differ-integral sedimentation equation was proven analytically.

Notably, the proposed approach can be applied not only for particle sedimentation and deposition problems but also for solving nutrient release problems from mineral fertilizers and the pneumatic classification of fine particles. Particularly, Equation (11) and its general solution Equation (27) allow us to clarify the concentration distribution for nutrients during their washing out from the fertilized soil [37]. The migration process of mineral components in the soil [38] can also be studied more precisely using the proposed fractional-calculus approach. The proposed methodology also extends the understanding of the distribution of fine particles' concentration and the height of the rhomb-shaped pneumatic classificators [39].

The fractional-order sedimentation equation can be further expanded to solve the problem of dispersed particle's deposition in separation channels [40] and for further precise modeling of the gas emission process [41] and designing gas-cleaning equipment [42].

#### **5. Conclusions**

In this article, the fractional origin for the Basset force was substantiated. As a result of modifying the particle motion equation moving with a time-varying velocity, the fractionalorder equation of particles' sedimentation was obtained. This equation differs from the traditional one by the semi-derivative of the particle velocity.

As a result of the integration of the proposed equation, its general solution was found analytically. This solution varies from the traditional one by the component based on the modified Mittag–Leffler function. Notably, for the particular case of zero Basset force, this equation reduces to the traditional one. Moreover, the belonging of this solution to the real range of values was theoretically proven.

In addition, the particle sedimentation's fractional-order differential equation was solved numerically using the analogous fractional-order integral equation. The numerical simulation was realized based on the S-approximation method using the block-pulse operational matrix. For normalizing the obtained solutions, the dimensionless particle velocity was introduced.

As a result, approximations of the numerically obtained case studies were presented graphically for a single variable dimensionless system parameter ranging from 0 to 1. The obtained analytical dependencies approximated the corresponding block-pulse curves. The adequacy of the proposed analytical approach for solving the fractional-order differential equation of particle sedimentation was proven analytically and numerically.

Notably, the Basset force significantly increases the sedimentation time. For example, if the Basset force is not considered, then the dimensional sedimentation time is equal to 3.0. However, when considering the Basset force with the dimensionless parameter equal to 0.05, this time increases by 43%.

Finally, the proposed mathematical model can be applied for modeling the processes of particle sedimentation in liquids, aerosol deposition in gas flows, and particle deposition in gas-dispersed systems as essential processes in the fields of mechanical and chemical engineering.

**Author Contributions:** Conceptualization and methodology—I.P.; software and validation—I.P., M.O., A.K. and R.O.; formal analysis—P.A.; investigation and resources—I.P., R.O. and B.M.; data curation—I.P. and M.O.; writing—original draft preparation—I.P.; writing—review and editing— M.O.; visualization—I.P. and A.K.; supervision—I.P. and M.O.; project administration—I.P. and P.A.; funding acquisition—M.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Education and Science of Poland: SBAD.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The authors acknowledge the research project of Sumy State University (state reg. no. 0120U102036) ordered by the Ministry of Education and Science of Ukraine and the International Association for Technological Development and Innovations, as well as Poland Ministry of Education and Science. The research results allow fulfilling the objectives of the perspective development plan of Sumy State University within the scientific direction "Technical Sciences".

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


Institute of Technology, New Delhi, India, 6–7 February 2020; Bansal, P., Tushir, M., Balas, V., Srivastava, R., Eds.; Advances in Intelligent Systems and Computing; Springer: Singapore, 2021; Volume 1164, pp. 93–101. [CrossRef]

