**Mathematical Analysis on an Asymmetrical Wavy Motion of Blood under the Influence Entropy Generation with Convective Boundary Conditions**

#### **Arshad Riaz 1, Muhammad Mubashir Bhatti 2, Rahmat Ellahi 3,4,***∗***, Ahmed Zeeshan <sup>3</sup> and Sadiq M. Sait <sup>5</sup>**


Received: 13 December 2019; Accepted: 3 January 2020; Published: 6 January 2020

**Abstract:** In this article, we discuss the entropy generation on the asymmetric peristaltic propulsion of non-Newtonian fluid with convective boundary conditions. The Williamson fluid model is considered for the analysis of flow properties. The current fluid model has the ability to reveal Newtonian and non-Newtonian behavior. The present model is formulated via momentum, entropy, and energy equations, under the approximation of small Reynolds number and long wavelength of the peristaltic wave. A regular perturbation scheme is employed to obtain the series solutions up to third-order approximation. All the leading parameters are discussed with the help of graphs for entropy and temperature profiles. The irreversibility process is also discussed with the help of Bejan number. Streamlines are plotted to examine the trapping phenomena. Results obtained provide an excellent benchmark for further study on the entropy production with mass transfer and peristaltic pumping mechanism.

**Keywords:** convection; entropy production; heat transfer engineering; blood flow

#### **1. Introduction**

In our daily life, living organisms require energy to do physical work and keep the body temperature under the influence of heat exchange to the environment, as well as to generate, replace, and propagate molecules to the relevant constituents. Such type of energy comes from the oxidation process of organic substances i.e., amino acids, fats, and carbohydrates fed to the organisms. As compared to the other heat engines (i.e., in which the chemical energy gets transformed to the thermal energy, and then is transformed to mechanical work), living organisms can transform the nutrient's chemical energy into work. It happens due to the oxidation of nutrients located internally in the organisms i.e., metabolism, pass through different steps, which helps to hold some energy from ATP (adenosine triphosphate). The ATP utilized entirely by all beings for the direct transformation of mechanical energy and also actively supports other biological reactions [1]. In recent years, various authors [2–6] have examined the heat production of mammals via calorimetry, and presented that for the given nutrients, both combustion and animal metabolism expends the same amount of oxygen. According to previous research [7], it is found that living things can produce thermal energy via fat, and combustion of carbohydrates in the living body, and is identical to the oxidation of heat of these elements. As a result, the amount of nutrients digested by a living being, and hence its input energy, can be determined by the chemical composition of food intake and the measurements of breathing i.e., CO2 and O2.

Hershey [8] and Hershey and Wang [9] examined the entropy production during the lifespan of a human being. They found that when the human body is in a state of rest, mostly the output energy due to the nutrient's metabolism occurs as a heat. They also reinforced the calorimeter to examine the heat transfer rate to the environment and verified that with the help of BMR (basal metabolic). According to their results, it was found that entropy generated over a lifespan was 10,678 KJ/kg·K for females, and 10,025 KJ/kg·K for males. Rahman [10] discussed the entropy generation for forced and free convection using a new mathematical model. He discussed forced and free convection at distinct mass influx, outflux (i.e., waste, air, food, and water, etc.), level of physical activity, clothing effects, and airspeeds. His results are similar to those from Hershey [8] and Aoki [11] but one order of magnitude higher in general. Annamalai and Puri [12] utilized the first law of thermodynamics to achieve the metabolic scaling for a biological system. They also used the second law of thermodynamics to determine the entropy generation in humans and prognosticate the lifespan of 77 years by assuming the maximal entropy generation as 10,000 KJ/kg·K. Bejan [13,14] introduced a constructal design principle and presented optimal geometric types scales to the power of their associated size and showed that different natural structures (i.e., lightening, river deltas, tree branches, and vascularized tissue) are periodic in nature. Rashidi et al. [15] discussed the entropy generation with magnetic effects and slip boundary conditions propagating among an infinite porous disk having variable features. According to their results, they observed that the disk is an essential root of entropy generation. Komurgoz et al. [16] examined the entropy generation through an inclined porous channel with magnetic effects. According to their results, they found that maximal entropy production can be gained in the absence of porosity and magnetic field.

Blood is an essential part of the human body, which comprises 7% of the total body weight. The leading role of blood is molecular oxygen for cellular metabolism and carry the nutrients as well as a significant role in thermoregulatory. Blood performs as a non-Newtonian fluid. The blood viscosity changes due to the shear rate. The viscosity of the blood can be analyzed by the hematocrit, plasma (constitute 54.3% of the whole blood) viscosity, and the mechanical features of red blood cells (constitute 45% of the whole blood). The human blood is a heterogeneous solution that contain multiple kinds of cells (known as corpuscles or formed elements), which consist of leukocytes, thrombocytes, and erythrocytes. In view of such importance, different authors discussed the entropy generation in blood. For instance, Rashidi et al. [17] discussed the magnetic effects on the blood flow propagating through a porous medium with a filtration and control process. Akbar et al. [18] examined the thermal conductivity on the peristaltic propulsion of H2+Cu nanofluids with entropy production. Rashidi et al. [19] obtained the series solution for the entropy generation of the blood flow of a nanofluid in the presence of a magnetic field. Endoscopic effect and entropy production on peristaltic nanofluid flow, having a thermal conductivity of 2 HO were investigated by Akbar et al. [20]. Abbas et al. [21] presented a detailed analysis of the peristaltic flow with nanofluids and entropy production through a finite channel with compliant walls. Bhatti et al. [22] considered the Casson blood flow to examine the entropy process with peristaltic movement under the uniform magnetic field. Ranjit and Shit [23] examined the entropy production on the electroosmotic flow under uniform magnetic field with peristaltic pumping. More studies on the blood flow and entropy generation can be found from the references [24–28].

According to the above survey, it is found that less attention has been given to the entropy production asymmetric peristaltic propulsion of blood flow with heat transfer. Therefore, in the present analysis, we discuss the entropy generation with convection on the asymmetric propulsion of the peristaltic blood of nonlinear Williamson fluid. An assumption of long peristaltic wavelength

is taken into account and Reynolds number is considered to be very small (Re ≈ 0). A regular perturbation method is used to obtain series solutions. The novelty of all the leading parameters is discussed and illustrated. The trapping mechanism is also examined to determine the nonlinear asymmetric peristaltic motion.

#### **2. Governing Equations**

In this section we analyze the incompressible peristaltic propulsion of Williamson fluid in a two-dimensional channel with a width *d*<sup>1</sup> + *d*2. The flow is initialized by a sinusoidal wave propagating with a constant speed *c* along the layout of channel (see Figure 1). The addition here is the extra equations of energy and entropy generation. It is assumed that the temperature at the upper wall of the channel is *T*<sup>1</sup> and lower wall has the temperature *T*<sup>0</sup> such that *T*<sup>0</sup> < *T*1. It depicts the physical reasoning that heat will transfer from lower to upper wall. The wall surfaces are suggested as:

$$\overline{Y} = \begin{cases} \begin{array}{c} H\_i = \overline{d\_i} + \overline{a\_i}\cos\left[2\pi\omega\omega\right], \text{ i } i = 1, \\ H\_j = -\overline{d\_j} - \overline{b\_i}\cos\left[2\pi\omega + \phi\right], \text{ } j = 2, \end{array} \tag{1}$$

and

$$
\omega = \frac{\overline{X} - c\overline{t}}{\lambda},
\tag{2}
$$

where *a*<sup>1</sup> and *b*<sup>1</sup> are the wave amplitudes, *λ* the wave length, *t* the time, *c* the velocity of the propagation, and *X* is the direction of wave propagation. The phase difference *φ* has the range 0 ≤ *φ* ≤ *π* i.e., waves out of phase *φ* = 0 associated to the symmetric channel, and *φ* = *π* associated to the waves are in phase. Moreover, *ai*, *bi*, *di*, *dj*, and *φ* satisfy the condition:

$$
\left(\overline{\pi\_i}^2 + \overline{b\_i}^2 + 2\overline{\pi\_i}\overline{b\_i}\cos\phi \le \left(\overline{d\_i} + \overline{d\_j}\right)^2. \tag{3}
$$

The equations of momentum in component forms are described as:

$$
\rho \left[ D\_{\overline{Y}} + \overline{U} D\_{\overline{X}} + \overline{V} D\_{\overline{Y}} \right] \overline{U} = -D\_{\overline{X}} P + D\_{\overline{X}} \overline{S}\_{\overline{X}\overline{X}} + D\_{\overline{Y}} \overline{S}\_{\overline{X}\overline{Y}} \tag{4}
$$

where *Dt* = *<sup>∂</sup> ∂t* , *DX* <sup>=</sup> *<sup>∂</sup> <sup>∂</sup><sup>X</sup>* , *DY* <sup>=</sup> *<sup>∂</sup> ∂Y* .

$$\rho \left[ D\_{\overline{Y}} + \overline{\mathcal{U}} D\_{\overline{X}} + \overline{\mathcal{V}} D\_{\overline{Y}} \right] \overline{\mathcal{V}} = D\_{\overline{Y}} P + D\_{\overline{X}} \overline{\mathcal{S}}\_{\overline{Y}\overline{X}} + D\_{\overline{Y}} \overline{\mathcal{S}}\_{\overline{Y}\overline{Y}}.\tag{5}$$

The stress tensor for the Williamson fluid model reads as:

$$\overline{\mathbf{S}} = \left[ \overline{\mu}\_{\infty} - (\overline{\mu}\_{\infty} - \overline{\mu}\_{0}) \left( 1 - \overline{\Gamma} \dot{\gamma} \right)^{-1} \right] \dot{\gamma} \tag{6}$$

where *<sup>μ</sup>*∞, *<sup>μ</sup>*<sup>0</sup> the infinite and zero shear rate viscosity, <sup>Γ</sup> the time constant, and . *γ* reads as:

$$\dot{\gamma} = \left(\frac{1}{2} \sum\_{m} \sum\_{n} \dot{\gamma}\_{mn} \dot{\gamma}\_{mn}\right)^{\frac{1}{2}} = \left(\frac{1}{2} S\_t\right)^{\frac{1}{2}},\tag{7}$$

where *St* is the second invariant strain tensor. For the present flow problem, we considered *μ*<sup>∞</sup> = 0 (the infinite shear rate viscosity is very small as compared to zero shear rate viscosity) and <sup>Γ</sup> . *γ* < *i* i.e. *i* = 1. Then, Equation (6) takes the following form:

$$\overline{\mathbf{S}} = \overline{\mu}\_0 \left( 1 - \overline{\Gamma} \dot{\gamma} \right)^{-1} \dot{\gamma}. \tag{8}$$

*Symmetry* **2020**, *12*, 102

The energy equation to represent the heat exchange in the channel is as stated below. The law of conservation of energy in the dimensional mathematical pattern is given by:

$$S\_h \left[ D\_{\overline{l}} + \overline{l} \overline{l} D\_{\overline{X}} + \overline{V} D\_{\overline{Y}} \right] T = \frac{K}{\rho} \left[ D\_{\overline{X}\overline{X}} T + D\_{\overline{Y}\overline{Y}} T \right] + \frac{\overline{S}\_{\overline{X}\overline{Y}}}{\rho} D\_{\overline{Y}} \overline{l} \overline{l}. \tag{9}$$

In the above equation, *Sh* is the specific heat coefficient, *K* the thermal conductivity, and *ρ* the density of the governing fluid.

Introducing wave frame coordinates transformations with propagation velocity *c* away from the fixed frame read as:

$$\{\overline{x} + c\overline{t}, \overline{u} + c, \overline{y}, \overline{v}, \overline{P}(\overline{x})\} = \{\overline{X}, \overline{U}, \overline{Y}, \overline{V}, \overline{P}\left(\overline{X}, t\right)\}\tag{10}$$

Defining the dimensionless quantities as:

$$x = \frac{\overline{x}}{\lambda}, y = \frac{\overline{y}}{\overline{d}\_i}, u = \frac{\overline{u}}{\lambda}, v = \frac{\overline{v}}{\lambda}, S\_{\overline{x}x} = \frac{\lambda}{\overline{\mu}\_0 c} \overline{S}\_{\overline{\nabla}\overline{x}}, S\_{xy} = \frac{\overline{d\_i}}{\overline{\mu}\_0 c} \overline{S}\_{\overline{\nabla}\overline{x}}, S\_{yy} = \frac{\overline{d\_i}}{\overline{\mu}\_0 c} \overline{S}\_{\overline{\nabla}\overline{x}},$$

$$\theta = \frac{T - T\_0}{T\_i - T\_0}, p = \frac{\overline{d\_i}^2}{c\lambda \overline{\mu}\_0} P, \dot{\gamma} = \frac{\overline{\dot{\gamma} d\_i}}{c}, a = \frac{\overline{a\_i}}{\overline{d\_i}}, b = \frac{\overline{b\_i}}{\overline{d\_i}}, d = \frac{\overline{d\_j}}{\overline{d\_i}}, h\_{i,j} = \frac{H\_{i,j}}{\overline{d\_i}},\tag{11}$$

where *θ* is the dimensionless temperature profile.

By invoking the above transformations in Equations (4)–(6), we arrive at (after ignoring the bars):

$$\operatorname{Re}\left[\delta u D\_{\mathbf{x}} u + v D\_{\mathbf{y}} u\right] = -D\_{\mathbf{x}} p + D\_{\mathbf{x}} \mathbf{S}\_{\mathbf{x}\mathbf{x}} + D\_{\mathbf{y}} \mathbf{S}\_{\mathbf{xy}},\tag{12}$$

$$\operatorname{Re}\delta\left[\delta u D\_{\mathbf{x}}\upsilon + \upsilon D\_{\mathbf{y}}\upsilon\right] = -D\_{\mathbf{y}}\upsilon + D\_{\mathbf{y}}S\_{\mathbf{y}\mathbf{y}} + \delta^{2}D\_{\mathbf{x}}S\_{\mathbf{y}\mathbf{x}},\tag{13}$$

$$P\_r \text{Re} \left[ \delta u D\_x \theta + \delta v D\_y \theta \right] = \left[ -\delta^2 D\_{xx} \theta + D\_{yy} \theta \right] + B\_r S\_{xy} D\_y u \theta,\tag{14}$$

and

$$S\_{xy} = -\left(1 + \mathcal{W}e\dot{\gamma}\right)\left(D\_y u + \delta D\_x v\right),\tag{15}$$

where,

$$\delta = \frac{\overline{d\_1}}{\lambda}, \text{Re} = \frac{\rho c \overline{d\_1}}{\mu\_0}, \text{We} = \frac{\Gamma c}{\overline{d\_i}}, B\_r = P\_r E\_\varepsilon, P\_r = \frac{v S\_{\hbar} \rho}{K}, E\_\varepsilon = \frac{c^2}{S\_{\hbar} \left(T\_i - T\_0\right)}. \tag{16}$$

In the above equation, *We* the Weissenberg number, *Ec* is the Eckert number, *Pr* the Prandlt number, Re the Reynolds number, and *Br* the Brinkman number. Under the assumptions of long wavelength and low Reynolds numbers (*δ* ≈ 1, Re ≈ 0), Equations (12)–(14) take the form:

$$D\_x p = D\_y \left[ \left( 1 + \mathcal{W} e D\_y u \right) D\_y u \right], \tag{17}$$

$$D\_{\Psi}p = 0,\tag{18}$$

$$D\_{yy}\theta = -B\_r\left[\left(D\_yu\right)^2 + \text{We}\left(D\_yu\right)^3\right].\tag{19}$$

This equation implies that *p* = *p*(*y*) so *∂p*/*∂x* can be written as *dp*/*dx*. At *We* = 0, the above equation turns into viscous fluid flow. The associated no slip and convective boundary conditions selected for the problem read as:

$$u = -1, \theta' - B\_i \theta = -B\_i \text{ at } y = h\_i \left(\mathbf{x}\right) = 1 + a \cos 2\mathbf{x} \,\pi,$$

$$u = -1, \theta = 0 \text{ at } y = h\_{\rangle}(\mathbf{x}) = -d - b \cos \left(\phi + 2\pi \mathbf{x}\right). \tag{20}$$

where *Bi* is the Biot number.

**Figure 1.** Flow structure.

#### **3. Entropy Generation Analysis**

According to the theory of thermodynamics, the physical process can be divided in to two types: Irreversible and reversible process. The characterization of such kind of procedures is associated with the change of entropy. Particularly, we say that the process is reversible if there is no change in the entropy, whereas, if the change occurs i.e., entropy is not zero, it shows that the process is irreversible. Therefore, the production of entropy is the measure of the irreversibility of a process. All the processes that arise in nature are irreversible and this reveals a significant obstacle in the study of that process.

The entropy generation in the dimensional form can be defined as:

$$S\_{\mathcal{S}^{\rm ren}}' = \frac{K}{T\_0^2} \left( D\_{\overline{\nabla}} T \right)^2 + \frac{\mathcal{S}\_{\overline{\nabla}\overline{\nabla}}}{T\_0} D\_{\overline{\nabla}} \overline{\mathcal{U}}.\tag{21}$$

Here we define some new dimensionless quantities in addition to those used above:

$$S\_{\mathcal{S}}' = \frac{K\left(T\_i - T\_0\right)}{T\_0^2 \overline{d\_i^2}}, \Delta = \frac{T\_0}{K\left(T\_i - T\_0\right)}.\tag{22}$$

Using Equation (22) in Equation (21), we get the dimensionless form of entropy generation:

$$N = \frac{S'\_{\text{gen}}}{S'\_{\mathcal{S}}} = \left(D\_y \theta\right)^2 + \Delta B\_{\mathcal{T}} \left[-1 + \mathcal{W}eD\_y u\right] \left(D\_y u\right)^2. \tag{23}$$

In the above expression, Δ shows the entropy production characteristics and temperature difference parameter. Equation (23) is divided into two parts. The first is due to the finite temperature difference whereas the second part defines the fluid frictional irreversibility.

The Bejan number is describe as the entropy production ratio because of heat transfer irreversibility to the total entropy production:

$$B\varepsilon = \frac{\left(D\_y\theta\right)^2}{\left(D\_y\theta\right)^2 + \Delta B\_r\left[-1 + \mathcal{W}\varepsilon D\_y u\right]\left(D\_y u\right)^2}.\tag{24}$$

Bejan number lies between 0 to 1. *Be* < 1 represents that the total entropy production dominates the total entropy production due to heat transfer. *Be* = 1 represents when the total entropy production is equal to entropy production due to heat transfer irreversibility.

#### **4. Series Solution**

Since Equation (17) is non linear, its exact solution may not be possible, therefore, we employ the regular perturbation method to find the solution. For perturbation solution, we expand *u*, *F* and *dp*/*dx* as:

$$\mu = \sum\_{n=0}^{\infty} \mathcal{W} e^n u\_{n\prime} \tag{25}$$

$$F = \sum\_{n=0}^{\infty} \mathcal{W} \mathcal{e}^n F\_{n\prime} \tag{26}$$

$$\frac{dp}{d\mathbf{x}} = \sum\_{n=0}^{\infty} \mathcal{W} \epsilon^n \frac{dp\_n}{d\mathbf{x}} \, , \tag{27}$$

Substituting above expression in Equation (17) and their boundary conditions in Equation (20) and comparing the coefficients of powers of *We* we get the zeroth and first order systems which can be manipulated easily by a mathematical computing tool Mathematica and are conclusively stated as:

$$\begin{split} u &= \frac{1}{2!} \left[ -2 + \mathbb{C}\_{1} h\_{1} h\_{2} - \mathbb{C}\_{1} \mathbb{C}\_{3} y + \mathbb{C}\_{1} y^{2} \right] + \frac{1}{3!} W \varepsilon \left[ \mathbb{C}\_{2} h\_{1} h\_{2} + \mathbb{C}\_{1}^{2} \mathbb{C}\_{3} h\_{1} h\_{2} \\ &- \mathbb{C}\_{2} \left( \mathbb{C}\_{3} + y \right) y + \mathbb{C}\_{1}^{2} \left( h\_{1}^{2} - 4 h\_{1} h\_{2} - h\_{2}^{2} \right) y + \mathbb{C}\_{1}^{2} (3 \mathbb{C}\_{3} - 2 y) y \right] + O(W \varepsilon^{2}), \end{split} \tag{28}$$

$$\frac{dp}{dx} = 12\text{C}\_1 + 12\text{WeC}\_2 + O(\mathcal{W}e^2),\tag{29}$$

where the constant are defined as:

$$\mathcal{C}\_1 = \frac{12\left(1 + d - h\_1 + h\_2 - Q\right)}{\left(h\_1 - h\_2\right)^3},\tag{30}$$

$$\mathcal{C}\_2 = \frac{3\hbar \left(1 + d - Q\right)}{\left(h\_1 - h\_2\right)^3},\tag{31}$$

$$C\_3 = h\_1 + h\_{2\prime} \tag{32}$$

$$\mathbf{C\_4} = \mathbf{C\_3^2} + h\_1 h\_2 \tag{33}$$

$$\mathcal{C}\_{\mathbb{S}} = 17\mathcal{C}\_{\mathbb{S}}^2 + 4h\_1h\_2 \tag{34}$$

$$\mathbf{C}\_6 = \mathbf{7}\mathbf{C}\_3^2 + 2h\_1h\_2. \tag{35}$$

The solution for velocity *u* obtained by above perturbation method can be used in Equation (19). The final solution for *θ* can be obtained by integrating Equation (19) along with their associated boundary conditions (See Equation (20)) and can be written as:

$$
\theta = \theta\_1 + \theta\_2 y + \theta\_3 y^2 + \theta\_4 y^3 + \theta\_5 y^4 + \theta\_6 y^5 + \theta\_7 y^6 + \theta\_8 y^7 + \theta\_9 y^8,\tag{36}
$$

where constants of integration *θ*<sup>1</sup> and *θ*<sup>2</sup> can be evaluated by using boundary conditions defined in Equation (20) and the expression obtained are very large and therefore are not presented here. The remaining constants are defined as:

$$\begin{split} \theta\_3 &= \frac{B\_r \sqrt{\mathsf{C}\_3}}{432} \left[ 3 \frac{dp\_0}{dx} + \left( \mathbb{C}\_3 \left( \frac{dp\_0}{dx} \right)^2 + 3 \frac{dp\_1}{dx} \right) \mathsf{W} \epsilon \right]^2 \\ &\times \left[ -6 + 3 \mathsf{C}\_3 \mathsf{W} \epsilon \frac{dp\_0}{dx} + \mathsf{C}\_3 \left( \mathsf{C}\_3 \left( \frac{dp\_0}{dx} \right)^2 + 3 \frac{dp\_1}{dx} \right) \mathsf{W} \epsilon^2 \right], \end{split} \tag{37}$$

$$\begin{split} \theta\_{4} &= -\frac{B\_{7}\mathbb{C}\_{3}}{72} \left( \frac{dp\_{0}}{dx} + \mathbb{C}\_{3} \mathcal{W} \mathcal{e} \left( \frac{dp\_{0}}{dx} \right)^{2} + \frac{dp\_{1}}{dx} \right) \\ &\times \left[ 3\frac{dp\_{0}}{dx} + \left( \mathbb{C}\_{3} \left( \frac{dp\_{0}}{dx} \right)^{2} + 3\frac{dp\_{1}}{dx} \right) \mathcal{W} \mathcal{e} \right] \\ &\times \left[ -4 + 3\mathbb{C}\_{3} \mathcal{W} \mathcal{e} \frac{dp\_{0}}{dx} + \mathbb{C}\_{3} \left( \mathbb{C}\_{3} \left( \frac{dp\_{0}}{dx} \right)^{2} + 3\frac{dp\_{1}}{dx} \right) \mathcal{W} \mathcal{e}^{2} \right], \end{split} \tag{38}$$

$$\begin{split} \theta\_{5} &= \frac{B\_{r}}{12} \left[ -12 \left( \frac{dp\_{0}}{dx} \right)^{2} - 6 \frac{dp\_{0}}{dx} \left( 3 \mathbf{C}\_{3} \left( \frac{dp\_{0}}{dx} \right)^{2} + 4 \frac{dp\_{1}}{dx} \right) \right. \\ &\left. + \left( (7h\_{1} + 5h\_{2})(5h\_{1} + 7h\_{2}) \left( \frac{dp\_{0}}{dx} \right)^{4} + 18 \mathbf{C}\_{3} \left( \frac{dp\_{0}}{dx} \right)^{2} \frac{dp\_{1}}{dx} - 12 \left( \frac{dp\_{1}}{dx} \right)^{2} \right) \mathcal{W} \mathbf{c}^{2} \\ &\left. + 6 \frac{dp\_{0}}{dx} \left( 6 \mathbf{C}\_{3} \mathbf{C}\_{4} \left( \frac{dp\_{0}}{dx} \right)^{4} + \mathcal{C}\_{5} \left( \frac{dp\_{0}}{dx} \right)^{2} \frac{dp\_{1}}{dx} + 9 \mathbf{C}\_{3} \left( \frac{dp\_{1}}{dx} \right)^{2} \right) \right. \\ &\left. + \mathcal{C}\_{3} \left( \mathcal{C}\_{3} \left( \frac{dp\_{0}}{dx} \right)^{2} + 3 \frac{dp\_{1}}{dx} \right) \left( \mathcal{C}\_{6} \left( \frac{dp\_{0}}{dx} \right)^{4} + 15 \mathcal{C}\_{3} \left( \frac{dp\_{0}}{dx} \right)^{2} \frac{dp\_{1}}{dx} + 6 \left( \frac{dp\_{1}}{dx} \right)^{2} \right) \mathcal{W} \mathbf{c}^{4} \right), \end{split}$$

$$\begin{split} \theta\_{6} &= -\frac{B\_{r}We}{20} \left( \frac{dp\_{0}}{dx} + \mathbb{C}\_{3} \left( \frac{dp\_{0}}{dx} \right)^{2} + \mathcal{W}\epsilon \frac{dp\_{1}}{dx} \right) \\ &\times \left[ -\left(\frac{dp\_{0}}{dx}\right)^{2} + \frac{dp\_{0}}{dx}\mathcal{W}\epsilon \left( 5\mathbb{C}\_{3} \left(\frac{dp\_{0}}{dx}\right)^{2} + 2\frac{dp\_{1}}{dx} \right) \mathcal{W}\epsilon \right] \\ &+ \left( 2\mathbb{C}\_{4} \left(\frac{dp\_{0}}{dx}\right)^{4} + 5\mathbb{C}\_{3} \left(\frac{dp\_{0}}{dx}\right)^{2} \frac{dp\_{1}}{dx} + \left(\frac{dp\_{1}}{dx}\right)^{2} \right) \mathcal{W}\epsilon^{2} \right], \tag{40} \end{split}$$

$$\theta = \frac{B\_r W \varepsilon^2}{60} \left(\frac{dp\_0}{dx}\right)^2 \left[ 4 \left(\frac{dp\_0}{dx}\right)^4 + 3 \frac{dp\_0}{dx} \left( 5 \text{C}\_3 \left(\frac{dp\_0}{dx}\right)^2 + 4 \frac{dp\_1}{dx} \right) \mathcal{W} \varepsilon \right]$$

$$+ \left( \mathcal{C}\_6 \left(\frac{dp\_0}{dx}\right)^4 + 15 \mathcal{C}\_3 \left(\frac{dp\_0}{dx}\right)^2 \frac{dp\_1}{dx} + 6 \left(\frac{dp\_1}{dx}\right)^2 \right) \mathcal{W} \varepsilon^2 \right],\tag{41}$$

$$\theta\_8 = -\frac{B\_r \mathcal{W} e^3}{14} \left(\frac{dp\_0}{dx}\right)^4 \left[\frac{dp\_0}{dx} + \mathcal{C}\_3 \mathcal{W} e \left(\frac{dp\_0}{dx}\right)^2 + \mathcal{W} e \frac{dp\_1}{dx}\right],\tag{42}$$

$$
\theta\_{\Theta} = \frac{1}{56} B\_r \mathcal{W} e^4 \left(\frac{dp\_0}{dx}\right)^6. \tag{43}
$$

The dimensionless mean flow reads as:

$$F = d + 1 - Q.\tag{44}$$

and

$$F = \int\_{h\_2}^{h\_1} \mathbf{u} \,\mathrm{d}y.\tag{45}$$

The expression for entropy generation and Bejan number can be easily obtained by incorporating value of *u* and *θ* in Equation (24).

#### **5. Discussion**

In this section, we present our results by varying the quantities under the variation of several factors. Figures of temperature profile *θ*, entropy generation coefficient *N*, and streamlines are illustrated below. Figures 2–5 reflect the behavior of *θ* for some useful parameters. Entropy generation graphs are given in Figures 6–11. The streamlines conducting the flow samples are depicted in Figures 12 and 13.

Figure 2 shows the impact of parameters *a* and *b* on temperature profile *θ*. It can be observed from this plot that temperature is getting increased for both parameters from the lower wall to the upper wall. Figure 3 shows the mechanism of the Biot number and Brinkman number. Biot number is an important mechanism to determine the heat transfer. It can be visualized from this figure that an enhancement in Biot number tends to boost the temperature profile while the contrary behavior has been observed with the Brinkman number. Brinkman number is the product of Eckert and Prandtl numbers *Br* = *PrEc*, or it is the ratio of the heat generated by viscous dissipation and propagation of heat by molecular conduction, such as, the ratio of the viscous heat production to extrinsic heating. Therefore, the enhancement of Brinkman's number tends to increase the temperature profile. It can be seen in Figure 4 that the volumetric flow rate significantly enhances the temperature profile. It can also be noticed that the temperature profile has a lower magnitude for smaller values of *d* whereas the behavior is converse for higher values. It can be viewed from Figure 5 that the Weissenberg number causes a remarkable resistance for higher values. By enhancing the Weissenberg number, the elastic forces are more dominant, which diminishes the temperature profile. However, the phase difference *φ* also produces a significant resistance in the temperature profile.

Figures 6–9 are presented for entropy profiles against the leading parameters. It can be viewed from Figure 6 that an increment in *a* and *b* tends to boost the entropy profile whereas the entropy profile is increasing along the whole channel. Figure 7 shows that by increasing the Brinkman number, the entropy profile rises, and it decreases by increasing the Weissenberg number. However, the entropy remains positive and growing along the entire channel. It is seen from Figure 8 that the Biot number enhances the entropy profile. It can be seen that at the lower wall, the entropy profile is maximum and minimum at the upper wall, whereas it is uniform in the middle of the channel. The entropy profile for various values of Δ is presented in Figure 9. It is noticed in this figure that the entropy profile is uniform, and no change occurs in the middle of the channel i.e., *y* ∈ (0, 0.5). Although it shows a decreasing pattern, but it rises along the upper wall of the channel and remains positive.

Figures 10 and 11 are plotted for the Bejan number profile against the governing parameters. It is observed from Figure 10 that the Bejan number profile diminishes for higher values of the Brinkman number and shows a converse behavior for the Weissenberg number. In Figure 11, we can see that the phase difference shows versatile behavior for higher values on the Bejan number profile. When Bejan number rises, then the phase difference's effects are negligible for the domain *y* ∈ (0, 1.3), while when the Bejan number is small, it decreases in a similar area.

The most interesting and useful phenomena of peristaltic motion are trapping, which is plotted in Figures 12 and 13 via streamlines. It was found that by enhancing the phase difference parameter, the effects are negligible on the trapping bolus despite the fact that an unusual movement in the magnitude of the bolus is noticed. Furthermore, we can see in Figure 13 that an increment in the Weissenberg number profile tends to diminish the width of the trapping bolus. The number of boluses disappeared more quickly in the lower region as compared with the upper one.

**Figure 2.** Temperature distribution for different values of *a* and *b*. Solid line: *a* = 0.1, dashed line: *a* = 0.15 and dot-dashed line: *a* = 0.2.

**Figure 3.** Temperature distribution for different values of *Bi* and *Br*. Solid line: *Bi* = 0.1, dashed line: *Bi* = 0.25 and dot-dashed line: *Bi* = 0.3.

**Figure 4.** Temperature distribution for different values of *Q* and *d*. Solid line: *Q* = 1.0, dashed line: *Q* = 1.2 and dot-dashed line: *Q* = 1.4.

**Figure 5.** Temperature distribution for different values of *φ* and *We*. Solid line: *φ* = 0.1, dashed line: *φ* = 0.5 and dot-dashed line: *φ* = 0.9.

**Figure 6.** Entropy profile for different values of *a* and *b*. Solid line: *a* = 0.1, dashed line: *a* = 0.15, and dot-dashed line: *a* = 0.2.

**Figure 7.** Entropy profile for different values of *We* and *Br*. Solid line: *Br* = 1.0, dashed line: *Br* = 1.2, and dot-dashed line: *Br* = 1.4.

**Figure 8.** Entropy profile for different values of *Bi* and *Br*. Solid line: *Bi* = 0.1, dashed line: *Bi* = 0.25, and dot-dashed line: *Bi* = 0.3.

**Figure 9.** Entropy profile for different values of Δ. Solid line: Δ = 0.1, dashed line: Δ = 0.2, dot-dashed line: Δ = 0.3, and dot-dashed line: Δ = 0.4.

**Figure 10.** Bejan number for different values of *We* and *Br*. Solid line: *Br* = 1.0, dashed line: *Br* = 1.2, and dot-dashed line: *Br* = 1.4.

**Figure 11.** Bejan number for different values of *φ* and *Bi*. Solid line: *φ* = 0.1, dashed line: *φ* = 0.5, and dot-dashed line: *φ* = 0.9.

**Figure 12.** Trapping mechanism for different values of *φ*.

**Figure 13.** Trapping mechanism for different values of *We*.

#### **6. Conclusions**

In this study, we analyzed the entropy generation on the asymmetric peristaltic propulsion of non-Newtonian fluid with convective boundary conditions. The Williamson fluid model was considered to examine the entropy profile. The mathematical modeling was performed under the approximation of small Reynolds number and long wavelength of the peristaltic wave. A regular perturbation method was employed to get the series solutions up to third-order approximation. The significant results of the governing flow problem are summarized below:


The present results provide an excellent benchmark for further study on the entropy production with mass transfer and peristaltic pumping mechanism. The mass transfer phenomena with magnetic and porosity effects that were not covered in this paper is a topic for future research.

**Author Contributions:** Investigation, A.R.; M.M.B., Methodology; Conceptualization, R.E.; Validation, A.Z.; Writing—review & editing, S.M.S. All authors have read and agreed to the published version of the manuscript.

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

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

#### **Nomenclature**



#### **Greek Symbol**


#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Serious Solutions for Unsteady Axisymmetric Flow over a Rotating Stretchable Disk with Deceleration**

#### **Muhammad Adil Sadiq**

Department of Mathematics, DCC-KFUPM, KFUPM Box 5084, Dhahran 31261, Saudi Arabia; adilsadiq@kfupm.edu.sa; Tel.: +966-13-868-3300

Received: 12 December 2019; Accepted: 2 January 2020; Published: 3 January 2020

**Abstract:** In this article, the author has examined the unsteady flow over a rotating stretchable disk with deceleration. The highly nonlinear partial differential equations of viscous fluid are simplified by existing similarity transformation. Reduced nonlinear ordinary differential equations are solved by homotopy analysis method (HAM). The convergence of HAM solutions is also obtained. A comparison table between analytical solutions and numerical solutions is also presented. Finally, the results for useful parameters, i.e., disk stretching parameters and unsteadiness parameters, are found.

**Keywords:** homotopy analysis method; rotating stretchable disk; newtonian fluid; axisymmetric flow

#### **1. Introduction**

Recently, due to the massive practical application in the scientific and technical field, the study of the rotating stretchable disk has become significant, such as thermal power generation system, medical equipment's, computer storage devices, rotating machinery, gas turbine routers, air cleaning machines, crystal growth process, and in aerodynamic applications [1]. Initially, von Kármán [2] conducted a study on rotating disk. Several researchers then illustrated the different aspects of this important analysis. Fang and Zhang [3] have highlighted the flow between two stretchable disks and found the exact solutions. The parameters analysis and optimization of entropy generation in unsteady magneto hydrodynamics flow over a rotating stretchable flow over a rotating disk using artificial neural network and practical swarm optimization algorithm was presented by Rashidi et al. [4]. Recently, Fang and Tao [5] wrote about the unsteady flow over a rotating stretchable disk with deceleration. After using the similarity analysis, they found the numerical solutions.

In many situations, exact solutions are very difficult and in most of the cases exact solutions are impossible. Therefore, series solutions are more useful if they satisfy the given initial and boundary value problems. Nevertheless, there are various analytical approaches, and each approach has certain limitations. However, homotopy analysis method (HAM) has many advantages over many analytical methods. Liao [6] introduced the idea of HAM, which is used by many researchers effectively. Some useful studies are cited in [7–15]. The purpose of this article is to illustrate the application of HAM for unsteady Newtonian fluid flow over a rotating stretchable disk with declaration. Tables provide a correlation between current HAM solution and Fang and Tao's [5] numerical solution.

#### **2. Formulation of the Problem**

Let us consider an incompressible, laminar, and unsteady flow of a viscous fluid or Newtonian fluid over a stretchable disk, which is rotating about the *z*-axis with time dependent angular velocity <sup>Ω</sup>(*t*) <sup>=</sup> <sup>Ω</sup> <sup>1</sup>−*bt*, where <sup>Ω</sup> is constant angular speed of the disk and 'b' is the measure of unsteadiness. Flow is due to the rotation of the stretchable disk and is axisymmetric about the *z*-axis. Figure 1 describe

the geometry of the proposed problem. The governing equations for an unsteady three-dimensional flow of viscous fluid in cylindrical coordinates are shown below.

$$\frac{1}{r}\frac{\partial}{\partial r}(ru) + \frac{\partial w}{\partial z} = 0,\tag{1}$$

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial r} + w \frac{\partial u}{\partial z} - \frac{v^2}{r} = -\frac{1}{\rho} \frac{\partial p}{\partial r} + v \left( \frac{1}{r} \frac{\partial}{\partial r} (r \tau\_{\mathcal{H}}) + \frac{\partial}{\partial z} (\tau\_{z\mathcal{T}}) - \frac{\tau\_{\partial \mathcal{O}}}{r} \right) \tag{2}$$

$$\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial r} + w\frac{\partial v}{\partial z} - \frac{uv}{r} = v\left(\frac{\partial}{\partial r}(\tau\_{r0}) + \frac{\partial}{\partial z}(\tau\_{z0}) + \frac{\tau\_{\partial r} + \tau\_{r0}}{r}\right) \tag{3}$$

$$\frac{\partial w}{\partial t} + u \frac{\partial w}{\partial r} + w \frac{\partial w}{\partial z} = -\frac{1}{\rho} \frac{\partial p}{\partial z} + v \left( \frac{1}{r} \frac{\partial}{\partial r} (r \tau\_{rz}) + \frac{\partial}{\partial z} (\tau\_{zz}) \right) \tag{4}$$

where *r* is along the radial direction, θ is along the azimuthal direction, and *z* is in normal direction to the axis. Here, Equation (1) is the continuity equation and Equations (2) and (4) represent the momentum equation for incompressible flow.

**Figure 1.** Geometry of the problem.

Where *u*, *v*, and *w* are the velocities along *r*, θ, and *z* directions, ρ is the density of fluid, *p* is the pressure, *v* is the kinematic viscosity, and τ*rr*, τ*r*θ, τ*zr*, τ*z*θ, τ*zz* are the stress which are defined as

$$
\begin{bmatrix}
\tau\_{rr} & \tau\_{r\theta} & \tau\_{rz} \\
\tau\_{\ell r} & \tau\_{\ell\theta} & \tau\_{\ell z} \\
\tau\_{zr} & \tau\_{z\theta} & \tau\_{zz}
\end{bmatrix} = \begin{bmatrix}
2\frac{\partial\mu}{\partial r} & \frac{\partial\mu}{\partial r} - \frac{\overline{v}}{r} & \frac{\partial\mu}{\partial z} + \frac{\partial\overline{w}}{\partial r} \\
\frac{\partial\overline{v}}{\partial r} - \frac{\overline{v}}{r} & \frac{2\mu}{r} & \frac{\partial\overline{w}}{\partial z} \\
\frac{\partial\mu}{\partial\mu} + \frac{\partial\overline{w}}{\partial r} & \frac{\partial\overline{v}}{\partial z} & 2\frac{\partial\overline{w}}{\partial z}
\end{bmatrix}.
\tag{5}
$$

The proposed boundary conditions are specified in accordance with the geometry of the problem as

$$\begin{aligned} u(r,\theta,0) &= \frac{a\Omega r}{1 - bt}, \ v(r,\theta,0) = \frac{\Omega r}{1 - bt}, \ w(r,\theta,0) = 0, \\ u(r,\theta,\infty) &= v(r,\theta,\infty) = 0, \end{aligned} \tag{6}$$

where 'a' is the disk stretching parameter.

Introducing, the similarity transformation used in [5] are

$$\begin{split} u &= \frac{\Omega r}{1 - bt} f'(\eta), \; v = \frac{\Omega r}{1 - bt} g(\eta), \; w = -2 \frac{\sqrt{\Omega v}}{\sqrt{1 - bt}} f(\eta), \\ p &= \frac{\rho v \Omega}{1 - bt} P(\eta), \; \text{and } \eta = \sqrt{\frac{\Omega}{v}} \frac{z}{\sqrt{1 - bt}}. \end{split} \tag{7}$$

*Symmetry* **2020**, *12*, 96

Applying these similarities into the above equations, following non-dimensional equations along with boundary conditions can be obtained as

$$f'''' + 2ff'' - f'^2 + g^2 = \mathcal{S} \{ \frac{\eta}{2} f'' + f' \}. \tag{8}$$

$$\mathbf{g}\,\mathbf{g}\,\mathbf{\prime}-2f\,\mathbf{g}-2fg\,\mathbf{\prime}=\mathbf{S}\Big(\frac{\eta}{2}\mathbf{g}\,\mathbf{\prime}+\mathbf{g}\Big).\tag{9}$$

$$P' = 2f'' + 4ff'' - \mathcal{S}(\eta f' + f),\tag{10}$$

$$f(0) = 0, \ f'(0) = a, \ g(0) = 1, \ f' \to 0 \text{ and } \ g \to 0, \text{ as } \eta \to \infty,\tag{11}$$

where *S* = *b*/Ω is the unsteadiness parameter.

#### **3. Homotopy Analysis Method**

Homotopy Analysis Method (HAM) [6–15] is used to find an analytical solution to Equations (8) and (11). The velocity distribution *f*(η) and *g*(η) can be expressed by a set of base functions

$$\left\{\eta^{n}\exp^{(-m\eta)}|m,n\rangle \ge 0\right\}\tag{12}$$

in the form

$$f(\eta) = a\_{0,0}^0 + \sum\_{n=0}^{\infty} \sum\_{k=0}^{\infty} a\_{m,n}^k \eta^n \exp^{(-m\eta)},\tag{13}$$

$$\log(\eta) = b\_{0,0}^0 + \sum\_{n=0}^{\infty} \sum\_{k=0}^{\infty} b\_{m,n}^k \eta^n \exp^{(-m\eta)}\tag{14}$$

in which *ak <sup>m</sup>*,*<sup>n</sup>* and *bk <sup>m</sup>*,*<sup>n</sup>* are the coefficients, the initial guesses *f*<sup>0</sup> and *g*<sup>0</sup> can be selected on the basis of the law of the solution expressions and of the boundary conditions:

$$f\_0(\eta) = a(1 - \exp^{-\eta}),\tag{15}$$

$$\mathfrak{g}\_0(\eta) = \exp^{-\eta}.\tag{16}$$

The auxiliary linear operators are

$$\mathcal{L}\_1 = \frac{d^3}{d\eta^3} + \frac{d^2}{d\eta^2},\tag{17}$$

$$\mathcal{L}\_2 = \frac{d^2}{d\eta^2} + \frac{d}{d\eta'} \tag{18}$$

which satisfy

$$\mathcal{L}\_1[\mathcal{C}\_1 + \mathcal{C}\_2 \exp^{-\eta} + \mathcal{C}\_3 \eta] = 0,\tag{19}$$

$$\mathcal{L}\!\_2[\mathbb{C}\_4 \exp^{-\eta} + \mathbb{C}\_5] = 0,\tag{20}$$

where *Ci*(*i* = 1 − 5) are integral constants.

#### *3.1. Zeroth-Order Deformation Equation*

If *<sup>q</sup>* <sup>∈</sup> [0, 1] denote an embedding parameter, *<sup>f</sup>* and *g* indicate the non zero auxiliary parameters for *f*(η) and *g*(η), the zeroth-order deformations for the given problem are

$$(1 - q)\_f \left[ \hat{f}(\eta; q) - f\_0(\eta) \right] = q \hbar\_f N\_f \left[ \hat{f}(\eta; q), \xi(\eta; q) \right],\tag{21}$$

$$(1 - q)\_{\mathcal{S}} [\pounds(\eta; q) - \pounds0(\eta)] = q \hbar\_{\mathcal{S}} N\_{\mathcal{S}} [\pounds(\eta; q), \hat{f}(\eta; q)].\tag{22}$$

*Symmetry* **2020**, *12*, 96

$$\hat{f}(0;q) = 0,\ \hat{f}^\flat(0;q) = a,\ \hat{\mathfrak{z}}(0;q) = 1\ \hat{f}^\flat\left(\infty;q\right) = \mathfrak{z}(\infty;q) = 0. \tag{23}$$

Defining the nonlinear operators for the above problem as

$$\begin{split} N\_f \left[ \hat{f}(\eta; q)\_\* \hat{g}(\eta; q) \right] &= \frac{\partial^3 \hat{f}(\eta; q)}{\partial \eta^3} + 2 \hat{f}(\eta; q) \frac{\partial^2 \hat{f}(\eta; q)}{\partial \eta^2} - \left( \frac{\partial \hat{f}(\eta; q)}{\partial \eta} \right)^2 \\ &+ \left( \hat{g}(\eta; q) \right)^2 - S \left( \frac{\eta}{2} \frac{\partial^2 \hat{f}(\eta; q)}{\partial \eta^2} + \frac{\partial \hat{f}(\eta; q)}{\partial \eta} \right) \end{split} \tag{24}$$

$$\begin{split} \mathcal{N}\_{\mathcal{S}} \Big[ \mathcal{G}(\eta;q), f(\eta;q) \Big] &= \frac{\partial^{2} \mathcal{G}(\eta;q)}{\partial \eta^{2}} + 2f(\eta;q) \frac{\partial \mathcal{G}(\eta;q)}{\partial \eta} \\ &- 2\mathcal{G}(\eta;q) \frac{\partial \hat{f}(\eta;q)}{\partial \eta} - s \Big( \frac{\eta}{2} \frac{\partial \hat{\xi}^{\circ}(\eta;q)}{\partial \eta} + \mathcal{G}(\eta;q) \Big). \end{split} \tag{25}$$

For *q* = 0 and *q* = 1, one can have

$$f(\eta;0) = f\_0(\eta), \; f(\eta:1) = f(\eta), \tag{26}$$

$$\mathfrak{g}(\eta;0) = \mathfrak{g}\_0(\eta), \mathfrak{g}(\eta:1) = \mathfrak{g}(\eta). \tag{27}$$

By Taylor's theorem

$$\hat{f}(\eta;q) = f\_0(\eta) + \sum\_{m=1}^{\infty} f\_m(\eta) q^m, \; f\_m(\eta) = \left. \frac{1}{m!} \frac{\partial^m f(\eta;q)}{\partial \eta^m} \right|\_{\eta=0\prime} \tag{28}$$

$$\hat{g}(\eta;q) = g\_0(\eta) + \sum\_{m=1}^{\infty} g\_m(\eta) q^m,\\ \ g\_m(\eta) = \left. \frac{1}{m!} \frac{\partial^m g(\eta;q)}{\partial \eta^m} \right|\_{\eta=0,} \tag{29}$$

and

$$f(\eta) = \left. f \right| \eta \left( \eta \right) + \sum\_{m=1}^{\infty} f\_m(\eta)\_{\prime} \tag{30}$$

$$\mathcal{g}(\eta) = \mathcal{g}\_0(\eta) + \sum\_{m=1}^{\infty} \mathcal{g}\_m(\eta). \tag{31}$$

#### *3.2. Mth-Order Deformation*

Differentiating the zeorth-order deformation Equations (21) and (23) with respect to *q*, then setting *q* = 0, and finally dividing them by *m*!, the mth-order deformation equations can be obtained as

$$
\mathcal{L}\_1[f\_m(\eta) - \chi\_m f\_{m-1}(\eta)] = \, ^\circ \hbar \prescript{f}{}{R}(\eta), \tag{32}
$$

$$
\mathcal{L}\_2[\mathcal{g}\_m(\eta) - \chi\_m \mathcal{g}\_{m-1}(\eta)] = \, ^\circ \hbar\_\S \prescript{\mathcal{G}}{}{\eta}\_m(\eta), \tag{33}
$$

$$f\_m(0) = f\_m'(0) = \ g\_m(0) = f\_m'(\infty) = \mathcal{g}\_m(\infty) = 0,\tag{34}$$

where

$$R\_m^f(\eta) = f\_{m-1}^{\prime\prime} + 2\sum\_{k=0}^{m-1} f\_k f\_{m-1-k}^{\prime\prime} - \sum\_{k=0}^{m-1} f\_k^{\prime} f\_{m-1-k}^{\prime} + \sum\_{k=0}^{m-1} g\_k g\_{m-1-k} - S\left(\frac{\eta}{2} f\_{m-1}^{\prime\prime} + f\_{m-1}^{\prime}\right) \tag{35}$$

$$R\_m^{\mathcal{S}}(\eta) = \mathcal{g}\_{m-1}^{\prime\prime} - 2\sum\_{k=0}^{m-1} g\_k f\_{m-1-k}^{\prime} + 2\sum\_{k=0}^{m-1} f\_k^{\prime} g\_{m-1-k}^{\prime} - S\left(\frac{\eta}{2}\mathcal{g}\_{m-1}^{\prime} + g\_{m-1}\right). \tag{36}$$

$$\chi\_m = \begin{cases} \begin{array}{cc} 0, & m \le 1, \\ 1, & m > 1, \end{array} \end{cases} \tag{37}$$

in which *fm*(η) and *gm*(η) denote the special solutions of Equations (32) and (33) and the *Ci*(*i* = 1 − 5) integral constants are calculated by the (34) boundary conditions. Equations (32) and (34) can be solved using Mathematica for *m* = 1, 2, 3 ... .

#### **4. Convergence of the HAM Solution**

The homotopy analysis method includes the regulating parameter *h*, which controls the region of convergence and HAM solution approximation. To ensure that the solutions converge within the admissible spectrum of auxiliary parameter values and *h f* and *hg*, *h* − *curves* were sketched for 15th-order approximation. The *h* − *curves* are plotted in Figures 2 and 3. The admissible ranges of values of *h f* and *hg* are −1.5 ≤ *hf* < −0.3 and − 1.5 ≤ *hg* < −0.3, these ranges vary with the change in parameters.

**Figure 3.** 15th-order *f* (0) for *a* = 1 and *S* = −1.

#### **5. Results and Discussion**

To solve Equations (8) and (9) homotopy analysis method (HAM) is applied as a subject to the boundary conditions (11). Homotopy analysis method is a strong analytical technique which is applied to obtain the convergent series solution of nonlinear differential equations. The convergence region for HAM through *h* − *curves* are sketched and analyzed in Figures 3 and 4. Homotopy analysis method provides great freedom to obtain the convergent result. The convergence region varies for different values of *a* and *S*. Tables 1–4 represent the convergence of solution for different values of parameters. The error analysis of the obtained approximated results is as follows.

$$E\_m = \int\_0^\infty e\_m^2(t)dt,\tag{38}$$

where *em*(*t*) is the residual error of Equations (8) and (9) at the mth-order approximation. It is observed that 10th-order approximation is in good agreement with the numerical result.

**Figure 4.** Comparison of convergence of *g*(η) when - = −0.333, *a* = 1 and *S* = −1 (line: 10th-order, dots: 5th-order).

**Table 1.** Comparison of the numerical result [5] with homotopy analysis method (HAM) convergent result when *<sup>a</sup>* <sup>=</sup> 1, *<sup>S</sup>* <sup>=</sup> <sup>−</sup>1 and *<sup>f</sup>* <sup>=</sup> <sup>−</sup>1/3, *<sup>g</sup>* = −1/4.


Numerical result [5] *f* (**0**) = −0.6520, *g* (0) = −1.2716.

**Table 2.** Comparison of the numerical result [5] with HAM convergent result when *a* = 1, S = −1/10 and *<sup>f</sup>* <sup>=</sup> <sup>−</sup>1/3, *<sup>g</sup>* = −1/4.


Numerical result [5] *f* (**0**) = −0.9189, *g* (0) = −1.4656.

**Table 3.** Comparison of the numerical result [5] with HAM convergent result when *a* = 2, S = −1/10 and *<sup>g</sup>* = *<sup>f</sup>* = −1/5.


Numerical result [5] *f* (**0**) = −3.1178, *g* (0) = −2.0530.

**Table 4.** Comparison of the numerical result [5] with HAM convergent result when *a* = 1, *S* = −1/2 and *<sup>g</sup>* = *<sup>f</sup>* = −1/4.


Numerical result [5] *f* (**0**) = −0.8007, *g* (0) = −1.3797.

The convergence control parameter plays an important role. In Tables 5 and 6, the effect of on convergence is shown. Tables 5 and 6 show that the convergence of the solution depends strongly on -. It can be seen easily that for one set of the convergence is faster than the other.

**Table 5.** The Convergence analysis of *f* (0) for different when *a* = 1, and *S* = −1/2.


Numerical result [5] *f* (**0**) = −0.8007.

**Table 6.** The Convergence analysis of *g* (0) for different when *a* = 1, and *S* = −1/2.


Numerical result [5] *g* (**0**) = 1.3797.

In Figures 5 and 6, the comparison of 5th-order approximation with 10th-order approximation is shown, which again provide the facts for convergence. The Mathematica software is used to compute the results for higher-order approximation. As the given problem is highly nonlinear, the computation time increases if higher-order approximation is computed or increases the value of the parameters. For *a* = 1 and *S* = 0, the given problem becomes a special case as mentioned in the numerical paper [5]. Table 7 provides the convergence result for this special case as well. The results obtained in the present research for this special case are also in very good agreement with the numerical result. This shows the strength of homotopy analysis methods. It is found that for small *S*, *f* (0) decreases with the increase of 'a' as shown in Tables 2 and 3. Figure 7 represents the velocity distribution for different values of *a*. It is observed that with the increase in disk stretching parameter the velocity decreases. Figures 8 and 9 show that with the decrease in the unsteadiness parameter, both tangential and radial velocities increase.

**Figure 5.** Comparison of convergence of *f*(η) when - = −0.333, *a* = 1 and *S* = −1 (line: 10th-order, dots: 5th-order).

**Figure 6.** For *S* = −1/2 solid line: *a* = 1, Dashed line: *a* = 2 10th order HAM approximation for *f* (η).

**Table 7.** Comparison of the numerical result [5] with HAM convergent result for special case when analysis = 1, *S* = 0, and *<sup>f</sup>* <sup>=</sup> <sup>−</sup>28/100, *g* = −1/3.


Numerical result [5] *f* (0) = −1.1737 , *g* (0) = −1.4541.

**Figure 7.** For *S* = −1/10 10th-order HAM approximation for *f* (η).

**Figure 8.** Variation of *f* (η) for different values of unsteadiness parameter for *a* = 1.

**Figure 9.** Variation of *g*(η) for different values of unsteadiness parameter for *a* = 1.

#### **6. Conclusions**

In this research, viscous axisymmetric flow is studied on a stretchable rotating disk with deceleration. It is found that Navier–Stokes equation admits a similarity solution, which depends on non-dimensionalized parameters *S* and *a* measuring unsteadiness and disk stretching, respectively. The resulting group of nonlinear ordinary differential equations is then solved analytically using homotopy analysis method (HAM). In numerical paper [5] it is mentioned that there are two solution branches. The upper solution branch is physically feasible, but the lower solution branch may not be practically possible. Here, the author has discussed and evaluated the outcome for a physical solution from the upper field.

The main results are summarized as


**Funding:** This research received funding from KFUPM.

**Acknowledgments:** The author wishes to express his thanks for the support received from King Fahd University of Petroleum and Minerals.

**Conflicts of Interest:** Author declare no conflicts of interest with any one.

#### **References**


*Symmetry* **2020**, *12*, 96

© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

*Article*
