A new time and spatial fractional heat conduction model for Maxwell nanofluid in porous medium

Computers and Mathematics with Applications 78 (2019) 1621–1636

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

A new time and spatial fractional heat conduction model for Maxwell nanofluid in porous medium ∗

Mengchen Zhang a , Ming Shen a , , Fawang Liu a,b , Hongmei Zhang a a b

College of Mathematics and Computer Science, Fuzhou University, Fuzhou, China School of Mathematical Sciences, Queensland University of Technology, Australia

article

info

Article history: Available online 25 February 2019 Keywords: Maxwell viscoelastic nanofluid Time and spatial fractional derivatives Brownian diffusion and thermophoresis Heat and mass transfer

a b s t r a c t A new time and spatial fractional heat conduction model with Brownian diffusion and thermophoresis is presented to investigate the heat and mass transfer of Maxwell viscoelastic nanofluid over a moving flat plate in porous medium. New dimensionless variables are introduced to nondimensionalize the governing equations which are solved numerically by L1 algorithm and shifted Grünwald formula. The main observations of the model are the effects of embedded time and space fractional parameters on velocity, temperature and concentration profiles which are analyzed graphically. Numerical computations for a comparison between the heat conduction model with the time and spatial fractional derivatives and the model with spatial fractional derivative are made. It is found that the heat transfer is enhanced by the time and spatial fractional heat conduction model. Moreover, the larger fractional derivatives refer to the stronger memory characteristic which exhibits the physical meanings of fractional derivative. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Fractional calculus [1] is ideal for characterizing materials and physical processes with memory and hereditary properties via fractional differential equations [2]. Compared with integer-order dynamic systems, the fractional-order differential equations can model complicated systems easily and define the physical meaning of some parameters exactly. It has attracted growing attention in recent decades as the application of the fractional calculus has penetrated into various fields of engineering technology and industrial production, especially in the study of rheology [3], mechanical systems [4], control systems [5], signal processing [6], anomalous diffusion [7,8] and viscoelastic materials [9]. Sun et al. concluded various nonlocal phenomena described by fractional calculus which were explored by distinguished researchers [10]. The investigations of fractional boundary layer flow and heat transfer have attracted tremendous attention of researchers [11–13]. Fractional derivative has been also found quite flexible in the description of viscoelastic constitutive relation. The non-local nature of fractional derivative operators can describe the long-term memory character of viscoelastic materials exactly [14–16]. The constitutive relationships for the viscoelastic fluids are modified from the Newtonian fluid models by replacing the time derivative of an integer order by fractional calculus operators. As one of the important viscoelastic fluids, Maxwell fluid has been investigated extensively by many researchers. Tripathi presented the flow of viscoelastic fluid with fractional Maxwell model through a channel [17]. Zhao et al. introduced the fractional Marangoni boundary condition to explore the unsteady heat transfer of Maxwell fluid with Cattaneo heat flux [18]. The flow of fluid through porous medium is a significant branch of various engineering and disciplines. From the perspective of disciplinary development, ∗ Corresponding author. E-mail address: [email protected] (M. Shen). https://doi.org/10.1016/j.camwa.2019.01.006 0898-1221/© 2019 Elsevier Ltd. All rights reserved.

1622

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

Nomenclature Da DB DT hp jp K k Le Nb Nt P Pr q Re r S t T u, v V x, y

Darcy number Brownian diffusion coefficient Thermophoresis diffusion Specific enthalpy Diffusion mass flux Permeability Thermal conductivity Lewis number Brownian motion parameter Thermophoresis parameter Pressure Generalized Prandtl number Energy flux Reynold numbers Darcy resistance Extra stress tensor (dimensionless) Time Maxwell nanofluid temperature (dimensionless) Velocity components Velocity vector (dimensionless) Coordinate components

Greek Symbols

Γ α, β αf γ δ ε θ λ µ νf ρf ρp (ρ c)f σ τ0 φ ϕ γ +1 ωl

Gamma function Time fractional derivative parameters Thermal diffusivity Spatial fractional derivative Weight coefficient Porosity Dimensionless temperature Relaxation time Dynamic viscosity Kinematic viscosity Density of the fluid Mass density of the nanoparticles Heat capacity of the fluid The coefficient to keep the dimension of constitutive equation balance Temperature relaxation time Concentration Dimensionless concentration Grünwald weight coefficient

Subscript f p

w ∞

Fluid Nanoparticles Wall condition Ambient condition

Superscript

∗

Dimensionless form

the heat and mass transfer of porous medium has penetrated into many disciplines and new technologies, including energy, biotechnology, environmental science, chemical engineering, medicine and agricultural engineering [19]. Recently,

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

1623

researchers have paid considerable attention to the study of flow of viscoelastic fluids in the porous medium. Tripathi and Bég [20] studied the oscillating peristaltic flow of generalized Maxwell fluids through a porous medium and employed a sinusoidal model to simulate the oscillating flow by solving the problem analytically and numerically. Nayak et al. [21] investigated the heat and mass transfer of an electrically conducting viscoelastic fluid through porous medium in the presence of chemical reaction. The thermal instability in variable gravitational field of a rotating Maxwell viscoelastic fluid in porous medium was investigated [22]. In recent years, the investigation of nanofluid has been a hot research field due to the fact that nanofluids exhibit superior heat transfer compared with the conventional heat transfer fluids [23]. Studies show that the thermal conductivity of nanofluids cannot be predicted by traditional heat transport correlations as the suspended particles remarkably increase thermal conductivity of nanofluids [24]. Due to many uncertain factors, there is still no unified framework to explain this abnormal property of nanofluids. Buongiorno [25] proposed seven factors and affirmed that only Brownian diffusion and thermophoresis are dominant slip mechanisms in nanofluids. In his model, the heat flux of nanofluid q = −k∇ T + hp jp , in which jp was written as the sum of Brownian diffusion and thermophoresis. It is significant and meaningful to apply the fractional calculus theory to describe anomalous heat transport due to the global dependency of fractional operators. Povstenko [26] discussed the space–time fractional heat conduction equation, and then he formulated the corresponding theory of thermal stresses by the generalized Cattaneo-type equations with Caputo’s time fractional derivatives [27]. Recently, Shen et al. [28] presented a renovated Buongiorno’s model with Caputo time fractional derivatives to investigate the heat and mass transfer characteristics of Sisko nanofluid over a continuously moving flat plate. In the model, the fractional Cattaneo heat flux was introduced to describe the anomalous heat transport of Sisko nanofluid considering the influences τ

β

∂β q

of Brownian diffusion and thermophoresis, and the heat flux q was written as q + β!0 ∂ t β = −k∇ T + hp jp with 0 ≤ β < 1. Studies showed that the renovated model with Caputo time fractional derivatives is more capable of explaining the abnormal thermal conductivity enhancement. On the other hand, it is worth noting that the space fractional derivative has been applied to model the non-local characteristics of the flow and heat transfer of boundary layer flow. Pan et al. [29] suggested a spatial fractional-order thermal transport equation to describe convective heat transfer of nanofluids within disordered porous media in boundary layer flow. Chen et al. [30] introduced spatial fractional derivative shear stress and Fourier’s law to study the MHD boundary layer equations of fractional viscoelastic fluid over a stretching sheet. A new generalization to the classical heat conduction problem was proposed by using the space–time Cattaneo heat conduction law [31]. Povstenko [32] considered the thermoelasticity based on space–time fractional heat conduction equation and he emphasized the significance of time and spatial fractional derivatives which can describe the non-local phenomena of thermoelasticity. Liu et al. [33] proposed a novel time and space fractional Cattaneo–Christov upper-convective derivative flux heat conduction model. But they only considered the energy equation and regarded the velocity as a constant, which means the natural conduction and velocity boundary layer are out of consideration. Therefore, the exploration of boundary layer flow with time and space fractional derivatives is still sparse, further research is needed to focus on this challenging problem. Motivated by the above discussions, we firstly investigate the flow as well as heat and mass transfer of Maxwell viscoelastic nanofluids with time and space fractional derivatives considering the influences of Brownian diffusion and thermophoresis. The corresponding governing equations are formulated and nondimensionalized by the new dimensionless variables. The time and space fractional derivatives are approximated by the L1-algorithm and the shifted Grünwald formulae [34], respectively. A comparison between the solutions of the problem with time and spatial fractional derivatives and the solutions with only spatial fractional derivative is presented. The influences of involved parameters including fractional derivatives on the velocity, temperature and concentration distributions are graphically analyzed in detail. 2. Mathematical formulation of the problem 2.1. Conservation equations for nanofluids with fractional calculus According to the force balance which acts on a volume unit of fluid, the momentum equation in a porous medium is formulated as [35]

ρf

dV dt

= −∇ p + ∇ · S + r,

(1)

where ρf is the density of the fluid, V is the velocity vector, d/dt is the material time derivative, p is the pressure, S is the extra shear stress and r is the Darcy resistance which is offered by the solid matrix. The modified fractional Darcy’s law of Maxwell model in a porous medium is given by [36] (1 + λα

∂α µε )r = − V, ∂tα K

(2)

where λ is the relaxation time, α (0 ≤ α < 1) is the fractional derivative parameter, µ is the dynamic viscosity, K and ε are the permeability and porosity of the porous medium. Eq. (2) is reduced to the Darcy’s law for the viscous fluid through a

1624

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636 ∂α ∂tα

porous medium when λ = 0. α

∂ f (t) 1 = α ∂t Γ (1 − α )

is the Caputo’s fractional derivative operator of order α which is defined as follows [1]

t

∫

(t − η)−α f ′ (η)dη,

(3)

0

where Γ (·) is the Gamma function. Substituting Eq. (1) into Eq. (2), the momentum equation can be obtained as follows (1 + λα

) µε ∂ α ( dV ) ρ + ∇ p − ∇ · S = − V. f α ∂t dt K

(4)

The energy equation for nanofluid can be written as [25] (ρ c)f

( ∂T

) + V · ∇ T = −∇ · q + hp ∇ · jp ,

∂t

(5)

where the diffusion mass flux for the nanoparticles jp can be written as the sum of Brownian diffusion and thermophoresis jp = jp,B + jp,T = −ρp DB ∇φ − ρp DT

∇T T∞

,

(6)

where T is the nanofluid temperature, (ρ c)f is the heat capacity of the fluid while ρp is the mass density of the nanoparticles, hp is the specific enthalpy of the nanoparticle material and ∇ hp has been set equal to cp ∇ T , q is the energy flux, DB and DT represent the Brownian diffusion coefficient and thermophoresis diffusion coefficient, respectively. In the renovated Buongiorno’s model [28], the Cattaneo model of heat conduction with time fractional derivative is introduced to describe the abnormal thermodynamic characteristics of Sisko nanofluid as the time fractional order derivative refers to memory character of heat transfer process. Moreover, the space fractional derivative can characterize non-local behavior which demonstrates that the state depends on not only the nearby positions but also the whole positions. Therefore, a new time and spatial fractional heat conduction model with Brownian diffusion and thermophoresis is proposed to formulate energy equation which is defined as β

q+

τ0 ∂ β q = −kσ γ −1 ∇ γ T + hp jp , (0 ≤ β < 1, 0 ≤ γ < 1) β! ∂ t β

(7)

where k is the thermal conductivity, σ is introduced to keep the dimension of constitutive equation balance and its dimension is meter, β is the time fractional derivative parameter, ∂ β /∂ t β is the Caputo’s fractional derivative of order β , γ is the spatial fractional derivative parameter, and when T = T (t , x, y, z) the ∇ γ T is defined as

∇ γ T = (δ

∂γ T ∂γ T ∂γ T ∂γ T ∂γ T ∂γ T − (1 − δ ) , δ γ − (1 − δ ) , δ γ − (1 − δ ) ) γ γ γ ∂x ∂ (−x) ∂y ∂ (−y) ∂z ∂ (−z)γ γ

with δ (0 ≤ δ ≤ 1) being the weight coefficient of forward versus backward transition probability. Here the symbols ∂∂ xγT γ and ∂ (∂−x)T γ are the left and right Riemann–Liouville fractional derivatives of order γ (n − 1 ≤ γ < n), the corresponding definitions on a finite domain [a, b] are given by [34]: 1 ∂n ∂γ T = ∂ xγ Γ (n − γ ) ∂ xn

∫

x

(x − ξ )−γ +n−1 T (t , ξ , y, z)dξ ,

(8)

a

∂γ T (−1)n ∂ n = γ ∂ (−x) Γ (n − γ ) ∂ xn

∫

b

(ξ − x)−γ +n−1 T (t , ξ , y, z)dξ ,

(9)

x

respectively. Substituting Eqs. (5) and (6) into Eq. (7), we get the fractional energy equation

(

(ρ c)f 1 +

β

β

) ) τ0 ∂ β )( ∂ T τ ∂β ( + V · ∇ T = k∇ · (σ γ −1 ∇ γ T ) − ∇ hp · jp + 0 β hp ∇ · jp , β β! ∂ t ∂t β! ∂ t

(10)

which reduces to the energy equation in the renovated Buongiorno’s model [28] when γ = 1. The detailed deducing process of Eq. (10) is given in Appendix A. The continuity equation for the nanoparticles in the absence of chemical reactions reads [25] 1 ∂φ + V · ∇φ = − ∇ · jp , ∂t ρp

(11)

where φ is the concentration of the fluid, jp is given by Eq. (6). 2.2. Problem statement Consider two-dimensional unsteady flow of Maxwell nanofluid in porous medium over a moving sheet. The X -axis is along the flat sheet while the Y -axis is perpendicular to the sheet. The sheet impulsively starts to move at time t = 0 with

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

1625

Fig. 1. The schematic diagram of the physical model.

a constant velocity. The surface of the moving sheet and ambient temperatures are assumed to be constants Tw and T∞ , α respectively. The constitutive relation of fraction Maxwell nanofluid satisfies [37] (1 + λα ∂∂t α )Sxy = µ ∂∂ uy with Sxy being the component of shear stress. It is assumed that ∂ p/∂ x = 0 as the flow along a horizontal sheet. The schematic diagram of the physical model is described in Fig. 1. Under the velocity and thermal boundary layer approximations, the governing equations of the two-dimensional unsteady Maxwell nanofluid flow with the new time and spatial fractional heat conduction read:

∂ u ∂v + = 0, ∂x ∂y ∂ 2 u νf ε ∂α ( ∂u ∂u ∂u ) (1 + λα α ) +u +v = νf 2 − u, ∂t ∂t ∂x ∂y ∂y K β τ ∂β ( ∂T ∂T ∂T ∂T ∂T ∂T ) +u +v + 0 β +u +v ∂t ∂x ∂y β! ∂ t ∂ t ∂x ∂y [ ∂γ T [ ∂φ ∂ T ∂γ T ] DT ( ∂ T )2 ] γ −1 ∂ = αf σ δ γ − (1 − δ ) + τ DB + γ ∂y ∂y ∂ (−y) ∂y ∂y T∞ ∂ y β τ0 ∂ β [ ∂ 2φ DT T ∂ 2 T ] −τ DB T + , β! ∂ t β ∂ y2 T∞ ∂ y 2 ∂φ ∂ 2φ ∂φ ∂φ DT ∂ 2 T +u +v = DB 2 + . ∂t ∂x ∂y ∂y T∞ ∂ y2

(12) (13)

(14) (15)

where νf = µ/ρf is the kinematic viscosity, K is the permeability of the medium, αf = k/(ρ c)f is the thermal diffusivity of the nanofluid, τ = (ρ c)p /(ρ c)f is the ratio of effective heat capacity of the nanoparticle material to the heat capacity of the fluid. The initial and boundary conditions for this problem are given by u(0, x, y) = 0, v (0, x, y) = 0, T (0, x, y) = T∞ , φ (0, x, y) = φ∞ ; 3

u(t , x, 0) = U0 Re− 2 , v (t , x, 0) = 0, T (t , x, 0) = Tw , φ (t , x, 0) = φw ; u(t , 0, y) = v (t , 0, y) = 0, T (t , 0, y) = T∞ , φ (t , 0, y) = φ∞ ; u(t , x, y) = v (t , x, y) = 0, T (t , x, y) = T∞ , φ (t , x, y) = φ∞ as y → ∞. Introduce the new dimensionless variables with σ in the following forms: 1

x∗ =

xRe 2

, y∗ =

y

, t∗ =

tU0

3

, u∗ =

uRe 2

, v∗ =

v Re

σ σ σ Re U0 U0 T − T∞ φ − φ∞ ∗ λU0 ∗ τ0 U0 θ= ,ϕ = ,λ = ,τ = , Tw − T∞ φw − φ∞ σ Re 0 σ Re

,

(16) (17)

in which Re = ρσ U0 /µ is the generalized Reynolds numbers. Noting that the boundary conditions of dimensionless temperature θ (t , x, 0) = 1 and θ (t , x, ∞) = 0, so we assume δ = 0 here. Omitting the dimensionless mart ∗ for brevity, the

1626

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

governing equations of continuity, momentum, energy and concentration can be written as

∂ u ∂v + = 0, ∂x ∂y ∂α ( ∂u ∂u ∂ u ) ∂ 2u ε (1 + λα α ) +u +v = 2 − u, ∂t ∂t ∂x ∂y ∂y Da β τ ∂ β ( ∂θ ∂θ ∂θ ∂θ ∂θ ) ∂θ +u +v + 0 β +u +v ∂t ∂x ∂y β! ∂ t ∂ t ∂x ∂y 1 ∂ γ +1 θ Nb ∂ϕ ∂θ Nt ( ∂θ )2 = + + Pr ∂ (−y)γ +1 Pr ∂ y ∂ y Pr ∂ y β τ0 ∂ β ( Nb ∂ 2 ϕ Nt ∂ 2 θ ) θ + θ 2 , − β! ∂ t β Pr ∂ y2 Pr ∂ y ∂ϕ ∂ϕ ∂ϕ 1 ∂ 2ϕ 1 Nt ∂ 2 θ +u +v = + , 2 ∂t ∂x ∂y PrLe ∂ y PrLe Nb ∂ y2

(18) (19)

(20) (21)

where Da is the Darcy number, Pr is the Prandtl number, Le denotes the Lewis number, Nb and Nt are the Brownian motion parameter and the thermophoresis parameter respectively, which are expressed as follows:

νf , σ α (ρ c)f τ DB (φw − φ∞ ) τ DT (Tw − T∞ ) α τ= , Nb = , Nt = , Le = , (ρ c)p α T∞ α DB Da =

K

2

, Pr =

The corresponding initial and boundary conditions become u(0, x, y) = 0, v (0, x, y) = 0, θ (0, x, y) = 0, ϕ (0, x, y) = 0;

(22)

u(t , x, 0) = 1, v (t , x, 0) = 0, θ (t , x, 0) = 1, ϕ (t , x, 0) = 1;

(23)

u(t , 0, y) = v (t , 0, y) = 0, θ (t , 0, y) = 0, ϕ (t , 0, y) = 0;

(24)

u(t , x, y) = v (t , x, y) = 0, θ (t , x, y) = 0, ϕ (t , x, y) = 0 as y → ∞.

(25)

3. Numerical technique 3.1. Discretization method Denote xi = i∆x(i = 0, 1, 2, . . . , M), yj = j∆y(j = 0, 1, 2, . . . , N) and tk = k∆t(k = 0, 1, 2, . . . , R), where ∆x = L/M and ∆y = Ymax /N are space steps, ∆t is time step. According to the definition of Caputo’s fractional derivative, the L1 algorithm [38] is introduced to discretize the time fractional derivative: k−1 ] ∂ β f (tk ) ∆t −β ∑ [ = bq f (tk−q ) − f (tk−q−1 ) + O(∆t 2−β ) β ∂t Γ (2 − β ) q=0

(26)

k−1

∑ [ ] ∆t f (tk ) − bk−1 f (t0 ) − (bq−1 − bq )f (tk−q ) + O(∆t 2−β ) = Γ (2 − β ) −β

q=1

where bq = (q + 1) − q , q = 0, 1, 2, . . . , R. Based on the shifted Grünwald formulae [34], the space fractional derivative of order 1 ≤ γ + 1 < 2 is approximated by: 1−β

1−β

N −j+1 ∑ γ +1 ∂ γ +1 θ (xi , yj , tk ) 1 ≈ ωl θ (xi , yj+l−1 , tk ), ∂ (−y)γ +1 ∆yγ +1

(27)

l=0

γ +1

where the symbol ωl

γ +1

refers to the Grünwald weight coefficient with ω0

γ +1

= 1 and ωl

( ) γ +1 = 1 − γ +l 2 ωl−1 . The detailed

difference forms and iterative equations are illustrated in Appendix B. 3.2. Solution procedure The values of u, v , θ and φ in the specific domain at t = 0 are obtained from the initial condition (22). The values of the variables of (k − 1)-level are regarded as constants during any one-time step. The iteration equations of velocity and concentration at each internal nodal point are formulated by tri-diagonal system of equations which can be solved by the

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

1627

Fig. 2. Grid independence test for different grid sizes.

Fig. 3. Velocity profiles for different values of α .

Thomas algorithm [39]. Moreover, due to the coefficient matrix of the iteration energy equation is a kind of nonsingular matrix, the solution of the energy equation is acquired by solving the inverse matrix. When the absolute values of the difference between temperature θ and velocity u at all the mesh nodes within two consecutive time steps are less than 10−4 , the iteration of time level stops and reaches to the steady state. Regarding the computational domain as a rectangle with the boundary of xmax = 1 and ymax = 10. Here ymax lies outside the porous medium and corresponds to y → ∞. Considering the accuracy of numerical solutions and the time of calculation, the spatial grid sizes and time steps are fixed as ∆x = 0.025, ∆y = 0.1 and ∆t = 0.1. The stability of the numerical solutions can be verified by the computation of dimensionless velocity u and temperature θ as depicted in Fig. 2. 4. Results and discussions The numerical solutions for governing equations are carried out with various values of involved physical parameters. In the following analyses, we focus on the effects of fractional derivatives α , β , γ , porosity ε and Darcy number Da on the dimensionless velocity u, temperature θ and concentration ϕ . It is evident that the time fractional derivative parameter α has no influence on temperature and concentration in view of Eqs. (18)–(21). So we only give the velocity profiles for different values of α depicted in Fig. 3. It is observed that the velocity profiles rise with the increment of the time fractional derivative parameter α monotonically. Meanwhile, the momentum boundary layer is thicker for larger α which indicates that viscoelasticity strengthens the flow resistance with the increase

1628

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

Fig. 4. Temperature profiles for different values of β .

Fig. 5. Concentration profiles for different values of β .

of the velocity fractional derivative parameter. The temperature and concentration distributions for various values of β are described in Figs. 4 and 5. It is shown in Fig. 4 that the larger fractional derivative parameter β leads to the higher temperature. The thermal boundary layer thickens as β rises. These results indicate that the heat transfer will reduce with increasing β which can be attributed to the heat transfer resistance. It is observed from Fig. 5 that the peak value of the concentration profiles declines when the parameter β increases. But the values of concentration increase with β when y > 2 which implies that the concentration profiles intersect with each other. The appearance of the intersection point suggests that the time fractional derivative exhibits a short memory of previous states. Porosity is a fraction of the volume of voids over the total volume. It can characterize the flow complexity of porous medium. It is shown in Fig. 6 that the velocity profiles decrease and the momentum boundary layer becomes thinner for larger values of ε . On the contrary, the larger ε is, the temperature and concentration rise as illustrated in Figs. 7 and 8. The higher values of temperature and concentration result from the fact that the capability of heat and mass transfer declines with the increase of porosity. These results indicate that the porosity promotes the convection flow but weakens the efficiency of heat and mass transfer for this problem. Darcy number is proportional to the permeability K which is a significant parameter in the theory of the porous medium. As is depicted in Fig. 9, the velocity and momentum boundary layer rise remarkably with the increment of Darcy number, which indicates that the increase of permeability can reduce Darcian resistance. Unlike the tendency of velocity profiles, the temperature and thermal boundary layer decrease slightly for the larger Darcy number as described in Fig. 10. Also, the concentration distributions decline dramatically as the values of

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

1629

Fig. 6. Velocity profiles for different values of ε .

Fig. 7. Temperature profiles for different values of ε .

Da go up. These observations indicate that the better permeability impedes the convection flow but enhances the efficiency of heat and mass transfer (see Fig. 11). A comparison between the heat conduction model with time and spacial fractional derivatives and the model with only spacial fractional derivative is presented in Fig. 12. When the relaxation parameter τ0 = 0, Eq. (20) reduces to the space fractional heat flux model. It is easy to observe that the temperature distributions rise with the increment of fractional derivative γ for both models. Furthermore, the thermal boundary layer of the new time and spatial fractional heat conduction model is thinner which suggests that this model leads to less heat conduction loss. The temperature of the space fractional model transports slower than another model which implies that the time and space fractional model leads to a higher efficiency of heat transfer. This result demonstrates the efficiency of heat transfer for the new heat conduction model with time and spacial fractional derivatives is better than that of the space fractional heat flux model. Brownian diffusion and thermophoresis are dominant slip mechanisms in nanofluids which can affect the efficiency of heat and mass transfer. The effects of Brownian diffusion and thermophoresis on temperature and concentration are illustrated in Figs. 13–16. It is observed from Figs. 13 and 14 that the temperature and concentration profiles decrease as Nb ascends. It is further found from Fig. 14 that the peak value of concentration profiles and concentration boundary layer decline significantly for larger Nb. These results indicate that the Brownian diffusion can enhance the efficiency of heat and mass transfer. The influences of thermophoresis on the temperature and concentration profiles are presented in Figs. 15– 16. The temperature profiles and thermal boundary layer increase slightly with Nt as depicted in Fig. 15. It is obvious that the thermophoresis can influence the concentration distributions considerably as described in Fig. 16. The peak value rises

1630

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

Fig. 8. Concentration profiles for different values of ε .

Fig. 9. Velocity profiles for different values of Da.

remarkably and the concentration boundary layer is thicker for larger Nt, which suggests that the thermophoresis plays an important role in the diffusion of nanoparticles. 5. Conclusions This paper deals with a study on two-dimensional unsteady boundary layer flow as well as heat and mass transfer of Maxwell viscoelastic nanofluids in a porous medium over a moving sheet. With the consideration of Brownian diffusion and thermophoresis, the new modified energy equation is obtained by replacing the time and space derivatives by the Caputo’s fractional derivative and the left and right Riemann–Liouville fractional derivatives, respectively. Based on the L1-algorithm and the shifted Grünwald formula, the governing equations are solved numerically. The influences of embedded fractional derivative parameters α , β , γ , the porosity, Darcy number, Brownian diffusion and thermophoresis on the velocity, temperature and concentration distributions are discussed by graphical illustrations in this paper. The porosity can enhance the convection flow but weaken the efficiency of heat and mass transfer, which is opposite to the effects of Darcy number on the velocity, temperature and concentration. Besides, the intersection point of the concentration profiles refers to the memory character of the heat conduction process. The thermal boundary layer thickness of the new time and spatial fractional heat conduction model is thinner which suggests that this model leads to less heat conduction loss. The Brownian diffusion can promote the efficiency of heat and mass transfer, which is opposite to the influence of thermophoresis. This study can provide a basis for further research on the application of fractional calculus in the field of viscoelastic nanofluid.

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

1631

Fig. 10. Temperature profiles for different values of Da.

Fig. 11. Concentration profiles for different values of Da.

Acknowledgments This research is supported by the Natural Science Foundation of Fujian Province, China (Grant No. 2017J01555). Author Liu wishes to acknowledge that this research was partially supported by Natural Science Foundation of China (Grant No. 11772046). Authors Liu and Turner also wish to acknowledge that this research was partially supported by the Australian Research Council (ARC), Australia via the Discovery Project (DP180103858). Appendix A The deducing process of the fractional energy equation is illustrated in detail. Eq. (5) can be rewritten as

∇ · q = −(ρ c)f

( ∂T ∂t

) + V · ∇ T + hp ∇ · jp ,

(28)

and Eq. (7) is rewritten in the following form

(

β

1+

τ0 ∂ β ) ∇ · q = −k∇ · (σ γ −1 ∇ γ T ) + ∇ · (hp jp ). β! ∂ t β

(29)

1632

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

Fig. 12. Temperature profiles for different values of γ .

Fig. 13. Temperature profiles for different values of Nb.

Substituting Eq. (28) into Eq. (29), we can obtain β β ( ) ( ) τ ∂ β )( ∂ T τ ∂ β )( −(ρ c)f 1 + 0 β + V · ∇ T + 1 + 0 β hp ∇ · j p β! ∂ t ∂t β! ∂ t = −k∇ · (σ γ −1 ∇ γ T ) + ∇ hp · jp + hp ∇ · jp ,

(30)

then Eq. (30) can be deduced to

(

(ρ c)f 1 +

β

β

) τ ∂β ( ) τ0 ∂ β )( ∂ T 0 + V · ∇ T − h ∇ · j = k∇ · (σ γ −1 ∇ γ T ) − ∇ hp · jp . p p β! ∂ t β ∂t β! ∂ t β

(31)

Finally, the fractional energy equation is obtained as follows

(

(ρ c)f 1 +

β

β

) ) τ ∂β ( τ0 ∂ β )( ∂ T + V · ∇ T = k∇ · (σ γ −1 ∇ γ T ) − ∇ hp · jp + 0 β hp ∇ · jp . β β! ∂ t ∂t β! ∂ t

(32)

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

1633

Fig. 14. Concentration profiles for different values of Nb.

Fig. 15. Temperature profiles for different values of Nt.

Appendix B In order to ensure diagonal dominance and avoid quadratic equation, the integer-order terms in the governing equations by linearization processing are discretized as

∂θ θ (xi , yj , tk ) − θ (xi , yj , tk−1 ) |t =tk = + O(∆t), ∂t ∆t ∂θ θ (xi , yj , tk ) − θ (xi−1 , yj , tk ) u |t =tk = u(xi , yj , tk−1 ) + O(∆t + ∆x), ∂x ∆x θ (xi , yj , tk ) − θ (xi , yj−1 , tk ) ∂θ v |t =tk = v (xi , yj , tk−1 ) + O(∆t + ∆y), ∂y ∆y ∂ 2θ θ (xi , yj+1 , tk ) − 2θ (xi , yj , tk ) + θ (xi , yj−1 , tk ) |t =tk = + O(∆y2 ), ∂ y2 ∆y2 ∂ϕ ∂θ ϕ (xi , yj+1 , tk−1 ) − ϕ (xi , yj , tk−1 ) θ (xi , yj+1 , tk ) − θ (xi , yj , tk ) |t =tk = + O(∆t + ∆y), ∂y ∂y ∆y ∆y

(33) (34) (35) (36) (37)

1634

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

Fig. 16. Concentration profiles for different values of Nt.

( ∂θ )2

θ (xi , yj+1 , tk−1 ) − θ (xi , yj , tk−1 ) θ (xi , yj+1 , tk ) − θ (xi , yj , tk ) |t =tk = + O(∆t + ∆y), ∂y ∆y ∆y ϕ (xi , yj+1 , tk ) − 2ϕ (xi , yj , tk ) + ϕ (xi , yj−1 , tk ) ∂ 2ϕ + O(∆t + ∆y2 ), θ 2 |t =tk = θ (xi , yj , tk−1 ) ∂y ∆y2 ∂ 2θ θ (xi , yj+1 , tk ) − 2θ (xi , yj , tk ) + θ (xi , yj−1 , tk ) θ 2 |t =tk = θ (xi , yj , tk−1 ) + O(∆t + ∆y2 ). ∂y ∆y 2

(38) (39) (40)

The fractional derivatives become:

∂ β+1 θik,j ∂ t β+1

k−1

=

)( )] ∑( ∆t −1−β [ k k−q k−q−1 θi,j − θik,j−1 − bq−1 − bq θi,j − θi,j , Γ (2 − β )

(41)

q=1

k−1

[ ( ) ∑( ) ( )] ∆t −β ∂ β ( ∂θ ) k−q−1 k−q k−q k−1 k k u = u θ − θ − b − b u θ − θ , q − 1 q i , j i − 1 , j i , j i , j i , j i − 1 , j ∂tβ ∂x ∆xΓ (2 − β )

(42)

q=1

where the truncation error is O(∆t 2−β + ∆x). k−1 ( [ ( ) ∑ ) ( )] ∆t −β ∂ β ( ∂θ ) k−q−1 k−q k−q k−1 k k v = v θ − θ − b − b v θ − θ , q − 1 q i , j i , j − 1 i,j i ,j i,j−1 ∂tβ ∂y ∆yΓ (2 − β ) i,j q=1

where the truncation error is O(∆t 2−β + ∆y). For simplicity, we note that β

τ0 ∆t −β λα ∆t −α ∆t ∆t , m2 = , m3 = , m4 = , 2 Γ (2 − α ) ∆y β!Γ (2 − β ) ∆yγ +1 k−1 ( k−1 ( )( ) ) ( ) ∑ ∑ A1 = αs−1 − αs uki,−j s − uki,−j s−1 , A2 = αs−1 − αs uki,−j s−1 uki,−j s − uki−−1s,j ,

m1 =

s=1

A3 =

s=1

k−1 ( ∑

)

(

)

αs−1 − αs vik,−j s−1 uki,−j s − uki,−j−s1 ,

C1 =

s=1

C2 =

)

( ) k−q−1 bq−1 − bq ui,j θik,j−q − θik−−1q,j ,

C3 =

q=1

E1 =

q=1

bq−1 − bq

)(

) θik,j−q − θik,j−q−1 ,

)

k−q−1

q=1

k−1 ( ∑

k−1 ( ∑

k−1 ( ∑

k−1 ( ∑

bq−1 − bq vi,j

(

) θik,j−q − θik,j−−q1 ,

q=1

) k−q−1 2 k−q bq−1 − bq θi,j δy ϕi,j ,

E2 =

k−1 ( ∑

)

k−q−1 2 k−q y i,j

bq−1 − bq θi,j

δ θ

,

q=1 k−q

where αs = (s + 1)1−α − s1−α (s = 0, 1, 2, . . . , R), δy2 ϕi,j

q k−q q 2 k−q = ϕik,− − ϕik,− = θik,j−+q1 + 2θik,j−q − θik,j−−q1 . j+1 + 2ϕi,j j−1 , δy θi,j

(43)

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

1635

Finally, the iteration equations are achieved as follows:

∆y k vik,j = vik,j−1 − (u − uki−1,j ), ∆x i,j [ ] ∆t ∆t − (1 + m1 ) uik,−j 1 uki−1,j − (1 + m1 ) vik,−j 1 + m2 uki,j−1 ∆x ∆y [ ( ∆t k−1 ) ∆t ε ∆t ] k + (1 + m1 ) 1 + ui,j + (1 + m1 ) vik,−j 1 + 2m2 + u − m2 uki,j+1 ∆x ∆y 2Da i,j ( ( ε ∆t ) k−1 ∆t ∆t ) = 1 + m1 − ui,j + m1 A1 + A2 + A3 , 2Da ∆x ∆y

(44)

(45)

[ ∆t ∆t Nt m4 γ +1 ] k −(1 + m3 ) uik,−j 1 θik−1,j − (1 + m3 ) vik,−j 1 − m2 m3 θik,j−1 + ω θi,j−1 ∆x ∆y Pr Pr 0 [ ( ) ∆t k−1 ∆t m4 γ +1 Nb 1 k−1 + (1 + m3 ) 1 + u + (1 + m3 ) vik,−j 1 − ω + m2 (ϕik,− j+1 − ϕi,j ) ∆x i,j ∆y Pr 1 Pr [m ] Nt Nt Nb 4 γ +1 1 k−1 + m2 (θik,j−+11 − θik,j−1 ) − 2m2 m3 θik,j−1 θik,j − ω2 + m2 (ϕik,− j+1 − ϕi,j ) + m2

Pr Nt Pr

(θik,j−+11 − θik,j−1 ) − m2 m3 Nb

[

Pr Nt Pr

Pr

Pr

N ] m4 ∑ γ +1 k θik,j−1 θik,j+1 − ωl−j+1 θi,l

Pr

]

l=j+2

= 1 + m3 − m2 m3 (ϕ − 2ϕ + ϕ θ Pr ( ∆t ∆t ) m2 m3 + m3 C1 + C2 + C3 + (NbE1 + NtE2 ), ∆x ∆y Pr ( ∆t k−1 1 ∆t ) k ∆t k−1 k − ui,j ϕi−1,j − vi,j + ϕ ∆x ∆y PrLe ∆y2 i,j−1 ( ) ∆t k−1 ∆t k−1 2m2 k m2 k + 1+ u v ϕ − ϕ + + ∆ x i ,j ∆y i,j PrLe i,j PrLe i,j+1 1 = ϕik,− + j

m2 Nt

PrLeNb

k i,j+1

k i,j

(θik,j+1 − 2θik,j + θik,j−1 ).

k i,j−1 )

k−1 i,j

(46)

(47)

References [1] I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. [2] S.G. Samko, A.A. Kilbas, O.I. Marichev, Fractional integrals and derivatives : theory and applications, Minsk Nauka Tekhn. (3) (1993) 397–414. [3] M.E. Reyes-Melo, V.A. González-González, C.A. Guerrero-Salazar, F. García-Cavazos, U. Ortiz-Méndez, Application of fractional calculus to the modeling of the complex rheological behavior of polymers: From the glass transition to flow behavior. i. the theoretical model, J. Appl. Polym. Sci. 108 (2) (2010) 731–737, http://dx.doi.org/10.1002/app.27435. [4] W. Grzesikiewicz, A. Wakulicz, A. Zbiciak, Non-linear problems of fractional calculus in modeling of mechanical systems, Int. J. Mech. Sci. 70 (5) (2013) 90–98, http://dx.doi.org/10.1016/j.ijmecsci.2013.02.007. [5] Y.Q. Chen, C. Ionescu, Special issue: applied fractional calculus in modelling, analysis and design of control systems, Internat. J. Control 90 (6) (2017) 1155–1156, http://dx.doi.org/10.1080/00207179.2017.1315242. [6] E.A. Gonzalez, I. Petráš, Advances in fractional calculus: Control and signal processing applications, in: Carpathian Control Conference, IEEE, 2015, pp. 147–152, http://dx.doi.org/10.1109/CarpathianCC.2015.7145064. [7] Y. Luchko, A. Nepomnyashchy, V. Volpert, A new fractional calculus model for the two-dimensional anomalous diffusion and its analysis, Math. Model. Nat. Phenom. 11 (3) (2016) 1–17, http://dx.doi.org/10.1051/mmnp/201611301. [8] F. Gao, General fractional calculus in nonsingular power-law kernel applied to model anomalous diffusion phenomena in heat-transfer problems, Therm. Sci. 21 (2017) http://dx.doi.org/10.2298/TSCI170310194G, 194–194. [9] R.L. Bagley, Power law and fractional calculus model of viscoelasticity, AIAA J. 27 (10) (2017) 1412–1417, http://dx.doi.org/10.2514/3.10279. [10] H.G. Sun, Z. Yong, D. Baleanu, C. Wen, Y.Q. Chen, A new collection of real world applications of fractional calculus in science and engineering, Commun. Nonlinear Sci. Numer. Simul. 64 (2018) 213–231, http://dx.doi.org/10.1016/j.cnsns.2018.04.019. [11] L. Liu, L. Zheng, Y. Chen, F. Liu, Fractional boundary layer flow and heat transfer over a stretching sheet with variable thickness, J. Heat Transfer 140 (9) (2018) 091701, http://dx.doi.org/10.1115/1.4039765. [12] L. Liu, L. Zheng, F. Liu, Research on macroscopic and microscopic heat transfer mechanisms based on non-fourier constitutive model, Int. J. Heat Mass Transfer 127 (2018) 165–172, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.06.011. [13] L. Liu, F. Liu, Boundary layer flow of fractional maxwell fluid over a stretching sheet with variable thickness, Appl. Math. Lett. 79 (2018) 92–99, http://dx.doi.org/10.1016/j.aml.2017.10.008. [14] M.A. Matlob, Y. Jamali, The concepts and applications of fractional order differential calculus in modelling of viscoelastic systems: A primer, Appl. Math. Lett. (2017) 1–36. [15] R. Metzler, T.F. Nonnenmacher, Fractional relaxation processes and fractional rheological models for the description of a class of viscoelastic materials, Int. J. Plast. 19 (7) (2003) 941–959, http://dx.doi.org/10.1016/s0749-6419(02)00087-6. [16] M. Du, Z. Wang, H. Hu, Measuring memory with the order of fractional derivative, Sci. Rep. 3 (7478) (2013) 3431, http://dx.doi.org/10.1038/srep03431. [17] D. Tripathi, S.K. Pandey, S. Das, Peristaltic flow of viscoelastic fluid with fractional maxwell model through a channel, Appl. Math. Comput. 215 (10) (2010) 3645–3654, http://dx.doi.org/10.1016/j.amc.2009.11.002. [18] J. Zhao, L. Zheng, X. Chen, X. Zhang, F. Liu, Unsteady marangoni convection heat transfer of fractional maxwell fluid with cattaneo heat flux, Appl. Math. Model. 44 (2017) 497–507, http://dx.doi.org/10.1016/j.apm.2017.02.021.

1636

M. Zhang, M. Shen, F. Liu et al. / Computers and Mathematics with Applications 78 (2019) 1621–1636

[19] A. Kasaeian, R. Daneshazarian, O. Mahian, L. Kolsi, A.J. Chamkha, S. Wongwises, I. Pop, Nanofluid flow and heat transfer in porous media: A review of the latest developments, Int. J. Heat Mass Transfer 107 (2017) 778–791, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.11.074. [20] D. Tripathi, O.A. Bég, A numerical study of oscillating peristaltic flow of generalized maxwell viscoelastic fluids through a porous medium, Transp. Porous Media 95 (2) (2012) 337–348, http://dx.doi.org/10.1007/s11242-012-0046-5. [21] M.K. Nayak, G.C. Dash, L.P. Singh, Heat and mass transfer effects on mhd viscoelastic fluid over a stretching sheet through porous medium in presence of chemical reaction, Propulsion Power Res. 5 (1) (2016) 70–80, http://dx.doi.org/10.1016/j.jppr.2016.01.006. [22] R. Chand, Thermal instability of rotating maxwell visco-elastic fluid with variable gravity in porous medium, J. Indian Math. Soc. 80 (1–2) (2013) 23–31. [23] Y. Xuan, L.I. Qiang, Heat transfer enhancement of nanofluids, J. Eng. Thermophys. 21 (1) (2000) 58–64, http://dx.doi.org/10.1016/S0142-727X(99) 00067-3. [24] M. Pan, L. Zheng, F. Liu, X. Zhang, Modeling heat transport in nanofluids with stagnation point flow using fractional calculus, Appl. Math. Model. 40 (21–22) (2016) 8974–8984, http://dx.doi.org/10.1016/j.apm.2016.05.044. [25] J. Buongiorno, Convective transport in nanofluids, J. Heat Transfer 128 (3) (2006) 240–250, http://dx.doi.org/10.1115/1.2150834. [26] Y.Z. Povstenko, Theory of thermoelasticity based on the space-time-fractional heat conduction equation, Phys. Scr. 2009 (T136) (2009) 014017, http://dx.doi.org/10.1088/0031-8949/2009/T136/014017. [27] Y.Z. Povstenko, Fractional cattaneo-type equations and generalized thermoelasticity, J. Therm. Stresses 34 (2) (2011) 97–114, http://dx.doi.org/10. 1080/01495739.2010.511931. [28] M. Shen, L. Chen, M. Zhang, F. Liu, A renovated buongiornos model for unsteady sisko nanofluid with fractional cattaneo heat flux, Int. J. Heat Mass Transfer 126 (2018) 277–286, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.05.131. [29] M. Pan, L. Zheng, F. Liu, C. Liu, X. Chen, A spatial-fractional thermal transport model for nanofluid in porous media, Appl. Math. Model. 53 (2018) 622–634, http://dx.doi.org/10.1016/j.apm.2017.08.026. [30] X. Chen, Y. Ye, X. Zhang, L. Zheng, Lie-group similarity solution and analysis for fractional viscoelastic mhd fluid over a stretching sheet, Comput. Math. Appl. 75 (8) (2018) http://dx.doi.org/10.1016/j.camwa.2018.01.028, S0898122118300440. [31] T. Atanacković, K. Sanja, O. Ljubica, Z. Dušan, The cattaneo type space-time fractional heat conduction equation, Contin. Mech. Thermodyn. 24 (4–6) (2012) 293–311, http://dx.doi.org/10.1007/s00161-011-0199-4. [32] Y. Povstenko, Thermoelasticity based on space-time-fractional heat conduction equation, Solid Mech. Appl. 219 (2015) 171–190, http://dx.doi.org/ 10.1007/978-3-319-15335-3_6. [33] L. Liu, L. Zheng, F. Liu, X. Zhang, Heat conduction with fractional cattaneo-christov upper-convective derivative flux model, Int. J. Therm. Sci. 112 (2017) 421–426, http://dx.doi.org/10.1016/j.ijthermalsci.2016.11.008. [34] F. Liu, P. Zhuang, Q. Liu, Numerical Methods of Fractional Partial Differential Equations and Applications, Science Press, Beijing, 2015. [35] W.C. Tan, T. Masuoka, Stokes First problem for an oldroyd-b fluid in a porous half space, Phys. Fluids 17 (2) (2005) 023101, http://dx.doi.org/10.1063/ 1.1850409. [36] M. Khan, T. Hayat, S. Asghar, Exact solution for mhd flow of a generalized oldroyd-b fluid with modified darcy’s law, Internat. J. Engrg. Sci. 44 (5) (2006) 333–339, http://dx.doi.org/10.1016/j.ijengsci.2005.12.004. [37] C. Friedrich, Relaxation and retardation functions of the maxwell model with fractional derivatives, Rheol. Acta 30 (2) (1991) 151–158, http: //dx.doi.org/10.1007/BF01134604. [38] F. Liu, P. Zhuang, V. Anh, I. Turner, K. Burrage, Stability and convergence next term of the difference methods for the space-time fractional advectiondiffusion equation, Appl. Math. Comput. 191 (1) (2007) 12–20, http://dx.doi.org/10.1016/j.amc.2006.08.162. [39] B. Carnahan, H. Luther, J.O. Wilkes, Applied Numerical Methods, John Wiley and Sons, New York, 1969, http://dx.doi.org/10.1002/aic.690160604.