Homotopy Simulation of Nonlinear Unsteady Rotating Nanofluid Flow from a Spinning Body

The development of new applications of nanofluids in chemical engineering and other technologies has stimulated significant interest in computational simulations. Motivated by coating applications of nanomaterials, we investigate the transient nanofluid flow from a time-dependent spinning sphere using laminar boundary layer theory. The free stream velocity varies continuously with time. The unsteady conservations equations are normalized with appropriate similarity transformations and rendered into a ninth-order system of nonlinear coupled, multidegree ordinary differential equations. The transformed nonlinear boundary value problem is solved using the homotopy analysis method (HAM), a semicomputational procedure achieving fast convergence. Computations are verified with an Adomian decomposition method (ADM). The influence of acceleration parameter, rotational body force parameter, Brownian motion number, thermophoresis number, Lewis number, and Prandtl number on surface shear stress, heat, and mass (nanoparticle volume fraction) transfer rates is evaluated. The influence on boundary layer behavior is also investigated. HAM demonstrates excellent stability and leads to highly accurate solutions.


Introduction
Nanofluids continue to stimulate significant interest in modern engineering and medical sciences [1]. These fluids offer significant thermal enhancement characteristics. Nanofluids contain suspended metallic nanoparticles, which increase the thermal conductivity of the base fluid by a substantial amount. Fabrication techniques for nanofluids are also being continuously refined [2]. The heat transfer coefficient of nanofluids increases with volume concentration and they offer other advantages. Nanofluids have been deployed in a tremendous spectrum of applications including sterilization of medical suspensions [3], nanomaterial processing [4,5], automotive coolants [6], microbial fuel cell technology [7], polymer coating [8], intelligent building design [9], microfluid delivery devices [10], and aerospace tribology [11]. Although initially nanofluid dynamics research concentrated on quantification of fundamental thermophysical properties of nanofluids (including thermal conductivity, density, viscosity, and heat transfer coefficient), over the past few years the focus has become increasingly orientated towards modelling and simulation. In this regard building on mathematical models of nanofluid transport, an extensive range of numerical approaches have been implemented to simulate increasingly more practical problems of nanofluids. The inherent nonlinearity of nanofluid flows which involve momentum, heat, and mass transfer necessitates very robust computational algorithms for their analysis. The approaches range from intensive molecular dynamics methods [12,13] which can model interfacial tension between nanoparticles, control volume methods [14], genetic algorithms [15], control volume finite elements [16], Nakamura finite difference codes [17], homotopy analysis techniques [18], ANSYS finite element commercial code [19], differential transform methods [20], MAPLE integration quadrature routines [21], Crank-Nicolson finite difference schemes [22], Keller box implicit methods [23], Mathematica integration subroutines [24], Blottner difference methods [25], and dual reciprocity boundary element methods [26]. These studies have generally examined Brownian motion and thermophoresis effects for various nanoparticle suspensions and considered two-dimensional boundary layer, channel, and cavity flows. They have ignored Coriolis body force effects which arise in rotational fluid mechanics. External flows from spinning bodies are significant in many branches of chemical and industrial engineering including electrolysis treatments [27] and polymer deposition on components [28]. In such systems the rotation strongly influences boundary layer growth and structure on the body periphery. This in turn controls heat and mass transfer rates whether flows are laminar, transitional, or indeed fully turbulent [29]. The bodies may be conical, spherical, elliptical, disk-shaped, and indeed concentric and eccentric systems. Many mathematical and numerical studies of such flows have been presented. Faltas and Saad [30] used a collocation method to analyze steady axisymmetric flow between two spinning eccentric spheres with a linear slip of Basset-type boundary condition at both surfaces. Andersson and Rousselet [31] studied partial momentum and thermal slip from a rotating disk with a Runge-Kutta method. Niazmand and Renksizbulut [32] employed a finite volume code to investigate unsteady heat transfer and thermal patterns around a rotating sphere (as a model of a particle) with surface blowing delineating three distinct wake regimes, namely, steady and axisymmetrical, steady but nonsymmetrical, and unsteady with vortex shedding. They found that although rotation strongly induces local modifications in flow patterns, the surface-averaged heat transfer rates were not altered markedly even at large rotational speeds. Roy and Anilkumar [33] used the Keller box method to simulate transient free and forced convection boundary layer flow from a rotating cone for the case when the free stream angular velocity and the angular velocity of the cone vary arbitrarily with the time. Subhashini et al. [34] used the Bellman-Kalaba quasilinearization method to address uniform slot injection/suction and nonuniform total enthalpy wall effects on steady nonsimilar laminar compressible boundary layer flow over a rotating sphere. They showed that greater rotation and total enthalpy at the wall encourages earlier flow separation whereas cooling delays this and furthermore that greater Mach number displaces the point of separation upstream as a result of the adverse pressure gradient.
The focus of the present work is to analyse rotating nanofluid boundary layer flows at the stagnation point on the external surface of a spinning sphere. Computational rotating nanofluid dynamics has recently attracted some interest, since the deployment of nanofluids in revolving chemical engineering devices offers significant improvements over existing designs. Rana et al. [35] used a variational finite element algorithm to study unsteady magnetonanofluid transport from a rotating stretching continuous sheet. They showed that greater rotational parameter reduces primary and secondary velocities, temperature, and nanoparticle concentration. They further showed that reduced Nusselt number (wall temperature gradient) was suppressed with both Brownian motion and thermophoresis effects whereas reduced Sherwood number (wall mass transfer gradient) was enhanced. Further studies of swirling nanofluid dynamics have been reported by Nadeem and Saleem [36] for a vertical cone and Malvandi [37] for stagnation point nanofluid flow from a spinning sphere. In the present study we use both a homotopy analysis method (HAM) and an Adomian decomposition method (ADM) to simulate stagnation point nanofluid flow from a rotating sphere. The influence of acceleration parameter, rotational body force parameter, Lewis number, Brownian motion number, and thermophoresis parameter on velocity, temperature, and nanoparticle distributions is examined. This work is relevant to coating applications in the polymeric industry.

Mathematical Model
The physical regime under investigation is illustrated in Figure 1, in ( , , ) coordinate system. Transient laminar boundary layer flow of an incompressible Newtonian nanofluid is studied, in the vicinity of the stagnation point region of an isothermal rotating sphere of radius, , rotating with angular velocity, Ω. Soret and Dufour effects are neglected. Following Anilkumar and Roy [38], the free stream and angular velocities depend on time in the form of ( , ) = / and Ω( ) = / , where and are arbitrary constants both greater than zero. The flow field is assumed to be axisymmetric and the fluid possesses constant thermophysical properties with the exception of those caused by density changes which generate the buoyancy force, under the Boussinesq approximation. The Buonjiornio nanofluid model is adopted prioritizing Brownian motion and thermophoresis effects [39]. In light of these approximations, the time-dependent conservation equations for mass, momentum, energy, and species (nanoparticle volume International Journal of Engineering Mathematics 3 fraction) may be presented, as follows, as documented by Malvandi [37]: The appropriate initial and boundary conditions take the form Here , V, and denote velocity components along the , , and coordinates where these coordinates are orientated, respectively, from the forward stagnation point along the surface (see Figure 1), normal to the surface and in the rotating directions, respectively, ( ) is the radial distance from a surface element to the axis of symmetry, is time, is initial condition, is the thermal diffusivity of the nanofluid, ] is kinematic viscosity of nanofluid, is the thermal conductivity of nanofluid, is the Brownian diffusion coefficient (a measure of the species diffusivity of nanoparticles), is the thermophoretic diffusion coefficient, = ( ) /( ) defines the ratio of effective heat capacity of the nanoparticles (e.g., titanium oxide) to the heat capacity of the fluid, is isothermal specific heat capacity and is the nanofluid temperature, ∞ is the free stream temperature (at the edge of the boundary layer), and ∞ is the free stream nanoparticle concentration. It is assumed that the base fluid and the nanoparticles are in thermal equilibrium and no slip occurs between them. We note that, in (4), which is a statement of Fick's law of mass (species) diffusion for nanoparticles, the first term on the left hand side is the transient concentration gradient, and the second and third terms are the convective mass transfer terms. The first term on the right hand side denotes the species diffusion and the last term is the relative contribution of thermophoresis to Brownian motion. These effects have also been considered in detail by Dib et al. [40] and Gupta and Saha Ray [41]. The boundary value problem defined in terms of primitive variables may be solved using numerical methods. However to provide a more amenable solution and one which enables the use of scaling parameters it is advantageous to introduce a set of similarity transformation variables. We define the following: Here is the acceleration (unsteadiness) parameter, is the transformed normal coordinate (perpendicular to sphere surface), is the normalized rotation parameter, is sphere radius, is dimensional stream function, is dimensionless stream function, is the dimensionless secondary velocity function, is the dimensionless temperature function, is dimensionless nanoparticle concentration function (volume fraction), Pr is Prandtl number, Le is Lewis number, Nb is Brownian motion parameter, and Nt is the thermophoresis parameter. Introducing these transformations into (1)-(5), the partial differential boundary layer equations contract to the following nonlinear coupled system of self-similar ordinary differential equations; namely, 1 Pr The corresponding transformed boundary conditions are specified as follows: As → ∞; In engineering simulations of nanofluid flows, not only the velocity, temperature, and nanoparticle volume fraction distributions but also primary ( -) skin friction and secondary ( -) skin friction coefficients, , and local Nusselt number function are important. We define these as Here Re = /] = 2 /] is the local Reynolds number. The set of ordinary differential equations (8)-(11) are highly nonlinear and purely analytical solutions are difficult if not intractable. An efficient homotopy analysis method (HAM) is therefore adopted. Solutions are validated with the Adomian decomposition method (ADM). Although full solutions were given based on a shooting algorithm in Malvandi [37], unfortunately the exact data needed for a comparison is not available in that work. Therefore another objective of the present study is to provide a dual approach for validated solutions both with HAM and ADM techniques, in order to document correct solutions for other researchers to utilize. This therefore allows future researchers who may wish to extend the model to, for example, magnetohydrodynamics, proper benchmark data with which to validate their own numerical methods.

Homotopy Analysis Method (HAM) Simulation
HAM has emerged as a significant alternative to conventional numerical methods for nonlinear systems of partial or ordinary differential equations. Liao [42] employed the basic ideas of homotopy in topology to propose an alternative and general analytical-numerical method for nonlinear problems, namely, the Homotopy Analysis Method (HAM). The validity of HAM is independent of whether or not there exist small parameters in the considered equation(s). Therefore HAM can overcome the foregoing restrictions of perturbation methods. In recent years, HAM has been successfully employed to solve many types of nonlinear problems in engineering sciences including thin membrane reactiondiffusion phenomena [43] and porous media convection [44]. We provided details of the application of HAM to the system of (8)- (11) in the next section.

Homotopy Analysis Method (HAM).
We write the initial guesses and linear operators as International Journal of Engineering Mathematics with the following properties: where ( = 1-9) are arbitrary constants. Let ∈ [0, 1] represent an embedding parameter and ℏ , ℏ , ℏ , ℏ denote the nonzero auxiliary linear operators and construct the following zeroth-order deformation equations: subject to the boundary conditionŝ The auxiliary parameters are properly chosen so that series (24) converge at = 1 and thus The resulting problems at the th order deformation are The general solution of (26) is where * ( ), * ( ), * ( ), and * ( ) are the particular solutions and the constants are to be determined by boundary conditions (27). (25) gives an analytical solution of the problem in series form. The convergence of the series solution given by HAM depends strongly upon auxiliary parameters ℏ , ℏ , ℏ , and ℏ . These parameters provide a convenient mechanism for adjusting and controlling the convergence region and convergence rate of the series solution. Therefore in order to select appropriate values for these auxiliary parameters, the so-called ℏ , ℏ , ℏ , and ℏ curves are displayed at 20th-order approximations as shown in Figure 2. This achieves excellent accuracy and is therefore adopted in all HAM numerical computations. HAM is executed in a symbolic code to investigate the influence of Further computations for primary skin friction, secondary skin friction, wall heat transfer, and wall mass transfer rate are presented in Tables 1-3, where a comparison is also given with the ADM algorithm [45], discussed in the next section.

Validation with Adomian Decomposition Method (ADM)
Validation of the HAM computations is achieved with ADM, a seminumerical technique which employs Adomian polynomials to achieve very accurate solutions which may be evaluated using symbolic packages such as Mathematica.  Introduced by American mathematician, George Adomian [46], it has been embraced extensively in computational engineering sciences over the past two decades. Interesting studies using ADM include enzyme kinetics in biological engineering [47], heat transfer [48], structural damping systems [49], non-Newtonian foam drainage problems [50],   and most recently magnetic biotribology [51] and nanofluid squeezing flows [40,41,52]. ADM [46] deploys an infinite series solution for the unknown functions, that is, , , , and , and utilizes recursive relations. The present ordinary differential nonlinear boundary value problem (BVP) is International Journal of Engineering Mathematics   rewritten using the standard operator, following Bég et al. [51]: where is the unknown function, is the highest-order derivative (assumed to be easily invertible), is a linear differential operator of order less than , designates the nonlinear terms, and is the source term. Applying the   inverse operator −1 to both sides of (34) and using the given conditions we obtain  where V represents the terms arising from integrating the source term and from the auxiliary conditions. ADM defines solution by the series

The solution for the nonlinear terms is
Here are the Adomian polynomials which are evaluated via the following relation [51]: If the nonlinear term is expressed as a nonlinear function ( ), the Adomian polynomials are arranged into the form Components 0 , 1 , 2 . . . are then determined recursively by using the relation where 0 is referred to as the zeroth component. Ancomponents truncated series solution is finally obtained as     Decomposition series (41) converges exceptionally fast, in particular on high memory dual processor machines [53]. The rapid convergence means that relatively few terms are required to obtain an approximate analytical solution. This is a considerable advantage of the ADM approach compared with other semianalytical methods such as perturbation expansions. The present HAM accuracy is compared with the ADM solutions in Tables 1-2. Excellent agreement guarantees confidence in the HAM computations. Further HAM computations are given in Table 3 for variation of all the control parameters. Figure 3 illustrates the influence of the Lewis number (Le) on the nanoparticle concentration distribution. Although the effects of this parameter were investigated on velocity functions and temperature, no tangible modifications were observed and therefore are not discussed further. Lewis number quantifies the ratio the thermal diffusivity to the mass diffusivity. An increase of Lewis number corresponds to a lower species diffusivity of the nanoparticles ( ) for a prescribed thermal diffusivity ( ). For this reason a rise in Le induces a significant reduction in the dimensionless nanoparticle volume fraction. For Le > 1 thermal diffusivity exceeds the species diffusivity and vice versa for Le < 1. For Le = 1 both thermal and species diffusivity will be the same, and thermal and nanoparticle concentration boundary layer thicknesses will be equal. Concentration boundary layer thickness for the nanoparticle species is significantly reduced with greater Lewis number. With greater Lewis number the decay in concentration profiles also progressively evolves from a linear descent (from the maximum at the sphere surface to zero in the free stream) to a more monotonic profile. In all cases asymptotically smooth profiles are achieved with HAM testifying to the prescription of a suitably large value for infinity, that is, 5. Figures 4-7 depict the influence of the unsteadiness parameter ( ) on flow characteristics. Primary velocity, ( ) (Figure 4), is observed to be strongly accelerated with greater values. The rotation of the sphere draws momentum from the -direction and redistributes this in the -direction. Secondary velocity ( Figure 5) is therefore strongly decreased with greater unsteadiness parameter. The primary and secondary profiles are also very different. Primary velocity grows with greater distance from the sphere surface attaining maxima in the free stream. Secondary velocity, ( ), however decays from a maximum at the sphere surface (wall) to vanish in the free stream. The rotation of the sphere acts like a fan drawing momentum from one direction and channeling it into another. Effectively as the -direction flow is accelerated the momentum boundary layer thickness is decreased. These computations concur well with the trends of Malvandi [37], although an erroneous interpretation is given in that paper. In Figures 6 and 7 both temperature and nanoparticle concentration are found to be strongly depressed with increasing acceleration parameter. Although arises multiple times in the primary and secondary momenta equations (8) and (9), it features only in a single term in each of the energy and concentration equations (10) and (11), specifically / in International Journal of Engineering Mathematics (10) and Le / in (11). These terms couple the thermal and species diffusion to the primary momentum field only.

Discussion and Interpretation of Results
The key influence from increasing unsteadiness is therefore an acceleration in primary flow which will counteract both heat and nanoparticle diffusion. Effectively as the primary momentum boundary layer is thickened, the decrease in temperatures will cool the boundary layer and reduce thermal boundary layer thickness. Species (nanoparticle) boundary layer thickness will also be reduced. Both temperature and concentration distributions exhibit a consistent descent from the wall (sphere surface) to the free stream. However the decay in temperatures is more gradual compared with concentrations which plummet more sharply. Generally the unsteadiness is found to induce a nontrivial influence on all flow characteristics.
In Figures 8-11, the effects of spin (rotation) parameter, (= Ω / ), on velocity functions, temperature, and nanoparticle concentration are depicted. This parameter embodies the influence of the secondary velocity field on the primary velocity field, that is, via the swirl effect. It is directly proportional to the rotational velocity of the sphere and arises in the coupling term 2 in (8). For the case of a stationary sphere, Ω → 0 and → 0, and for this scenario the primary flow (Figure 8) is weakest and the secondary flow ( Figure 9) is strongest. As is increased, the rotation becomes more intense and this boosts primary momentum leading to escalation in / values. The converse response is computed for secondary flow which is suppressed with greater values. The reduction in secondary flow however is weaker than the growth in the primary flow. Temperatures ( Figure 10) are found to be weakly reduced with greater rotation effect, implying a slight thinning in thermal boundary layers. Similarly nanoparticle concentration ( Figure 11) is also marginally decreased with increasing values. Heat and mass transfer are therefore weakly resisted with greater rotation. They are maximized for the stationary sphere case. Better control of thermal and species diffusion is achieved with rotation of the sphere. This may be beneficial therefore in spin coating operations employing nanomaterials. Figures 12 and 13 illustrate the response of temperature ( ) and species concentration ( ) to a change in Brownian motion parameter (Nb). Temperature is slightly increased as Nb is increased. The reverse trend is noticed in the case of concentration. Physically smaller nanoparticles yield higher Nb values, which assist in thermal diffusion in the boundary layer via increased thermal conduction. On the contrary, larger nanoparticles show lower Nb values and this depresses thermal conduction. Higher Nb values will conversely stifle the diffusion of nanoparticles away from the surface into the fluid regime which will manifest in a decrease in nanoparticle concentration values in the boundary layer. The distribution of nanoparticles in the boundary layer regime can therefore be regulated via the Brownian motion mechanism (higher Nb values) and cooling of the regime can also be achieved via larger Nb values. Heat transfer from the fluid to the sphere surface (wall) is promoted with higher Nb values. Thicker thermal boundary layers are produced with higher Nb values whereas larger concentration boundary layer thickness is associated with lower Nb values. The influence of Brownian motion on the velocity fields was found to be inconsequential and these plots are therefore excluded here.
Finally Figures 14 and 15 illustrate the effects of thermophoresis parameter (Nt) on temperature and nanoparticle concentration distributions. Increasing thermophoresis effect (greater Nt values) slightly elevates nanofluid temperatures ( Figure 14). Higher Nt values also increase nanoparticle concentrations, since lesser particle deposition will occur at the wall and greater migration of nanoparticles from the wall to the fluid regime will result. Thermal boundary layer thickness is slightly increased with thermophoresis whereas concentration boundary layer thickness is more significantly enhanced. It is further noted that the strongest influence of thermophoresis on nanoparticle distribution is at intermediate distances from the sphere, transverse to the sphere surface. Table 3 documents the influence of many parameters on the skin friction components, heat and mass transfer rates. With greater rotation effects ( ), primary skin friction ( (0)) is strongly elevated whereas secondary skin friction (− (0)) is weakly elevated. There is also a weak increase in the surface heat (− (0)) and mass transfer (− (0)) rates. With an increase in Prandtl number (Pr), skin friction components are unaffected whereas heat transfer rate is strongly increased and mass transfer rate (nanoparticle diffusion rate at the sphere surface) is decreased. Cooling is therefore achieved successfully in the rotating boundary layer regime with larger Prandtl number, Pr (decreasing nanofluid thermal conductivity), since more heat is conducted away from the fluid to the sphere. This is one of the main attractions of nanofluids. Greater thermophoresis (Nt) boosts the heat transfer rate whereas it decreases the mass transfer rate. It exerts no tangible influence on the skin friction magnitudes. Greater Brownian motion effect (Nb) decreases wall heat transfer rate but elevates the mass transfer rate. Increasing unsteadiness parameter ( ) enhances both primary and secondary friction and furthermore increases both heat and mass transfer rates. Greater Lewis number (Le) results in a reduction in the surface heat transfer rate and increase in the surface mass transfer rate but does not alter the primary or secondary skin friction components.

Conclusions
Computational algorithms have been developed to study the transient nanofluid flow in the stagnation region from a spinning spherical body. The Buongiorno model has been employed to simulate nanoparticle Brownian motion and thermophoresis effects, for the case of dilute nanofluids. The nonlinear boundary value problem has been solved with HAM. ADM has also been used to verify the HAM solutions. The computations have shown that, with greater rotation effect, the primary flow is enhanced whereas the secondary flow is weakened. With increasing unsteadiness, both primary and secondary velocity fields are aided, as are the wall heat and mass transfer rates. An increase in nanoscale parameters (Brownian motion and thermophoresis) is found to mainly influence the temperature and nanoparticle distributions, although a slight alteration is computed in surface skin friction components. Thermophoresis tends to enhance the wall heat transfer rate and reduces the mass transfer rate. Brownian motion exerts the opposite influence to thermophoresis. The current study is relevant to nanotechnological coating applications in the polymer industry. In this study we have employed a Newtonian nanofluid model. Future investigations will use non-Newtonian nanofluid models (e.g., micropolar theory) [8] and will be communicated imminently. Furthermore, the current study, it is envisaged, has demonstrated the advantage of HAM in being able to achieve very high order approximations in symbolic packages. It is a computer-extended series expansion method, a modern analogy to Van Dykes asymptotic expansion/perturbation series method of the 1970s (which was used in inviscid and viscous supersonic flows). The popularity of this method among Eastern researchers is immense. However very few British researchers have explored this technique. Although HAM is algebraically laborious it is nevertheless an elegant approach and avoids the traditional pitfall of other numerical schemes, namely, the time-consuming nature of discretization processes. We hope that the present paper will further popularize the scheme with British researchers who may not have encountered it thus far.