The slope and transformation of an elliptical trochoidal path was previously studied [30]. In this research, the elliptical trochoid was graphically analysed, as seen in Figures 9 and 21.

**Figure 21.** Graphical analysis obtaining of *aeff* in an elliptical trochoidal path.

Figure 22 shows the variation of *aeff* and the milling displacement angle with increments of 5◦.

**Figure 22.** *aeff* for the circular and elliptical trochoidal paths.

The 20◦ elliptical path keeps the value of *aeff* more stable, being necessary to optimize its parameters [30,31]. However, the circular trochoidal path was tested, with the objective of studying the feed per tooth variation to maintain the chip thickness constant.

Figure 23 shows the path of the milling tool centre, combining Figures 15 and 16, thus extending the trochoidal arch and rounding entry and exit paths.

During the tests, the milling tool "marena 965" (Figure 24) was used [32], with z = 5 teeth, ∅ = 12 mm, maximum cutting length = 36 mm. The cutting speed was 90 m/min and the feed was 0.08 mm/tooth for peripheral milling, with *ae* = 0.5 mm and *p* = 35 mm. With these conditions, the spiral path was milled several times with three similar machines (Figure 17).

**Figure 23.** Trochoidal roughing adapted to a spiral groove.

**Figure 24.** Milling tool used in the experiments.

Figure 25 shows the path in the NC display (FAGOR 8060TM).

**Figure 25.** Trochoidal-spiral path.

Figure 26 shows the real spiral after milling it (left), and the initial (a), intermediate (b) and final (c) Titanium chips.

**Figure 26.** One of the milling processes and types of chips (**a**–**c**).

Finally, Figures 27 and 28 show the tool wear obtained with constant *F* and with variable *F*, respectively.

**Figure 27.** Tool wear for constant *F*.

**Figure 28.** Tool wear for variable *F*.

#### **4. Discussion**

As shown in Figure 22, the 20◦ elliptical path can be interesting because *aeff* stays more stable, although its value is higher, decreasing the feed per tooth. In contrast, the 0◦ elliptical path starts with a high value of *aeff*, making it necessary to decrease the feed per tooth, although it could be increased at the exit.

The circular trochoidal path was chosen because its initial *aeff* is not high, obtaining its maximum value before tracing the centre of the semitrochoid. The feed per tooth has a high value at the entry and decreases until *aeffmax* and, after it, its value increases to the exit. From the exit to the following entry, the feed value can be significantly increased.

After analysing references [30] and [31], it would be interesting to test these new paths by varying the *F* in further research.

The tool holders must be hydraulic, particularly when the loads are high or the feeds vary a lot during the milling process, because other kinds of tool holders (e.g.,: Weldon or clamps) were not able to maintain the position of the tool.

In both milling processes, a frontal tool tooth was broken (Figures 27 and 28). The tool wear was slightly lighter with variable *F* and the tool coating was not peeled in any case. The diameter deviation between teeth, before and after the milling, was below 5 micrometres.

The milling process with variable *F*, even adapting S [28,32], suggests improvements which should be tested in further research.

An additional possibility would be elevating the milling tool in the exit points (helicoidally), descending in the entry points (Figures 15 and 17). This possibility was not implemented, according to the tool manufacturer indications (to avoid stress in the *Z* spindle and the entrapment of chip residues).

#### **5. Conclusions**

The machining time with the conditions shown in Figure 19, was 1 h 17 min with constant *F*, and 51 min with variable *F*. When the cutting speed was raised to 100 m/min and the feed per tooth was 0.07 mm, the machining times were 1 h 36 min and 59 min, respectively.

The wear of the tool edge decreases with the variable *F* strategy.

Finally, it was concluded that the variable *F* strategy was advantageous, but it can only be applied with milling machines that allow high acceleration of their axes (*X*, *Y*).

It is possible to reduce the energy consumption as well as the time required, increasing tool life and improving the sustainability of machining processes, thanks to several aspects. These aspects include the improved features of current machine tools, new machining strategies and advances in tool cooling, e.g., using cryogenic products [33] with phase change. It is also possible to use cryogenic cooling to improve the effectiveness of lubrication [34,35], significantly reducing the quantity required. Another interesting option is based on adding nanoparticles [36] to reduce tool wear or by locally applying a process of laser preheating [37].

**Author Contributions:** Conceptualization, C.G.-H., P.U.-A., J.-L.H.-T. and P.K.; methodology, C.G.- H., J.-J.G.-B., J.-J.V.-S., P.U.-A., J.-L.H.-T. and P.K.; software, J.-L.H.-T.; validation, C.G.-H., J.-J.G.-B., J.-A.B.-P., B.C.-Á., M.-Á.A.-S., F.V.-C., R.P.-C. and J.-L.H.-T.; formal analysis, C.G.-H. and J.-L.H.-T.; investigation, C.G.-H., P.U.-A., J.-A.B.-P., B.C.-Á., M.-Á.A.-S., F.V.-C., J.-L.H.-T. and P.K.; resources, J.-J.G.-B.; data curation, C.G.-H. and J.-L.H.-T.; writing—original draft preparation, C.G.-H., J.-J.G.-B., J.-J.V.-S. and J.-L.H.-T.; writing—review and editing, C.G.-H., P.U.-A., R.P.-C., J.-L.H.-T. and P.K.; supervision, C.G.-H., J.-J.G.-B., J.-J.V.-S., P.U.-A., J.-A.B.-P., B.C.-Á., M.-Á.A.-S., F.V.-C., R.P.-C. and J.-L.H.-T.; project administration, J.-J.G.-B. and F.V.-C.; funding acquisition, J.-J.G.-B. and F.V.-C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Spanish "Secretaría de Estado de Educación y Formación Profesional" grant for applied innovation and knowledge transfer in vocational training: AFP 19/00017, 16 September 2019.

**Data Availability Statement:** All data generated or analysed during this study are included in this published article.

**Acknowledgments:** The authors want to thank the following: Pixel Sistemas, for the CAM code in Solid Edge CAM Pro; Hexagon, for the CAM code in EDGECAM; Students of the NC and CAM modules, for simulating the NC and CAM code during the remote (due to COVID-19) practical sessions; Kiomac S.L., for the technical advice; Ysamco S.L., for preparing the Ti6V4Al parts; Marena S.L., for facilitating the access to all the required milling tools and for the technical advice.

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

#### **References**

## *Article* **Numerical Solution to Anomalous Diffusion Equations for Levy Walks**

**Viacheslav V. Saenko 1,2,\*, Vladislav N. Kovalnogov 1, Ruslan V. Fedorov 1,\* and Yuri E. Chamchiyan 1**

**Abstract:** The process of Levy random walks is considered in view of the constant velocity of a particle. A kinetic equation is obtained that describes the process of walks, and fractional differential equations are obtained that describe the asymptotic behavior of the process. It is shown that, in the case of finite and infinite mathematical expectation of paths, these equations have a completely different form. To solve the obtained equations, the method of local estimation of the Monte Carlo method is described. The solution algorithm is described and the advantages and disadvantages of the considered method are indicated.

**Keywords:** Levy walks; anomalous diffusion; fractional material derivative; combustion process; local estimate; Monte Carlo method

**MSC:** 58J65; 60J60; 60J70; 65C05

#### **1. Introduction**

At present, the theory of anomalous diffusion is rarely used to describe the combustion processes of a substance, although there are all the prerequisites for this. For example, when a substance is burning, various kinds of turbulence are developing. As a result, the diffusion packet width Δ(*t*) stops obeying the law Δ(*t*) ∝ *t γ* with an exponent *γ* = 1/2 and starts growing with time by the law with an exponent *γ* = 1/2, which testifies to the appearance of anomalous diffusion. Signs of the appearance of anomalous diffusion at thermal transport in a low-dimensional system are indicated in the papers [1,2]. This assumption is also supported by experimental data. For example, in the work [3], the thermal radiation in the combustion chamber during the combustion of natural gas and acetylene was studied, depending on the level of fuel enrichment with oxygen. Time series analysis showed that the combustion process at any oxygen concentration is subdiffusive in nature.

The assumption about the formation of anomalous diffusion during the combustion of a substance allows us to introduce into consideration the fractional differential equations of anomalous diffusion. An effective coefficient of heat conductivity for the Levy–Fokker– Plank equation was obtained in the papers [4,5]. In the papers [6–8], to describe the combustion process, it is proposed to use the fractional differential equation of anomalous diffusion:

$$\_0D\_t^\beta u = \partial\_{xx}u + f(u), \quad t > 0, \quad 0 < x < L,\tag{1}$$

where 0*Dβ t* is the fractional Riemann–Liouville derivative [9] of the order 0 < *β* < 1 by time and *∂xx* is the classical particular derivative of the second order by coordinate. In the paper [10], a two-dimensional combustion model with a fractional time derivative was studied. To solve the obtained equation of diffusion, the authors develop an adaptive

**Citation:** Saenko, V.V.; Kovalnogov, V.N.; Fedorov, R.V.; Chamchiyan, Y.E. Numerical Solution to Anomalous Diffusion Equations for Levy Walks. *Mathematics* **2021**, *9*, 3219. https:// doi.org/10.3390/math9243219

Academic Editors: Camelia Petrescu and Valeriu David

Received: 20 November 2021 Accepted: 7 December 2021 Published: 13 December 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/).

finite-difference discontinuous Galerkin method. A modification of the Equation (1) in the case of the dependence of the diffusivity on temperature and relaxation time has been considered in the article [11]. In the papers [8,12], a fractional-differential combustion model with the first derivative with respect to time and a fractional derivative with respect to the spatial variable is considered:

$$
\partial\_t u = \frac{\partial^\mathbf{a}}{\partial |\mathbf{x}|^4} u + f(u), \quad u(\mathbf{x}, 0) = u\_0(\mathbf{x}), \quad u(a, t) = u(b, t) = 0. \tag{2}
$$

Here, *∂t* is the partial time derivative, and *∂α ∂*|*x*|*α* is the fractional-differential Riesz operator [9]. To solve this equation, a finite-difference scheme using an adaptive strategy is described in the papers [8,12]. In the paper [13], a fractional-differential generalization of the kinetic equation was obtained that describes the relationship between the radius of the ball and time in the model of the combustion of a fireball, theoretically predicted by the Soviet physicist Ya.B. Zeldovich [14].

However, when using the Equations (1) and (2), it should be kept in mind that they are based on the process of Continuous Time Random Walk (CTRW) [15–22]. This process assumes that a particle instantly moves from one point in space to another at a random distance *Ri*, *i* = 1, 2, 3, ... , and then rests at this point for a random time *Ti*, *i* = 1, 2, 3, ... All these random values *Ri* and *Ti* are independent of each other and between one another and they are distributed by the laws *pR*(*x*) ∝ *x*−*α*−1, *x* → ∞ and *qT*(*t*) ∝ *t* −*β*−1, *t* → ∞, respectively.

Depending on the value *α*, the path distribution has different properties. If *α* = 2, then the distribution has finite mathematical expectation and variance, at 1 < *α* < 2, mathematical expectation is finite and variance is infinite, at 0 < *α* - 1 both mathematical expectation and variance are infinite. The instantaneousness of jumps means that, for an arbitrarily small time interval from the initial one, the particle can be at an arbitrarily large distance from the source. In some cases, this non-physical behavior of a particle does not lead to any contradiction with the experiment. For example, in the case of the normal diffusion (*α* = 2, *β* = 1), when the distribution of paths and rest times have a finite variance, the process is described by the classical diffusion equation, the solution to which is expressed in terms of the normal distribution. It is well known that the normal distribution is nonzero on the entire number axis, which indicates the infinite velocity of the particle in the inherent walk model. However, in view of the finite variance of the distribution of paths, the infinite velocity of motion is compensated for by a small value of the paths. The situation is completely different in the case of anomalous diffusion. As the exponent *α* decreases, the probability of greater paths increases and at values *α* < 1 this probability turns out to be significant. Therefore, it is necessary to use the anomalous diffusion model with a certain degree of caution, especially when considering problems with limited spatial geometry or processes limited in time, for example, to describe combustion processes in furnaces.

One of the ways to eliminate the difficulty described above is to introduce a constant final velocity of the particle. One of the first works in which a constant velocity of particle motion was introduced into the model of anomalous diffusion is the work [23]. In this paper, the authors called this model Levy walks. Further study of this model was carried out in the works [24–26]. The work [27] is devoted to the study of Levy walks in bounded and semi-bounded spaces. In the work [28], kinetic equations of anomalous diffusion with a finite velocity are obtained, the root-mean-square deviation is investigated, and an exact solution to the kinetic equation in the one-dimensional case is obtained. The work [29] is devoted to the study of statistical moments for the case of Levy random walks with a random finite velocity without traps, and the case of multidimensional walks with traps of arbitrary type with constant velocity is considered in the works and the case of multidimensional walks with traps of arbitrary type with constant velocity is considered in the works [30–33]. The work [34] examines the influence of the final velocity on the spatial distribution of particles in Levy walks with exponential traps. In this work, it was found that, in the case 1 < *α* - 2, taking account of the constant velocity of particle motion is reduced to a decrease in the diffusion coefficient in the equation of anomalous diffusion (2). In the case *α* < 1 distributions take *W*-like or *U*-like form and cannot be described by the Equation (2). The papers [35–37] study the process of Levy random walks without traps and it was shown that in the case of infinite mathematical expectation of the distribution of paths, the asymptotic distributions have *U*-like and *W*-like form. In the paper [38] an expression for current within framework of the Levy walks model was obtained. The obtained expression generalize the Fourier's law to the case of anomalous thermal transport for the Levy diffusion model.

The work [39] succeeded in addressing the problem of describing Levy random walks in the case *α* < 1. In this work, it was shown that, in the case of an infinite mathematical expectation of the distribution of paths, to take account of the finite velocity, it is necessary to replace the fractional Laplacian in the equation of anomalous diffusion by a material derivative of a fractional order. Later, the authors of the works [40–44] come to the same conclusion. The solution to the equation of anomalous diffusion with a material derivative of the fractional order was obtained in the works [44–46]. In these works, it was shown that the solution to this equation is expressed in terms of the Lamperti distribution.

In this paper, we consider a method for the numerical solution to the equations of anomalous diffusion taking account of the constant finite velocity of the particle motion between collisions. The work is organized as follows. In Section 2, a kinetic equation describing the considered process of walks is derived. In Section 3, fractional differential equations and the solutions of these equations are obtained, describing the asymptotic (at *t* → ∞ and *x* → ∞) distribution of particles. It was shown that, in the case of an infinite mathematical expectation of the distribution of paths (0 < *α* < 1) and finite mathematical expectation 1 < *α* < 2, the behavior of the process is described by completely different equations. In Section 4, a numerical method for solving the kinetic equation based on the method of local density estimation is considered.

#### **2. Kinetic Equation of the Random Walk Process**

To obtain the kinetic equation, we will use the approach proposed in the paper [28], which was developed in the paper [47]. We will consider the density of collisions *f*(**r**, **p**, *t*), where **r** is the radius-vector of the particle, **p**—is the particle momentum, *t*—is time. The value *f*(**r**, **p**, *t*)*d***r** *d***p** *dt* is the number of collisions in the volume element *d***r** of the vicinity of the point **r** for the interval of time *dt*, at which the momentum of the particle takes on a value from **p** to **p** + *d***p**. We will consider the nonrelativistic case **p** = *m***v**. Without loss of generality, we assume that *m* = 1. The paper [47] shows that, with the presence *n* of discrete states, the value *f*(**r**, **v**, *t*) can be represented in the form:

$$f(\mathbf{r}, \mathbf{v}, t) = \sum\_{j=1}^{n} f\_j(\mathbf{r}, \mathbf{v}, t)\_{\prime} \tag{3}$$

where

$$f\_{\vec{j}}(\mathbf{r}, \mathbf{v}, t) = s\_{\vec{j}}(\mathbf{r}, \mathbf{v}, t) + \sum\_{i=1}^{n} c\_{ij} \int\_{0}^{t} k\_{i}(\mathbf{r}) d\tau \int \mathcal{W}\_{\vec{ij}}(\mathbf{\Omega}', \mathbf{\Omega}) d\mathbf{\Omega}' \times$$

$$\int f\_{i}(\mathbf{r} - \mathbf{v}' \mathbf{\Omega}' \tau, \mathbf{v}' \mathbf{\Omega}', t - \tau) h\_{\vec{ij}}(\mathbf{v}', \mathbf{v}) d\mathbf{v}'. \tag{4}$$

Here, *ki*(*τ*) is the probability density distribution of the residence time in the state *i*, *cij*—the probabilities of transition from the state *i* into the state *j*, *Wij*(**Ω** , **Ω**) is the probability density that before collision the velocity had the direction **Ω** , after collision the direction took the value **Ω**, *hij*(*v* , *v*) is the density of the probability of the change in velocity from the value *v* to *v*, *sj*(**r**, **v**, *t*) is the density of new particle sources in the state *j*, **v** = *v***Ω**, *v* = |**v**|, *d***v** = *dv d***Ω**, the summation is carried out over all possible previous states. The values *cij*, *Wij*(**Ω** , **Ω**) and *hij*(*v* , *v*) are normalized:

$$\sum\_{j=1}^{n} \mathfrak{c}\_{ij} = 1, \quad \int \mathcal{W}\_{ij}(\Omega', \Omega) d\Omega = 1, \quad \int h\_{ij}(v', v) dv = 1. \tag{5}$$

The transition from the density of collisions *f*(**r**, **v**, *t*) to the phase density *ψ*(**r**, **v**, *t*) is carried out with the help of the integral,

$$
\psi(\mathbf{r}, \mathbf{v}, t) = \int\_0^t K(\tau) f(\mathbf{r} - v\Omega\tau, v\Omega, t - \tau) d\tau,\tag{6}
$$

where *K*(*t*) = & ∞ *t k*(*τ*)*dτ*. Substituting (3) in (6) we get that the phase density has the form of the sum,

$$
\psi(\mathbf{r}, \mathbf{v}, t) = \sum\_{j=1}^{n} \psi\_{j}(\mathbf{r}, \mathbf{v}, t),
\tag{7}
$$

where

$$\psi\_{\dot{\jmath}}(\mathbf{r}, \mathbf{v}, t) = \int\_0^t K\_{\dot{\jmath}}(\tau) f\_{\dot{\jmath}}(\mathbf{r} - v\Omega\tau, v\Omega, t - \tau) d\tau,\tag{8}$$

where *Kj*(*t*) = & ∞ *t kj*(*τ*)*dτ*. The physical interpretation of the last expression is simple. To detect the particle in the state *j* of the vicinity *d***r** of the point **r** with the velocity in the interval from **v** to **v** + *d***v** at the moment of time from *t* to *t* + *dt* the particle must pass to this state at the point **r** − *v***Ω***τ* at the moment of time *t* − *τ* and stay in this state during the time greater than *τ*. Transition to the density of particles *ρ*(**r**, *t*) is carried out with the help of the integral,

$$
\rho(\mathbf{r},t) = \int \psi(\mathbf{r}, \mathbf{v}, t) d\mathbf{v}.\tag{9}
$$

The system of Equations (4), (7) and (8), together with conditions (5) describes practically any process of random walks with *n* discrete states under fairly general assumptions about the scattering indicatrix *Wij*(**Ω** , **Ω**) and the law of redistribution of speed *hij*(*v* , *v*). In this work, using these equations, we obtain kinetic equations describing Levy walks with a constant velocity of motion between two successive scatterings of a particle.

We define the process of walks as follows. There is only this state—the state of motion (*n* = 1). A particle moves at a constant velocity *v* between two successive collisions. After the collision, the particle changes its direction, which is determined by the scattering indicatrix *W*(**Ω**). After which, the particle continues to move in a new direction with the same constant velocity *v*. Random times between two successive collisions of a particle *Ti*, *i* = 1, 2, 3, ... are independent random values. Since the motion occurs with a finite speed, then for times *Ti* the particle covers the path *Ri* = *vTi*, *i* = 1, 2, 3, . . . . The values *Ri* are the paths of particles.

Since there is only one state, then in Equations (4), (7) and (8), we need to put *n* = 1, *c*11 = 1. We will also assume that the source can be represented in the form *s*(**r**, **v**, *t*) = *s*(**r**, *t*)*hs*(*v*)*Ws*(**Ω**) and, for brevity, we omit the subscript indicating the status number. As a result, we obtain:

$$\psi(\mathbf{r}, \mathbf{v}, t) = \int\_0^t K(\tau) f(\mathbf{r} - v\Omega \tau, v\Omega, t - \tau) d\tau,\tag{10}$$

$$\begin{split} f(\mathbf{r}, \mathbf{v}, t) &= \mathbf{s}(\mathbf{r}, t) h\_{\mathbf{s}}(\mathbf{v}) \mathcal{W}\_{\mathbf{s}}(\boldsymbol{\Omega}) \\ &+ \int\_{0}^{t} k(\tau) d\tau \int \mathcal{W}(\boldsymbol{\Omega}', \boldsymbol{\Omega}) d\boldsymbol{\Omega}' \int f(\mathbf{r} - \mathbf{v}' \boldsymbol{\Omega}' \tau, \mathbf{v}' \boldsymbol{\Omega}', t - \tau) h(\mathbf{v}', \mathbf{v}) d\boldsymbol{v}'. \end{split} \tag{11}$$

The motion velocity between two successive collisions will be denoted by *v*0. Let us also assume that the direction of motion after each collision does not depend on the previous direction of motion. In view of the foregoing, we get:

$$h(v',v) = h\_s(v) = \delta(v - v\_0), \quad \mathcal{W}(\Omega',\Omega) = \mathcal{W}\_s(\Omega) = \mathcal{W}(\Omega).$$

Now, substituting these expressions in Equation (11), we obtain:

$$f(\mathbf{r}, \mathbf{v}, t) = \mathcal{W}(\Omega) \delta(\mathbf{v} - \mathbf{v}\_0) \left( s(\mathbf{r}, t) + \int\_0^t k(\tau) d\tau \int d\Omega' \int f\left(\mathbf{r} - v' \Omega' \tau, v' \Omega', t - \tau\right) d\upsilon' \right). \tag{12}$$

From this relation, it is clear that, in cases when the densities of transition probabilities *h*(*v* , *v*) and *Wij*(**Ω** , **Ω**) do not depend on the previous value *v* and **Ω** , then the density of the collision can be represented in the form of the product *f*(**r**, *v***Ω**, *t*) = *W*(**Ω**)*h*(*v*)*F*(**r**, **t**). Thus, we obtain:

$$f(\mathbf{r}, v\mathbf{\Omega}, t) = \mathcal{W}(\mathbf{\Omega})\delta(v - v\_0)F(\mathbf{r}, \mathbf{t}).\tag{13}$$

Now substituting this relation in (12) and by integrating over *dv d***Ω** we get the equation for *F*(**r**, *t*):

$$F(\mathbf{r},t) = \mathbf{s}(\mathbf{r},t) + \int\_0^t k(\tau)d\tau \int \mathcal{W}(\Omega')d\Omega' \int F(\mathbf{r} - \mathbf{v}'\Omega'\mathbf{r}, t - \tau)\delta(\mathbf{v}' - \mathbf{v}\_0)d\mathbf{v}'$$

$$= \mathbf{s}(\mathbf{r},t) + \int\_0^t k(\tau)d\tau \int F(\mathbf{r} - \mathbf{v}\_0\Omega'\mathbf{r}, t - \tau)\mathcal{W}(\Omega')d\Omega'. \tag{14}$$

Here, we used the normalization condition & *W*(**Ω** )*d***Ω** = 1. The physical meaning of the quantity *F*(**r**, *t*) is quite simple. This is the density of collisions in the volume element *d***r** of the vicinity of the point **r**.

Let us now pass from the collision density to the phase density *ψ*(**r**, **v**, *t*), and then to the density of particles *ρ*(**r**, *t*). To this end, we put the expression (13) in (10). As result, we obtain:

$$
\psi(\mathbf{r}, \mathbf{v}, t) = \int\_0^t K(\tau) W(\Omega) \delta(v - v\_0) F(\mathbf{r} - v \Omega \tau, t - \tau) d\tau.
$$

Now integrating this expression over *dv d***Ω** and taking account of the ratio (9), we get the equation for the density of particles:

$$\rho(\mathbf{r},t) = \int\_0^t \mathcal{K}(\tau)d\tau \int \mathcal{W}(\Omega)F(\mathbf{r}-\upsilon\eta\Omega\tau, t-\tau)d\Omega. \tag{15}$$

For the further solution to the obtained equations, it turns out to be convenient to pass from the time *τ* to the particle path *ξ*. By substituting the integration variable *τ* = *ξ*/*v*0, in Equations (14) and (15), we get:

$$\rho(\mathbf{r},t) = \frac{1}{v\_0} \int\_0^{v\_0t} P(\boldsymbol{\xi}) d\boldsymbol{\xi} \int F(\mathbf{r} - \boldsymbol{\xi}\mathbf{\Omega}, t - \boldsymbol{\xi}/v\_0) W(\mathbf{\Omega}) d\mathbf{\Omega},\tag{16}$$

$$F(\mathbf{r},t) = s(\mathbf{r},t) + \int\_0^{v\_0t} p(\xi)d\xi \int F(\mathbf{r} - \xi\mathbf{\Omega}', t - \xi'/v\_0) \mathcal{W}(\mathbf{\Omega}')d\mathbf{\Omega}'... \tag{17}$$

Here, the following notation was introduced: *p*(*ξ*) = 1 *v*0 *k*(*ξ*/*v*0) is the density of probability of the path distribution and *K*( *ξ v*0 ) ≡ *P*(*ξ*) = & ∞ *ξ p*(*y*)*dy*.

From this system of equations, it is possible to exclude the equation for *F*(**r**, *t*). To do this, we will substitute Equation (17) in (16) and change the order of integration in the second summand. As a result, we obtain the equation for density,

$$\begin{split} \rho(\mathbf{r},t) = \frac{1}{v\_0} \int\_0^{v\_0t} P(\boldsymbol{\xi}) d\boldsymbol{\xi} \int \mathbf{s}(\mathbf{r} - \boldsymbol{\xi}\boldsymbol{\Omega}, t - \boldsymbol{\xi}/v\_0) \mathcal{W}(\boldsymbol{\Omega}) d\boldsymbol{\Omega} \\ &+ \int\_0^{v\_0t} p(\boldsymbol{\xi}) d\boldsymbol{\xi} \int \rho(\mathbf{r} - \boldsymbol{\xi}\boldsymbol{\Omega}', t - \boldsymbol{\xi}/v\_0) \mathcal{W}(\boldsymbol{\Omega}') d\boldsymbol{\Omega}'. \end{split} \tag{18}$$

The first summand in this equation describes unscattered radiation. The second summand describes multiply scattered radiation. This equation describes the random walk of a particle with constant velocity in three-dimensional space with an arbitrary distribution of paths.

Let us simplify the problem and consider one-dimensional particle walks. Let the random walk process occur along the axis *x*. In this case, the function *W*(**Ω**) takes the form:

$$\mathcal{W}(\Omega) = \mathcal{W}(\theta, \varphi) = \frac{1}{\sin \theta} (\omega\_1 \delta(\varphi) + \omega\_2 \delta(\varphi - \pi)) \delta(\theta - \pi/2), \tag{19}$$

where *ω*1 and *ω*2 are the probabilities of motion in the positive and negative directions of the axis *Ox* respectively and *ω*1 + *ω*2 = 1. Now substituting (19) in Equation (18) and considering that **Ω** = (sin *θ* cos *ϕ*, sin *θ* sin *ϕ*, cos *θ*), *ρ*(**r**, *t*) = *ρ*(*x*, *y*, *z*, *t*), *s*(**r**, *t*) = *s*(*x*, *y*, *z*, *t*), *d***Ω** = sin *θdθdϕ* and integrating the resulting equation over the angular variables, we obtain:

$$\begin{split} \rho(\mathbf{x},y,z,t) &= \frac{1}{\upsilon\_{0}} \int\_{0}^{\upsilon\_{0}t} P(\xi) (\omega\_{1}\mathbf{s}(\mathbf{x}-\boldsymbol{\xi},\mathbf{y},z,t-\boldsymbol{\xi}/\upsilon\_{0}) + \omega\_{2}\mathbf{s}(\mathbf{x}+\boldsymbol{\xi},\mathbf{y},z,t-\boldsymbol{\xi}/\upsilon\_{0})) d\mathbf{f} \\ &+ \int\_{0}^{\upsilon\_{0}t} p(\boldsymbol{\xi}) (\omega\_{1}\rho(\mathbf{x}-\boldsymbol{\xi},\mathbf{y},z,t-\boldsymbol{\xi}/\upsilon\_{0}) + \omega\_{2}\rho(\mathbf{x}+\boldsymbol{\xi},\mathbf{y},z,t-\boldsymbol{\xi}/\upsilon\_{0})) d\boldsymbol{\xi}. \end{split}$$

Since random walks along the axis *Ox* are considered, then,

$$
\rho(\mathbf{x}, y, z, t) = \rho(\mathbf{x}, t)\delta(y)\delta(z), \quad \mathbf{s}(\mathbf{x}, y, z, t) = \mathbf{s}(\mathbf{x}, t)\delta(y)\delta(z).
$$

Now substituting these expressions into the previous equation and integrating over the variables *y* and *z*, we finally obtain:

$$\rho(\mathbf{x},t) = \frac{1}{\upsilon\_0} \int\_0^{\upsilon\_0 t} P(\xi) (\omega\_1 \mathbf{s}(\mathbf{x} - \boldsymbol{\xi}, t - \boldsymbol{\xi}/\upsilon\_0) + \omega\_2 \mathbf{s}(\mathbf{x} + \boldsymbol{\xi}, t - \boldsymbol{\xi}/\upsilon\_0)) d\xi$$

$$+ \int\_0^{\upsilon\_0 t} p(\xi) (\omega\_1 \rho(\mathbf{x} - \boldsymbol{\xi}, t - \boldsymbol{\xi}/\upsilon\_0) + \omega\_2 \rho(\mathbf{x} + \boldsymbol{\xi}, t - \boldsymbol{\xi}/\upsilon\_0)) d\xi. \tag{20}$$

The first component in this equation describes unscattered particles that, after escaping from the source, did not have a single collision and move in positive and negative directions, respectively. The second component describes multiply scattered particles, which at the moment of time *t* − *ξ*/*v*0 had a collision and after that they began their motion in positive and negative directions, respectively.

#### **3. Asymptotic Solution to a Kinetic Equation**

An asymptotic solution to this equation can be found. To do this, we perform the Fourier–Laplace transform:

$$\rho(k,\lambda) = \int\_0^\infty dt \int\_{-\infty}^\infty e^{ikx - \lambda t} \rho(\mathbf{x}, t) d\mathbf{x} \tag{21}$$

of the Equation (20). As a result, we obtain:

$$\begin{split} \boldsymbol{\hat{\rho}}(k,\lambda) = \frac{1}{v\_0} \boldsymbol{\hat{s}}\_2(k,\lambda) \Big( \boldsymbol{\omega}\_1 \boldsymbol{\hat{P}}(\lambda/v\_0 - ik) + \boldsymbol{\omega}\_2 \boldsymbol{\hat{P}}(\lambda/v\_0 + ik) \Big) \\ &+ \boldsymbol{\hat{\rho}}(k,\lambda) (\boldsymbol{\omega}\_1 \boldsymbol{\hat{p}}(\lambda/v\_0 - ik) + \boldsymbol{\omega}\_2 \boldsymbol{\hat{p}}(\lambda/v\_0 + ik)). \end{split}$$

As a result, we obtained a simple algebraic equation, the solution to this equation has the form:

$$
\hat{\rho}(k,\lambda) = \frac{\frac{1}{v\_0}\\$(k,\lambda)\left(\omega\_1\not{\mathcal{P}}(\lambda/v\_0 - ik) + \omega\_2\not{\mathcal{P}}(\lambda/v\_0 + ik)\right)}{1 - \omega\_1\not{\mathcal{P}}(\lambda/v\_0 - ik) - \omega\_2\not{\mathcal{P}}(\lambda/v\_0 + ik)}.
$$

Let there be a point instantaneous source *s*(*x*, *t*) = *δ*(*x*)*δ*(*t*). This means that *s*ˆ(*k*, *λ*) = 1. Considering that

$$
\hat{P}(\lambda) = \int\_0^\infty e^{-\lambda \frac{\pi}{4}} P(\xi) d\xi = \frac{1 - \hat{p}(\lambda)}{\lambda},
$$

we obtain

$$\rho(k,\lambda) = \frac{1/v\_0}{1 - \mathcal{W}(k, v\_0, \lambda)} \left( \omega\_1 \frac{1 - \mathcal{P}\left(\frac{\lambda}{v\_0} - ik\right)}{\frac{\lambda}{v\_0} - ik} + \omega\_2 \frac{1 - \mathcal{P}\left(\frac{\lambda}{v\_0} + ik\right)}{\frac{\lambda}{v\_0} + ik} \right),\tag{22}$$

where *W*(*k*, *v*, *λ*) = *ω*1 *p*ˆ(*λ*/*v* − *ik*) − *ω*2 *p*ˆ(*λ*/*v* + *ik*). This solution describes the spatial distribution of particles with random walks of a particle at a constant velocity, with an arbitrary distribution of paths. This solution is not new and was obtained earlier (see, for example, [40,42,48–52]). In this paper, we will consider the asymptotic solution to this equation in the case when the distribution of paths has asymptotics of the form:

$$p(\mathbf{x}) \propto a \mathbf{x}\_0^{\mathbf{a}} \mathbf{x}^{-a-1}, \quad 0 < \mathbf{a} < \mathbf{2}, \quad \mathbf{x} \to \infty,\tag{23}$$

where *x*0 = (Γ(1 − *α*) sin( *π* 2 (1 − *α*)))−1/*α*. The cases of other distributions of paths are considered in the work [52].

Let us consider the case 0 < *α* < 1. In this case, the Laplace transform of the density (23) has the form (see [33]):

$$\mathcal{P}(\lambda) \approx 1 - (\lambda \ge\_0)^{\alpha} \Gamma(1-\alpha), \quad \lambda \to 0.$$

Now substituting this expression into the solution (22), we obtain:

$$\rho(k,\lambda) = \frac{\omega\_1(\lambda - ikv\_0)^{a-1} + \omega\_2(\lambda + ikv\_0)^{a-1}}{\omega\_1(\lambda - ikv\_0)^a + \omega\_2(\lambda + ikv\_0)^a}.\tag{24}$$

This expression completely coincides with the result obtained in [37,49,53]. Taking account of the fact that the multiplier (*λ* ± *ikv*)*α* is the Fourier–Laplace transform of the fractional material derivative [40,54],

$$(\lambda \pm \imath kv)^{a} f(k, \lambda) = \int\_{0}^{\infty} dt \int\_{-\infty}^{\infty} e^{\imath kx - \lambda t} f(x, t) \left(\frac{\partial}{\partial t} \pm v \frac{\partial}{\partial x}\right)^{a} dx \,\omega$$

we obtain that the process is described by the fractional differential equation [44,53–55]:

$$\begin{aligned} \left[\omega\_1 \left(\frac{\partial}{\partial t} - v\_0 \frac{\partial}{\partial \mathbf{x}}\right)^a + \omega\_2 \left(\frac{\partial}{\partial t} + v\_0 \frac{\partial}{\partial \mathbf{x}}\right)^a\right] \rho(\mathbf{x}, t) &= \\ &= \left[\omega\_1 \left(\frac{\partial}{\partial t} - v\_0 \frac{\partial}{\partial \mathbf{x}}\right)^{a-1} + \omega\_2 \left(\frac{\partial}{\partial t} + v\_0 \frac{\partial}{\partial \mathbf{x}}\right)^{a-1}\right] \delta(\mathbf{x}) \delta(t). \end{aligned} \tag{25}$$

By changing the variable *s* = *ik* in (24), one can perform the inverse transformation using the method described in the work [56] (see also [37]). As result, we have:

$$\rho(\mathbf{x},t) = \frac{2\sin(\pi a)}{\pi v\_0 t} \frac{\omega\_1 \omega\_2 \left(1 - (\mathbf{x}/v\_0 t)^2\right)^{a-1}}{\omega\_1^2 (1 - \mathbf{x}/v\_0 t)^{2a} + \omega\_2^2 (1 + \mathbf{x}/v\_0 t)^{2a} + 2\omega\_1 \omega\_2 (1 - (\mathbf{x}/v\_0 t)^2) \cos(\pi a)}. \tag{26}$$

In the case 1 < *α* < 2, the distribution (23) has a mathematical expectation. In view of the fact that the Laplace transform of density (23) takes the form (see [33]):

$$
\mathfrak{p}(\lambda) \approx 1 - m\_1 \lambda + m\_2 \lambda^\alpha, \quad \lambda \to 0,
$$

where *m*1 = *αx*0 *α*−1 , *m*2 = *xα* 0Γ(2−*α*) *α*−1 . Now substituting this expression in (22) and considering asymptotics *λ* → 0, *k* → 0, |*λ*/*kv*| → 0, we get:

$$\rho(k,\lambda) \approx \frac{1}{\lambda + ikv\_0(\omega\_2 - \omega\_1) - \frac{m\_2}{m\_1}(\omega\_1(-ikv\_0)^a + \omega\_2(ikv\_0)^a)}.\tag{27}$$

To obtain an equation describing the process of random walks, we rewrite this expression in the form:

$$
\lambda \dot{\rho}(k,\lambda) - 1 = \dot{\rho}(k,\lambda) \left( ikv\_0(\omega\_1 - \omega\_2) + \frac{m\_2}{m\_1} (\omega\_1(-ikv\_0)^a + \omega\_2(ikv\_0)^a) \right).
$$

Let us simplify the problem and consider symmetric random walks *ω*1 = *ω*2 = 1/2, we obtain:

$$
\lambda \rho(k, \lambda) - 1 = \rho(k, \lambda) \frac{m\_2}{2m\_1} v\_0^{\alpha} |k|^{\alpha} \left( (-i)^{\alpha} + i^{\alpha} \right).
$$

In view of the fact that 1 2 ((−*i*)*α* + *i α*) = cos( *π* 2 *α*) we get:

$$
\lambda \not{\rho}(k,\lambda) - 1 = -D\_a \not{\rho}(k,\lambda)|k|^a. \tag{28}
$$

Here, *Dα* = *m*2 *m*1 *vα* 0 sin *π* 2 (*α* − 1) . To perform the inverse Fourier–Laplace transform of this equation, we need the Riesz fractional differentiation operator:

$$\frac{\partial^{\nu}}{\partial |\boldsymbol{x}|^{\boldsymbol{\nu}}} = \frac{-1}{2\cos(\pi\nu/2)}( \_{-\infty}D\_{\boldsymbol{x}}^{\nu} + \prescript{\infty}{}{D\_{\boldsymbol{x}}^{\nu}}) .$$

The Fourier transform of this operator has the form (see, for example, [9]):

$$\int\_{-\infty}^{\infty} e^{ikx} \frac{\partial^{\nu} f(x)}{\partial |x|^{\nu}} dx = -|k|^{\nu} \hat{f}(k).$$

Taking this relation into account, it is possible to perform the inverse Fourier–Laplace transform of Equation (28). As a result, we obtain the equation of anomalous diffusion:

$$\frac{\partial \rho(\mathbf{x}, \mathbf{t})}{\partial t} = D\_a \frac{\partial^a \rho(\mathbf{x}, \mathbf{t})}{\partial |\mathbf{x}|^a},\tag{29}$$

with the initial condition *ρ*(*x*, 0) = *δ*(*x*)*δ*(*t*).

The solution to this equation can be obtained by performing the inverse Fourier– Laplace transform of the solution (27). It has already been done in the works [52,53]. As a result, the inverse Fourier–Laplace transform has the form:

$$\rho(\mathbf{x},t) = \frac{1}{(D\_{\mathbf{a}}t)^{1/a}} \mathbf{g}\left((\mathbf{x} - \omega vt)(D\_{\mathbf{a}}t)^{-1/a}; \mathbf{a}, \omega, \mathbf{0}, 1\right),$$

where *g*(*x*; *α*, *ω*, 0, 1) is the density of the standard stable law and *ω* = *ω*1 − *ω*2. In the case of symmetric walks *ω*1 = *ω*2 = 1/2, this expression takes the form:

$$\rho(\mathbf{x},t) = \frac{1}{(D\_{\mathbf{a}}t)^{1/a}} \mathbf{g}\left(\mathbf{x}(D\_{\mathbf{a}}t)^{-1/a}; \mathbf{a}, \mathbf{0}, \mathbf{0}, \mathbf{1}\right). \tag{30}$$

As we can see, the equation solution (29) is expressed in terms of a symmetrical stable law [20,57,58].

Thus, taking account of the constant velocity of the particle in the Levy walk model leads to the need to consider two cases 0 < *α* < 1 and 1 < *α* < 2. The principal difference between these two cases is that, at 0 < *α* < 1, the mathematical expectation of the distribution (23) is equal to infinity and in the case 1 < *α* < 2, the mathematical expectation is finite. It is this fact that leads to different asymptotic distributions of particles. In the case of an infinite mathematical 0 < *α* < 1 the probability of the appearance of large paths turns out to be significant, and leads to the fact that the overwhelming majority of particles move along the front of the distribution *x* = ±*vt*. This is especially clear at values *α* < 0.6. In this case, the *U*—like form of the asymptotic distribution of particles is formed. As the parameter *α* increases, the share of large paths decreases. This leads to the appearance of multiply scattered particles, which form a hump in the central part of the distributions. As a result, the *W*—like asymptotic distribution of particles is formed. These solutions are not new and were obtained earlier when considering random walks with constant velocity [49,52,53], and when considering similar models of random walks [37,42,48,59,60] (see also the overview work [50]). In the case of the finite mathematical expectation (1 < *α* < 2), the process of random walks of a particle with a constant velocity falls under the action of the generalized central limit theorem. As a result, the asymptotic distribution of the particle coordinate is described by a stable law. In this case, the influence of the final velocity is reduced to a decrease in the diffusion coefficient *D* → *Dα* [34].

As we can see, in the case of the finite mathematical expectation and infinite mathematical expectation, not only are the asymptotic distributions different, but so are the equations describing the process of random walks. In the case of the finite mathematical expectation, the random walk process is described by the anomalous diffusion equation with the first time derivative (29), and the influence of the finite velocity is reduced to replacing the diffusion coefficient in this equation. In the case of the infinite mathematical expectation for taking account of the finite velocity of motion, as noted in the work [34], it is no longer sufficient to simply replace the diffusion coefficient. For this, it is necessary to replace the operator of the fractional Laplacian with the material derivative of the fractional order, which was first noted in the work [39]. Later, the authors of the works [40–44] come to the conclusion like this. Thus, in the case 0 < *α* < 1 in asymptotics *t* → ∞ the random walk process is described by Equation (25).

#### **4. Numerical Solution to a Kinetic Equation**

It should be noted that the numerical methods for solving Equation (29) are well investigated (see, for example, [8,12]). For a more detailed familiarization with the methods of numerical solutions of equations with fractional derivatives we refer the reader to the reviews [61,62]. However, there are no numerical methods for solving Equation (25) or similar equations with the material derivative of the fractional order. At least, the authors are not familiar with works devoted to methods of the numerical solution of such equations.

We consider a method for a numerical solution to Equations (25) and (29), which is based on the Monte Carlo method. From Section 2, it is clear that the random walk process of a particle with a constant velocity is described by an integral transport equation, which in the one-dimensional case takes the form (20). This allows the use of the Monte Carlo method to find a solution to this equation. The advantage of Monte Carlo methods is that they allow one to find a solution in multidimensional problems, as well as for various boundary and initial conditions.

From Section 2, the simplest method of numerical solution immediately follows based on modeling trajectories and histogram estimates of the particle distribution density. Each trajectory starts at a moment in time *t* = 0 from the origin of coordinates *x* = 0. We will consider one-dimensional symmetric walks *ω*1 = *ω*2 = 1/2. Therefore, from the origin of coordinates, a particle with equal probability moves to the right or to the left at a random distance *R*1, spending time *T*1 = *R*1/*v*0 on it. After that, the direction of motion is again modeled and with equal probability the particle continues to move to the right or to the left, moving a random distance *R*2, spending time *T*2 = *R*2/*v*0 on it. After that, the process continues in the same way. The construction of the trajectory continues as long as the condition is met: *N*(*T*∗)

$$\sum\_{k=1}^{T(T^\*)} T\_k \lessapprox T^\* \,. \tag{31}$$

where *T*∗ is a given moment in time at which it is necessary to find a solution. As soon as this condition is no longer met, the trajectory is terminated and a new trajectory begins. Random paths *Rk*, *k* = 1, 2, . . . are distributed with a density,

$$p(\mathbf{x}) = \begin{cases} \ \ \mathfrak{a}\mathfrak{x}\_0^{\mathfrak{a}}\mathfrak{x}^{-\mathfrak{a}-1} & \mathfrak{x} \ \geqslant \mathfrak{x}\_{0}, \\ 0 & \mathfrak{x} < \mathfrak{x}\_{0}. \end{cases}$$

Thus, random quantities *Rk*, *k* = 1, 2, ... are modeled according to the formula *R* = *x*0*ζ*−1/*α*, where *ζ* is a uniformly distributed random variable on a segment (0, 1].

To construct the simplest histogram estimate for the solution to the kinetic equation, the entire area of interest Δ = [*a*, *b*] is broken down into non-overlapping intervals Δ*i* = (*xi*, *xi*+1], *i* = 1, 2, ... *M* − 1, *x*1 = *a*, *xM* = *b*. To construct a histogram, the trajectory is modeled until the condition is met (31). As soon as this condition ceases to be met, the trajectory is terminated and the contribution from this trajectory is calculated:

$$h\_{\vec{j}}(\Delta\_{\vec{i}}) = \frac{I(\Delta\_{\vec{i}})}{\Delta\_{\vec{i}}},$$

where *I*(Δ*i*) is the interval indicator Δ*i*,

$$I(\Delta\_i) = \begin{cases} \mathbf{1}\_{\prime} & \mathbf{x}^\* \in \Delta\_{i\prime} \\ \mathbf{0}\_{\prime} & \mathbf{x}^\* \notin \Delta\_{i\prime} \end{cases}$$

where *x*∗ is the coordinate of a particle at the moment of time *T*∗. As a result, the estimate of the density for the interval Δ*i* is given by the expression:

$$\tilde{\rho}(\Delta\_i, t) = \frac{1}{N} \sum\_{j=1}^{N} h\_j(\Delta\_i)\_{\prime} \tag{32}$$

where summation is performed over the ensemble *N* of independent trajectories.

Despite the simplicity of this estimate, it has several disadvantages. Firstly, the estimate of the solution *ρ*˜(Δ*i*, *t*) is sought for the interval Δ*i*. This is the source of the systematic (horizontal) component of the error *δx*. Secondly, this estimate also contains the statistical component of the error ˆ *δ*, which decreases like *N*−1/2 at *N* → ∞. It is impossible to eliminate these errors completely; one can only reduce their value. However, a decrease in one of these values leads to an increase in the other value or to an increase in the calculation time.

It is possible to get rid of the systematic component of the error completely *δx*, if to consider one of the varieties of a local estimate. As in the case of the histogram estimate, the problem is to estimate the probability density of detecting a particle at the point *x*∗ at the moment of time *T*∗. The main element of solving the problem of transport theory by the Monte Carlo method—trajectory modeling remains unchanged. The difference lies in

the estimation method. In the case of a local estimate, the probability of a particle reaching a point is calculated (*x*∗, *T*∗), provided that the next collision is the final one.

To find this probability, let us consider the process of a random walk. Suppose that the density needs to be estimated at the point *x*∗ at the moment of time *T*∗ (see Figure 1). We will consider symmetrical random walks *ω*1 = *ω*2 = 1/2. The particle trajectory begins at the point *x* = 0 at the moment of time *t* = 0. Let the particle, as a result of a random walk at the moment of time *t* , reach the point *x* . To estimate the density at the point *x*∗, it is necessary to calculate the probability of transition from the point (*x* , *t* ) to the point (*x*∗, *T*∗). It is obvious that, to get from the point (*x* , *t* ) to the point (*x*∗, *T*∗), there are two possible alternatives of the trajectory continuation (see, Figure 1). The first option is to cover the distance *ξ*1 in the positive direction of the axis *Ox*. Then, to scatter at the point (*x* 1, *t* 1) and to cover the distance greater than *ξ*2 in the negative direction of the axis *Ox*. The second option is to cover the distance *ξ*2 first in the negative direction of the axis *Ox*. Then, to scatter at the point (*x* 2, *t* 2) and to cover the distance greater than *ξ*1 in the positive direction of the axis *Ox*. Since the particle is moving with a constant velocity *v*0 all the time, then the paths *ξ*1 and *ξ*2 will take time *τ*1 = *ξ*1/*v*0 and *τ*2 = *ξ*2/*v*0. The values of the paths *ξ*1 and *ξ*2 are determined by the formulas:

**Figure 1.** Possible continuation of the particle trajectory for constructing a local estimate.

Thus, the probability of transition from the point (*x* , *t* ) to the point (*x*∗, *T*∗) is determined by the formula:

$$
\Psi(\mathbf{x}^\*, \mathbf{x}', T^\*, t') = \frac{1}{4} (p(\mathfrak{J}\_1) P(\mathfrak{J}\_2) + p(\mathfrak{J}\_2) P(\mathfrak{J}\_1)),
\tag{33}
$$

)).

where *P*(*ξ*) = & ∞ *ξ p*(*ξ* )*dξ* , and the multiplier 1/4 appears as a result of the fact that, on the considered section of the trajectory, the particle changes its direction of motion twice. This transition probability is calculated after each particle scattering.

Theoretically, the contribution to the sought density can be made at the point *x*∗ by those particles that can get to the point *x*∗ from the current point without having a single scatter in the remaining time. Suppose that, as a result of random walks, the particle is scattered at the point *x* 1 at the moment of time *t* 1 (see Figure 1). Thus, staying at the point *x* 1 at the moment of time *t* 1, the particle can be found at the point *x*∗ at the moment of time *T*∗ if its path *ξ* satisfies the inequality *ξ* |*x*∗ − *x* 1|. In this case, the contribution to the density estimate has the form,

$$\psi\_{umc}(\mathfrak{x}^\*, \mathfrak{x}', T^\*, t') = \frac{1}{2}P(|\mathfrak{x}^\* - \mathfrak{x}'|),$$

where we need to substitute *x* = *x* 1, *t* = *t* 1. However, it should be noted that, for such a situation to be realized, during the random walk, the particle must undergo scattering at some point exactly lying on the segments *OA* or *OB*. Considering that the probability that the coordinate of a particle will take exactly a given value is equal to zero, then the contribution to the sought density from such particles will be zero. Nevertheless, unscattered particles make a significant contribution to the points lying on the half lines *x* = ±*vt*. These half lines determine the front of the distribution. Beyond this front, the probability of detecting a particle is zero. In view of the fact that the source is at the point *x* = 0, then it is obvious that the probability to detect a particle at the point *x* = 0 at the moment of time *t* = 0 is equal to one. Therefore, if after its appearance, the particle did not have any scattering and its path *ξ v*0*T*∗, then such particles reach the points *x* = ±*v*0*T*∗ and form the front of the distribution. The contribution of such particles at the points *x* = ±*v*0*T*∗ is calculated by the formula:

$$
\psi\_{unsc} \left( \pm \upsilon\_0 T^\* , 0 , T^\* , 0 \right) = \frac{1}{2} P(\upsilon\_0 T^\*) .
$$

The particles will also form the front of the distribution which, after their appearance during scattering, did not change their original direction of motion. Density contribution at points *x* = ±*v*0*T*∗ from such particles after each collision is calculated by a similar formula:

$$
\psi(\pm v\_0 T^\*, \mathbf{x'}, T^\*, t') = \frac{1}{2} P(|\!\!\!\! \pm v\_0 T^\* - \mathbf{x'}|). \tag{34}
$$

Here, *x* , *t* is the coordinate and moment of time at which the particle scattered. It should be noted that if, after its appearance, the particle began its motion in the positive direction, then unscattered or multiply scattered in one direction particles contribute only at the point *v*0*T*∗. Particles that began their motion in the negative direction—to the point −*v*0*T*∗. As a result, the contribution to the density estimate from each individual trajectory has the form:

$$h\_j(\mathbf{x}^\*, T^\*) = \sum\_{i=1}^{K(T^\*)} \Psi(\mathbf{x}^\*, \mathbf{x}\_{i\prime} T^\*, t\_i)\_{\prime\prime}$$

where *K*(*T*∗) is the number of scatterings that occurred during the time interval (0, *T*∗). To estimate the density at the points *x* = ±*v*0*T*∗, the function *ψ*(*x*∗, *xi*, *T*∗, *ti*) has the form (33), to estimate the density at the points *x* = ±*v*0*T*∗ the function *ψ*(*x*∗, *xi*, *T*∗, *ti*) takes the form (34). We finally obtain that the density at the point *x*∗ at the moment of time *T*∗, is estimated by the formula:

$$\bar{\rho}(\mathbf{x}^\*, T^\*) = \frac{1}{N} \sum\_{j=1}^N h\_j(\mathbf{x}^\*, T^\*)\_{\prime} \tag{35}$$

where the summation is performed over an ensemble of *N* independent trajectories.

In Figure 2, the results of the numerical solution to the Equation (25) are given for the values *α* = 0.5, *v* = 1, *T*∗ = 1000. Circles show the results of the histogram's estimation of the Monte Carlo method (32), asterisks are the results of the local estimate (35) and the solid line is the solution (26). The calculation results are transformed for the variable −1 *ζ* - 1 with the help of transformation *ρ*(*ζ*, *t*) = *vtρ*(*x*/*vt*, *t*), where −*vt x vt*. As we can see from the figure, at the value *α* = 0.5, the results of all three solutions coincide. The figure also shows the advantages of the considered solution method. The results of the local estimate (35) do not contain the horizontal component of the error, which is present in the histogram estimate of the Monte Carlo method (32) and is connected with a finite quantity of the interval Δ*i*. In view of the fact that the contribution from one trajectory to the point *x*∗ is calculated *K*(*T*∗) times, then this leads to a decrease in the statistical (vertical) error. In the figure, the magnitude of the statistical error does not exceed the size of the symbol. It should be noted that the solution (26) is an asymptotic solution and describes the distribution of particles at *t* → ∞. As we can see from the figure, in the case *α* = 0.5 and time *T*∗ = 1000, the results of all three methods of solution coincide. This means that

the process of the random walk of a particle has already entered the asymptotic regime. However, the moment when the random walk process enters the asymptotic regime for different values *α* is different.

**Figure 2.** Particle distribution density at *α* = 0.5, *v* = 1.0, *t* = 1000. Number of trajectories *N* = 105.

Figure 3 contains the solution results of a kinetic Equation (20) using the method of the histogram estimation (circles), the method of the local estimation (asterisks). Solution results are given for three times *T*∗ = 102, 103, 104. The figure also shows the solution to the Equation (25) (solid curve), which describes the asymptotic distribution of particles. As we can see in the figure, for the value *α* = 0.8 at the indicated times, the random walk process does not yet reach the asymptotic regime. However, when time *T*∗ increases, the solutions to the kinetic equation gradually approach the asymptotic distribution. It can also be seen from the figure that, at times *T*∗ = 103 and *T*∗ = 104, the solutions to the histogram estimate of the Monte Carlo method and the local estimate coincide, which confirms the validity of the results of the local estimation method. It can be seen from this figure that, at large times *T*∗ = 103 and *T*∗ = 104, there is an increase in the magnitude of the statistical error in the results of the local estimation. This increase in the calculation error is due to the fact that, at greater times, the value of the contribution (33) turns out to be a small magnitude. As a result, the calculation error can already affect it.

**Figure 3.** The density of particle distribution for the value *α* = 0.8 and specified values *T*∗. Number of trajectories *N* = 106.

The solution to the kinetic equation for the case *α* > 1 and specified values *T*∗ is given in Figure 4. In this figure, the points are the results of solving the kinetic Equation (20) with the use of the local estimation method, the solid curve is the solution (30). As in the previous cases, the calculation results are given for the variable *ζ* = *x*/*vt*, where −*vt x vt*. The figure shows that, at times *T*∗ = 100 and *T*∗ = 1000, the random walk process has not yet reached the asymptotic regime. The asymptotic regime of random walks is reached at *T*∗ = 104. As one can see, for a given time, the asymptotic solution (30) and the local estimation results are in good agreement with each other.

**Figure 4.** Density of the particle distribution for the value *α* = 1.3 and specified values *T*∗. Number of trajectories *N* = 105.

One obvious shortcoming of the local estimation method should also be pointed out. In the considered method, it is necessary to calculate the contribution (33) after each scattering of a particle. This leads to an increase in the computational operations of the processor, which in turn leads to an increase in the calculation time in comparison with the histogram method for estimating the density. However, this increase in calculation time is a necessary price to pay for the clear advantages of the local estimation method. As noted earlier, the main benefit of the local estimation (35) is the possibility to find a solution at a given point *x*∗. This means that the results of the local estimation do not contain a systematic (horizontal) component of the error. If we give a set of points *x*∗ *i* , *i* = 1, ... , *n*, then one trajectory will contribute at once to all points of the set *x*∗ *i* . This allows the solution to be constructed as a smooth function of the coordinate *x*. If we also take into account that the contribution (33) is calculated after each particle scattering, then the contribution from one trajectory is calculated *K*(*T*∗) times, where *K*(*T*∗) is the number of scatterings of a particle in the interval [0, *T*∗]. This leads to a decrease in the statistical component of the error.

#### **5. Conclusions**

The use of the theory of anomalous diffusion to describe combustion processes is only at the initial stage of development. At present, there are a few works devoted to the studies undertaken in this area of research [3,6–8,10,12]. In all the works, to describe combustion processes, the use of the equations of anomalous diffusion (1) or (2) is proposed. As is known, these equations describe the asymptotic (*t* → ∞, *x* → ∞) distribution of particles in the CTRW process, which is based on the assumption of the instant travel of particles from one point of space to another. This non-physical behavior of the particle leads to the fact that, at a time instant that is arbitrarily close to the initial one, the particle can be at an

arbitrarily large distance from the source. Therefore, it is necessary to use these equations to describe the processes occurring in a limited area of space and develop them in time with a certain degree of caution.

An obvious way to eliminate the instantaneous movement of a particle is to introduce into consideration a constant final velocity of the movement of a particle, which was done in this work. Taking account of the constant velocity of motion shows that, depending on whether the mathematical expectation of the travel is finite or infinite, the asymptotic distribution of particles is described by completely different equations. In the case of the finite mathematical expectation of the path value (1 < *α* < 2), the asymptotic process is described by the anomalous diffusion Equation (29), and the consideration of the finite velocity is reduced to the substitution of the anomalous diffusion coefficient *D* → *Dα*. In the case of the infinite mathematical expectation (0 < *α* < 1), the consideration of the finite velocity leads to a completely different Equation (25) containing the fractional-order material derivative operator.

The main difficulty in using equations in fractional derivatives is to find solutions to these equations. Analytical methods for solving these equations are only at the stage of development. Therefore, the main method of solution is to apply numerical methods. In this paper, a numerical solution method is considered which is based on a local estimate of the Monte Carlo method. This method is based on the idea of modeling the trajectory of a particle's random walk. The idea of the proposed method consists of the following. Now, after each scattering of a particle, one should calculate the probability of transition from the scattering point to a given point at which it is necessary to estimate the particle density. By giving a set of points, the transition probability should be calculated for each point from the given set. The advantages of this method over the standard histogram estimation of the Monte Carlo method are obvious. Firstly, since the solution is estimated at a specified point, then the estimation results do not contain a systematic (horizontal) component of the error. Secondly, each trajectory will contribute to all points of the given set at once. Taking account of the fact that, before the termination of the simulation of the trajectory, the particle undergoes *K*(*T*∗) of scatterings, then one trajectory will contribute *K*(*T*∗) times to each point of the given set. This leads to a decrease in the statistical error. However, since now the transition probability is calculated after each particle scattering, this leads to an increase in the arithmetic operations of the processor. As a result, in comparison with the histogram estimate of the Monte Carlo method, more time is spent on modeling the same number of trajectories. However, an increase in the calculation time should be considered a necessary payment for the complete absence of horizontal error.

The calculations made show that the results of the local estimation of the Monte Carlo method are in good agreement with both the results of the histogram estimation of the Monte Carlo method and the results of solving Equations (25) and (29), describing the asymptotic behavior of the process. From the calculations presented it also follows that, at different values of the parameter *α*, the process of random walks reaches the asymptotic regime at different times. This indicates another advantage of using the Monte Carlo method. In fact, Equations (25) and (29) describe the asymptotic behavior of the process at *t* → ∞, while the Monte Carlo method allows us to find a solution to the kinetic equation at any arbitrary moment in time and, thus, trace the evolution of the distribution of particles in time.

**Author Contributions:** Conceptualization, V.V.S.; methodology, V.V.S.; software, V.V.S.; validation, V.V.S. and Y.E.C.; investigation, V.V.S.; writing—original draft preparation, V.V.S.; writing—review and editing, V.V.S., V.N.K. and R.V.F.; visualization, V.V.S. and Y.E.C.; supervision, V.N.K.; project administration, R.V.F.; funding acquisition, V.N.K. and R.V.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Megagrant of the Government of the Russian Federation in the framework of the federal project No. 075-15-2021-584.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**

## *Article* **Application of a Pattern-Recognition Neural Network for Detecting Analog Electronic Circuit Faults**

**M. Isabel Dieste-Velasco**

Electromechanical Engineering Department, Higher Polytechnic School, University of Burgos, 09006 Burgos, Spain<EMAIL_ADDRESS>
**Abstract:** In this study, machine learning techniques based on the development of a pattern– recognition neural network were used for fault diagnosis in an analog electronic circuit to detect the individual hard faults (open circuits and short circuits) that may arise in a circuit. The ability to determine faults in the circuit was analyzed through the availability of a small number of measurements in the circuit, as test points are generally not accessible for verifying the behavior of all the components of an electronic circuit. It was shown that, despite the existence of a small number of measurements in the circuit that characterize the existing faults, the network based on pattern-recognition functioned adequately for the detection and classification of the hard faults. In addition, once the neural network has been trained, it can be used to analyze the behavior of the circuit versus variations in its components, with a wider range than that used to develop the neural network, in order to analyze the ability of the ANN to predict situations different from those used to train the ANN and to extract valuable information that may explain the behavior of the circuit.

**Keywords:** modeling; analog circuits; fault diagnosis; neural networks

**Citation:** Dieste-Velasco, M.I. Application of a Pattern-Recognition Neural Network for Detecting Analog Electronic Circuit Faults. *Mathematics* **2021**, *9*, 3247. https:// doi.org/10.3390/math9243247

Academic Editors: Camelia Petrescu and Valeriu David

Received: 3 November 2021 Accepted: 12 December 2021 Published: 15 December 2021

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

**Copyright:** © 2021 by the author. 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/).

#### **1. Introduction**

In analog electronic circuits, the limited access to measurement points makes determining faulty components a very complex task. On the other hand, when defining a set of measurement variables to characterize faults, many of the states that are generated by faults in the circuit are equivalent from the point of view of the values of the measured inputs, because the test points are generally not accessible to verify the behavior of all the components of the electronic circuit. In addition, performing measurements in each component of the circuit is not feasible from a practical point of view.

The present study deals with an application of supervised learning, based on the use of a pattern-recognition artificial neural network (ANN), for the detection of the individual hard faults (open circuits and short circuits) that may arise in an analog electronic circuit. The fact that the test points cannot be placed at all locations may cause several equivalent states to exist, depending on the points chosen to monitor the behavior of the circuit. This makes the detection of existing faults in an analog circuit a very complex task and much less developed than the same task in digital electronic circuits. In order to detect the hard faults that may arise in an electronic circuit, measurements are to be taken at accessible points in the circuit. Specifically, for the analysis to be carried out in this study, measurements of DC voltage and voltage gain were considered as input values so that it was possible to monitor the circuit and to determine, from these easily obtainable measures, whether the circuit was in a hard fault situation (open circuit or short circuit).

The first circuit under test (CUT) used in the present study was a single-stage smallsignal BJT amplifier, in which it is difficult to detect the hard faults that may arise because some faults lead to an equivalent state, from the point of view of the inputs used to monitor the behavior of the circuit, and later a more complex CUT was also studied. First, in the present study, the outputs of the CUTs versus variations that may arise from the tolerances of the passive elements of the circuit were obtained through a Monte Carlo analysis by

using Cadence® OrCAD® (CA, USA) design electronic simulation software. The values thus obtained were then used to train the ANN applied to predict the faulty components of the circuit. Moreover, a dataset obtained from the simulation software was used to validate and test the obtained results. In addition, once the pattern-classification neural network had been obtained, it was used to predict the behavior of the circuit subject to variations in the faulty components at wider ranges than those used to develop the neural network. This was carried out to determine the ranges of the parameters from which it is possible to detect hard faults in the CUTs.

Nowadays, determining faults in analog electronic circuits is being deeply studied by several research studies. For example, as shown in a review of Binu and Kariyappa [1], fault diagnosis in electronic circuits has been extensively researched in the last few years, for which machine learning approaches have been widely applied for fault detection. As shown by Binu and Kariyappa [1], open circuits and short circuits are some of the main failure sources in analog electronic circuits, and these hard faults can be modeled by including a 1 Ω parallel resistance with the component in a short-circuit situation and a 1 MΩ series resistance with the component in an open-circuit situation.

As previously mentioned, fault diagnosis in analog electronic circuits is a very complex task and is much less developed than the equivalent task in digital electronic circuits. The methods for analyzing faults in analog electronic circuits may be classified, roughly speaking, into two main categories: simulation before test (SBT) and simulation after test (SAT), as shown in the research study of Aizenberg et al. [2].

In the SBT approach, the development of a fault dictionary is very useful for detecting the faults in a circuit. In that way, the main faults that may arise in the circuit are simulated along with the nominal behavior of the circuit. In addition, in order to detect the faults that can occur in the analog circuit, it is important to consider both ambiguity groups, that is, the set of components of the electronic circuit that do not provide a unique solution if considered as a potential fault, and the canonical ambiguity groups, where a canonical ambiguity group is a group that does not contain other ambiguity groups [2–4], because it is very difficult to determine which component is faulty within one of these ambiguity groups.

Over the last few years, soft computing techniques for modeling and analyzing the behavior of electronic devices, as well as other kind of devices, have been widely used. As a consequence, several research studies dealing with this subject have been developed, as can be observed, for example, in [5–9], among many others research studies. With regard to the application of ANNs for detecting faults in analog electronic circuits, the study of Gao et al. [10] could be mentioned, where a dual-input fault diagnosis model based on convolutional neural networks, gated recurrent unit networks, and a softmax classifier was proposed. Likewise, Zhang et al. [11] used a convolutional neuronal network and backward difference for soft fault diagnosis in analog circuits, where the circuits being tested were the Sallen–Key band-pass filter and a four-opamp biquad high-pass filter.

Another studyworth mentioning is that of Wang et al. [12], which used a long shortterm memory neural network for fault detection and classification in modular multilevel converters in high-voltage direct current systems.

On the other hand, Xiao and Feng [13] used Monte Carlo analysis and SPICE simulation along with particle swarm optimization to tune the neural networks for analog fault diagnosis. Likewise, Aizenberg et al. [2] presented a method for detecting single parametric faults in analog circuits. They used a multi-valued neuron-based multilayer neural network (MLMVN) as a classifier, and a comparison with support vector machines (SVMs) was also presented in their study. These authors found that the MLMVN was highly accurate for classifying the fault class (FC) in the circuits under analysis in their study. Likewise, in the research of Kalpana et al. [14], Monte Carlo analysis was combined with machine learning techniques for fault diagnosis in analog circuits based on SBT.

Neural networks and genetic algorithms were also used in Tan et al. [15] for analog fault diagnosis, in which PSPICE simulations were used, and three different circuits were analyzed. These authors applied back propagation neural networks with 28–36 hidden layers, depending on the CUT, and with a binary coding scheme for the outputs, where the open-circuit faults were modeled with 1 × 106 times the nominal parameters and the short circuit as 1 × 10−6 the nominal values of each element. On the other hand, Viveros-Wacher et al. [16] used a CMOS RF negative feedback amplifier as the CUT for diagnosing faults using ANNs.

Some other studies on diagnosing analog circuit faults used neural networks and fuzzy logic, as shown by Bo et al. [17], who used a negative feedback amplifier as the CUT. Simulation and deep learning were also used by Pawlowski et al. [18] for identifying circuit faults in post-market circuit boards. In other studies, Li et al. [19] used a radial basis function (RFB) neural network and a back propagation algorithm for fault detection in a differential amplifier circuit. An RBF and back propagation were also used by Wuming and Peiliang [20], who employed a particle swarm optimization algorithm to adjust the neural network. SPICE and a quantum Hopfield neural network were employed by Li et al. [21] for fault analysis in a Sallen–Key band pass filter. Likewise, Monte Carlo analysis combined with deep learning and convolutional neuronal networks were used in Moezi and Kargar [22] for fault detection in analog circuits. In another study, Mosin [23] applied a three-layer feedforward neural network for fault diagnosis, where a tan-sigmoid function was used as the transfer function for the input and intermediate layers, and a log-sigmoid function was employed for the output layer, with a Sallen–Key bandpass filter being the CUT.

Further studies are that of Grasso et al. [24], which applied a procedure based on multifrequency fault diagnosis, where the CUT was a two-stage CE audio amplifier, and that of Li and Xie [25], which used a method based on the cross-entropy between a circuit under nominal behavior and one with faults, where the CUT was analyzed by Monte Carlo simulation. Some other studies are that of Sheikhan and Sha'bani [26], which used a modular neural model for soft fault diagnosis in analog circuits; that of Liang et al. [27], which applied a support vector machine classifier and fuzzy feature selection for analog circuit fault diagnosis; and that of Wang et al. [28], which used a semi-supervised algorithm for parametric fault diagnosis in analog circuits, among many others.

The remainder of this article is structured as follows: In Section 2, the methodology used to develop the ANN used to detect circuit faults is shown. In Section 3, the results are presented. A discussion of these results is provided in Section 4. Finally, the main conclusions of this study are outlined in Section 5.

#### **2. Fault Diagnosis Method**

As previously mentioned, this study analyzed the application of a pattern-classification ANN to detect hard faults in two analog circuits in which the faults that could arise were difficult to diagnose because several faults could provide similar results, from the point of view of the selected test points used to monitor the behavior of the circuit, because the test points should be selected in accessible points of the circuit and cannot simply be located anywhere due to practical considerations. Therefore, to detect the faults that may arise in the circuit, three measurements of DC voltage and the gain voltage were considered as input variables in the first CUT and six measurements of DC voltage and the gain voltage were considered as input variables in the second one. In order to develop the ANN used in this study, the software Cadence® OrCAD® Design Systems was first used in order to carry out a Monte Carlo analysis of the tolerances of the passive components of the circuit. The first CUT is shown in Figure 1, for which it is assumed that only the DC voltage in the transistor and the gain voltage are available. In addition, a simulation was first used to determine the failure situations that presented ambiguity because it was not possible to determine precisely which was the faulty component.

Figure 1 shows the first CUT used in the present study, which is a single-stage smallsignal BJT amplifier, similar to that shown in [29]. Likewise, Figure 1 shows the test points in this study. As can be observed, these test points are easily accessible. The nominal

values of the circuit's components are shown in Tables 1 and 2. First, the application of the pattern-recognition ANN to the CUT shown in Figure 1 is analyzed, and later, a more complex circuit that incorporates two amplification stages is analyzed in order to show that the ANN developed is capable of adequately predicting fault situations, as well as nominal behavior, in the second CUT.

**Figure 1.** Electrical diagram of the first CUT (a single-stage small-signal BJT amplifier).

**Table 1.** Nominal values and tolerances of the passive components of the circuit.

**Table 2.** Values of the load resistance and voltage sources of the circuit.

Table 3 shows the possible individual faults that may arise in the first CUT obtained when the hard faults (a short circuit (sc) or an open circuit (oc)) arise in the passive components. From Table 3, it is possible to see that, in this CUT, there are 14 individual hard faults, as well as the nominal behavior of the circuit {Nominal, R1oc, R1sc, R2oc, R2sc, R3oc, R3sc, R4oc,R4sc, C1oc, C1sc, C2oc, C2sc, C3oc, C3sc}, which were coded as {F01, F02, F03, F04, F05, F06, F07, F08, F09, F10, F11, F12, F13, F14, F15 }. Therefore, these were the working modes that were analyzed. As previously mentioned, in order to characterize the behavior of the CUT, an electronic simulation was carried out by using Cadence® OrCAD® design electronic simulation software for each of the failure modes shown in Table 3. As can be observed, the faults were grouped into ambiguity groups, from the point of view of the inputs considered to diagnose the circuit's behavior, where a hard fault in a component of the circuit due to an open circuit (oc) was simulated by placing a resistance (RFault = 10 MΩ) in series with the component, and a hard fault due to a short circuit (sc) was simulated by placing a resistance (RFault = 1 Ω) in parallel with the component. The ambiguity groups were determined from the values of the inputs, which were obtained from an electronic simulation. These ambiguity groups (Mj classes) were coded as {M01, M02, M03, M04, M05, M06, M07, M08, M09, M10, M11}. It should be mentioned that there are some fault events, such as those obtained, for example, in the M04 class, which include hard faults {F10, F12}, in the event of which it would not be possible to determine the faulty component. Moreover,

in case of a situation due to a catastrophic fault leading to an actual open circuit or a short circuit, the DC voltages and the gain voltage (Av) could be different from those obtained by the model employed in this study. These situations were obtained from the simulation when a 10 MΩ resistance was placed in series with the faulty component to simulate the open circuit (oc), and when a value of 1 Ω resistance was placed in parallel with the faulty component. Therefore, to consider the actual catastrophic fault, the values obtained in the test points were also obtained from the simulation and considered as additional inputs to those provided by the Monte Carlo analysis in order to train the ANN.

**Table 3.** Nominal behavior and hard faults grouped by ambiguity groups (Mj classes).

A Monte Carlo analysis considering the tolerances of the passive elements of the circuit shown in Table 1 was first carried out for each of the hard faults (open circuits and short circuits) in order to train the ANN, and 64 results were generated for each fault (63 results from the Monte Carlo analysis and 1 additional result from the actual catastrophic fault). Likewise, 64 results were obtained for the nominal behavior. These results were then used to train the pattern-recognition ANN considered in this study. Figure 2 shows the ANN applied, which was trained to detect the nominal behavior and the individual faults shown in Table 3. The hard faults that may arise in the circuit shown in Figure 1, as well as the nominal behavior, were characterized from the outputs of the ANN as shown in Equation (1):

$$S\_{\circ} = \operatorname{colum} m\_{\circ -th} \{ I \} \tag{1}$$

where *Sj* corresponds to the ANN outputs, so that the *j*-*th* output class corresponds to the *jth* column of the identity matrix (*I*). The nominal behavior corresponds to *M*01(*F*1) and the remaining classes shown in Table 3 correspond to the short-circuit and open-circuit faults, where the hard faults were grouped by the ambiguity groups obtained from the inputs used to characterize the behavior of the circuit. Therefore, the coding used to characterize a fault should provide a "1" in the position of the fault and "0" in the rest of the outputs, and hence, all outputs will have a "0" value except the *j-th* class (the fault class to be identified), which will have a "1" value. The same is applicable for the nominal value.

As can be observed, in Figure 2, the first ANN used in the present study was made up of an input layer that has four inputs (VB, VC, VE, AV), which correspond to the DC voltages in the base, collector, and emitter of the BJT transistor and to the gain voltage (AV), respectively, as well as a single hidden layer (with two neurons and a hyperbolic tangent as the transfer function) and one output layer with a softmax transfer function, which is commonly used in pattern-recognition neural networks. As can be noted, the output layer has 11 outputs, which correspond to the 10 fault classes identified in the electronic circuit and to the nominal working mode.

As shown later in this study, with the configuration given in Figure 2, it is possible to have high accuracy in the ANN for detecting both the hard faults of the circuit and the nominal behavior. It should be mentioned that different ANN topologies were analyzed with one and two hidden layers and by using different training algorithms to adjust the ANN parameters. Finally, a Levenberg–Marquardt back propagation algorithm was selected to update the weights and biases of the ANN by using the Deep Learning Toolbox™ of MATLABTM 2020a [30]. The ANN shown in Figure 2 was used, since it was able to provide accurate results without having to increase the number of neurons or the number of hidden layers. The metric used to test the models was the mean squared error (MSE). Different transfer functions were also analyzed in the hidden layer but, finally, a hyperbolic tangent was used in this study. On the other hand, the Levenberg–Marquardt algorithm was able to provide, in this case, more accurate results than the others analyzed, such as

the Scaled Conjugate Gradient. Therefore, the topology was that shown in Figure 2, where *W*1 and *b*1 are the weights and bias of the hidden layer, and *W*2 and *b*2 are those of the output layer. As previously mentioned, a hyperbolic tangent (*f*1) was used as the transfer function in the hidden layer and a softmax transfer function (*f*2) was used in the output layer.

$$Outputs = f\_2(\mathcal{W}\_2 \* f\_1(\mathcal{W}\_1 \* Inputs + b\_1) + b\_2) \tag{2}$$

As Figure 2 shows, the number of outputs was 11, where each output corresponds to the class identified (Mj); one of them represents the nominal behavior and the remaining classes representing the ambiguity groups, where the outputs of the ANN can be obtained from Equation (2). In order to obtain the results shown in this study, the Deep Learning ToolboxTM of MATLABTM R2020a [30] was used.

**Figure 2.** Pattern-recognition ANN applied in this study.

#### **3. Results**

After training the ANN shown in Figure 2 with the data obtained from the Monte Carlo simulations and following the procedure shown in the previous section, it was possible to obtain the confusion matrices shown in Figures 3–6, for training, validation, testing, and all data, respectively, where 70% of data were used for training, 15% for testing and 15% for the validation. As can be observed, in Figures 3–6, a perfect classification of the results was obtained with this ANN comprising a single hidden layer that contains two neurons.

**Figure 3.** Confusion matrix obtained with the ANN (training data).

**Figure 4.** Confusion matrix obtained with the ANN (validation data).

Figure 3 shows the results obtained in the confusion matrix when 70% of the data from the Monte Carlo analysis were employed to train the ANN. It can be seen that there are fault classes that present a larger amount of data due to the fact that they agglutinate fault configurations that belong to the same ambiguity group. As can be seen, 100% of the data are classified correctly.

Figure 4 shows the results obtained in the confusion matrix when 15% of the Monte Carlo data were used for validation of the ANN, and Figure 5 shows the results for the test case. Similar to the results obtained during training, the ANN was able to diagnose 100% of the working modes correctly (hard faults and nominal behavior).

**Figure 5.** Confusion matrix obtained with the ANN (test data).

**Figure 6.** Confusion matrix obtained with the ANN (all data).

Figure 6 shows the results of the confusion matrix for each of the hard faults and for the nominal behavior considering all the Monte Carlo data. As previously mentioned, 64 values were used for each fault and for the nominal behavior in the Monte Carlo analysis by considering the tolerances of the components.

It can be noted in Figure 6 that, in the confusion matrix generated from all data, there are classes with a greater number of elements because the ambiguous failure modes were grouped into failure classes. Thus, for example, class M03 has 192 elements, since it encompasses three failure modes (R1oc, R2sc, C1sc). In addition, it can be seen that 100% of the data were classified correctly.

Figure 7a shows the mean squared error (MSE) obtained with the ANN and Figure 7b shows the error histograms for training, validation, and testing.

**Figure 7.** (**a**) Validation performance (MSE) and (**b)** error histogram.

As Figure 8 shows, the ROC (receiver operating characteristic) curves have an area under the curve (AUC) of 1, which demonstrates that the pattern-recognition ANN developed was able to diagnose the working modes of the BJT amplifier once they were classified into the 11 classes shown in Table 3. In order to show that the pattern-recognition ANN is capable of predicting the behavior of other circuits, a two-stage small-signal BJT amplifier, such as that shown in Figure 9, was also analyzed in the present study.

**Figure 8.** ROC curves.

The nominal values of the circuit's components are shown in Tables 4–6.

**Table 4.** Nominal values and tolerances of the passive components of the circuit (R).

**Table 5.** Nominal values and tolerances of the passive components of the circuit (C).

**Table 6.** Values of the load resistance and the voltage sources of the circuit.

**Figure 9.** Electrical diagram of the second CUT (a two-stage small-signal BJT amplifier).

Table 7 shows the possible individual faults that may arise in the second CUT obtained when the hard faults (short circuit (sc) or open circuit (oc)) arise in the passive components. From Table 7, it is possible to see that, in this second CUT, there are 27 individual hard faults, as well as the nominal behavior of the circuit {nominal, C1oc, C1sc, C2oc, C2sc, C3oc, C3sc, C4oc, C4sc, C5oc, C5sc, R1oc, R1sc, R2oc, R2sc, R3oc, R3sc, R4oc, R4sc, R5oc, R5sc, R6oc, R6sc, R7oc, R7sc, R8oc, R8sc }, which were coded as {F01, F02, ... , F26, F27}, where F01 corresponds to the nominal behavior. Therefore, these are the working modes that are analyzed in the second case. As can be observed, the faults were grouped into ambiguity groups, from the point of view of the inputs considered to diagnose the circuit's behavior, following the previously mentioned procedure, where a hard fault in a component of the circuit due to an open circuit (oc) was simulated by placing a resistance (RFault = 10 MΩ) in series with the component, and a hard fault due to a short circuit (sc) was simulated by placing a resistance (RFault = 1 Ω) in parallel with the component.

**Table 7.** Nominal behavior and hard faults in the second CUT grouped by ambiguity groups (Mj classes).

As in the previous case, the ambiguity groups were determined from the values of the inputs, which were obtained from an electronic simulation. These ambiguity groups (Mj classes) were coded as {M01, M02, ... , M19, M20} because, in this second case, 20 classes were detected. It should be mentioned that there are some fault events, such as those obtained, for example, in the M03 class, that include hard faults {R1oc, R2sc, C1sc} for which it would not be possible to determine the faulty component. Moreover, in case of a situation due to a catastrophic fault leading to an actual open circuit or a short circuit, the DC voltages and the gain voltage (Av) could be different to those obtained by the model employed in this study. Therefore, to consider an actual catastrophic fault, as in the previous case, the values obtained in the test points were also obtained from the simulation and considered as additional inputs to those provided by the Monte Carlo analysis in order to train the ANN.

Figure 10 shows the ANN for the second CUT, which is shown in Figure 9. As can be seen in this case, the number of inputs is seven, which correspond to the voltages at the base, emitter and collector of both transistors as well as the gain voltage (VB1, VC1, VE1, VB2, VC2, VE2, AV), and the outputs are 20, corresponding to the detected fault classes and the nominal behavior. Similar to the previous case, the same network topology is used, although, in this case, there are four neurons in the hidden layer. As was done with the ANN developed for the first CUT, a Levenberg–Marquardt back propagation algorithm was selected to update the weights and biases of the ANN by using the Deep Learning Toolbox™ in MATLABTM 2020a [30].

**Figure 10.** Pattern-recognition ANN applied to the second CUT.

Figure 11a shows the MSE values obtained versus the number of epochs, and Figure 11b shows the training state values for the ANN employed to analyze the faults in the second CUT.

**Figure 11.** (**a**) Validation performance (MSE) and (**b**) training state of the ANN for the second CUT.

Figure 12 shows the results obtained in the confusion matrix when 70% of the data from the Monte Carlo analysis were employed to train the ANN shown in Figure 10, which was employed to model the behavior of the second CUT. It can be seen that there are fault classes that present a larger amount of data due to the fact that they agglutinate fault configurations that belong to the same ambiguity group. As can be seen, 100% of the data are classified correctly in the second case, similar to the previous one. As can be observed, the fault classes do not have the same number of elements because the data used to train, validate, and test the ANN were randomly selected.

**Figure 12.** Confusion matrix obtained with the ANN (training data) for the second CUT.

Figure 13 shows the results obtained in the confusion matrix when 15% of the Monte Carlo data were used for validation of the ANN, and Figure 14 shows the results for the test case. Similar to the results obtained during training, the ANN was able to diagnose 100% of the working modes correctly (hard faults and nominal behavior). As can be observed, the same results as those obtained with the first CUT were obtained with the second CUT.

**Figure 13.** Confusion matrix obtained with the ANN (validation data) for the second CUT.

**Figure 14.** Confusion matrix obtained with the ANN (test data) for the second CUT.

Finally, Figure 15 shows the confusion chart for all data, and Figure 16 shows the ROC curve (all data) for the second CUT. In this curve, the true positive rate (TPR) versus the false positive rate (FPR) was plotted at different threshold settings. The ANN developed in this study is a perfect classifier for the electronic faults in the second CUT because it is perfectly able to distinguish each fault class for any FPR. Similar to the results obtained with the first CUT, it can be seen that 100% of the data were classified correctly.

**Figure 15.** Confusion matrix obtained with the ANN (all data) for the second CUT.

**Figure 16.** ROC curve (all data) for the second CUT.

#### **4. Discussion**

As Figures 8 and 16 show, the pattern-recognition ANN can correctly diagnose both the nominal behavior and the fault classes of the CUTs considered in the present study, with 100% of the data correctly classified when considering 64 Monte Carlo points for each fault and for the nominal behavior. However, in order to analyze the ability of the ANN to explain situations different from those used to train the ANN and to extract valuable information that may explain the behavior of the circuit, wider ranges of the fault resistances placed in series and in parallel with the components to simulate the hard faults in the CUT were used. These values were chosen in order to generate different fault scenarios to determine the ability of the developed ANN to diagnose possible fault situations before a hard fault occurs. To test the ANN with these fault resistances, a new Monte Carlo analysis was performed. In this later case, the number of runs generated for the fault resistance was 1024, for each fault, instead of the 64 runs used to train the ANN, following a uniform distribution, as shown in Figure 17a, for the case of a fault resistance in series with the faulty component to simulate an open circuit, and in Figure 17b for the case of a fault resistance in parallel with the faulty component to simulate a short circuit. On the other hand, the rest of the components of the CUT were allowed to vary within the specified tolerances.

**Figure 17.** Histograms showing the 1024 values of the (**a**) series and (**b**) parallel fault resistances used in the Monte Carlo analysis (grouped into 25 bins).

That is, to train the neural network, an initial Monte Carlo analysis was performed for the CUT in which 64 runs were used for each fault event as well as for the nominal value. This results in a total of 960 input vectors (VB, VC, VE, AV) in the case of the first CUT and a total of 1728 input vectors (VB1, VC1, VE1, VB2, VC2, VE2, AV) in the case of the second CUT. Therefore, this study employed a supervised learning technique in order to develop a pattern-recognition neural network. It should be noted that, in the case of the first CUT, 70% of these 960 data, obtained by Monte Carlo analysis, was used to train the neural network (i.e., 672 data). The remaining 15% of the data was used for validation and the other 15% for testing. The same procedure was followed for the second CUT (70% train, 15% test, 15% validation). Once the neural network was developed, the values predicted by the network for the different modes of operation were analyzed. This first Monte Carlo analysis was generated from the tolerances of the circuit components, which were considered commercial and standardized values with tolerances of 10% for resistors and 20% for capacitors. As shown in the present study, the proposed ANN is a perfect classifier since it is able to discriminate 100% of the data, not only with those used for training, but also with those used for validation and testing, in both CUTs. This can be observed in the ROC curves shown in Figure 8 (for the first CUT) and Figure 16 (for the second CUT). Once the network was developed, another Monte Carlo analysis was carried out to analyze how the ANN is able to predict other fault events, where the resistances used were different from those used to develop the ANN. This was done by varying the fault resistances (which are placed in series and in parallel with the potentially faulty components) with values of 10 MΩ ± 99.9% to simulate the open circuit and values of 1 Ω ± 99.9% for the short circuit. These values of the resistors are shown in Figure 17 and were chosen in order to generate different fault scenarios to determine the ability of the ANN to diagnose possible fault situations before a hard fault occurs. In the latter case, the Monte Carlo analysis was carried out using 1024 values for each fault event. From this, it was possible to obtain the outputs of the ANN for these fault events and to determine the thresholds from which the fault will be detected in each component. In the case of the nominal behavior, it was also considered that the tolerance of the components was increased by 50% relative to the nominal values, as shown in Figure 18, so that, in this case, the resistance tolerances were increased to 15% and up to a value of 30% in the case of the capacitors.

**Figure 18.** Histograms showing the 1024 values used in the Monte Carlo analysis (grouped into 25 bins) for each component in the case of the nominal behavior.

Figure 19 shows the results of the confusion matrix obtained in the case of a wider range of variation in the fault resistance, for the first CUT. As can be noted, the ANN was

able to correctly diagnose most of the fault classes that may arise in the first CUT, as well as the nominal behavior. A similar analysis could be carried out with the second CUT.

**Figure 19.** Predicted faults in the first CUT when the fault resistances had 99.9% variation.

Specifically, it can be observed in Figure 19 that, when the first CUT works with the nominal values of the passive components, with their tolerances increased by 50%, the ANN predicted nominal behavior in all cases (100%), which is logical, since the BJT amplifier considered as the first CUT in this study was robust to variations in the tolerances of the passive components, so it was not greatly affected by the fact that these tolerances were increased by 50% with respect to the design values, as can be observed in Figure 20.

On the other hand, regarding the M02 and M03 classes, these were correctly diagnosed. In the case of M04, there were some faults that were classified as M02 and M05 classes, which, at first, may seem like a detection failure by the ANN, but may have actually been

caused by the fact that varying the resistance in series with C1 and C2 within the range of values analyzed (by setting 99.9% variation in the series fault resistance) can lead to a similar configuration from the point of view of the DC voltages of the transistor. In any case, 93.6% of the cases analyzed were correctly detected. Additionally, in the case of M05, M06, and M07, 100% of the cases were detected correctly. Moreover, regarding the M08 class the network predicted 94.6% of the faults. For the rest of the classes (M09, M10, and M11), the ANN detected 100% of the faults.

Therefore, the ANN developed in this study could accurately predict the behavior of the first CUT when faced with variations in the fault resistance. Figures 21–26 show the values predicted by the ANN versus the values of the fault resistance.

**Figure 21.** Predicted values (ANN) vs. fault resistance (R) in parallel with C2.

Figure 21 shows the values predicted by the ANN for the fault class M02. It can be noted that the ANN detected all the faults in the circuit for the values of parallel resistance considered.

Figure 22 shows the values predicted by the ANN for fault class M03. It can be noted that the ANN detected all the faults in the circuit for the values of serial and parallel resistances (R) considered. Likewise, Figure 23a shows that, in the case of C1 for low values of resistance in series (R) with the faulty component, some of these situations could be detected as M02{C2sc} and M05{C3oc} since the values of the fault resistance in series with C1 presented a minimum value of 64 kΩ, which was obtained in this study through Monte Carlo analysis with 1024 runs. The same behavior was obtained in the case of C2, although for different thresholds of resistance (R), as can be observed in Figure 23b.

**Figure 22.** Predicted values (ANN) vs. fault resistance (R) (**a**) in series with R1 and (**b**) in parallel with R2 or C1.

**Figure 23.** Predicted values (ANN) vs. fault resistance (R) (**a**) in series with C1 and (**b**) in series with C2.

Figure 24a shows the values predicted by the ANN for the hard faults of class M05 {C3oc}, and Figure 24b shows those for the M06 {R4sc, C3sc} and M07 {R3sc} fault classes. As can be observed, the ANN detected all the faults in the circuit.

Figure 25 shows the results predicted by the ANN for the M08 fault class (94.6% faults were detected). Finally, Figure 26 shows the results predicted for the remaining fault classes. As can be observed, 100% of fault data were correctly diagnosed in the case of R1sc (M09), R2oc (M10), and R3oc (M11).

**Figure 24.** Predicted values (ANN) vs. fault resistance (R) (**a**) in series with C3 and (**b**) in parallel with R4, C3, or R3.

**Figure 25.** Predicted values (ANN) vs. fault resistance (R) in series with R4.

**Figure 26.** Predicted values (ANN) vs. fault resistance (R) (**a**) in parallel with R1 and (**b**) in series with R2 or R3.

#### **5. Conclusions**

In the present study, a pattern-recognition neural network with a hyperbolic tangent as the transfer function in the hidden layer, along with a softmax transfer function in the output layer, was used to diagnose individual hard faults in two CUTs. First, a single-stage small-signal BJT amplifier was studied, followed by a two-stage small-signal BJT amplifier. It was shown that the ANN was able to predict the hard faults accurately in both CUTs considered in this work.

It was shown that a pattern-recognition ANN such as the one considered in this study can be used to model hard faults in the CUTs by training the ANN with a reduced number of measurements that have been taken at accessible points of the circuit and by using only one hidden layer with a reduced number of neurons.

Moreover, in case of applying the neural network to situations different from those used to train the neural network, where these situations have been modeled by expanding the ranges of failure resistance, it was shown that the ANN developed had high precision in diagnosing the failures in the first CUT, and it was able to explain situations different from those used to train the ANN and to extract valuable information that may explain the behavior of the circuit.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**

## *Article* **Rotating Flow in a Nanofluid with CNT Nanoparticles over a Stretching/Shrinking Surface**

**Nor Azizah Yacob 1, Nor Fadhilah Dzulkifli 1, Siti Nur Alwani Salleh 2, Anuar Ishak 3,\* and Ioan Pop 4**

	- Merbok 08400, Kedah, Malaysia<EMAIL_ADDRESS>
**Abstract:** The steady three-dimensional rotating flow past a stretching/shrinking surface in water and kerosene-based nanofluids containing single and multi-walled carbon nanotubes (CNTs) is investigated. The governing equations are converted to similarity equations, and then numerically solved using MATLAB software. The impacts of rotational, suction, and nanoparticle volume fraction on the flow and the thermal fields, as well as velocity and temperature gradients at the surface, are represented graphically and are analyzed. Further, the friction factor and the heat transfer rate for different parameters are presented in tables. It is found that the heat transfer rate increases with increasing nanoparticle volume fraction as well as suction parameter in water and kerosenebased nanofluids of single and multi-walled CNTs. However, the increment in the rotating flow parameter decreases the rate of heat transfer. Multi-walled carbon nanotubes and kerosene-based nanofluid contribute to heat transfer rates better than single-walled carbon nanotubes and waterbased nanofluid, respectively. A unique solution exists for the stretching surface, while two solutions are obtained for the shrinking surface. Further analysis of their stabilities shows that only one of them is stable over time.

**Keywords:** carbon nanotubes; heat transfer; nanofluid; rotating; stretching/shrinking

#### **1. Introduction**

Researchers have begun looking for fluids with high thermal conductivity in the last few decades since traditional fluids have low thermal conductivity. In 1873, Maxwell studied the dispersion of solid particles with a size of millimeters and micrometers in traditional fluids. However, the existence of these particles caused a drastic drop in pressure, clogged flow channels, caused sedimentation of the particles to occur faster and erosion on the channel, and certain particles are too large for the microsystem [1]. Choi and Eastman [1] proposed a new fluid called nanofluid, in which solid particles with sizes ranging from 1–100 nm are dispersed in traditional fluids such as water, oil, and ethylene glycol. This fluid is significant because of its unique chemical and physical properties, as well as its high thermal conductivity, which can improve the heat transfer rate [2]. Waqas et al. [3] investigated the impact of nanoparticle shape's factor in porous media and discovered that platelet and cylindrical shapes had the highest thermal conductivity compared to brick and spherical shapes. A comprehensive collection of nanofluids studies can be found in [4–9].

Recently, many studies have considered carbon nanotubes as nanoparticles since their thermal conductivity is better and higher than those of traditional fluids. Carbon nanotubes (CNTs) are a sheet of graphene rolled into a cylinder introduced by Iigima in 1991. Carbon nanotubes are classified as single-walled (SWCNT) and multi-walled (MWCNT). They have

**Citation:** Yacob, N.A.; Dzulkifli, N.F.; Salleh, S.N.A.; Ishak, A.; Pop, I. Rotating Flow in a Nanofluid with CNT Nanoparticles over a Stretching/Shrinking Surface. *Mathematics* **2022**, *10*, 7. https:// doi.org/10.3390/math10010007

Academic Editors: Camelia Petrescu and Valeriu David

Received: 17 November 2021 Accepted: 15 December 2021 Published: 21 December 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/).

a nanometer diameter size which is much smaller than other advanced semiconductor devices, and the length of nanotubes can exceed 1 micrometer. In addition, carbon nanotubes also offer unique electrical characteristics that other semiconductors lack. Further studies on unique features and review of carbon nanotubes can be found through investigation by Popov [10], Patel et al. [11], and Taherian et al. [12].

Due to their performance in optimizing thermal conductivity, researchers carried out many works to study the influence of carbon nanotubes availability and other effects in the fluids on the flow behavior. Khan et al. [13] investigated the fluid flow and heat transmission of CNTs nanofluid over a flat plate using three types of base fluids: water, engine oil, and kerosene. In comparison to water and kerosene-CNTs, they discovered that engine oil-CNTs have the highest heat transfer rate. Using a homogeneous model, Akbar et al. [14] investigated the stagnation-point flow of carbon nanotubes towards a stretching sheet with a magnetic field, velocity slip effect, and convective boundary condition. It was discovered that a higher magnetic field reduces the heat transfer rate at the surface and that single-walled carbon nanotubes have a higher heat transfer rate than multi-walled carbon nanotubes. Later, using similar CNTs nanofluid, Hayat et al. [15] explored the stagnationpoint flow over a nonlinear stretching surface with homogeneous–heterogeneous reactions. Hussain et al. [16] studied a three-dimensional flow of carbon nanotubes nanofluids over a nonlinearly stretching sheet with heat generation/absorption effect in a porous medium. Sreedevi et al. [17] investigated the heat and mass transfer characteristics of nanofluids containing CNTs towards a vertical cone in a porous medium by considering the magnetic field, chemical reaction, and suction/injection parameters. The findings revealed that increasing the nanoparticle volume fraction improves the heat transfer rate, while employing multi-walled CNT. as nanoparticles enhance the heat transfer over single-walled CNT. The latter conclusion was similar to Anuar et al. [18], who analyzed the mixed convection flow towards a moving plate. However, Anuar et al. [18] obtained a different result when the nanoparticle volume fraction increased. Reddy and Sreedevi [19] studied the effects of thermal radiation and single-walled carbon nanotubes inside a square chamber. From their observation, increasing these two factors enhances the heat transfer rate. Ramzan et al. [20] applied the Cattaneo–Christov heat flux model to analyze the impact of nonlinear thermal radiation and thermal stratification on the unsteady magnetohydrodynamics (MHD) flow between two stretching rotatory discs in nanofluid with carbon nanotubes. They discovered that nonlinear thermal radiation gives rise to the temperature of both SWCNT and MWCNT. Relevant studies on CNT nanofluids can be found in [21–23].

Investigation of the rotating flow towards a stretching surface is vital due to its application in a variety of scientific, technical, and product applications, designing and modeling jet engines, pumps, and vacuum cleaners, as well as geophysical flows [24]. Numerical studies on rotating flow have been explored by many researchers such as Singh and Sathi [25], who obtained an exact solution using Laplace transform for an unsteady rotating flow over an infinite plate. Wang [26] was the first to study the rotating flow over a stretching sheet using a perturbation method to identify the velocity profiles of small parameter values. Later, Nazar et al. [27] extended the work of Wang [26] by considering an unsteady flow. This study was expanded by Ali et al. [28] to the shrinking surface, and they obtained dual solutions. Further, Sreelakshmi et al. [29] included the magnetic field and thermal radiation effects in their study. They reported that the magnetic field, thermal radiation parameter, and rotation parameter increase the heat transfer rate. Rosali et al. [30] investigated the rotating flow over an exponentially permeable shrinking surface and found dual solutions. Rana et al. [31] implemented the Boungiorno model to investigate the MHD flow across a stretching surface in a nanofluid and numerically solved using a finite element method. It was revealed that the local Nusselt number and the local Sherwood number drop with the growth of the rotational parameter. By combining the Boungiorno and Tiwari-Das models, Bakar et al. [32] solved a rotating flow towards a shrinking surface with three nanoparticles types: copper, alumina, and titania with water as the base fluid using a shooting method. Different from Rana et al. [31], they found dual solutions, and the

local Nusselt number and the local Sherwood number increase with rotational parameter. Further, Krishna [33] explored the MHD effect in a Jeffrey fluid on a porous surface with hall and ion slip.

Nadeem et al. [34] explored a rotating flow towards a stretching surface utilizing copper and titania as the nanoparticles with water as the base fluid by applying the Tiwari-Das model. They discovered that increasing the volume fraction of nanoparticles boosted the heat transfer rate while decreasing it with the rotation parameter. Hayat et al. [15] reported the analytical results of rotating flow with carbon nanotubes over a stretching porosity surface. Dzulkifli et al. [35] continued the study of Ali et al. [28] by considering the unsteady rotation flow on a stretched surface in a nanofluid with copper as nanoparticles. According to the findings, as the rotating parameter and nanoparticle volume percentage grew, the heat transfer rate dropped. Later, Salleh et al. [36] expanded the study by Nadeem et al. [34] to a permeable shrinking surface and obtained dual solutions. A stability analysis was carried out to inspect the stability of the solutions in the long run. As stated by Merkin [37] and Harris et al. [38], the upper branch is always stable, and this discovery is similar to the works made by Salleh et al. [36], Mustafa et al. [39], and Tshivhi and Makinde [40]. Hafeez et al. [41] explored the MHD rotating flow in a hybrid nanofluid using a finite element method. Mehdipour et al. [42] designed a rotating liquid sheet (RLS) contactor to explore carbon dioxide absorption from a gas stream using water-based SiO2 and ZnO nanofluids. CO2 absorption is more significant in the presence of ZnO than in the presence of SiO2, and tube rotation can improve absorption performance.

Acharya et al. [43] conducted a comparative investigation of SWCNTs and MWCNTs in a rotating regime over a stretching sheet in the presence of a magnetic field using the RK-4 method. They noticed that increasing the volume fraction of SWCNTs and MWCNTs improved the nanofluid temperature. Later, Shah et al. [44] analyzed a similar investigation over a stretching sheet using an analytical method. Noranuar et al. [45] studied a noncoaxial rotation flow towards a permeable moving disk in an MHD Casson nanofluid. Hussain et al. [46] extended the work of Nadeem et al. [34] by considering engine oil as the base fluid with copper and aluminum dioxide as nanoparticles over a stretching surface with slip condition. They found that copper nanoparticles react as heat carriers better than aluminum oxide. Manjunatha et al. [47] recently performed a numerical analysis to study the effect of a magnetic field on the fluid flow over a rotating disc in a nanofluid using single and multi-walled carbon nanotubes as nanoparticles and water as the base fluid. They reported that a unique solution exists, and the increase in the magnetic field caused an increase in the radial velocity and temperature.

Based on the previous studies mentioned above, many researchers focused on a rotational flow/disk over a stretching surface rather than a permeable stretching/shrinking surface and only obtained a unique solution. A permeable stretching/shrinking surface is vital in industrial applications since the surface may not be static such as polymer extrusion process, liquid crystal, and ceramic production [48]. Therefore, the present study aims to extend the work of Nadeem et al. [34] and Salleh et al. [36] by exploring the impact of rotation flow parameter, suction and nanoparticle volume fraction on the heat and mass transfer characteristics of carbon nanotubes (SWCNTs and MWCNTs) over a permeable stretching/shrinking surface in a nanofluid. Since no studies have considered this mathematical model with carbon nanotubes nanoparticles, the results obtained are genuine and can contribute to the development of this field. The numerical results are obtained using the bvp4c solver in MATLAB software, and a stability analysis is performed to investigate the stability of the solutions as time evolves.

#### **2. Mathematical Model Description**

A steady 3D rotating boundary layer flow of a carbon nanotube past a stretching or shrinking surface is examined, and its physical model is illustrated in Figure 1. In Figure 1, *x*, *y*, and *z* are the Cartesian coordinates with *x*- and *y*-axes measured in the plane *z* = 0 where nanofluid is in the region *z* ≥ 0. The boundary *z* = 0 is elastic and linearly stretched

or shrunk by two equal and opposite forces in the direction of the *x*-axis with the velocity *Uw* = *ax*. The fluid is rotating at an angular velocity Ω in the *z*-direction. It is assumed that the nanofluid has a uniform wall and ambient temperatures, *Tw* and *T*∞, respectively. Based on the above assumptions, the governing equations of mass, momentum, and energy are (see Nadeem et al. [34] and Salleh et al. [36]):

$$
\frac{\partial \mu}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0 \tag{1}
$$

$$u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z} - 2\Omega v = \frac{\mu\_{nf}}{\rho\_{nf}}\frac{\partial^2 u}{\partial z^2} \tag{2}$$

$$u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + w\frac{\partial v}{\partial z} + 2\Omega u = \frac{\mu\_{nf}}{\rho\_{nf}}\frac{\partial^2 v}{\partial z^2} \tag{3}$$

**Figure 1.** Physical model of the problem.

$$
\mu \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y} + w \frac{\partial T}{\partial z} = \alpha\_{nf} \frac{\partial^2 T}{\partial z^2} \tag{4}
$$

The boundary conditions are

$$\begin{aligned} \mu &= \varepsilon l I\_{w\prime} \ v = 0, \ w = -w\_0, \ T = T\_w \text{ at } z = 0, \\\ u &\to 0, \ v \to 0, \ T \to T\_\infty \text{ as } z \to \infty \end{aligned} \tag{5}$$

where *u*, *v*, and *w* represent the velocity components in the *x*-, *y*-, and *z*-directions, respectively, *T* represents the fluid temperature, *w*0 is the constant mass flux with *w*0 > 0 indicates injection, while *w*0 < 0 is for suction. Further, *ε* is the stretching or shrinking parameter with *ε* > 0 for stretching, *ε* < 0 for shrinking and *ε* = 0 for the static surface. In this study, we adopt the theoretical model proposed by Xue [49] as follows:

$$\mu\_{nf} = \frac{\mu\_f}{(1-\varrho)^{25}}, \; a\_{nf} = \frac{k\_{nf}}{\left(\rho \mathbb{C}\_p\right)\_{nf}}, \; \frac{k\_{nf}}{k\_f} = \frac{1-\varrho+2\varrho\frac{k\_{cnt}}{k\_{nt}-k\_f}\ln\frac{k\_{cnt}+k\_f}{2k\_f}}{1-\varrho+2\varrho\frac{k\_f}{k\_{nt}-k\_f}\ln\frac{k\_{cnt}+k\_f}{2k\_f}}\tag{6}$$

$$\rho\_{nf} = (1-\varrho)\rho\_f + \varrho\rho\_{cnt}, \; \left(\rho\mathbb{C}\_p\right)\_{nf} = (1-\varrho)\left(\rho\mathbb{C}\_p\right)\_f + \varrho\left(\rho\mathbb{C}\_p\right)\_{cnt}.$$

In Equation (6), *μ* is dynamic viscosity, *ρ* is density, *α* is thermal diffusivity, *k* is thermal conductivity, *Cp* is the specific heat capacity at constant pressure, *ϕ* is nanoparticle volume fraction parameter in which the subscripts *nf*, *f*, and *cnt*, respectively, refer to nanofluid, fluid, and carbon nanotube.

The thermophysical properties of SWCNT and MWCNT for different base fluids namely, water (Pr = 6.2) and kerosene (Pr = 21) at the temperature 20–25 ◦C are given in Table 1 (Khan et al. [13]). Water and kerosene are chosen as the base fluids for comparison with the previous studies. Furthermore, water has the highest specific heat that can stabilize the temperature and is used as a cooling agent, while oil-based nanofluids are often utilized in industrial operations for lubrication or high-temperature applications [50]. The following similarity transformations (7) are introduced for solving Equations (2)–(4) together with the boundary conditions (5) that satisfy the continuity Equation (1).

$$u = a\mathbf{x}f'(\eta), \ v = a\mathbf{x}h(\eta), \ w = -\sqrt{a\nu\_f}f(\eta), \ \eta = \sqrt{\frac{a}{\nu\_f}}z, \ \theta(\eta) = \frac{T - T\_{\infty}}{T\_w - T\_{\infty}}\tag{7}$$

**Table 1.** Thermophysical properties of base fluids and carbon nanotubes.

From Equation (7), the velocity mass flux can be defined as

$$w\_0 = -\sqrt{a\nu\_f} \mathcal{S} \tag{8}$$

where *S* = *f*(0), *S* is the suction parameter in which *S* < 0 for injection and *S* > 0 for suction. Substituting Equations (6) and (7) into Equations (2)–(5) yields

$$\frac{1}{(1-q)^{2.5}\left[1-q+\varphi\left(\rho\_{cnt}/\rho\_f\right)\right]}f'''+f f''-f'^2+2\omega h = 0\tag{9}$$

$$\frac{1}{(1-\varrho)^{2.5}\left[1-\varrho+\varrho\left(\rho\_{\rm cut}/\rho\_f\right)\right]}h''+fh'-f'h-2\omega f'=0\tag{10}$$

$$\frac{k\_{nf}/k\_f}{\Pr\left[1-\varrho+\varrho(\rho\mathbb{C}\_p)\_{\text{cut}}/\left(\rho\mathbb{C}\_p\right)\_f\right]}\theta''+f\theta'=0\tag{11}$$

The BCs (5) become

$$\begin{aligned} (0) &= S, \ h(0) = 0, \ f'(0) = \varepsilon, \ \theta(0) = 1, \\ h(\eta) &\to 0, \ f'(\eta) \to 0, \ \theta(\eta) \to 0 \text{ as } \eta \to \infty. \end{aligned} \tag{12}$$

In Equations (9)–(12), prime denotes the differentiation w.r.t. similarity variable *η*, *ω* = Ω/*a* is the rotation parameter, and Pr is the Prandtl number.

The quantities of physical interest are the skin friction coefficients and the local Nusselt number (heat transfer rate at the surface). The skin friction coefficients in the *x*- and *y*directions, respectively, are given by:

$$\mathcal{C}\_{fx} = \frac{\tau\_{xz}}{\rho\_f \mathcal{U}\_w^2}, \mathcal{C}\_{fy} = \frac{\tau\_{yz}}{\rho\_f \mathcal{U}\_w^2} \tag{13}$$

The local Nusselt number *Nux* is defined as:

$$Nu\_{\lambda} = \frac{xq\_w}{k\_f(T\_w - T\_\infty)}\tag{14}$$

The shear stresses *τxz*, *τyz* and the heat flux *qw* at the surface are given by:

$$\tau\_{xz} = \mu\_{nf} \left(\frac{\partial u}{\partial z}\right)\_{z=0}, \tau\_{yz} = \mu\_{nf} \left(\frac{\partial v}{\partial z}\right)\_{z=0}, q\_{w} = -k\_{nf} \left(\frac{\partial T}{\partial z}\right)\_{z=0} \tag{15}$$

Using Equation (15), Equations (13) and (14) become

$$\begin{aligned} \operatorname{Re}\_x^{\frac{1}{2}} \mathbb{C}\_{fx} &= \frac{1}{(1-q)^{25}} f''(0), \; \operatorname{Re}\_x^{\frac{1}{2}} \mathbb{C}\_{fy} = \frac{1}{(1-q)^{25}} h'(0),\\ \operatorname{Re}\_x^{-1/2} \operatorname{Nu}\_x &= -\frac{k\_{uf}}{k\_f} \theta'(0). \end{aligned} \tag{16}$$

where *Rex* = *Ux*/*νf* is the local Reynold number.

#### **3. Stability Analysis**

To investigate the stability of the solutions of Equations (9)–(12) over time, the unsteady case of Equations (2)–(4) are considered as [37,51]

$$u\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z} - 2\Omega v = \frac{\mu\_{nf}}{\rho\_{nf}}\frac{\partial^2 u}{\partial z^2} \tag{17}$$

$$u\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + w\frac{\partial v}{\partial z} + 2\Omega u = \frac{\mu\_{nf}}{\rho\_{nf}}\frac{\partial^2 v}{\partial z^2} \tag{18}$$

$$
\frac{\partial T}{\partial t} + \mu \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y} + w \frac{\partial T}{\partial z} = \alpha\_{nf} \frac{\partial^2 T}{\partial z^2} \tag{19}
$$

where *t* refers to time. Next, new dimensionless variables are introduced as follows:

$$\tau = \frac{\text{lft}}{\text{x}}, \; u = a \text{x} \frac{\partial f}{\partial \eta}(\eta, \tau), \; v = a \text{x} \mathbf{l}(\eta, \tau), \; w = -\sqrt{\text{n} \mathbf{v}\_f} f(\eta, \tau), \; \eta = \sqrt{\frac{a}{\text{v}\_f}} z,\tag{20}$$

$$\theta(\eta, \tau) = \frac{T - T\_{\infty}}{\text{T}\_w - \text{T}\_{\infty}}$$

By utilizing Equation (20), Equations (17)–(19) become:

$$\frac{1}{(1-\varrho)^{2.5}\left[1-\varrho+\varrho\left(\rho\_{\rm cut}/\rho\_f\right)\right]}\frac{\partial^3 f}{\partial \eta^3} + f \frac{\partial^2 f}{\partial \eta^2} - \left(\frac{\partial f}{\partial \eta}\right)^2 + 2\omega h - \frac{\partial^2 f}{\partial \eta \partial \tau} = 0. \tag{21}$$

$$\frac{1}{(1-\varrho)^{2.5}\left[1-\varrho+\varrho\left(\rho\_{cnt}/\rho\_f\right)\right]}\frac{\partial^2 h}{\partial \eta^2} + f \frac{\partial h}{\partial \eta} - h \frac{\partial f}{\partial \eta} - 2\omega \frac{\partial f}{\partial \eta} - \frac{\partial h}{\partial \tau} = 0. \tag{22}$$

$$\frac{k\_{nf}/k\_f}{\Pr\left[1-\varrho+\varrho(\rho\mathbb{C}\_p)\_{cnt}/\left(\rho\mathbb{C}\_p\right)\_f\right]}\frac{\partial^2\theta}{\partial\eta^2} + f\frac{\partial\theta}{\partial\eta} - \frac{\partial\theta}{\partial\tau} = 0\tag{23}$$

where their respective boundary conditions are

$$\begin{aligned} f(0, \tau) &= S, \; h(0, \tau) = 0, \; \frac{\partial f}{\partial \eta}(0, \tau) = \varepsilon, \; \theta(0, \tau) = 1, \\\ h(\eta, \tau) &\to 0, \; \frac{\partial f}{\partial \eta}(\eta, \tau) \to 0, \; \theta(\eta, \tau) \to 0 \text{ as } \eta \to \infty. \end{aligned} \tag{24}$$

The stability of the steady flow solutions *f*(*η*) = *f*0(*η*), *h*(*η*) = *h*0(*η*), and *θ*(*η*) = *θ*0(*η*) can be determined by writing *f*(*η*, *τ*), *h*(*η*, *τ*) and *θ*(*η*, *τ*) as follows [51]:

$$f(\eta, \tau) = f\_0(\eta) + e^{-\gamma \tau} F(\eta, \tau), \; h(\eta, \tau) = h\_0(\eta) + e^{-\gamma \tau} H(\eta, \tau), \tag{25}$$

$$\theta(\eta, \tau) = \theta\_0(\eta) + e^{-\gamma \tau} G(\eta, \tau). \tag{25}$$

In Equation (25), *F*(*η*, *τ*), *H*(*η*, *τ*), and *G*(*η*, *τ*) are smaller relative to *f* = *f*0(*η*), *h* = *h*0(*η*) and *θ* = *θ*0(*η*), respectively, and *γ* is the eigenvalue. Substituting Equation (25) into Equations (21)–(24), the following system of linear eigenvalue problem is obtained:

$$\begin{split} \frac{1}{(1-\varrho)^{2J}\left[1-\varrho+\varrho(\rho\_{cnt}/\rho\_f)\right]} \frac{\partial^3 F}{\partial \eta^3} + f\_0 \frac{\partial^2 F}{\partial \eta^2} + F \frac{\partial^2 f\_0}{\partial \eta^2} - 2 \frac{\partial f\_0}{\partial \eta} \frac{\partial F}{\partial \eta} + 2\omega H \\ + \gamma \frac{\partial F}{\partial \eta} - \frac{\partial^2 F}{\partial \eta \partial \tau} = 0 \end{split} \tag{26}$$

$$\begin{split} \frac{1}{(1-\varrho)^{2.5}\left[1-\varrho+\varrho\left(\frac{\rho\_{\rm cut}}{\vartheta\_f}\right)\right]} \frac{\partial^2 H}{\partial \eta^2} + f\_0 \frac{\partial H}{\partial \eta} + F \frac{\partial h\_0}{\partial \eta} - H \frac{\partial f\_0}{\partial \eta} - h\_0 \frac{\partial F}{\partial \eta} + \gamma H \\ -\frac{\partial H}{\partial \tau} - 2\omega \frac{\partial F}{\partial \eta} = 0 \end{split} \tag{27}$$

$$\frac{k\_{nf}/k\_f}{\Pr\left[1-\varrho+\varrho(\rho\mathbb{C}\_{\mathbb{P}})\_{cnt}/(\rho\mathbb{C}\_{\mathbb{P}})\_{f}\right]}\frac{\partial^2 G}{\partial\eta^2} + f\_0 \frac{\partial G}{\partial\eta} + F \frac{\partial\theta\_0}{\partial\eta} + \gamma G - \frac{\partial G}{\partial\tau} = 0 \tag{28}$$

subject to the boundary conditions

$$\begin{aligned} F(0, \tau) &= 0, \ H(0, \tau) = 0, \ \frac{\partial F}{\partial \eta}(0, \tau) = 0, \ G(0, \tau) = 0 \\\ H(\eta, \tau) &\to 0, \ \frac{\partial F}{\partial \eta}(\eta, \tau) \to 0, \ G(\eta, \tau) \to 0 \text{ as } \eta \to \infty. \end{aligned} \tag{29}$$

Following Weidman et al. [51], we take *τ* = 0 to determine the initial decay or growth of the solutions in Equation (25) and functions *F*, *H*, and *G* can be written as *F*0(*η*), *H*0(*η*), and *G*0(*η*), respectively. Thus, Equations (26)–(29) can be simplified as

$$\frac{1}{(1-\varrho)^{2.5}\left[1-\varrho+\varrho\left(\rho\_{cnt}/\rho\_f\right)\right]}F\_0''' + f\_0F\_0'' + F\_0f\_0'' - 2f\_0'F\_0' + 2\omega H\_0 + \gamma F\_0' = 0 \tag{30}$$

$$\frac{1}{(1-\varrho)^{2.5}\left[1-\varrho+\varrho\left(\frac{\rho\_{cut}}{\rho\_f}\right)\right]}H\_0'' + f\_0H\_0' + F\_0h\_0' - H\_0f\_0' - h\_0F\_0' + \gamma H\_0 - 2\omega F\_0' = 0 \quad \text{(31)}$$