**1. Introduction**

The very first consideration of the problem involving fractional differential equation is traced back to the very end of the seventeenth century. It was Leibniz who introduced the notation *dn*/*dx<sup>n</sup>* and L'Hospital who asked Leibniz: "What if *n* be 1/2?". Leibniz response to this question has laid the foundations of what is today known as fractional calculus (FC) [1]. Despite these early attempts made by outstanding minds who were definitely ahead of their time, further development towards practical applications appeared to be somewhat slow. An interested reader can refer to an excellent review of Ross for details on early stages of FC development [2].

The fundamentals of FC along with a brief applications overview can be found in classical books by Miller and Ross [3], Podlubny [4], or more recent one by Mainardi [5]. In the following decades it became obvious that with the incredible potential of FC to tackle problems in absolutely different fields, each of them should be addressed separately. By these means the general approaches in solving fractional differential equations could be tailored to reflect the specifics of a given field, optimize solution process, and formulate reliable constraints and ranges for, sometimes, purely mathematical parameters. Along these lines, readers with a background in physics can refer to [6–9], while those from the engineering field can refer to [10,11].

The comprehensive overview summarizing state-of-the-art practical applications of FC has been recently published by The Royal Society Publishing. The sixteen-paper issue entitled "Advanced materials modeling via fractional calculus: challenges and perspectives" [12–27] covers applications of constant-order (CO) and variable-order (VO) fractional differential operators to several fundamental phenomena. These include anomalous diffusion, [13,16] heat conduction [14,27], fractional viscoelasticity of fluids [19], and materials [12,18,22]. The approach to model viscoelastic properties of materials with VO FC operators is undoubtedly among the most promising ones, as it allows for the consideration of fractional order dynamics with respect to time, space, and material variables [22].

Viscoelasticity has been in the scope of researchers' interest since the nineteenth century and has gone through gradual development. With mathematical apparatus available at the time, several mechanical models have been proposed, developed, and generalized to address the problem of a more accurate description of material properties observed. Those classical models named after researchers who made a significant contribution to the field include but are not limited to Maxwell, Kelvin, Voigt, Zener, and Burgers. Based on two principal elements, spring and dashpot, connected in series and/or in parallel, those models met the needs of adequate material response description towards stresses and strains applied. However, as it often happened, experimental results started to accumulate, which illustrated the behaviors beyond the above-mentioned models. This issue was addressed from scratch, namely via rethinking of a basic element by Scott-Blair in 1947 [28]. He came up with an idea of considering a single element capable of a simultaneous description of viscous and elastic properties of the material. This approach capitalized on a notion of fractional derivative and allowed for the alteration of the balance between viscous and elastic properties without any additional complexity, but covering a wide range of materials. The systematic study of fractional calculus applications to viscoelasticity was made by Bagley and Torvik [29], who laid theoretical foundations of this approach. Since that time the field had emerged, especially when biomedical applications were outlined. Historically, several approaches were developed to implement balancing between viscous and elastic properties of materials. Initially, ordinary first-order constitutive equations were straightforwardly replaced with fractional counterparts [30–34]. An obvious drawback of this approach was its purely phenomenological character. An alternative approach implied physical representation of fractional constitutive equations via hierarchical combinations of dashpots and springs [35–38]. This concept (so-called ladder models) has been successfully implemented in further studies by Schiessel et al. [39] and Friedrich et al. [40]. The authors have considered generalized viscoelastic models, replaced them with fractional ones, and obtained their analytical solutions in terms of relaxation modulus and creep compliance.

Numerous studies have demonstrated the superiority of fractional calculus approach compared to classical one in predicting viscoelastic properties of materials and flows in various geometries. Hernandez et al. [41] studied the behavior of relaxation modulus for polymethyl methacrylate (PMMA) and polytetrafluoroethelene (PTFE) and demonstrated much more accurate fitting of experimental data using fractional Maxwell model compared to integer-order one. Markis et al. [42] proposed and experimentally verified the generalized fractional Maxwell model in the design of damper systems for seismic and vibration isolation. Zhang and coworkers [43] have studied the stress-relaxation behavior of fabrics coated with PTFE under various temperatures. The authors demonstrated the superiority of fractional Maxwell model in predicting stress-relaxation behavior. This behavior was experimentally proven to be nonlinear, while predicted to be linear by the classic Maxwell model. Moreover, compared to the generalized Maxwell model, the fractional one appeared to be much easier, as it did not require a large number of structural units to increase accuracy. The similar results for fractional Maxwell and fractional Zener models were obtained for elastomers (carbon-black filled resins) [44–46], polymers and rocks [47], and biological materials [48]. Fractional Kelvin–Voigt model was found to be efficient in predicting viscoelastic behavior of sludge [49].

Fractional calculus approach to viscoelasticity found applications not only for solids, as shown above, but also for various fluids. A brief historic retrospective reveals that both classic and fractional viscoelastic fluid flows in constrained geometries, especially pipes and ducts, turned out to be of a particular researchers' interest for about a century. This can be attributed to the fact that such geometries are widely used for practical applications and also admit relatively easy analytical solutions. Oscillatory pipe flows were among thoroughly investigated. Its classical consideration was traced back to 1920–1930s [50,51]. These seminal studies were further extended and generalized by Wornersley, Uchida, and others [52–55]. Later researchers considered dynamics of oscillatory pipe flows at various Reynolds numbers [56,57], both experimentally and theoretically, and specifically investigated transition to turbulence [58,59]. The foundations of fractional derivatives towards fluids viscoelastic behavior were laid by Bagley and Torvik [29]. Wood [60] studied viscoelastic transient flows in cylindrical pipes and annulus. Yin et al. [61] provided a theoretical framework for oscillating viscoelastic pipe flow using fractional Maxwell model. The authors demonstrated a drastic difference in velocity profile compared to the integer-order Maxwell model. Viscoelastic start-up flow with the fractional Maxwell model was considered by Yang et al. [62]. A similar problem in the annular pipe using fractional Burgers model was solved by Shah et al. [63]. An unsteady viscoelastic flow in a cylinder [64] and rectangular duct [65] (using fractional Maxwell model) were also considered. The exact solutions of unsteady flow in cylindrical domains with Maxwell fractional model were derived by Khandelwal and Mathur [66,67]. Maqbool and coworkers [68] considered a flow of generalized fractional Burgers fluid in inclined tube. Tang et al. [69] studied nonlinear free vibrations of a pipe conveying fractional viscoelastic fluid. They demonstrated decreasing mode amplitudes with increasing fractional order. Wang and Chen [70] considered a similar problem for the pipeline conveying fractional fluid more accurately employing Legendre polynomials. Javadi et al. [71] investigated the effect of gravity on fractional viscoelastic fluid flow in a pipe and addressed the problem of stability for it. This selected list of fractional calculus applications is far from complete. An interested reader can refer to a recent review of Sun et al. [72]. A big picture of various applications of FC is also given in [17,23].

However, despite abundant theoretical, numerical, and experimental results on fractional viscoelastic flow in pipes, the general, unifying theoretical approach to tackle flow dynamics is still missing. To fill in this research gap, this paper provides exact analytical solutions for velocity profiles and shear stresses. We demonstrate that the same solution form is applicable for different viscoelastic models including Maxwell, Kelvin–Voigt, Zener, Poynting–Thomson, and Burgers. Velocity profiles and shear stresses are studied parametrically with respect to fractional order and elastic properties (Young's modulus ratio, starting from 3-element model). This paper is organized as follows. We first provide general problem formulation for fluid flow along with domain definition (Section 2.1). In Section 2.2 we introduce the notion of the fractional element and its governing equation. Sections 2.3–2.7 provide constitutive equations for the most common fractional viscoelastic models. The main results of this study appear in Section 3. In particular, Section 3.1 describes the general approach of seeking the solution along with brief revisiting of applicable transformations. Sections 3.2–3.6 present and discuss centerline velocity and shear stress profiles dynamics with varying fractional order and elastic parameters. Two-, three-, and four-component fractional models are considered. Finally, the results of this study are summarized in Section 4.

#### **2. Problem Formulation**

#### *2.1. Domain Definition*

Introduce cylindrical coordinate system (*r*, *θ*, *z*) and consider laminar flow of incompressible viscoelastic fluid along *z* axis of infinitely long pipe of radius *R* with circular cross-section. The domain of interest is shown in Figure 1. The governing equations are momentum and continuity that read as:

$$
\rho \frac{d\mathbf{u}}{dt} = -\nabla p + \text{div}\sigma \tag{1}
$$

$$\text{div}\!\!u = 0,\tag{2}$$

where *<sup>ρ</sup>*, *u*, *<sup>p</sup>*, *<sup>σ</sup>* are fluid mass density, flow velocity, pressure, and stress, respectively and *<sup>d</sup> dt* stands for total derivative. Provided *<sup>u</sup>* <sup>=</sup> *uz*(*r*, *<sup>t</sup>*)*ez*, where *ez* is a unit vector in *<sup>z</sup>*-direction and all but *<sup>σ</sup>rz* stress tensor components are set to 0, the momentum equation reduces to:

$$
\rho \frac{\partial u\_z}{\partial t} = -\frac{\partial p}{\partial z} + \frac{1}{r} \frac{\partial (r v\_{rz})}{\partial r},
\tag{3}
$$

where *uz* stands for *z*-component of flow velocity and all other quantities are defined above. Additionally, we apply nonslip boundary condition at the wall.

**Figure 1.** Domain definition for the circular pipe.

#### *2.2. Fractional Element*

The idea of introducing basic element accounting for both viscous and elastic properties of the material belongs to Blair [28]. Here we follow a conventional approach of constructing more complex mechanical models based on this element. For such an element, the stress-strain relation (constitutive equation) reads as:

$$
\sigma(t) = E \pi^{\mathfrak{a}}{}\_{\mathfrak{a}} D^{\mathfrak{a}}\_{t} \varepsilon(t), \quad 0 \le \mathfrak{a} \le 1,\tag{4}
$$

where *σ*(*t*),*ε*(*t*), *E*, *τ*, *α* stand for stress, strain, Young's modulus, relaxation time and fractional order, respectively. The limiting cases of *α* = 0 and *α* = 1 represent spring and dashpot, respectively. The schematic of the element is given in Figure 2a. Here we use a conventional notation of Riemann–Liouville fractional derivative of a smooth function, *f*(*t*), given by:

$$\_aD\_t^a f(t) = \frac{1}{\Gamma(k-a)} \frac{d^k}{dt^k} \int\_a^t (t-t\_0)^{k-a-1} f(t\_0) dt\_0. \tag{5}$$

where *k* is the integer, *α* is the fractional order, and Γ(·) stands for Gamma-function defined as:

$$\Gamma(\mathbf{x}) = \int\_0^\infty \mathbf{e}^{-t\_0} t\_0^{\mathbf{x}-1} dt\_0 \tag{6}$$

In the following subsections we will provide constitutive equations for two-, three-, and four-component fractional viscoelastic models.

**Figure 2.** Fractional models: (**a**) basic element; (**b**) Maxwell; (**c**) Kelvin–Voigt; (**d**) Zener; (**e**) Poynting–Thomson; (**f**) Burgers.

#### *2.3. Fractional Maxwell Model*

We start our consideration from the simplest yet commonly used two-component fractional model with basic elements connected in series (fractional Maxwell model, Figure 2b). The constitutive equation for this model reads as [39]:

$$
\sigma(t) + \tau^{a-\beta} \, \_d D\_t^{a-\beta} \sigma(t) = E \tau^a \, \_d D\_t^a \varepsilon(t), \quad 0 \le a \le 1, \quad 0 \le \beta \le 1, \quad a \ge \beta,\tag{7}
$$

where *α* and *β* stand for the fractional orders, *τ*, *E* are defined as:

$$\tau = \left(\frac{E\_1 \tau\_1^a}{E\_2 \tau\_2^b}\right)^{\frac{1}{a-b}}, \quad E = E\_1 \left(\frac{\tau\_1}{\tau}\right)^a, \quad E = E\_2 \left(\frac{\tau\_2}{\tau}\right)^b,\tag{8}$$

and *E*1, *τ*1, *E*2, *τ*<sup>2</sup> are the Young's moduli, relaxation times for elements "1" and "2", respectively. Constitutive equation of fractional Maxwell model (7) can be reduced to classical one, if *α* = 1 and *β* = 0. Moreover, if we set *α* = 1 and *β* = 1, Newtonian fluid can be obtained. The ranges set for fractional orders *α* and *β* have both mathematical and physical meaning. Non-negative values of fractional orders and their difference reflect the fact that the dynamics of processes considered are described by fractional derivatives not fractional integrals. The upper limit of fractional orders range is set to get the corresponding classical viscoelastic model as a limiting case of a fractional one. For *τ* as it was defined above, *α* = *β* represents a special case. Indeed, it corresponds to relaxation time blowing up. Physically it means that a fluid becomes a critical gel. This special case is, however, beyond the scope of the current study. It is worth mentioning that relaxation time blows up only provided that the base of *τ* is greater than unity. If, however, that is not the case, i.e., *τ* base is less than unity, *τ* itself remains bound.

#### *2.4. Fractional Kelvin–Voigt Model*

Now let us connect two basic fractional elements in parallel. The dynamics of corresponding model (Kelvin–Voigt fractional model) are described as follows [39]:

$$
\sigma(t) = E \tau^a{}\_a D\_t^a \varepsilon(t) + E \tau^\beta{}\_a D\_t^\beta \varepsilon(t), \quad E = E\_2 \left(\frac{\tau\_2}{\tau}\right)^\beta,\tag{9}
$$

where all the notations are similar to those for fractional Maxwell model. The schematic of the model is given in Figure 2c. If *α* = 1 and *β* = 0, we appear at classical Kelvin–Voigt model.

#### *2.5. Fractional Zener Model*

Let us now consider 3-element model referred to as fractional Zener model (Figure 2d). The corresponding constitutive equation reads as [39]:

$$\sigma(t) + \tau^{a-\beta} \, \_a D\_t^{a-\beta} \sigma(t) = E \tau^a \, \_a D\_t^a \varepsilon(t) + E\_0 \tau^\gamma \, \_a D\_t^\gamma \varepsilon(t) + E\_0 \tau^{\gamma+a-\beta} \, \_a D\_t^{\gamma+a-\beta} \varepsilon(t),$$

$$0 \le a \le 1, \quad 0 \le \beta \le 1, \quad 0 \le \gamma \le 1, \quad a - \beta \ge 0, \quad \gamma + a - \beta \ge 0, \quad \text{(10)}$$

where *γ* stands for the fractional order of element "3" and *E*<sup>0</sup> is given by:

$$E\_0 = E\_3 \left(\frac{\tau\_3}{\tau}\right)^\gamma,\tag{11}$$

with *E*3, *τ*<sup>3</sup> being Young's modulus and relaxation time for element "3" and all other notations being the same as above. Assuming *α* = 1, *β* = 0, and *γ* = 0, classical Zener model can be restored.

#### *2.6. Fractional Poynting–Thomson Model*

For fractional Poynting–Thomson model (Figure 2e) the constitutive equation is [39]:

$$\sigma(t) + \frac{E}{E\_0} \tau^{a-\gamma} {}\_a D\_t^{a-\gamma} \sigma(t) + \frac{E}{E\_0} \tau^{\beta-\gamma} {}\_a D\_t^{\beta-\gamma} \sigma(t) = E \tau^a {}\_a D\_t^a \varepsilon(t) + E \tau^{\beta} {}\_a D\_t^{\beta} \varepsilon(t),$$

$$0 \le a \le 1, \quad 0 \le \beta \le 1, \quad 0 \le \gamma \le 1, \quad a - \gamma \ge 0, \quad \beta - \gamma \ge 0, \quad \text{(12)}$$

with all the notations similar to above, but different limitations applied to fractional orders *α*, *β* and *γ*. Provided *α* = 1, *β* = 0, and *γ* = 0, we get classical Poynting–Thomson model.

#### *2.7. Fractional Burgers Model*

Finally, for the most complex, 4-element model considered in this study (fractional Burgers model, Figure 2f) the constitutive equation reads as [5]:

$$\sigma(t) + \frac{\mathbb{E}}{\mathbb{E}\_0} \pi\_a^{a-\gamma} D\_t^{a-\gamma} \sigma(t) + \frac{\mathbb{E}}{\mathbb{E}\_0} \pi^{a-\delta} \,\_d D\_t^{a-\delta} \sigma(t) + \frac{\mathbb{E}}{\mathbb{E}\_0} \pi\_a^{\delta-\gamma} D\_t^{\delta-\gamma} \sigma(t) + \frac{\mathbb{E}}{\mathbb{E}\_0} \pi\_a^{\delta-\delta} D\_t^{\delta-\delta} \sigma(t) = $$
 
$$= E \pi^a \,\_d D\_t^a \varepsilon(t) + E \pi\_d^\delta D\_t^\beta \varepsilon(t), \tag{13}$$
 
$$0 \le a \le 1, \quad 1 \le \beta \le 2, \quad 0 \le \gamma \le 1, \quad a - \gamma \ge 0, \quad \beta - \gamma \ge 0, \quad a - \delta \ge 0, \quad \beta - \delta \ge 0$$

where *δ* stands for the fractional order of element "4", (11) is still valid and additionally:

$$E\_0 = E\_4 \left(\frac{\tau\_4}{\tau}\right)^\delta,\tag{14}$$

with *E*4, *τ*<sup>4</sup> being Young's modulus and relaxation time for element "4" and all other notations being the same as above. If we set *α* = 1, *β* = 2, *γ* = 0, and *δ* = 0, classical Burgers model can be restored. Here *β* − *γ* ≥ 0 and *β* − *δ* ≥ 0 are satisfied automatically, given the range for *β*, *γ*, and *δ*.

#### **3. Results and Discussion**

#### *3.1. General Solution*

Here we follow the approach proposed earlier by Yin [61] for Maxwell model and then generalize it for all other models. Consider a specific type of pressure gradient given as follows:

$$\frac{\partial p}{\partial z} = P\_0 \mathbf{e}^{i\omega t} \tag{15}$$

Introduce Fourier and inverse Fourier transforms of flow velocity and pressure:

$$\mathrm{d}I\_{z}(r,\omega) = \int\_{-\infty}^{+\infty} u\_{z}(r,t) \mathrm{e}^{-i\omega t} dt, \quad u\_{z}(r,t) = \frac{1}{2\pi} \int\_{-\infty}^{+\infty} \mathrm{d}I\_{z}(r,\omega) \mathrm{e}^{i\omega t} d\omega \tag{16}$$

$$P(r,\omega) = \int\_{-\infty}^{+\infty} p(r,t) \mathbf{e}^{-i\omega t} dt, \quad p(r,t) = \frac{1}{2\pi} \int\_{-\infty}^{+\infty} P(r,\omega) \mathbf{e}^{i\omega t} d\omega \tag{17}$$

Recall Fourier transform rule for fractional derivative:

$$\mathcal{F}\{\_{a}D\_{t}^{a}\mu\_{z}\} = (i\omega)^{a}\mathcal{U}\_{z}(r,\omega), \quad \mathcal{F}\{\_{a}D\_{t}^{a-\beta}\frac{\partial p}{\partial z}\} = (i\omega)^{a-\beta}\frac{\partial P(r,\omega)}{\partial z}, \tag{18}$$

where appropriate integration limits are implied. In particular, lower terminal value is set to: *a* = −∞. Then the general algorithm to find flow velocity is the following: (1) express stress and its derivative from the constitutive eq-n; (2) plug stress and its derivative in momentum eq-n; (3) eliminate stress and its derivative and get modified momentum eq-n; (4) perform Fourier transform of modified momentum eq-n; (5) solve corresponding ODE for Fourier transform of velocity (*Uz*(*r*, *ω*)); (6) change variables with *ξ* for simplicity; (7) perform inverse Fourier transform of *Uz*(*r*, *ω*) to get *uz*(*r*, *t*). Regardless of the model considered, *z*-component of the velocity reads as:

$$u\_z(r,t) = \frac{i}{\rho \omega} \frac{\partial p(z,t)}{\partial z} \left[1 - \frac{f\_0(\xi r)}{f\_0(\xi R)}\right],\tag{19}$$

where *Jn*(·) is the Bessel function of the first kind of order *n*. The only difference that defines behavior of the system is hidden in *ξ*. At the same time, the form of the solution reproduces a well-known classical result (see, for example Reference [73]). Thus it represents a natural extension of integer-order models using apparatus of FC. As long as velocity profile was defined, we could also get expression for the stress. The direct integration of (3) accounting for (15) results in:

$$
\sigma\_{rz} = \frac{\partial p}{\partial z} \frac{J\_1(\not\xi r)}{\not\xi J\_0(\not\xi R)},\tag{20}
$$

where again the individual properties of the model are hidden in parameter *ξ*. Physical meaning of *ξ* is worth consideration. First, of all, [*ξ*]=[1/m], so that it can be considered as an "effective" wave number. Moreover, as shown below, *<sup>ξ</sup>* is related to complex compliance. Introducing *<sup>σ</sup>*<sup>0</sup> = *<sup>∂</sup><sup>p</sup> <sup>∂</sup><sup>z</sup> R*, nondimensional stress *σ*˜ = *σrz*/*σ*<sup>0</sup> is given by:

$$\mathcal{O} = \frac{J\_1(\mathcal{J}r)}{\mathcal{J}\mathcal{R}J\_0(\mathcal{J}\mathcal{R})} \tag{21}$$

*Appl. Sci.* **2020**, *10*, 9093

For practical purposes and analysis simplification, it is more convenient to introduce the following nondimensional quantities and relationships:

$$
\vec{r} = \frac{r}{\sqrt{\nu \tau}}, \quad \vec{a} = \frac{R}{\sqrt{\nu \tau}}, \quad \vec{\omega} = \omega \tau, \quad \Phi = \frac{\vec{r}}{\vec{a}}, \quad u\_0 = \frac{1}{4\rho\omega} \frac{\partial p}{\partial z}, \quad r^2 \xi^2 = \bar{\omega} \vec{r}^2,\tag{22}
$$

where *ν* is a kinematic viscosity and all other quantities were defined earlier. Then velocity and stress profiles are given by:

$$\vec{u} = \left| \frac{u\_z}{u\_0} \right| = \left( 1 - \frac{J\_0(\sqrt{\vec{\alpha}} \vec{a} \phi)}{J\_0(\sqrt{\vec{\alpha}} \vec{a})} \right), \quad \vec{v} = \left| \frac{\sigma\_{rz}}{\sigma\_0} \right| = \frac{J\_1(\sqrt{\vec{\alpha}} \vec{a} \phi)}{\sqrt{\vec{\alpha}} \vec{a} J\_0(\sqrt{\vec{\alpha}} \vec{a})} \tag{23}$$

As velocity profiles and shear stresses have been defined, we can now proceed with considering specific fractional viscoelastic models and outlining their specifics. Six parameters will be considered: fractional order (*α*, *β*, *γ*, *δ*), nondimensional radius, *a*˜, and elastic ratio, *ψ*.

#### *3.2. Fractional Maxwell Model*

Let us start from the simplest two-parameter fractional Maxwell model. As outlined above, the key quantity that defines the behavior of a model is *ξ*. For this model *ξ* is given by:

$$\zeta^2 = \frac{\rho \omega^2}{E} \left[ \frac{1}{(i\omega \tau)^a} + \frac{1}{(i\omega \tau)^\beta} \right],\tag{24}$$

or alternatively in nondimensional form:

$$
\omega = \tilde{\omega}^2 \left[ \frac{1}{(i\tilde{\omega})^a} + \frac{1}{(i\tilde{\omega})^\beta} \right] \tag{25}
$$

Centerline velocity profiles for fractional Maxwell model are shown in Figure 3. For all the plots *α* = 0.5. The value of *a*˜ ranges from *a*˜ = 0.01 (Figure 3a) to *a*˜ = 0.1 (Figure 3c). For this model an additional condition is imposed: *α* ≥ *β*, thus only a half of *β* range is considered given the value of *α* set. As can be seen from this figure, the behavior of centerline velocity changes dramatically from monotonically increasing to resonant or oscillatory. Centerline velocity profile also changes with increasing fractional order (*β*). For lower values of *a*˜, an increase of *β* results in decreasing of centerline velocity up to the constant value (as *β* → 0.5). For higher values of *a*˜ (Figure 3b,c) the system exhibits aperiodic oscillations at low *β* and demonstrates a switch from resonant behavior to monotonically increasing as *β* increases. It is also worth mentioning that peak values decrease with increasing *β*. What is the reason for such a dramatic change of a velocity profile dynamics with varying *a*˜? In fact, *a*˜<sup>2</sup> turns out to be nothing else but a Reynolds number: *a*˜<sup>2</sup> = *Re*, one of the key flow parameters.

**Figure 3.** Centerline velocity profile for fractional Maxwell model: (**a**) *a*˜ = 0.01; (**b**) *a*˜ = 0.05; (**c**) *a*˜ = 0.1. For all plots *α* = 0.5. Subscript "c" stands for "centerline".

Now let us have a look at the shear stress dynamics (at the wall) shown in Figure 4. Here again we fix *α* = 0.5 and vary *β* (*α* ≥ *β*). Similar to velocity profile, shear stress increases monotonically with *ω*˜ for lower values of *a*˜. The trend, however, changes to almost constant with increasing fractional order. As *a*˜ increases, the system exhibits more complex behaviors. In particular, there are aperiodic oscillations observed for lower values of fractional order that switch to resonant ones and monotonic decay as the fractional order starts to increase. Moreover, peak values shift to higher values of *ω*˜ with increasing fractional order and the peaks themselves become more dispersive. At the same time, an increase of *a*˜ results in resonant peaks becoming sharper and shifting to lower values of *ω*˜ .

**Figure 4.** Shear stress profile at the pipe wall for fractional Maxwell model: (**a**) *a*˜ = 0.01; (**b**) *a*˜ = 0.05; (**c**) *a*˜ = 0.1. For all plots *α* = 0.5. Subscript "w" stands for "wall".

Finally, for relatively high values of *ω*˜ and *a*˜, shear stress plots converge and become nearly independent from the fractional order *β*. To better understand physical meaning behind fractional orders *α* and *β*, we have examined its influence on centerline velocity and shear stress dynamics in a wider range of *ω*˜ (0 ≤ *ω*˜ ≤ 1000) with the fixed value of *a*˜ (*a*˜ = 0.1) as shown in Figure 5. We have first fixed *α* = 0.5 and varied *β*: 0.1 ≤ *β* ≤ 0.5. For both centerline velocity (Figure 5a) and shear stress (Figure 5b) low values of *β* resulted in oscillatory profiles that in turn reflects the fact that the fractional order *β* describes elastic properties of the fluid considered. In contrast, if we fix *β* = 0.5 and vary *α*: 0.6 ≤ *α* ≤ 1, the oscillations are damped much quicker for both centerline velocity and shear stress (Figure 5c,d). Finally, let us set *α* = 1 and vary *β*: 0.25 ≤ *β* ≤ 1. Corresponding plots for centerline velocity and shear stress are shown in Figure 5e,f, respectively. When *α* = 1 and *β* = 1, Newtonian fluid is obtained with centerline velocity monotonically reaching constant value and shear stress monotonically decaying.

**Figure 5.** *Cont.*

**Figure 5.** Fractional Maxwell model: (**a**) centerline velocity, *α* = 0.5; (**b**) shear stress, *α* = 0.5; (**c**) centerline velocity, *β* = 0.5; (**d**) shear stress, *β* = 0.5; (**e**) centerline velocity, *α* = 1; (**f**) shear stress, *α* = 1. For all plots *a*˜ = 0.1. Subscripts "c" and "w" stand for "centerline" and "wall", respectively.

#### *3.3. Fractional Kelvin–Voigt Model*

Next, we consider another commonly used two-parameter fractional model referred to as Kelvin–Voigt. We start from parameter *ξ*, responsible for system behavior. It reads as:

$$\zeta^2 = \frac{\rho \omega^2}{E} \left[ \frac{1}{(i\omega \tau)^a + (i\omega \tau)^\beta} \right],\tag{26}$$

or alternatively in nondimensional form:

$$
\bar{\omega} = \frac{\bar{\omega}^2}{\left[ (i\bar{\omega})^a + (i\bar{\omega})^\beta \right]} \tag{27}
$$

In contrast to fractional Maxwell model, no additional relations between *α* and *β* are implied here. Thus, we have more freedom to set the values of these fractional orders within the entire range.

Centerline velocity profiles and shear stresses for fractional Kelvin–Voigt model are given in Figure 6. Here again we fix *α* = 0.5 and vary *β*. Different from fractional Maxwell model, centerline velocity profiles (Figure 6a,c,e) do not exhibit aperiodic oscillations with varying both fractional order and *a*˜. Centerline velocity amplitudes increase with increasing *a*˜ and decrease with increasing fractional order *β*. Shear stresses' dynamics for the same model are presented in Figure 6b,d,f. Here again all three types of behavior encountered for fractional Maxwell model are present. More specifically, as *a*˜ increases, shear stresses experience three different types of behavior: from monotonically increasing through resonant to oscillatory. Relative stress amplitudes decrease with increasing fractional order *β*.

**Figure 6.** Fractional Kelvin–Voigt model: (**a**) centerline velocity, *a*˜ = 0.01; (**b**) shear stress, *a*˜ = 0.01; (**c**) centerline velocity, *a*˜ = 0.05; (**d**) shear stress, *a*˜ = 0.05; (**e**) centerline velocity, *a*˜ = 0.1; (**f**) shear stress, *a*˜ = 0.1. For all plots *α* = 0.5. Subscripts "c" and "w" stand for "centerline" and "wall", respectively.

#### *3.4. Fractional Zener Model*

Now examine the behavior of three-element fractional Zener model. The key model variable, *ξ*, in this case becomes:

$$\mathcal{J}^2 = \frac{\rho \omega^2}{E} \left[ \frac{1}{(E\_0/E)(i\omega \tau)^\gamma + \frac{(i\omega \tau)^{a+\beta}}{(i\omega \tau)^a + (i\omega \tau)^\beta}} \right],\tag{28}$$

or alternatively in nondimensional form:

$$
\bar{\omega} = \frac{\bar{\omega}^2}{\psi (i\bar{\omega})^\gamma + \frac{(i\bar{\omega})^{a+\beta}}{(i\bar{\omega})^a + (i\bar{\omega})^\beta}} \,\,\,\tag{29}
$$

where *ψ* = *<sup>E</sup>*<sup>0</sup> *<sup>E</sup>* is an elastic parameter. By definition elastic parameter *ψ* can be rewritten as:

$$
\psi = \frac{E\_3 \pi\_3^{\gamma}}{E\_1 \pi\_1^{\alpha}} \pi^{\gamma - \alpha} \tag{30}
$$

For classical Zener model (*α* = 1, *β* = 0, and *γ* = 0), *ψ* reduces to: *ψ* = *E*3/*E*2, i.e., The ratio of Young's moduli, thus justifying the notion introduced above. Similar to fractional Maxwell model, we impose additional conditions for fractional orders *α*, *β*, and *γ*. From two inequalities: *α* − *β* ≥ 0 and *γ* + *α* − *β* ≥ 0, we end up with the first one. The second one is satisfied automatically, provided the first one is and *γ* being non-negative.

We first examine the dynamics of centerline velocity. Different from two-element models, where we considered the influence of three parameters (*α*, *β*, and *a*˜) on the system behavior, here we need to account for five parameters (*a*˜, *α*, *β*, *γ* and *ψ*). Let us fix *α* = 0.5 (as previously), *γ* = 0.5 and *ψ* = 1. Since *α* − *β* ≥ 0, we only consider 0 ≤ *β* ≤ 0.5. Corresponding centerline velocity profiles are given in Figure 7a,d,g. Centerline velocity changes its behavior from monotonically increasing to resonant with increasing *a*˜. Moreover, relative velocity amplitudes increase with increasing *a*˜. In addition to it, an interesting phenomenon is observed in Figure 7g. Relative velocity amplitudes tend to decrease with increasing fractional-order *β* up to a certain value (*ω*˜ ∼ 150). At this point, centerline velocities for all fractional orders become close, and then the trend inverses, namely velocity amplitudes start to increase with increasing fractional-order *β*.

**Figure 7.** Centerline velocity for fractional Zener model: (**a**) *a*˜ = 0.01, *α* = *γ* = 0.2, *ψ* = 1; (**b**) *a*˜ = 0.01, *α* = *β* = 0.5, *ψ* = 1; (**c**) *a*˜ = 0.01, *α* = *β* = *γ* = 0.5; (**d**) *a*˜ = 0.05, *α* = *γ* = 0.2, *ψ* = 1; (**e**) *a*˜ = 0.05, *α* = *β* = 0.5, *ψ* = 1; (**f**) *a*˜ = 0.05, *α* = *β* = *γ* = 0.5; (**g**) *a*˜ = 0.1, *α* = *γ* = 0.2, *ψ* = 1; (**h**) *a*˜ = 0.1, *α* = *β* = 0.5, *ψ* = 1; (**i**) *a*˜ = 0.1, *α* = *β* = *γ* = 0.5. Subscript "c" stands for "centerline".

Next, let us fix *α* = *β* = 0.5, *ψ* = 1 varying *γ* and *a*˜. The corresponding plots for centerline velocity are shown in Figure 7b,e,h. Here we can observe switching from monotonically increasing to resonant behavior of the system with *a*˜ increasing and relative velocity amplitude decrease with increasing fractional order *γ*. Moreover, for larger values of *a*˜ centerline velocity behavior changes from resonant to almost linearly increasing with increasing fractional order *γ*.

Finally, we fixed *α* = *β* = *γ* = 0.5 and varied *a*˜ and *ψ* as shown in Figure 7c,f,i. As can be seen, the system behavior changes from monotonically increasing to resonant with increasing *a*˜. Moreover, resonant curves become less dispersive, with increasing *a*˜ and more dispersive with increasing *ψ*. Relative amplitudes of centerline velocity also decrease with increasing *ψ*.

The dynamics of shear stresses are given in Figure 8. The first column in Figure 8 (Figure 8a,d,g) corresponds to *α* = *γ* = 0.5, *ψ* = 1 and varying *a*˜ and *β*. In Figure 8g there exist two values of *ω*˜ where shear stresses become very close for different values of fractional order *β*. Upon reaching the one corresponding to lower value of *ω*˜ , the trend reverses with relative stress amplitude decreasing for increasing fractional order, *β*. For large values of *ω*˜ , shear stress becomes almost independent of the fractional order, *β*. For the second column of Figure 8 (Figure 8b,e,h), we set *α* = *β* = 0.5, *ψ* = 1 and varied *a*˜ along with *γ*. The new phenomenon observed can be seen in Figure 8e,h. That is, the resonant behavior changes to monotonically increasing and finally to constant with increasing fractional order *γ*. The similar trend is observed in case of varying *ψ* as shown in Figure 8c,f,i (*α* = *β* = *γ* = 0.5).

**Figure 8.** Shear stress at the pipe wall for fractional Zener model: (**a**) *a*˜ = 0.01, *α* = *γ* = 0.2, *ψ* = 1; (**b**) *a*˜ = 0.01, *α* = *β* = 0.5, *ψ* = 1; (**c**) *a*˜ = 0.01, *α* = *β* = *γ* = 0.5; (**d**) *a*˜ = 0.05, *α* = *γ* = 0.2, *ψ* = 1; (**e**) *a*˜ = 0.05, *α* = *β* = 0.5, *ψ* = 1; (**f**) *a*˜ = 0.05, *α* = *β* = *γ* = 0.5; (**g**) *a*˜ = 0.1, *α* = *γ* = 0.2, *ψ* = 1; (**h**) *a*˜ = 0.1, *α* = *β* = 0.5, *ψ* = 1; (**i**) *a*˜ = 0.1, *α* = *β* = *γ* = 0.5. Subscript "w" stands for "wall".

### *3.5. Fractional Poynting–Thomson Model*

The next model to be considered is fractional Poynting–Thomson one. Here *ξ* reads as:

$$\zeta^2 = \frac{\rho \omega^2}{E} \left[ \frac{1}{(E\_0/E)(i\omega \tau)^\gamma} + \frac{1}{\left((i\omega \tau)^a + (i\omega \tau)^\beta\right)} \right],\tag{31}$$

or alternatively in nondimensional form:

$$
\bar{\omega} = \bar{\omega}^2 \left( \frac{1}{\psi (i\bar{\omega})^\gamma} + \frac{1}{(i\bar{\omega})^\alpha + (i\tilde{\omega})^\beta} \right) \tag{32}
$$

Here we impose two additional conditions for fractional orders *α*, *β*, and *γ*: *α* − *γ* ≥ 0 and *β* − *γ* ≥ 0. As previously, we first examine the behavior of centerline velocity (Figure 9). Fixing *α* = *γ* = 0.2, *ψ* = 1, we vary *a*˜ and *β*. The trends are somewhat repeatable compared to fractional Maxwell model (Figure 3a,d,g) and demonstrate elastic behavior for the entire range of *β*. Next, we fixed *α* = *β* = 0.2, *ψ* = 1 and varied *a*˜ along with fractional order *γ*. The corresponding plots

are presented in Figure 9b,e,h. Here aperiodic oscillations for lower values of fractional order *γ*, oscillatory and monotonically increasing velocity profiles with increasing *γ* are observed. Finally, we set *α* = *β* = *γ* = 0.5 and varied *a*˜ and *ψ* (Figure 9c,f,i). Here, resonant behavior is observed even for lower values of *a*˜ (Figure 9c). Moreover, for *ψ* < 1 aperiodic oscillations of velocity are observed followed by it becoming independent from *ω*˜ . The trend, however, changes to resonant for *ψ* > 1. Overall, the velocity oscillations are damped slower with *ψ* increasing.

**Figure 9.** Centerline velocity for fractional Poynting-Thomson model: (**a**) *a*˜ = 0.01, *α* = *γ* = 0.2, *ψ* = 1; (**b**) *a*˜ = 0.01, *α* = *β* = 0.5, *ψ* = 1; (**c**) *a*˜ = 0.01, *α* = *β* = *γ* = 0.5; (**d**) *a*˜ = 0.05, *α* = *γ* = 0.2, *ψ* = 1; (**e**) *a*˜ = 0.05, *α* = *β* = 0.5, *ψ* = 1; (**f**) *a*˜ = 0.05, *α* = *β* = *γ* = 0.5; (**g**) *a*˜ = 0.1, *α* = *γ* = 0.2, *ψ* = 1; (**h**) *a*˜ = 0.1, *α* = *β* = 0.5, *ψ* = 1; (**i**) *a*˜ = 0.1, *α* = *β* = *γ* = 0.5. Subscript "c" stands for "centerline".

The shear stress dynamics at the wall are shown in Figure 10. As previously, we first fixed *α* = *γ* = 0.2, *ψ* = 1 and varied *a*˜ and *β* (Figure 10a,d,g). After major peak, shear stress drops dramatically with oscillations being damped for the entire range of *β*, except for Figure 10a, where it increases monotonically. The value of the major peak itself oscillates with increasing *β*. Next, as shown in Figure 10b,e,h, system parameters were set at *α* = *β* = 0.2, *ψ* = 1 with varying *a*˜ and *γ*. Monotonically increasing trend for lower values of *a*˜ changes to aperiodic oscillations and resonant for higher ones. Moreover, the value of the major peak decreases monotonically with increasing fractional order *γ*. Finally, we set *α* = *β* = *γ* = 0.5 and varied *a*˜ and *ψ*. The corresponding plots are shown in Figure 10c,f,i. As can be seen from it, resonant behavior is present for all values of *a*˜. However, it switches to monotonically increasing with increasing *ψ*. Moreover, the switch happens at higher values of *ψ* as *a*˜ increases and is absent for *a*˜ = 0.1 (Figure 10i).

**Figure 10.** Shear stress at the pipe wall for fractional Poynting-Thomson model: (**a**) *a*˜ = 0.01, *α* = *γ* = 0.2, *ψ* = 1; (**b**) *a*˜ = 0.01, *α* = *β* = 0.5, *ψ* = 1; (**c**) *a*˜ = 0.01, *α* = *β* = *γ* = 0.5; (**d**) *a*˜ = 0.05, *α* = *γ* = 0.2, *ψ* = 1; (**e**) *a*˜ = 0.05, *α* = *β* = 0.5, *ψ* = 1; (**f**) *a*˜ = 0.05, *α* = *β* = *γ* = 0.5; (**g**) *a*˜ = 0.1, *α* = *γ* = 0.2, *ψ* = 1; (**h**) *a*˜ = 0.1, *α* = *β* = 0.5, *ψ* = 1; (**i**) *a*˜ = 0.1, *α* = *β* = *γ* = 0.5. Subscript "w" stands for "wall".

#### *3.6. Fractional Burgers Model*

Here we discuss the most complex fractional viscoelastic model, Burgers model. For this model *ξ* is given by:

$$\mathcal{Z}^2 = \frac{\rho \omega^2}{E} \left[ \frac{1}{(E\_0/E)\left( (i\omega \tau)^\gamma + (i\omega \tau)^\delta \right)} + \frac{1}{\left( (i\omega \tau)^a + (i\omega \tau)^\beta \right)} \right],\tag{33}$$

or alternatively in nondimensional form:

$$
\omega = \tilde{\omega}^2 \left( \frac{1}{\psi[(i\tilde{\omega})^\gamma + (i\tilde{\omega})^\delta]} + \frac{1}{(i\tilde{\omega})^a + (i\tilde{\omega})^\delta} \right) \tag{34}
$$

It is worth noting that for fractional Burgers model the range of fractional orders differs from all other models considered above. In particular, while the range for orders *α*, *γ*, and *δ* is still from 0 to 1, 1 ≤ *β* ≤ 2. This change affects the dynamics of the centerline velocity profile dramatically. Centerline velocity profiles are shown in Figure 11. We first fixed *α* = *γ* = *δ* = 0.5, *ψ* = 1 and varied *a*˜ and *β*. The corresponding plots are shown in Figure 11a,d,g. The profile appears to be independent from the fractional order *β*. The situation changes when we set *β* = 1.5 and vary *δ* (Figure 11b,e,h). Velocity profile is very sensitive to the changes in fractional order *δ* when it has relatively low values. Aperiodic oscillations with reducing peak values for increasing values of *δ* are observed. For varying *ψ* (Figure 11 c,f,i), trends are overall similar to that for the fractional Poynting–Thomson model with the peaks for the fractional Burgers model being slightly sharper. This similarity can be attributed to the specific values of the fractional orders *α*, *γ* and *δ* set in Figure 11c,f,i (*α* = *γ* = *δ* = 0.5). This specific setting has consequences for constitutive equations in both models. That is, the terms containing fractional derivatives of the orders *α* − *γ* and *α* − *δ* become of the order zero.

**Figure 11.** Centerline velocity for fractional Burgers model: (**a**) *a*˜ = 0.01, *α* = *γ* = 0.5, *ψ* = 1; (**b**) *a*˜ = 0.01, *α* = *β* = 0.5, *ψ* = 1; (**c**) *a*˜ = 0.01, *α* = *β* = *γ* = 0.5; (**d**) *a*˜ = 0.05, *α* = *γ* = 0.5, *ψ* = 1; (**e**) *a*˜ = 0.05, *α* = *β* = 0.5, *ψ* = 1; (**f**) *a*˜ = 0.05, *α* = *β* = *γ* = 0.5; (**g**) *a*˜ = 0.1, *α* = *γ* = 0.5, *ψ* = 1; (**h**) *a*˜ = 0.1, *α* = *β* = 0.5, *ψ* = 1; (**i**) *a*˜ = 0.1, *α* = *β* = *γ* = 0.5. Subscript "c" stands for "centerline".

Finally, let us have a look at the shear stresses shown in Figure 12. As can be seen from it, shear stress, too, is independent from the fractional order *β* (Figure 12a,d,g). Stress dynamics change in Figure 12b,e,h. Quickly decaying aperiodic oscillations are observed. In particular, major peak values decrease monotonically with increasing value of fractional order *δ*. Moreover, these peaks become sharper and shift to lower values of *ω*˜ with increasing *a*˜. Next, we fixed *α* = *γ* = *δ* = 0.5 and *β* = 1.5. Here we observed resonant behavior for all values of *a*˜ and quickly damped oscillations for low values of *ψ* (Figure 12c,f,i).

A closer look at expressions for *ξ* corresponding to different models reveals their similarity. To illustrate this statement, let us introduce complex compliance, *J*∗(*ω*). For all the models considered, it has the following functional form: [...]/*E*, where an expression in square brackets is model-specific. Both *ξ* and *ω*¯ can be expressed in terms of complex compliance:

$$\xi^2 = \rho \omega^2 I^\*(\omega), \quad \bar{\omega} = \tilde{\omega}^2 E f^\*(\omega), \tag{35}$$

as it was defined in [39]. Thus, generally speaking, *ξ* is a complex quantity. Moreover, as we got complex compliance, other physical quantities including creep compliance, complex modulus, and relaxation modulus can be restored. For instance, *J*∗(*ω*) = 1/*G*∗(*ω*), where *G*∗(*ω*) stands for complex modulus. Then for fractional Maxwell model the loss modulus, *G*"(*ω*) = Im(*G*∗(*ω*)), reads as [39]:

$$G''(\omega) = E \frac{(\omega \tau)^a \sin(\pi a/2) + (\omega \tau)^{2a - \beta} \sin(\pi \beta/2)}{1 + (\omega \tau)^{2(a - \beta)} + 2(\omega \tau)^{a - \beta} \cos(\pi (a - \beta)/2)} \tag{36}$$

Another insight can be obtained upon examination of Figure 2b–f and corresponding constitutive equations. The general observation is that more complex models can be obtained via consecutive combination of simpler ones with fractional elements. In particular, by adding fractional element to fractional Maxwell model in parallel, we arrive at fractional Zener model. Alternatively, starting from fractional Kelvin–Voigt model and adding fractional element in series, we get fractional Poynting–Thomson model. Adding one more fractional element in series, we end up with fractional Burgers model. The same "additive" behavior is observed for complex compliances.

**Figure 12.** Shear strees at the pipe wall for fractional Burgers model: (**a**) *a*˜ = 0.01, *α* = *γ* = 0.5, *ψ* = 1; (**b**) *a*˜ = 0.01, *α* = *β* = 0.5, *ψ* = 1; (**c**) *a*˜ = 0.01, *α* = *β* = *γ* = 0.5; (**d**) *a*˜ = 0.05, *α* = *γ* = 0.5, *ψ* = 1; (**e**) *a*˜ = 0.05, *α* = *β* = 0.5, *ψ* = 1; (**f**) *a*˜ = 0.05, *α* = *β* = *γ* = 0.5; (**g**) *a*˜ = 0.1, *α* = *γ* = 0.5, *ψ* = 1; (**h**) *a*˜ = 0.1, *α* = *β* = 0.5, *ψ* = 1; (**i**) *a*˜ = 0.1, *α* = *β* = *γ* = 0.5. Subscript "w" stands for "wall".

It is worth emphasizing the generality of the approach developed. Provided the flow is axisymmetric, the pipe having a circular cross-section, and no-slip condition at the wall is satisfied, the solutions derived are suitable for various viscoelastic models commonly used for practical applications. For the same form of the solution of velocity and stress profiles, the prediction of flow dynamics for completely different fluids can be obtained via deriving a single model parameter, *ξ*. Moreover, the functional form of *ξ* itself is also identical for all the models considered. What makes a given model specific is the expression for the complex compliance. This powerful approach simplifies significantly the implementation of other possible viscoelastic models within a given framework. The classical viscoelastic and Newtonian fluid models present the limiting cases of the theoretical approach proposed in this study. Velocity profiles' dynamics for fractional Maxwell and Newtonian fluids in a circular pipe are shown in Supplementary Materials 1. These clearly demonstrate the complexity of behavior for fractional Maxwell model.

#### **4. Conclusions**

We have obtained exact analytical solutions for velocity profiles and shear stresses. Fractional Maxwell, Kelvin-Voigt, Zener, Poynting–Thomson, and Burgers models were considered. We demonstrated that the same form of the solution is applicable to all the models considered, thus generalizing prior studies. For both centerline velocity and shear profiles, three types of behavior (monotonically increasing, resonant, oscillating aperiodic) have been identified. In addition to it, switching between trends of relative velocity amplitudes have been predicted. Monotonically decreasing trends were found to be

typical for relatively low values of normalized pipe radius with more complex behaviors taking place for higher values of the normalized radius. In addition to this, centerline velocity profiles featured almost constant plateaus. The proposed models cover the widest range of viscoelastic materials both in terms of the balance between viscous and elastic properties (via fractional-order) and the ratio of elastic properties for complex materials.

The approach we developed features an advantage that is worth mentioning. It allows the wide range of fractional viscoelastic models to be represented and solved in the same functional form. In fact, the entire system behavior can be described by a single function *ξ* (or *ω*¯ ). It, in turn, is related to a complex compliance with the functional form identical for all five models considered in this study. Knowing complex compliance, complex and relaxation moduli and creep compliance are readily obtained. These quantities can be directly measured thus providing immediate practical applications for materials analysis. However, for some specific fields (e.g., polymers), experimental creep compliance data for short times has not been well-established so far. This problem was outlined in recent studies (see, for example, Reference [74]). It definitely deserves attention of the research community and will hopefully be addressed in future studies.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-3417/10/24/9093/ s1 .

**Author Contributions:** Conceptualization, D.G.; methodology, D.G.; software, D.G.; validation, R.P.; formal analysis, D.G.; investigation, D.G.; resources, R.P.; writing—original draft preparation, D.G.; writing—review and editing, R.P.; funding acquisition, R.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Argonne National Laboratory through grant #ANL 4J-303061-0030A "Multiscale modeling of complex flows".

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