Flow and Heat Transfer of Bingham Plastic Fluid over a Rotating Disk with Variable Thickness

1 Introduction

Bingham plastic fluid is featured by the existence of the threshold stress level (known as yield stress proposed by Bingham and Green [1]) which leads to the fluid holding at rest like elastic solid for the extra stress below this stress level. Upon the force exerted on the fluid exceeding the value of the yield stress, the fluid will exhibit normal shear flow. The Bingham plastic fluid has received considerable attention for their extensive applications in industry and engineering such as polymer solutions, paints, clay coatings and cosmetic products, etc. The flow, heat and mass transfer of Bingham fluid in laminar boundary layer over a steady rotating disk were investigated by Rashaida [2] and Rashaida et al. [3]. Then, Osalusi et al. [4], [5] considered the joint effects of viscous dissipation and Joule heating on the steady flow and heat transfer of magnetohydrodynamics (MHD) type Bingham fluid over a porous rotating disk under Hall and ion-slip currents. Nirmalkar et al. [6] presented the 2-D creeping flow of Bingham plastic fluids past a cylinder of square cross-section by finite element method. Later, they further studied the flow of Bingham plastic fluids past an elliptical cylinder [7], [8]. As a matter of fact, no study has yet considered the physical model of Bingham plastic fluid over a rotating disk with variable thickness radially.

The steady flow of Newtonian fluid over a rotating disk was first discussed by von Kármán in 1921 [9] with introducing the well-known von Kármán similar transformation, which can reduce complicated partial differential equations into ordinary differential equations. From then on, many studies had considered rotating disks with uniform thickness such as Turkyilmazoglu [10], [11], [12], [13], Ming et al. [14], Yin et al. [15], etc., while in 2002, Eraslan and Orcan [16] and Eraslan and Argeso [17] studied the elastic-plastic stress distribution in rotating variable thickness solid disks with limit angular velocities. The disk is symmetric with respect to the mid-plane, and its profile is assumed to vary radially in exponential form h(r)=h0e−n(rb)k or power form h(r)=(1−nrb)k. Actually, in fluid mechanics, some authors had also developed models for these problems of considering variable thickness boundary. Lee [18] first derived the governing equations of incompressible fluid flowing axially over a thin paraboloid. Due to the special shape of the boundary, the similarity solution exists, and thereby it provides the opportunity to investigate the problem in detail. Recently, Salahuddin et al. [19] presented a numerical investigation of MHD flow of Cattanneo–Christov heat flux model for Williamson fluid over a stretching sheet with variable thickness subjected to power formula. The objective of this paper is to study the flow and heat transfer of non-Newtonian Bingham fluid over a variable thickness rotating disk. To the best of our knowledge, the model of variable thickness disk with the exponential form has never been considered before.

Homotopy analysis method (HAM), introduced by Liao in 1992 [20], is one of the powerful methods for solving nonlinear equations, and the effectiveness of the HAM has been validated by himself [21], [22], [23] and many other researchers [24], [25], [26]. In this paper, we employ the HAM to solve the semi-similarity ordinary differential equations. The paper is organized as follows: in Section 2, the mathematical formulation of the model is proposed. The detailed similarity reduction procedures for the governing equations are presented in Section 3. In Section 4, the set of equations is solved by using the HAM. The analyses of results and discussions are given in Section 5, followed by conclusions in Section 6.

2 Mathematical Formulation of the Physical Model

Let us consider the motion of viscous incompressible Bingham plastic fluid over a rotating disk with variable thickness radially. The cylindrical coordinate system is adopted to formulate the physical model properly. The flow of Bingham plastic fluid is characterized by two-layered structure consisting of the plug layer and the shear layer [see Fig. 1, the cylindrical coordinate system (R, Φ, Z)]. The continuity equation, momentum equation, and energy equation may be written as follows

Figure 1:

The physical model of a rotating disk with variable thickness.

Andersson et al. [27] obtained the result of ∂P/∂Z=0 in the boundary layer. The motion of the fluid is considered to be fast while maintaining laminar flow, and therefore the viscous effects will be confined within a thin layer near the disk [28]:

1. The component of velocity W is much smaller in magnitude than either of the other two components.

2. The rate of change of any physical quantity in the direction normal to the disk is much greater than its variation in the radial and tangential direction.

3. The only significant fluid stress components are τ_φz and τ_rz.

4. The pressure depends only on the axial distance from the axis of rotation.

5. The fluid has constant properties with negligible dissipation.

Based on the above analysis, the governing equations (1)–(5) are reduced into

where U, V, and W are velocity components in the directions of R, Φ, and Z, respectively. T denotes the temperature of the fluid, ρ represents the density of the fluid, and α is the thermal diffusivity of fluid.

The boundary conditions are

Here a is the variable thickness positive parameter of disk, and m is the shape parameter of disk. R₀ is a characteristic radius which is convenient for the following dimensionless analysis and similarity transformation products. Ω is the angular velocity of the rotating disk, T_w is the temperature of the disk, and T_∞ is the temperature of the ambient fluid.

Generally speaking, the Bingham constitutive equation adequately describes the stress-deformation behavior of materials demonstrating the yield stress. This model relates the rate of deformation tensor e_ij which is defined in terms of the velocity components [2]

to the deviatoric stress tensor τ_ij using the relation below:

while e_ij=0, τ<τ_y. When the magnitude of the extra stress tensor |τ| is greater than the yield stress τ_y, the apparent viscosity is assumed as

where μ_p is the viscosity of the Bingham fluid, referred to as the plastic viscosity. The extra stress tensor and deformation rate are adopted as [3]

where the summation convention denotes repeated indices. Using the approximations (I–V) and the rotational symmetry assumption, one then obtains

where the conventional index notations represent the individual components.

Equation (13) implies that the rheological behavior of Bingham fluid is characterized by two different flow regimes: if |τ|<τ_y, the materials behave as a rigid solid; if |τ|≥τ_y, with the assumptions in boundary layer and the rotational symmetry, the material flows with an apparent viscosity η˜ taking the form [2], [3]

For the cylindrical coordinate, the two stress components in the plastic region are

Due to the assumptions ∂/∂Z≫∂/∂Φ and ∂/∂Z≫∂/∂R, the stress components (τ_φz and τ_rz) become

Substituting (21) into (24) and (25) gives

They are the components of stress required to the momentum equations given by (7) and (8). An important parameter of Bingham fluid is the “Bingham number” B_y, which is the ratio of the yield stress and the viscous stress. It is used to describe the viscoplastic character of the flow and is defined as

which is expressed by the following relation [28]

where ν=μ_p/ρ is the kinematic plastic viscosity of the fluid, L is a characteristic length scale, and R indicates that this is a local Bingham number.

The following dimensionless variables are introduced

Equations (6)–(9) then become

and the boundary conditions (10)–(11) are converted into

where Pr=ν/α is Prandtl number, Rev=ΩR02/ν is Reynolds number in the velocity direction, and β=a/R₀ is the variable thickness parameter.

3 Nonlinear Boundary Value Problems

Through the following similarity transformation, the generalized von Kármán transformation

where η is similarity variable, one can obtain the semi-similarity equations [29], [30] as follows

The corresponding boundary conditions become

where the primes denote the differentiation with respect to the similarity variable η.

Employing the following translation transform

Equations (38)–(41) yield

and the boundary conditions (42)–(43) are

where the prime denotes the differentiation with respect to the similarity variable ξ.

Now we consider the skin friction coefficient and the local Nusselt number. The radial and transversal stress components τ_rz and τ_φz, respectively, at the surface of the disk are given by

The coefficient of skin friction is expressed as

The Fouriers law is applied to compute the rate of heat transfer

The local Nusselt number Nu_r is obtained as

It is worthwhile to note that in Bingham plastic fluid the viscosity varies during the flow. Therefore, an effective viscosity expressed as η˜=μp+τy/|γ˙| might be more representative of the viscous stress within the flow than the constant plastic viscosity μ_p. Therefore, the Reynolds number and Bingham number can be defined more appropriately by using η˜ instead of μ_p. The following modified definitions of these parameters are [31] Rev∗=Rev/(1+By), By∗=By/(1+By). Correspondingly, the skin friction coefficient Cfr∗ can be more effectively defined as Cfr∗=emr.[(f′(0))2+(g′(0))2]12/(rRev∗)+2By∗/(rRev∗).

4 Homotopy Analysis Method

In this section, the nonlinear values (45)–(48) and the boundary conditions (49)–(50) are solved by using the HAM. The main procedures of solving the boundary value problems are presented as follows: The velocity and temperature distributions f(ξ), g(ξ), h(ξ), and s(ξ) can be expressed by the set of base functions

Based on the rule of solution expressions (56) and the boundary conditions (49)–(50), the initial guesses f₀, g₀, h₀, and s₀ of f(ξ), g(ξ), h(ξ), and s(ξ) are selected as follows

and the auxiliary linear operators are defined as

satisfying the following properties

where C_i(i=1, 2, …, 7) are the arbitrary constants. Upon making use of the above definitions, we construct the zero-order deformation problems

with the boundary conditions

Based on (45)–(48), the nonlinear operators N_f, N_g, N_h, and N_s are given by

where q∈[0, 1] is an embedding parameter. ℏ_f, ℏ_g, ℏ_h, and ℏ_s are the non-zero auxiliary parameters. H˜f(ξ),H˜g(ξ),H˜h(ξ) , and H˜s(ξ) are the auxiliary functions. As q increases from 0 to 1, F(ξ; q), G(ξ; q), H(ξ; q), and S(ξ; q) vary from f₀(ξ), g₀(ξ), h₀(ξ), and s₀(ξ) to f(ξ), g(ξ), h(ξ), and s(ξ), respectively. Using Taylor’s theorem, we then write F(ξ; q), G(ξ; q), H(ξ; q), and S(ξ; q) into the Taylor’s series with respect to q=0 as

If ℏ_f, ℏ_g, ℏ_h, and ℏ_s are chosen properly, the series (73)–(76) are convergent at q=1, and therefore, we obtain the solution series

Furthermore, the high-order deformation equations are written as

subject to boundary conditions

where

Under considerations of the rules of solution expression and existence [21], we choose the auxiliary function as

5 Results and Discussions

5.1 Convergence of the Series Solutions

The nonlinear ordinary differential equations (45)–(48) under the boundary conditions (49)–(50) are solved by using HAM with the symbolic computation software Mathematica. During the computation process, we set ℏ_f=ℏ_g=ℏ_h=ℏ_s=ℏ and calculate the fifth-order approximation solutions. The proper values of ℏ , determined by the ℏ curves and the averaged residual error curves, are chosen to guarantee the convergence of the series solution and the minimized error. Taking the set of parameters r=1/2, m=1, β=1/2, Pr=5, Re_v=10, and B_y=1 as an example, the ℏ curves and the averaged residual error curves are plotted in Figures 2 and 3.

Figure 2:

The ℏ curves of f′(0), g′(0), h′(0), and s′(0) for the fifth-order approximation solutions.

Figure 3:

The averaged residual error curves of f(ξ), g(ξ), h(ξ), and s(ξ) for the fifth-order approximation solutions.

5.2 Analysis of the Shape Parameter m and the Thickness Parameter β

We first consider the shape parameter m which dominates the radial direction decreasing way of the thickness of the disk from the center to the edge. Figures 4–6 illustrate the effects of the shape parameter m on the dimensionless radial f(ξ), tangential g(ξ), and axial velocities h(ξ). With the increasing m, the value of the peak velocity and the thickness of the radial velocity boundary layer both decrease. Especially the peak of the radial velocity moves away from the disk when the shape parameter m is increased. The influence of m on the component of tangential velocity is depicted in Figure 5. One can observe that the tangential velocity g(ξ) decreases uniformly in the boundary layer for each of the shape parameter m. The tangential velocity and associated boundary layer thickness increase when m becomes larger. Figure 6 presents the variations of the axial velocity h(ξ) for different m. It is shown that the axial velocity −h(ξ) and the boundary layer thickness both decrease for the increasing m. From the physical point of view, this is because the decrease in the velocity in the radial direction slows down the vertically downward flow which compensates for the outflow fluid. The influence of the shape parameter m on the dimensionless temperature distribution s(ξ) is presented in Figure 7. The results show that in the thermal boundary layer the temperature increases at each position when m is increased. The thickness of the thermal boundary layer also increases with the increasing m which implies that the thermal resistance is strengthened.

Figure 4:

Changes of radial velocity profile f(ξ) for different values of m.

Figure 5:

Changes of tangential velocity profile g(ξ) for different values of m.

Figure 6:

Changes of axial velocity profile −h(ξ) for different values of m.

Figure 7:

Changes of temperature profile s(ξ) for different values of m.

The analytical solutions show that the thickness parameter β has very little effect on the velocity components of the radial f(ξ), the tangential g(ξ), and the temperature distribution s(ξ) but not the same situation as the axial velocity h(ξ). Figure 8 illustrates the effects of the thickness parameter β on the dimensionless axial velocity h(ξ). The axial velocity −h(ξ) increases uniformly in the boundary layer for each β. It is further seen that the axial velocity −h(ξ) increases when the thickness parameter β is increased, and consequently, the momentum boundary layer thickness becomes thinner.

Figure 8:

Changes of axial velocity profile −h(ξ) for different values of β.

5.3 Analysis of the Reynolds Number Re_v and the Prandtl Number Pr

Figures 9–11 depict the effects of the Reynolds number Re_v on the dimensionless velocity of the radial f(ξ), tangential g(ξ), and axial h(ξ). In Figure 9, the values of the maximum velocity increase with the increasing Re_v, and its location moves slightly towards the disk. The variations of the tangential velocity for different Reynolds number Re_v are shown in Figure 10. The results show that the tangential velocity g(ξ) decreases uniformly in the boundary layer for each Reynolds number Re_v. Figure 11 presents the effects of the Reynolds number Re_v for the axial velocity −h(ξ). It is observed that the axial velocity −h(ξ) increases with the increasing Re_v, and the boundary layer thickness decreases correspondingly. Thus, it can be deduced that the faster the speed of the rotating disk, the faster the motion of the fluid moving towards the disk.

Figure 9:

Changes of radial profile f(ξ) for different values of Re_v.

Figure 10:

Changes of tangential profile g(ξ) for different values of Re_v.

Figure 11:

Changes of axial profile −h(ξ) for different values of Re_v.

Figure 12 shows the effects of Re_v on the dimensionless temperature distribution s(ξ). The effects of Re_v on the temperature fields are much smaller compared with the effects on velocity distribution. When investigated precisely, it is found that the temperature s(ξ) decreases with the increasing Re_v. Therefore, the thickness of the temperature boundary layer becomes thinner and the ability of heat transfer is enhanced. The effects of Prandtl number Pr on the temperature distribution are demonstrated in Figure 13. For the temperature field, it is showed that the temperature distribution decreases with the increasing of the Prandtl number Pr. Accordingly, the thickness of the thermal boundary layer becomes thinner. Hence, the thermal resistance is weakened for the thermal boundary layer and the efficiency of heat transport is improved.

Figure 12:

Changes of temperature profile s(ξ) for different values of Re_v.

Figure 13:

Changes of temperature profile s(ξ) for different values of Pr.

5.4 Analysis of the Bingham Number B_y

Figures 14–17 show the differences of the velocity and temperature distributions between the Reynolds number and the modified Reynolds number for different Bingham number. The Bingham number B_y is the characteristic material parameter to describe and measure Bingham fluid. We can observe that the changes of the velocity and temperature profiles under the modified Reynolds number are kept consistent with the conventional Reynolds number. In Figure 14, the results show that with the increase of the Bingham number B_y the values of the maximum velocity decrease and the momentum boundary layer thickness becomes thinner.

Figure 14:

Comparisons between radial velocity profiles f(ξ) of Bingham fluid with the Reynolds number and the modified Reynolds number.

Figure 15:

Comparisons between tangential velocity profiles g(ξ) of Bingham fluid with the Reynolds number and the modified Reynolds number.

Figure 16:

Comparisons between axial velocity profiles −h(ξ) of Bingham fluid with the Reynolds number and the modified Reynolds number.

Figure 17:

Comparisons between temperature profiles s(ξ) of Bingham fluid with the Reynolds number and the modified Reynolds number.

Figure 15 illustrates the variations of the tangential velocity for different Bingham number B_y. With the increase of the Bingham number B_y, the boundary layer thickness of tangential velocity becomes thicker. Analogous to the previously analysis, the reduction of the radial speed results in the reduction of the axial speed of the fluid because the radial outflow needs to carry away the incoming axial flow. Therefore, the axial velocity −h(ξ) and the boundary layer thickness both decrease with the increasing B_y as presented in Figure 16. Figure 17 shows the effects of the Bingham number B_y on the dimensionless temperature distribution s(ξ). The effects of the B_y on the temperature fields are relatively smaller compared with the effects on velocity distribution. The temperature boundary layer becomes thicker, and consequently the thermal resistance is strengthened in the thermal boundary layer.

5.5 Analysis of the Skin Friction Coefficient C_fr and Nu_r

The important physical quantities are the skin friction coefficient and the local Nusselt number. Figure 18a and b illustrates the variations of the skin friction coefficient C_fr and the local Nusselt number Nu_r with the shape parameter m. One can see that both of these two parameters increase as the value of m increases. The growth of the skin friction coefficient C_fr is nonlinear with respect to the shape parameter m, but the growth of the local Nu_r is almost linear with m. The increasing rate of Nu_r at position r=1/2 is larger than the case of r=2/5 for each m. It can be deduced that the heat transfer efficiency becomes better gradually when m increases.

Figure 18:

The skin friction coefficient and local Nusselt number at different positions of r for variation of m with β=1/2, Pr=5, Re_v=10, and B_y=1.

Figure 19 displays the effects of the modified Reynolds number and Bingham number and shows the results of the modified skin friction coefficient. Figure 19a shows the modified Reynolds number results in the increase of the skin friction coefficient, while the modified Bingham number reduces the skin friction coefficient as showed in Figure 19b. For the skin friction coefficient showed in Figure 19c, the modified skin friction coefficient is less than the original skin friction coefficient, but the growth rate of the modified skin friction coefficient is greater than the original skin friction coefficient for the shape parameter m.

Figure 19:

Comparisons for different definitions of the Reynolds number, Bingham number, and skin friction coefficient with m at β=1/2, Pr=5, Re_v=10, and B_y=1.

6 Conclusions

This paper investigates the flow and heat transfer characteristics of Bingham plastic fluid over a rotating disk with variable thickness radially. The HAM is exploited and analyzed thoroughly. Some important findings of the paper are the following:

1. The radial velocity depends on the shape parameter m and the Bingham number B_y strongly. The maximal value of the radial velocity moves away from the disk with the increasing m and backward with the increasing B_y.

2. The shape parameter m, the Reynolds number Re_v, and the Bingham number B_y have obvious influences on the tangential velocity, but it is not the same case for the thickness parameter β. When the shape parameter m or the Bingham number B_y increases, the thickness of the boundary layer and the absolute value of the tangential velocity increase. However, with the increase of the Reynolds number Re_v, the associated boundary layer decreases.

3. The variations of radial velocity exhibit acceleration or opposition to the motion of axial velocity. When the shape parameter m or the Bingham number B_y increases, the axial velocity h(ξ) decreases and obtains the external field velocity rapidly.

4. The thermal boundary layer increases when the shape parameter m or the Bingham number B_y is increased. The increase of Prandtl number Pr decreases the thickness of the thermal boundary layer, and hence the thermal resistance is weakened.

5. The variations of the skin friction coefficient C_fr are nonlinearly dependent on the shape parameter m, but the growth of the local Nu_r exhibits an almost linearly dependent relation with m. The efficiency of heat transfer improves when m is increased.

Nomenclature

R

R axis

Φ

Φ axis

Z

Z axis

U

Velocity in R axis direction

V

Velocity in Φ axis direction

W

Velocity in Z axis direction

T

Temperature of fluid

ν

Kinematic viscosity of fluid

α

Thermal diffusivity of fluid

P

Pressure

R₀

Feature radius

Ω

Angular velocity of rotating disk

T_w

Surface temperature of disk

T_∞

Temperature of audient fluid

r

Non-dimensional variable of R

z

Non-dimensional variable of Z

u

Non-dimensional variable of U

v

Non-dimensional variable of V

w

Non-dimensional variable of W

θ

Non-dimensional variable of T

Rev∗

Modified Reynolds number

By∗

Modified Bingham number

Cfr∗

Modified coefficient of skin friction

q_w

Heat flux

τ_rz

Radial stress tensor

τ_φz

Tangential stress tensor

ρ

Density of fluid

μ_p

Viscosity of fluid

e_ij

Rate-of-deformation tensor

ϕ

Tangential coordinate

L

Characteristic length scale

τ

Extra stress

τ_y

Yield stress

τ_ij

Deviatoric stress tensor

γ˙

Shear rate

m

Shape parameter of disk

β

Variable thickness parameter of disk

η˜

Apparent viscosity

Pr

Prandtl number

Nu_r

Local Nusselt number

λ

Thermal conductivity

c_p

Specific heat capacity

Re_v

Reynolds number

B_y

Bingham number

C_fr

Coefficient of skin friction