Hydromagnetic Flow and Heat Transfer of a Jeffrey Fluid over an Oscillatory Stretching Surface

1 Introduction

The phenomenon of flow and heat transfer of non-Newtonian fluids over continuous stretching sheet plays a vital role in many industrial processes. The most valuable applications of flow of non-Newtonian fluids over stretching sheet include paper production, glass blowing, wire drawing, and aerodynamic extrusion of plastic sheet. The fundamental work on flow over continuous moving sheet was provided by Sakiadis [1]. Crane [2] extended the work of Sakiadis [1] for a stretching sheet. The study of flows of non-Newtonian fluids has attracted the attention of many researchers because of its wide range of applications in technology and industry. The dynamics of non-Newtonian fluids is entirely different from that of Newtonian fluids. This is because non-Newtonian fluids exhibit many typical characteristics such as shear-thinning or thickening, yield stress, the Weissenberg effects, fluid memory, and die swelling, which are absent in Newtonian fluids. Therefore, these fluids do not obey the well-known Newton’s law of viscosity, which is valid for Newtonian fluids. Due to the diversity of non-Newtonian fluids, many constitutive relations are proposed in the literature. Fluid models in which viscosity depends on the instantaneous rate of deformation are known as generalised Newtonian fluid (GNF) models. These models are capable of describing shear thinning, shear thickening, and yield shear stress effects. The stretching sheet problems using GNF models were analysed by Anderson et al. [3], Hassanien [4], Chen [5], Nadeem et al. [6, 7], Mukhopadhyay [8], etc. The Weissenberg or rod-climbing effect is associated with nonlinear effects and normal stresses. Such an effect is predicted by those models for which normal stress differences are non-zero. Second-grade fluid is one of them. Studies pertaining to the flow of second-grade fluid over a stretching sheet can be found in [9–20]. To account for viscoelastic effects like fluid memory or die swelling, models like linear/convected Maxwell and Jeffrey models have been proposed in the literature. Some studies regarding flow over a stretching sheet using these models were also carried out in the literature [21–28]. Fluid models that can simultaneously predict shear thinning/shear thickening and normal stress effects have also been used to study flow over stretching sheets. For instance, Sajid and Hayat [29] discussed the boundary layer of a third-grade fluid over a stretching sheet. In another study, they carried out flow and heat transfer analysis of a third-grade fluid over a stretching surface [30].

The literature reveals that very little attention is paid to the flow and heat transfer analysis due to the oscillatory stretching sheet. The viscous flow analysis over the oscillatory stretching sheet was initially performed by Wang [31]. Later, Abbas et al. [32] extended Wang’s problem by including wall slip and performing heat transfer analysis. Furthermore they [33], complemented Wang’s problem by carrying out a study of the unsteady magnetohydrodynamic (MHD) boundary flow of a second-grade fluid over an oscillatory stretching surface. Some recent attempts regarding flow due to oscillatory stretching sheet can be found in [34–36]. Despite recent developments in non-Newtonian fluids, it is still of interest to analyse stretching flows involving non-Newtonian fluids. For instance, no attempts are available in the literature that deal with the flow and heat transfer of Jeffrey fluid over an oscillatory stretching surface. The purpose of this article is to provide such an analysis. The Jeffrey fluid model is a fairly simple linear model using a time derivative instead of a convected derivative as used in the Oldroyd-B model. This model has three material parameters: the viscosity, the ratio of relaxation to retardation time, and retardation time. The Jeffrey model is capable of predicting elastic and memory effects exhibited by dilute polymer solutions and biological liquids. Relaxation time is a time constant associated with viscoelastic fluids. For viscoelastic fluids, the internal molecular configuration of the fluid can sustain stress for some time. This time, called relaxation time, varies widely among materials. In fact, the Jeffrey model is a simple extension of the Newtonian model to include elastic and memory effects. By choosing the Jeffrey fluid model, it becomes possible to study the oscillatory stretching sheet problem for various Newtonian and non-Newtonian fluids. However, this model is not capable of predicting shear-thinning and shear-thickening effects and thus cannot be used to discuss the fluids that are driven to flow at high shear rates.

The problem presented here has many useful engineering and industrial applications. In fact, the production of sheeting material is involved in many manufacturing processes and includes both metal and polymer sheets. The rate of heat transfer over a surface has to play a useful role in the quality of the final product. Industrial applications of stretching sheet include hot rolling, fibers spinning, manufacturing of plastic and rubber sheet, continuous casting, and glass blowing. The present study also describes MHD effects that are very prominent in MHD power-generating systems, plasma studies, cooling of nuclear reactors, geothermal energy extraction, and many others.

The structure of the article is as follows. The problem is formulated with underlying assumptions in section 2. A solution using the homotopy analysis method (HAM) is provided in Section 3. A scheme of the numerical solution is illustrated in Section 4. The convergence of the HAM solution and its comparison with a numerical solution are presented in Section 5. A discussion of the results obtained is given in Section 6. The main conclusions of the study are summarised in Section 7.

2 Flow Analysis

We consider the unsteady, two-dimensional, and MHDs flow of an incompressible Jeffrey fluid over an oscillatory stretching sheet coinciding with the plane y̅=0. The elastic sheet is stretched back and forth periodically with velocity u_ω=bx̅sinωt (b is the maximum stretching rate, x̅ is the coordinate along the sheet, and ω represents the frequency). A constant magnetic field of strength B₀ is applied perpendicular to the stretching surface. The magnetic Reynolds number is assumed to be small so that the induced magnetic field can be neglected. Let T_w be the temperature of the surface of the sheet and T_∞ is the temperature of the fluid far away from the sheet, where T_w> T_∞. A schematic diagram of the flow geometry is illustrated in Figure 1. In addition to these assumptions, the boundary-layer approximations are used and viscous dissipation effects are neglected. Thus, the governing equations for the flow under consideration in the presence of magnetic field are [28]

(1)∂u∂x¯ + ∂v∂y¯ = 0, (1)

(2)∂u∂t + u∂u∂x¯ + v∂u∂y¯ = ν1 + λ[∂2u∂y¯2 + λ1(∂3u∂y¯2∂t + u∂3u∂x¯∂y¯2− ∂u∂x¯∂2u∂y¯2 + ∂u∂y¯∂2u∂x¯∂y¯ + v∂3u∂y¯3)] − σB02ρu, (2)

(3)ρcp(∂T∂t + u∂T∂x¯ + v∂T∂y¯) = k∂2T∂y¯2, (3)

Figure 1:

Geometry of the problem.

where u and v are velocity components along x̅ and y̅ directions, respectively. In the above equations, ν represents the kinematic viscosity, which is the ratio of the dynamic viscosity of a fluid to its density; λ is the ratio of relaxation to retardation time (ratio parameter); λ₁ is the retardation parameter; T represents the temperature; c_p is the specific heat, which shows the relationship between heat and temperature changes; and k is the thermal conductivity, which primarily appears in Fourier’s law for heat conduction. The thermal conductivity is the property of a material’s ability to conduct heat.

The boundary conditions to be satisfied by the velocity components and temperature are

(4)u = uω = bx¯sinωt,  v = 0,  T = Tw  at  y¯ = 0, t > 0, (4)

(5)u → 0,  ∂u∂y¯ → 0,  T → T∞  as  y¯ → ∞. (5)

The second boundary condition in (5) is the augmented boundary condition [33] and is imposed due to the consideration of the flow in the unbounded domain. To nondimensionalise the flow problem, we use the following dimensionless variables [31, 32]:

(6)u = bx¯fy(y, τ), v = −νbf(y, τ), y = bνy¯ and τ = tω, (6)

(7)θ(y, τ) = T − T∞Tw − T∞. (7)

With the help of (6) and (7), the continuity equation is identically satisfied, and (2) and (3) take the form

(8)(1 + λ)(Sfyτ + fy2 − ffyy + M2fy) = fyyy + β(Sfyyyτ + fyy2 − ffyyyy), (8)

(9)θyy + Pr(fθy − Sθτ) = 0, (9)

where subscripts denote differentiation. It is pointed out that a reduction in the number of independent variables by one is also achieved as an implication of transforming (2) and (3) in the dimensionless form. The boundary conditions (4) and (5) becomes

(10)fy(y, τ) = sinτ,  f(y, τ) = 0,  θ(y, τ) = 1 at y = 0, (10)

(11)fy(y, τ) → 0, fyy(y, τ) → 0, θ(y, τ) → 0 as y → ∞. (11)

Here, S≡ ω/b is the ratio of oscillation frequency of the sheet to its stretching rate, β=λ₁b represents the Deborah number, Pr=μc_p/k is the Prandtl number, and M = σB02/ρb represents the Hartmann number.

We define the skin friction coefficient C_f and the local Nusselt number Nu_x as

(12)Cf = τwρuw2, Nux = x¯qwk(Tw − T∞), (12)

where τ_w and q_w are the shear stress and heat flux at wall, respectively, which are defined as [28]

(13)τw = μ1 + λ[∂u∂y¯ + λ1(∂2u∂t∂y¯ + u∂2u∂y¯∂x¯ + v∂2u∂y¯2)]y¯ = 0, qw = −k(∂T∂y¯)y¯ = 0. (13)

In view of (6) and (7), (12) gives

(14)Rex1/2Cf = 11 + λ[fyy + β(Sfyyτ + fyfyy − ffyyy)]y¯ = 0, Rex−1/2Nux = −θy(0, τ), (14)

where Rex = uwx¯/ν denotes the local Reynolds number.

3 Solution of the Problem

We have employed two methods for obtaining the solution of (8) and (9) subject to boundary conditions (10) and (11). Firstly, we have used HAM to solve nonlinear partial differential equations. The procedural details of this method can be found in the latest book by Liao [37]. Secondly, a numerical method based on an implicit finite difference scheme is used. The details of numerical method are explained in the next section.

4 Direct Numerical Solution of the Problem

The set of nonlinear partial differential equations (8) and (9) with boundary conditions (10) and (11) are solved numerically by using the finite difference scheme. To transform the semi-infinite physical domain y∈ [0, ∞) to a finite calculation domain η∈ [0, 1], a coordinate transformation η=1/(y + 1) is employed. With the help of this coordinate transformation, (8) and (9) become

(15)S((1 + λ) − 6βη2)∂2f∂τ∂η − Sβη4∂4f∂η3∂τ − 6Sβη3∂3f∂η2∂τ = ((1+λ)η2−4βη4)(∂f∂η)2+ (6η2 + 24βη3f − (1 + λ)M2 − 2(1 + λ)ηf)∂f∂η + η4∂3f∂η3+ (6η3 − (1 + λ)η2f + 36βη4f)∂2f∂η2− 4βη5∂f∂η∂2f∂η2 − βη6(∂2f∂η2)2 + 12βη5f∂3f∂η3 + βη6f∂4f∂η4, (15)

(16)η4∂2θ∂η2 + 2η3∂θ∂η − Pr(fη2∂θ∂η + S∂θ∂τ) = 0, (16)

(17)fη = 0, fηη = 0, θ = 0 at η = 0, (17)

(18)f = 0,  fη = −sinτ,  θ = 1  at  η = 1. (18)

Equations (15) and (16) are discretised for L uniformly distributed discrete points in η=(η₁, η₂, ……, η_{L})∈ (0, 1) with a space grid size of Δη=1/(L + 1) and time level t=(t¹, t², …). Hence, the discrete values (f1n, f2n, …, fLn) and (θ1n, θ2n, …, θLn) at these grid points for time levels tⁿ=nΔt (Δt is the time step size) can be numerically obtained by solving a system of linear equations. Such a system of linear equations can be obtained by constructing a semi-implicit time difference for f and θ and assuming that only linear equation for the new time step (n + 1) needs to be solved. The following semi-implicit time difference for f and θ is constructed:

(19)S((1 + λ) − 6βη2)1Δt(∂f(n + 1)∂η − ∂f(n)∂η) − Sβη41Δt(∂3f(n + 1)∂η3 − ∂3f(n)∂η3)− 6Sβη31Δt(∂2f(n + 1)∂η2 − ∂2f(n)∂η2) = ((1 + λ)η2 − 4βη4)(∂f(n)∂η)2 + (6η2 − (1 + λ)M2)∂f(n + 1)∂η  + (24βη3 − 21 + λ)η)f(n)∂f(n)∂η + 6η3∂2f(n + 1)∂η2 + (36βη4 − (1 + λ)η2)f(n)∂2f(n)∂η2  − 4βη5∂f(n)∂η∂2f(n)∂η2 − βη6(∂2f(n)∂η2)2 + η4∂3f(n + 1)∂η3 + 12βη5f(n)∂3f(n)∂η3 + βη6f(n)∂4f(n)∂η4. (19)

(20)SPr(θ(n + 1) − θ(n))Δt = η4∂2θ(n + 1)∂η2 + 2η3∂θ(n + 1)∂η − Prf(n)η2∂θ(n + 1)∂η. (20)

Equations (19) and (20) are converted into a system of linear equations by means of the finite difference method, which can be solved by any method, e.g. Gaussian elimination. The simulations are started from motionless velocity field and a uniform temperature distribution equal to temperature at infinity as

(21)f(η, τ = 0) = 0 and θ(η, τ = 0) = 0. (21)

5 Convergence of HAM Solution and its Comparison with a Numerical Solution

It is a well-established fact that the convergence of the HAM solution is largely dependent on the proper choice of the auxiliary parameters h_f and h_θ. For a particular set of parameters, the convergence region can be obtained by plotting h-curves. Figure 2a and b presents two such curves showing the plausible values of h_f and h_θ for a given set of parameters. This figure indicates for convergent solution –1.1 ≤ h_f< 0 and –1.7 ≤ h_θ< –0.2. Table 1 presents values of f″(0, τ) at different order of approximations for a specific set of parameters. It is observed from this table that the convergent values of f″(0, τ) are achieved at the 10^th-order of approximation. A comparison of the HAM solution with a numerical solution is also presented in this section in Figures 3 and 4. These figures show an excellent agreement between both the solutions as the order of approximation increases. This is because at higher order of approximation, the HAM solution becomes closer to the exact solution. The values of f″(0, τ) obtained by HAM and numerical solution for a Newtonian fluid are also compared with the corresponding values in [33, 35]. Table 2 shows the comparison of such values. It can be seen that our results are in excellent agreement with those presented in [33, 35].

Figure 2:

The h-curves at 10^th-order of approximation (a) for velocity and (b) temperature profile.

Figure 3:

Comparison of f′(y, τ) obtained from the HAM solution (solid lines) and the numerical solution (open circles).

Figure 4:

Comparison of θ(y, τ) obtained from the HAM solution (solid lines) and the numerical solution (open circles).

6 Results and Discussion

In this section, the effects of pertinent parameters of interest on the velocity component f′, skin friction coefficient, temperature profile θ, and the Nusselt number are illustrated through various plots. The effects of the magnetic field on various flow characteristics are not shown to reduce the length of the article. The interested reader may consult [33, 35] for illustration of such effects in similar situations.

Figure 5 shows the effects of the Deborah number β and ratio of relaxation to retardation parameter λ on dimensionless velocity at a fixed distance y=0.25 from the sheet as a function of time. In Figure 5a, the effects of the Deborah number β are illustrated. This figure shows an increase in the amplitude of velocity f′ by increasing the Deborah number β. Moreover, a phase difference occurs for the non-zero value of β. Thus, a viscoelastic fluid with larger retardation time oscillates in time at a fixed location from the sheet with a different amplitude and phase than that of a Newtonian fluid. Figure 5b shows the effects of λ on the velocity profile f′. It is observed that the amplitude of flow oscillations decreases by increasing the ratio parameter.

Figure 5:

The velocity profile f′ as a function of time in the first five periods τ ∈ [0, 10π] at a fixed distance from the sheet, y=0.25: (a) effects of β and (b) effects of λ.

The transverse profiles of velocity f′ for various values of β and at four different time instants in the fifth period are shown in Figure 6a–d. Figure 6a shows that at τ=8.5π, the velocity f′ increases by increasing β. Moreover, at this time instant, there are no oscillations in the velocity profile f′. At time instant τ=9π (Fig. 6b), the velocity f′ oscillates near the wall and amplitude of oscillation increases by increasing β. Figure 6c elucidates that at τ=9.5π, velocity f′ increases from –1 at the sheet to zero at infinity. Moreover, f′ decreases by increasing β at this time instant. When τ=10π, again oscillation in f′ is observed near the oscillating sheet and its amplitude increases with β.

Figure 6:

Transverse profiles of velocity field f′ for different values of β in the fifth period τ ∈ [8π, 10π] with S=1, M=12, and λ=0.1.

Figure 7 shows the effects of ratio of relaxation to retardation time, λ on the transverse profiles of velocity f′ at four different instants τ=8.5π, τ=9π, τ=9.5π, and τ=10π. It is observed from this figure that, for τ=8.5π, velocity f′ decreases from unity at the wall to zero far away from the wall. Moreover, it decreases by increasing λ. At time instant τ=9π, f′ has zero value both at the wall and far away from the wall. Furthermore, it oscillates near the wall and amplitude of oscillation decreases by increasing λ. Similar results can be observed at τ=10π. However, at τ=9.5π, f′ decreases from –1 at the wall to zero far away from the wall. Further, its magnitude decreases by increasing λ.

Figure 7:

Transverse profiles of velocity field f′ for different values λ in the fifth period τ ∈ [8π, 10π] with S=0.5, M=12, and β=0.1.

The variation of skin friction coefficient Rex1/2Cf with time for different values of β and λ is depicted in Figure 8. We observe from this figure that, like velocity profile f′, skin friction coefficient also oscillates with time and amplitude of oscillation increases by increasing β and λ. Further, like velocity profile f′, the oscillations in the skin friction for β=0 are not in phase with the oscillations for non-zero of β. In fact, a phase lag is observed in skin friction for the non-zero value of β. The profile of temperature θ for different values of the Prandtl number Pr, β, M, and λ are shown in Figure 9. It is observed through this figure that the temperature θ of the of the fluid decreases by increasing the Prandtl number and β. However, it increases by increasing M and λ. The variation of the local Nusselt number with time for different values of Pr, β, M, and λ is displayed in Figure 10a–d. Figure 10a illustrates the effects of the Prandtl number on the local Nusselt number Rex−1/2Nux. It is noted that the local Nusselt number increases with the increase of Pr. In Figure 10b, the effects of the Deborah number β are illustrated. This figure predicts an increase in the local Nusselt number by increasing the Deborah number β. The profiles of the local Nusselt number follow a reverse trend in Figure 10c and d. This figure indicates that the local Nusselt number decreases by increasing ratio parameter λ and the Hartmann number M.

Figure 8:

The skin friction coefficient Rex1/2Cf in first five periods τ ∈ [0, 10π]: (a) effects of β and (b) effects of λ.

Figure 9:

Transverse profiles of the temperature field θ at the time point τ=8π (a) effects of Pr, (b) effects of β, (c) effects of M, and (d) effects of λ.

Figure 10:

Variation of the local Nusselt number Rex−1/2Nux in the first five periods τ ∈ [0, 10]: (a) effects of Pr, (b) effects of β, (c) effects of M, and (d) effects of λ.

7 Concluding Remarks

In this article, we have studied the unsteady flow and heat transfer of Jeffrey fluid over an oscillatory stretching surface, which is maintained at a constant temperature. The resulting nonlinear partial differential equations are solved by using HAM and a finite difference scheme. The main conclusions of the present study are as follows:

• The amplitude of flow velocity increases by increasing the Deborah number β, while it decreases with ratio of relaxation to retardation time λ. Moreover, a phase shift is observed for non-zero values of the Deborah number β.

• The amplitude of skin friction coefficient increases by increasing the Deborah number β and ratio of relaxation to retardation time λ.

• An increase in temperature is found by increasing the Hartmann number M and ratio of relaxation to retardation time λ. However, it decreases by increasing the Deborah number and the Prandtl number Pr. Moreover, the rate of heat transfer increases by increasing the Prandtl number Pr and the Deborah number β.

• The solution presented in this article is more general and includes the solution corresponding to radiative problem as a special case. This conclusion is also supported by recent papers by Magyari and Pantokratoras [38], Fetecau et al. [39], and Vieru et al. [40], where it is explicitly shown that an evaluation of the effect of thermal radiation in the linearized Rosseland approximation does not require any additional research effort. Once the problem has been solved for a comprehensive set of values of the Prandtl number in the absence of radiation, it has been automatically solved for the radiative case too. In view of above argument, the solution of problem with radiative effects can be obtained by replacing the Prandtl number in our solution by the effective Prandtl number Pr_eff=Pr/(1 + N_r), where N_r is the radiation parameter.