Fractional Burgers Fluid Flow Due to Metachronal Ciliary Motion in an Inclined Tube

Cilia-induced flow of fractional Burgers fluid is studied in an inclined tube for both symplectic and antiplectic wave patterns. The solution of the problem is persued under the long wave length limitation. The fractional Adomian decomposition method is employed to evaluate the pressure gradient. Mathematical expressions for the axial velocity, frictional force, pressure gradient, and stream function are obtained and the influence of the main operating parameters is discussed in detail. It is noted that the velocity profile is more dominant in the case of antiplectic metachronal waves compared to symplectic ones, which confirms former results on the better capability of antiplectic waves to transport mucus, obtained with more complex numerical solvers.


INTRODUCTION
Dutch light microscopist Antoni was the first to discover cilia in 1675 and Sharpey was the first to discuss cilia in English language in 1835. Exhaustive studies on the ciliary structures have been carried out during the nineteenth century (Sleigh et al., 1988). Cilia and flagella oscillate in a waving fashion during motion to transport fluids and propel cells (Vélez-Cordero and Lauga, 2013). Cilia motion constitutes a pivotal role in a wide variety of physiological processes, such as alimentation, circulation, locomotion, respiration, and reproduction (Maiti and Pandey, 2017). Ciliated cells can indeed be found in many human organs e.g., photoreceptor cells in retina, in hair bundles on ear, epithelial cells in the respiratory tract, in Fallopian tubes (Ashraf et al., 2018), in kidney (Guirao and Joanny, 2007), in the ependymal cells of the brain that generate cerebrospinal flow among few examples (Lieberstein, 1975). Malfunctioning of the ciliary activity may be responsible of many respiratory diseases (Rubin, 2014), like severe asthma.
As the amalgamated motion of cilia occurs, the cilia's upper layer can be seen to have a metachronal wave generated as a result of a small phase lag between neighboring cilia. This collective motion of cilia supports many physiological processes (Lieberstein, 1975;Murakami and Takahashi, 1975;Takahashi and Shingyoji, 2002). Different types of metachronal waves are classified depending upon the dynamics and strokes of the moving cilia. Symplectic beat patterns are produced if the directions of the propagative metachronal waves and the main flow are the same and antiplectic patterns are recognized when the directions of the wave propagation and the flow are opposite (Knight-Jones, 1954;Blake, 1972).
To date, many efforts have been made mathematically to understand the involvement of ciliary motion related to different biofluids, such as bronchial mucus, semen . . . in humans which are represented by Newtonian and non-Newtonian fluid models. Maiti and Pandey (2017) applied a numerical approach to study the flow of a power-law fluid representing semen in an axisymmetric tube (representing the efferent ducts of the male human reproductive channel) due to cilia motion and concluded that, more than only the ciliary activity, there are several other factors like the smoothness of muscles, the constant fluid secretion or the vacuum created during ejaculation, which may be responsible for the semen flow. Siddiqui et al. (2010Siddiqui et al. ( , 2014 obtained exact solutions of cilia induced flow problems of viscous fluid and power-law fluid models representing semen in a cylindrical tube and infinite channel under the long wavelength approximation. Mucus plays a significant role in individual's health, therefore many researchers discussed the bronchial mucus transportation due to ciliary activity (see for examples in Barton and Raynor, 1967;Ross and Corrsin, 1974;Fulford and Blake, 1986;Maqbool et al., 2016). Norton et al. (2011) developed a transport model, where the mucus, considered successively as a Doi-Edwards, Jeffrey, and Maxwell fluid, is transported as a rigid body and the metachronal wave exhibits a symplectic behavior. Vélez-Cordero and Lauga (2013) applied the regular perturbation method to solve a problem in which the tracheobronchial mucus is considered as a Carreau fluid. This solution only portrays the Newtonian effects when second-order perturbations are considered and non-Newtonian effects are captured when the perturbation analysis is pushed up to the fourth-order. Maqbool et al. (2016) considered the geometry of Siddiqui et al. (2014) and studied the mucus flow treating mucus as a Jeffrey fluid. Smith et al. (2008) noted that the most advanced model to investigate the mucociliary clearance process is the Maxwell viscoelastic model and proposed a fluid-structure interaction model to examine the complex fluid flow problem arising due to ciliary activity.
More recently, the immersed boundary (IB) method (Hao and Zhu, 2010) has been extensively applied to study biofluid flows representative of humans/animals. Dillon et al. (2007) used the IB method to examine the two dimensional flow of three cilia in a mucus layer such that the mucus layer is treated as an elastic solid instead of a viscoelastic fluid. Dauptain et al. (2008) used the IB method to examine the motion of fluid due to one row of cilia on a ctenophore Pleurobrachia pileus, which is commonly known as a sea gooseberry for Reynolds numbers Re within the range [50 − 200]. They found that as the beating of cilia increases, it spreads more power to the interacting fluid and this work may be considered as a guideline for solving the fluid-structure interaction problem. Very recently, Chatelin and Poncet (2016) investigated by 3D simulations the influence of the mucus viscosity, fluid height, cilia length, and beating frequency on the mucociliary process in a two-phase environment. Chateau et al. (2018) performed 3D simulations of the transport and mixing induced by beating cilia at Re up to 20 in a two-phase environment composed of Newtonian fluids using a coupled IB/lattice Boltzmann method.
All reported studies (Barton and Raynor, 1967;Ross and Corrsin, 1974;Fulford and Blake, 1986;Dillon et al., 2007;Smith et al., 2008;Vélez-Cordero and Lauga, 2013;Maqbool et al., 2016) confirmed the fact that more efforts are still required in terms of scientific research to better understand the internal flow structure due to cilia motion and their interaction with the surrounding fluid. In the present work, one will focus on the mucociliary clearance process. It is well-known now that bronchial mucus exhibits complex rheological properties: stress relaxation, tensile stresses, shear thinning, yielding stress, and thixotropic behavior (see in Lafforgue et al., 2017Lafforgue et al., , 2018. Though the advanced numerical solvers developed by Chatelin and Poncet (2016) and Chateau et al. (2018) considered the two-phase character of the problem and the behavior of each single cilium, such approaches do not account for the rheology of mucus. Moreover, simulating the mucociliary clearance process in human airways remains very challenging for such methods due to the multiscale character of the problem: from the micrometer scale when considering each individual cilium to the decimeter scale when looking at the main air flow within the trachea. So analytical approaches or simplified models like the envelope model are still deemed necessary if one wants to simulate the complete problem. The present paper is an attempt to demonstrate that analytical solutions obtained by the Adomian decomposition method can provide useful informations regarding the mucociliary clearance process. Momani and Odibat (2006) used successfully the Adomian decomposition method to solve a time-fractional Navier-Stokes equation in a tube and demonstrated both the efficiency and simplicity of its method. In this paper, bronchial mucus is considered as a fractional Burgers fluid. Motion is generated by linear pressure produced by the tips of the moving cilia under the long wavelength and low Reynolds number approximations (Shapiro et al., 1969). Various illustrations highlighting the effects of the most important parameters are also sketched. Another motivation of the present paper is that the literature is scarce on fractional fluid models (see the monograph of Oldham and Spanier, 1974 for example).

MATHEMATICAL MODEL
The fluid motion characteristics of an incompressible fractional Burgers fluid in a ciliary tube having an inclination angle θ are considered. The metachronal wave and inclined tube move with the same speed c to the right as shown in Figure 1.
The cilia tips follow an elliptical path as suggested by Maqbool et al. (2016) and Siddiqui et al. (2010Siddiqui et al. ( , 2014, which can be represented by where Equations (1, 2) are the parametric equations representing the cilia motion in which a is the mean tube radius, c is the wave speed, t is time, Z 0 is the reference position of cilia, α * is the eccentricity of ellipse, ǫ is a dimensionless varying parameter and η is the wavelength of the metachronal wave. One assumes that the fluid and adjacent cilia tips have the same velocity resulting in a no-slip condition. Therefore, the axial and radial velocities can be written as and Introducing Equations (1, 2) into Equations (3, 4) yields and where U and W are the radial and axial velocity components. The wave and fixed frames are related through the transformations and P( Z, R, T) = p( z, r), where P p , R ( r), and Z ( z) are the pressure, radius, and axial position in the fixed (wave) frame. The extra stress tensor for fractional Burgers' fluid is given as where Ŵ is the shear stress tensor, · ̹ is the strain rate, λ i (i = 1, 2, 3) are constitutive parameters, α and β are the fractional derivative and integral defined as Hilfer (2000) D (10) and with 0 ≤ α ≤ 1. The fractional Burgers fluid model reduces to fractional Oldroyd B model for λ 2 = 0 and the classical viscous model can be obtained if λ i = 0 (i = 1, 2, 3) . The velocity components for the fractional Burgers' fluid model for an inclined tubular flow should satisfy the following equations where g is the gravitational acceleration, u and w are the radial and axial velocity components, τ ij are the shear stress tensor components and ρ is the fluid density. To solve the problem, the following non-dimensional parameters are introduced where η, a, c, and µ denote the wavelength, tube radius, wave speed, and dynamic viscosity, respectively, Re and Fr are the Reynolds and Froude numbers while β * is the wave number. Under the long wavelength and low Reynolds number approximations, the flow may be considered as a Stokes flow. Thus, during the non-dimensionalizing of Equations (12, 13) with the help of Equation (15), we have ignored the terms involving β * , β * 2 , β * 3 , and Re, Re 2 , Re 3 . . . but terms involving Re/Fr are retained as the orders of Re and Fr numbers are the same. Equations (8-14) (after dropping hats) simplify to the following non-dimensional forms with appropriate boundary conditions given as where The following two equations are used to determine the initial guesses required by the Adomian decomposition method Integrating Equation (17) with respect to r and using the boundary condition (Equation 21), one gets Further integrating Equation (24) and applying the boundary conditions (Equations 18, 19) will yield The results for a Newtonian fluid in a horizontal tube can be deduced using the limits λ i → 0 for i = 1, 2, 3 and θ → 0, therefore one gets The volume flow rate is defined as which, in light of Equation (24), becomes so that q the dimensional volume rate and Q the dimensionless volume flow rate in the fixed frame are related as The mean volume flow rate Q can be calculated using the time period T in Equation (30) Equation (26) in view of Equation (29) gives The stream function ψ in the wave frame is computed with the help of Equations (25, 28, 31) as

SOLUTION METHODOLOGY
Equation (31) can be simplified as where l(z, t) = ∂p ∂z with the initial conditions and so that Equation (32) will take the form An infinite series solution for l (z, t) by using the Adomian decomposition method (Adomian, 1994;Babolian and Biazar, 2002) is given by where l 0 , l 1 , l 2 , l 3 ,...l n+1 are determined as From l n (z, t) (n 0), the other components can also be obtained.
Finally an approximate solution of Equation (37) by truncating the series can be written as where φ n (z, t) = N−1 n=0 l n (z, t).
The pressure difference p and friction force F across one wavelength are given by

RESULTS AND DISCUSSION
In order to evaluate the pressure rise p, pressure gradient dp/dz, frictional force F, and streamlines ψ, the software Mathematica 8.0 has been used. The flow characteristics of Burgers fluids  in an inclined ciliated tube are presented by controlling the fractional parameters α and β, the cilia length ǫ, and the angle of inclination θ . The Reynolds number Re and the Froude number Fr have been both fixed to 0.1. The Reynolds number based on the cilia tip speed is usually around 10 −5 in human airways but Chateau et al. (2018) recently demonstrated that there is no significant influence of Re as long as it remains lower than 1.
The variations of p withQ are examined in Figure 2 for different values of α, β, ǫ, and θ . It is found that the pressure rise p decreases with an increase of the fractional parameter α, cilia length ǫ in the pumping expanse region ( p > 0) of the tube and an opposite trend is noted in the copumping expanse region ( p < 0) of the tube. Also in the pumping (resp. copumping) expanse region, the pressure rise increases (resp. decreases) for increasing values of the fractional parameter β. As expected increasing the cilia length or decreasing the tube radius provide similar results in terms of pressure variations but the effect is more pronounced by changing the tube radius. The ratio β/ǫ may be the relevant geometric parameter governing the transport of fluid by the mucociliary clearance process. The influence of the inclination angle θ on p is more straightforward as shown in Figure 2D. p increases by increasing θ in the pumping expanse region and an opposite trend is reported in the copumping expanse region. Larger pressure differences are observed for a vertical tube.
The variations of the pressure gradient dp dz are examined in Figure 3 for different values of α, β, ǫ, and θ . It is noted that dp dz remains small and the fluid can flow smoothly without the application of a large pressure gradient in the expanse regions 0 ≤ z ≤ 0.2 and 0.8 ≤ z ≤ 1. On the other hand, for 0.2 < z < 0.8, a large amount of pressure gradient is required to maintain the flow. The magnitude of the pressure gradient increases by increasing the parameters β, ǫ, and θ . These parameters provide the resistive force to the flow thus a larger value of dp dz is required to maintain the fluid flow, whereas the parameter α provides the deriving force to the flow thus a smaller value of dp dz is required to maintain the fluid flow.
It is observed through Figure 4 that the frictional force F varies linearly by increasing α, β, and ǫ and it tends to increase in magnitude by increasing β and ǫ. The magnitude of the resistive force decreases by increasing α, in the expanse regionQ > 0.5. On the contrary, its magnitude increases in the same region when increasing β and/or ǫ. The magnitude of the frictional force increases in the expanseQ < 1 by increasing the tube inclination θ but an opposite trend is noted forQ > 1. The inclination of the tube provides then a large magnitude of F to oppose the flow. Figure 5 displays the streamline patterns and trapping for α * = β * = 0.7 andQ = 1.5 by increasing the dimensionless cilia length ǫ. The center line symmetry bifurcates the boluses of the fluid particles circulating along the closed stream lines. The boluses are confined and move with the velocity of the metachronal waves. In Figures 5A-C, keeping all the other parameters constant, the number and size of the boluses increase by increasing ǫ. Thus, the cilia length significantly affects the flow dynamics by generating boluses in the inclined tube.   Figure 6 displays a 3D view of the axial and radial velocity components for ǫ = 0.2, α = 0.7, β = 0.7, andQ = 1.5. It is observed that the waves are elliptic waves and the fluid velocity is in both forward and backward directions with a symmetric behavior due to metachronism of the ciliary motion. Figure 7 displays the radial distributions of the axial and radial velocity components for four flow ratesQ and both antiplectic and symplectic waves formed by the cilia tips. As expected, both components increase by increasing the flow rate. More interestingly, at any value ofQ, the axial velocity produced by antiplectic waves is larger than the one due to symplectic waves. The increase in the velocity profiles is also faster for antiplectic waves compared to symplectic waves. Regarding the radial velocity component, the two types of metachronal waves provide symmetric profiles in the radial direction. All in all, it confirms the former results of Chateau et al. (2018) obtained using a more advanced numerical model but for a two Newtonian fluid configuration. For a single layer of fractional Burgers fluid, antiplectic waves are also more efficient to transport fluid than symplectic waves.

CONCLUSIONS
In this paper, fractional Burgers' fluid flow in an inclined ciliated tube is examined. Using the long wavelength approximation, a semi-analytic solution is developed. Frictional force, pressure rise, pressure gradient and streamlines are plotted for different values of the main operating parameters and the main results can be summarized as follows: • p decreases by increasing α and ǫ in the pumping region and an opposite trend is noted in the copumping region of the tube. • p increases in the pumping region and decreases in the copumping region when increasing β. • dp dz is insignificant for 0 ≤ z ≤ 0.2 and 0.8 ≤ z ≤ 1. On the contrary, for 0.2 < z < 0.8, a higher value of dp dz is required to maintain the flux.
• The magnitude of dp dz increases by increasing β, ǫ, and θ .
• It is noticed that the parameters β, ǫ, and θ greatly influence the pressure gradient whereas the parameter α provides a smaller amount of the pressure gradient to the fluid flow. • The frictional force varies linearly when increasing α, β, and ǫ and its magnitude increases by increasing β and ǫ but it decreases in magnitude by increasing α and θ in the pumping region. • The magnitude of the frictional force increases in the region Q < 1 by increasing the tube inclination θ but an opposite trend is noted in the regionQ > 1.
• The magnitude of the pressure difference is larger for fractional generalized Burgers model in comparison to the generalized Burgers model forQ < 0.6 and a reverse trend is observed forQ > 0.6. • A greater magnitude of the pressure gradient is needed to the fluid flow for fractional generalized Burgers fluid as compared with the generalized Burgers fluid to pump the same amount of fluid. • The fractional generalized Burgers fluid provides a greater amount of frictional force in comparison to the generalized Burgers fluid. • The number and size of boluses increase by increasing ǫ.
• The magnitude of the velocity profile increases significantly for antiplectic waves as compared to symplectic waves confirming that they are more efficient to transport mucus.
The fractional Adomian decomposition method appears then as a valuable tool to get useful and realistic results for the mucociliary transport process. Such method provides much faster results compared to more complex solvers (Chatelin and Poncet, 2016;Chateau et al., 2018) and could be associated to 3D flow solvers to solve the multiscale problem of the mucus clearance in the airways.

AUTHOR CONTRIBUTIONS
This work is the product of intellectual effort of all authors. They read and approved the final draft of the manuscript. KM and AM suggested and formulated the problem. KM, SS, and AM solved the problem jointly. SS wrote the paper. SP helped in data analysis and reviewed the manuscript.

FUNDING
SP acknowledges the Natural Sciences and Engineering Research Council of Canada for its financial support through a discovery grant (RGPIN-2015-06512).