Effects of Shear Dependent Viscosity and Variable Thermal Conductivity on the Flow and Heat Transfer in a Slurry

In this paper we study the effects of variable viscosity and thermal conductivity on the heat transfer in the pressure-driven fully developed flow of a slurry (suspension) between two horizontal flat plates. The fluid is assumed to be described by a constitutive relation for a generalized second grade fluid where the shear viscosity is a function of the shear rate, temperature and concentration. The heat flux vector for the slurry is assumed to follow a generalized form of the Fourier's equation where the thermal conductivity k depends on the temperature as well as the shear rate. We numerically solve the governing equations of motion in the non-dimensional form and perform a parametric study to see the effects of various dimensionless numbers on the velocity, volume fraction and temperature profiles. The different cases of shear thinning and thickening, and the effect of the exponent in the Reynolds viscosity model, for the temperature variation in viscosity, are also considered. The results indicate that the variable thermal conductivity can play an important role in controlling the temperature variation in the flow.


Introduction
In fossil fuel combustion applications, coal particles are generally transported to the combustion chamber via either pneumatic or hydraulic transport [1,2].In certain cases, coal-water slurry is even heated prior to testing and use [3].Coal slurries exhibit non-Newtonian flow characteristics.Theoretical modeling of these processes requires a complete understanding of all the phenomena involved, requiring constitutive modeling of the various components along with the usual balance laws.Transport properties such as the viscosity and the thermal conductivity for these complex fluids are important parts of the formulation.To model the behavior of such non-linear fluids, one can use either statistical theories or continuum theories, in addition to the phenomenological/experimental approaches.Within the continuum mechanics approach, there are two distinct ways of studying these flows.In certain applications such as fluidization or gasification, one may represent the mixture as a mixture of two or more interacting continua; this is the multi-phase approach [2,4,5].In this approach, more specific information such as particle distribution, particle velocity and temperature can be obtained.Alternatively one can model the complex mixture, using a single-phase approach, with a constitutive relation for the suspension as a whole.In this approach one only works with the equations of motion for the suspension and obtains information about the global characteristics of the suspension.This is the approach that we take in this paper.
Many industrial fluids exhibit non-Newtonian behavior.These non-linear fluids, such as polymers, molten plastics, suspensions and slurries have been widely utilized in many manufacturing and processing industries and understanding their friction and heat transfer characteristics is becoming more important.Much work has been devoted to the study of non-Newtonian fluids for the best performance in practical applications [6][7][8][9].Considerable efforts have been directed toward the measurement of viscous properties in isothermal systems.It has been established that the thermal conductivity of Newtonian fluids is shear rate independent while there could be a dependency on temperature.Unlike Newtonian fluids, however, it has been observed that many non-Newtonian fluids are affected by the temperature, shear rate, volume fraction, etc.Many of the experimental studies on the thermal conductivity of non-Newtonian fluids have been performed under static conditions [10,11]; very few studies have been conducted over a range of shear rates.Among these cases where the influence of shear rate on the thermal conductivity of viscoelastic polymer solutions is considered one can name Cocci and Picot [12], Chitrangad and Picot [13], Wallace et al. [14], Loulou et al. [15], Chaliche et al. [16], Shin [17], Lee and Irvine [1], and Kostic and Tong [18].The results, however, have been contradictory.
The way the thermal conductivity depends on temperature, shear rate, viscosity and volume concentration, is of practical significance in mathematical formulation and the numerical modeling.Chaliche et al. [16] showed that by using cone-and-plane systems, that thermal conductivity increases gradually with increasing shear rate, depending either on temperature or on polymer concentrations.A higher temperature or a lower concentration results in a higher thermal conductivity.Sohn and Chen [19] analytically studied the influence of shear-rate-dependent thermal conductivity upon the heat transfer characteristics of suspensions.Their theoretical study explained the increase of thermal conductivity by introducing the concept of micro-scale convection, which is induced by the rotation of particles in a shear flow.This concept is useful in understanding the mechanism of increasing thermal conductivity in suspensions.Charunyakorn et al. [20] numerically studied the heat transfer in a slurry using a micro-encapsulated phase-change material, with a shear-rate-dependent thermal conductivity.
Lin et al. [21] reported the measurements of thermal conductivity for two concentrated fruit juices as well as some of their rheological properties; their results indicate that the thermal conductivity increases as the shear rate increases for the shear-thinning juice concentrates, which also agrees with the results of previous research for polymer.They stated that the fluid structures become more aligned along the streamlines and thus more orderly owing to the shearing.Such orientation of the structures should necessarily affect the thermal conductivity in a similar manner through the impact of structural arrangement on the effective diffusivity [22].Shin and Lee [23] proposed a new correlation for the shear-rate-dependent thermal conductivity of suspensions based on their experimental investigations.Their study shows that the thermal conductivity is also strongly related to the size of the dispersed particles.Many recent works have been devoted to investigating thermal conductivity of nanofluids [24][25][26].
The objective of the present paper is to study the heat transfer in the fully developed flow of a slurry, and explore the effects of shear-rate-dependent viscosity and thermal conductivity.The arrangement of the paper is as follows.In the next section, the governing equations of motion and heat transfer are provided.Section 3 focuses on the constitutive relations for the stress tensor and the heat flux vector.In Section 4, we describe the geometry of the problem and provide the derivation for the one-dimensional form of the governing equations, as well as the boundary conditions.In Section 5, we outline the numerical scheme we have applied.We also perform a mesh-independence study, and present a comparison between the numerical results and analytical results.In Section 6, the numerical results are presented through a parametric study by varying the dimensionless numbers.

Governing Equations of Motion and Heat Transfer
If the slurry is treated as a single component (phase) material (If the slurry is modeled as a multi-component material, then the governing equations should be given for all the components, and a multi-phase approach should be taken.This requires not only constitutive relations for each component, but also for the interactions among the components), then, in the absence of any electro-magnetic effects, the governing equations of motion are the conservation of mass, linear momentum, and energy equations: Conservation of mass: where  is the density of the fluid; /t is the partial derivative with respect to time; div is the divergence operator; and u is the velocity vector.For an isochoric motion,  =0.
Conservation of linear momentum: where b is the body force vector; T is the Cauchy stress tensor; and d/dt is the total time derivative, given by where  is the specific internal energy; L is the gradient of velocity; q is the heat flux vector; and r is the specific radiant energy.Thermodynamical considerations require the application of the second law of thermodynamics or the entropy inequality.The local form of the entropy inequality is given by (see [27], p. 130): where (, ) is the specific entropy density; (, ) is the entropy flux; and s is the entropy supply density due to external sources, and the dot denotes the material time derivative., where θ is the absolute temperature, then Equation ( 4) reduces to the Clausius-Duhem inequality: Even though we do not consider the effects of the Clausius-Duhem inequality in our problem, for a complete thermo-mechanical study, the Second Law of Thermodynamics must be considered [27][28][29][30].To achieve "closure" for these equations, constitutive relations are needed for T, r, and q.In the next section we will discuss briefly the constitutive modeling issues.

Constitutive Relations
For a thermo-mechanical description of a problem, in addition to the balance laws, one requires constitutive relations for the stress tensor, the heat flux vector, the internal energy, and the radiation.The effects of radiation will be ignored in the present study.We will next briefly discuss the constitutive relations which we will use for the stress tensor and the heat flux vector.

Stress Tensor
Some studies indicate that the viscosity of coal-water mixtures depends not only on the volume fraction of the coal particles, and the mean size and the size distribution of the coal, but also on the shear rate, since the slurry behaves as a shear-rate dependent fluid [31].There are also studies which indicate that preheating the fuel (slurry) results in better performance [32,33]; as a result of such heating, the viscosity of the slurry changes.Constitutive modeling of these non-linear fluids, commonly referred to as non-Newtonian fluids, has received much attention [34].
In general, non-Newtonian fluids differ from Newtonian fluids in at least two ways: (1) they exhibit normal stress effects, such as rod-climbing and die-swell; and (2) they exhibit shear-thinning or shear-thickening which is the decrease or increase in viscosity with increasing shear rate, respectively.Both these phenomena introduce non-linearities into the equations.Among other differences one can name viscoelastic effects such as creep and relaxation or the presence of yield stress for some non-linear fluids.Perhaps the simplest model which can capture the normal stress effects (which could lead to phenomena such as 'die-swell' and 'rod-climbing', which are manifestations of the stresses that develop orthogonal to planes of shear) is the second grade fluid, or the Rivlin-Ericksen fluid of grade two [30,35].We will use a modified form of the second grade fluid (See [36] for a review).
In this paper we use the constitutive model given in [37] and advocated by Miao et al. [38,39], where: where: is the second invariant of the symmetric part of the velocity gradient,  is the identity tensor, and m is a material parameter.When m = 0 and there are no temperature or volume fraction effects, we recover the standard second grade fluid model, and when m < 0, the fluid is shear-thinning, and if m > 0, the fluid is shear-thickening [40].Notice that m = n − 1, where n is the usual power-law exponent for the power-law fluid model.Also  is the volume fraction: the function  is an independent kinematical variable called the volume distribution or volume fraction function (related to concentration) having the property 0 ≤ (x, ) ≤   < 1.The function  is represented as a continuous function of position and time; in reality,  in such a system is either one or zero at any position and time, depending upon whether one is pointing to a particle or to the void space (fluid) at that location.Now,  is related to   (density of pure fluid) and  through  = (1 − )  .  is the maximum crystal fraction in which flow can occur,  is the temperature,  0 and  0 are reference values,  is a constant and  1 and  2 are material moduli, which are commonly referred to as the normal stress coefficients.In an important paper, Fosdick and Rajagopal [41] show that irrespective of whether  1 +  2 is positive, the fluid is unsuitable if  1 is negative.In particular, they showed that if it is assumed that  > 0,  1 < 0,  1 +  2 ≠ 0 (*), which as many experiments have reported to be the case "for those fluids which the experimentalists assume to be constitutively determined by (4), at least sufficiently well as a second order approximation" ( [41], p. 147), then certain anomalous results follow.Fosdick and Rajagopal [41] proved a theorem which indicates that if (*)2,3 hold, then an unusual behavioral property, not to be expected for any rheological fluid occurs, namely, "that the larger the viscosity, keeping everything else fixed, the faster that initial data is amplified in motions which take place under strict isolation."For further details on this and other relevant issues in fluids of differential type, we refer the reader to the review article by Dunn and Rajagopal [42].The kinematical tensors A 1 and A 2 are defined through: Equation ( 6) uses some of the ideas proposed by Gupta and Massoudi [31] who also provided a generalized form of a constitutive model for the second grade fluid by allowing the shear viscosity to be a function of temperature where () was given by the Reynolds viscosity model: where  = ( 1 −  0 ) .This viscosity law was first proposed by Reynolds [43] in his theory of lubrication, where he deduced an empirical formula based on the experimental results of the viscosity of olive oil at different temperatures; this expression is oftentimes sued in lubrication processes.
In the model described in Equation ( 6), the apparent viscosity is assumed to be a function of temperature and volume fraction, following the Einstein-Roscoe relation [44,45] Note that in Equation ( 6) we interpret p as the pressure of the fluid suspension which will be a function of the volume fraction and thus in a sense the model is a generalization of the second grade compressible-type fluid.That is, in this model, we are assuming that the presence of the particles affects the viscosity and the pressure of the fluid.In this paper, we will study the slurry flow in a channel, based on the assumption that the stress tensor for the fluid can be described by the constitutive Equation ( 6).This model incorporates into the second grade fluid model, a viscosity function which depends on the shear rate, temperature and volume fraction.
The thermodynamics and stability of second grade fluids have been studied in detail by Dunn and Fosdick [46].They show that if the fluid is to be thermodynamically consistent in the sense that all motions of the fluid meet the Clausius-Duhem inequality and that the specific Helmholtz free energy of the fluid is a minimum in equilibrium, then

Heat Flux Vector
The classical theory of heat conduction, first proposed by Fourier [47,48] and later generalized by Duhamel (see [49]) assumes that the constitutive relation for the heating flux h is a linear function of the temperature gradient  = (, )∇, where  is the temperature, F is the deformation gradient, and K is the thermal conductivity tensor.Fourier was the first person [50] to state that heat conduction depends on the temperature gradient and not on the temperature difference between the two adjacent parts of a solid body, where: where for more complex materials k is considered as an effective or modified form of the thermal conductivity [51].In general, k can also depend on concentration, temperature, etc., and in fact, for anisotropic material, k becomes a second order tensor.For complex materials, such as polymers and granular materials, whether stationary or in motion, the thermal conductivity of the material is assumed to depend on parameters such as volume fraction, particle size, shear rate, etc. [1,52,53].
In heat transfer studies of suspensions, composites, porous and granular materials one needs to know a priori whether the material is a porous one with a second phase interspersed in it, or whether the material is an assembly of distinct particles with the interstices filled with a gas or a liquid [54,55].In general, as Batchelor [56] points out: "Most of the bulk properties of two-phase systems about which one would like to be able to make predictions are 'transport' properties; that is to say, they represent the ability of the material as a whole to transfer some conservable quantity, such as heat, mass, electric, or momentum, in response to an imposed gradient of intensity of that conservable quantity."From the continuum mechanics perspective, these transport properties are modeled as parts of the constitutive relations for the stress tensor, heat flux vector, etc.
Recently, Massoudi [57,58] assumed that the heat flux vector q for a density-gradient-dependent fluid can be given by:  = (, , , , , grad ) (12) where u is the velocity;  is the temperature;  is the density [  = (1 − )  ];  = grad and  = grad .Then frame-indifference implies  = (, , , , ).This is a generalization of a case presented in Bowen ([59], p. 161) for a compressible heat conducting viscous fluid.That is, the density gradient (or volume fraction gradient) plays a role, not only in the distribution of the materials, but also in the way in which it influences the heat conduction.Massoudi [57,58] then showed that a general representation for the heat flux vector given by Equation ( 12) is: where a1 -a6 are scalar functions of the appropriate invariants of the tensor and vector terms.It can be seen that (i) when  2 =  3 =  4 =  5 =  6 = 0 , and  1 =  = − , then we recover the standard Fourier's Law,  = −∇, and (ii) when  1 =  = −, and  2 = , we have Interestingly, Soto et al. [60] showed that based on molecular dynamics (MD) simulations of inelastic hard spheres (IHS), the basic Fourier's law has to be modified for the case of fluidized granular media, similar to Equation (14).It is noted that Wang [61] also derived a general expression for the heat flux vector for a fluid where heat convection is also important; he assumed that  = (, ∇, , , ) where f is a vector-valued function,  the temperature, ∇ is the gradient of temperature, v the velocity vector, L its gradient, and X designates other scalar-valued thermophysical parameters.A simplified form of Equation ( 13) was used in a recent study by Yang et al. [62].In general, this method is very difficult (If we assume that q can depend explicitly on the temperature gradient, concentration gradient, etc., then clearly the problem would become more non-linear.An example of this type of constitutive relation is a power-law type model suggested by Rodrigues and Urbano [63], where  = −∇ = −|∇| −2 ∇ where 1 <  < ∞ .When  = 2 , the above equation reduces to the usual Fourier's law) as it requires knowledge of the various coefficients in Equation ( 14) and in the absence of many physical experiments, the only option is to do a parametric/numerical study.This was the case in [62].In the present study, we take a different, and perhaps a more practical approach, namely we use the results of experiments to obtain correlations for the shear-dependent thermal conductivity.
Based on the experimental data of Lee and Irvine [1] we obtain an expression for the thermal conductivity which depends on the shear rate and temperature: When the values of the constants are chosen as follows: In this paper, we will use Equations ( 11) and ( 15) for the thermal conductivity.In the next section, we will discuss a simple boundary value problem.

Geometry and the Kinematics of the Flow
We consider the fully developed flow of a slurry (a suspension of solid particles in a fluid) characterized by Equation ( 6) in a heated channel driven by a pressure gradient along the channel (see Figure 1).The two plates are at y = 0 and y = H.We assume that the constitutive equations for the stress tensor and the heat flux vector are given by Equations ( 6) and (11) Using Equation ( 16), the conservation of mass is automatically satisfied.Substituting Equation ( 6) into Equation ( 2), the momentum equations are obtained in x and y directions: The pressure p is assumed to depend on the volume fraction (For example, as an engineering approximation, it is sometime assumed that the overall pressure drop in a pneumatic conveying system has two components: one is the pressure drop due to the gas alone and another one is the pressure drop due to the presence of the solid particles (see [46], p. 172).That is, ∆   = ∆  + ∆  where it is assumed that the ∆  is related to the gas pressure through a parameter  where  is based on material properties, such as density, and systems parameters such as the loading ratio, etc.): where we assume () =  1 ; where  1 and  2 are constants; and  1 ,  2 < 0. With these assumptions Equations ( 17) and ( 18) become: For the balance of energy (Equation ( 3)), the specific internal energy  is related to the specific Helmholtz free energy through (see [46]): where  is specific Helmholtz free energy;  is the specific entropy.With Equation ( 16) it follows that: Also, the specific radiant energy, r, is assumed to be negligible.The energy equation becomes: where the first term is a measure of viscous dissipation (see [64]).Substituting Equation ( 15) into Equation ( 11) and we have: To solve the governing equations, we need two boundary conditions for the velocity field and two boundary conditions for the temperature, and one boundary condition for the volume fraction.At the two walls we apply the no-slip boundary condition for the velocity [5], and we impose constant values for the temperature: where  1 >  0 .We also assume an integral condition for the volume fraction, where the integral value of the volume fraction over the cross section is constrained to be a constant: where   is a measure of the amount of particles in the system.
In order to non-dimensionalize the equations, we introduce the following dimensionless parameters: where H is the distance between the two walls; V is a reference velocity;  0 and  1 are the temperature of the lower and upper walls respectively.The dimensionless forms of the governing equations are: where:  1 ( ̅ ) =  1 ( ̅ + ) 3 Now,  1 is a measure of the effect of the pressure gradient along the flow direction compared to the viscous effects;  2 represents the ratio of the normal stress differences to the force caused by the pressure in the y direction;  3 is a measure of viscous dissipation;  4 is the ratio of the gravity to the pressure gradient;  is a reference temperature representing the ratio of the upper wall temperature to the lower wall temperature; and  1 ~7 represent the dimensionless forms of the coefficients in Equation (15).
The dimensionless forms of the boundary conditions are: where N is an average value of  ̅ integrated over the cross section (an input to the problem).

Numerical Scheme
Numerical results are obtained by solving the governing Equations ( 30)-( 32) together with the boundary conditions Equations ( 36)-(38) using finite difference method.The three equations need to be integrated simutaneously for the three unknowns: velocity, volume fraction, and temperature.We apply the Newton Raphson's method.Multiplying both sides of Equation ( 30) by  ̅  ̅ , and then substituting Equation ( 31) into (30) we obtain: where We first assume initial guesses for the velocity, volume fraction and temperature to start the iteration.The value of the volume fraction at point j,  ̅  , can be obtained using the Newton's method through the following scheme: where: and  ̅  is replaced with  ̅ _ in the next step until it satisfies the convergence criterion: With the new values of volume fraction at each point, the velocity is simply calculated through: The solution of temperature follows a similar procedure as the volume fraction while: where: With all the values of velocity, volume fraction and temperature obtained through Newton iteration, we then check if the results (in the current step) have converged compared to the previous step by: where: To apply the integral condition given by Equation ( 38), the volume fraction values are integrated using the composite Simpson's rule: where ℎ = 1  with the domain split up in n subintervals equally.This estimated value is then used to correct the initial guess of the volume fraction for the next iteration by  ̅ _ =  ̅  + ∆ ̅  , where: until the iteration finally satisfies the criterion: We also performed a mesh-independence study to make sure the results are obtained independent of the mesh size.The tolerance values also change as the mesh is refined: We can see from Figure 2 that the numerical scheme can be considered to be independent of the mesh size.Considering Figure 2a which shows some minor deviations in the velocity profiles when ℎ = 0.1, 0.05, we choose ℎ = 0.02 in the current study.Also, in order to verify the accuracy of the numerical scheme we attepmt to find an analytical solution for the governing equations by testing a limiting case (an idelaized or special case) where  = 0,  = 0,  1 = −2,  2 = 0,  3 = 50,  4 = 0,  = 1.5,  = 0.3,  1 ~7 = 0.The solution is: where . The analytical solutions are compared with the numerical solution with the same values of the dimensionless numbers and mesh size of 0.02.The two solutions are plotted together in Figure 3.It is clear that the numerical solution coincide with the analytical solution, verifying that the numerical scheme with the mesh size of 0.02 is capable of providing good results.

Results and Discussion
We will now look at the effects of m, M,  1 ,  2 ,  3 ,  4 by changing their values within a designated range as shown in Table 1.
Table 1.Designated values of the dimensionless numbers and parameters.The representative values of other dimensionless numbers and parameters will be shown below the following figures when we discuss the effects of each parameter.

Effect of m
Figure 4 shows the effect of the power-law exponent m on velocity, volume fraction and temperature.It can be seen that the index m has major effects on the velocity distribution.It is evident from Figure 4a that the velocity increases as the exponent m becomes larger.When the fluid is shear-thickening (m > 0), it moves faster and produces a larger velocity gradient in the flow domain.This behavior is a common feature of shear-thickening fluids as shown for example in Figure 3.3 of Chhabra and Richardson [65].In these situations, the center-line velocity, when compared to the center-line velocity for a Newtonian fluid (m = 0), is higher for a shear-thickening fluid ( > 0) than for a shear-thinning fluid ( < 0).We also observe a similar phenomenon in Figure 4a,b depicts the volume fraction distribution, which is almost not affected by m.The exponent m also has some minor influence on the temperature profile, which can be seen in Figure 4c.The slope of the temperature profile slightly changes towards a constant value as the value of m becomes more negative.

Effect of M
Figure 5a displays the effect of M, the exponent in the Reynolds viscosity model, on the velocity profile.While M affects the velocity profilel significantly, it seems it has little influence on the temperature as shown in Figure 5c, where it is noticed that the viscous term (1 −  ̅ ) −2.5  − ̅  /2 decreases when the value of M increases.This conforms to the results in Figure 5a where velocity becomes larger as M increases, due to a smaller viscosity.Figure 5b also shows that a larger M leads to a greater gradient of volume fraction only in the region near the upper wall.The greater velocity gradient causes the volume fraction to vary in a broader range.Such difference in the four curves of volume fraction is more evident at the upper wall where the temperature is higher than most of the flow domain.When M approaches zero, the temperature is changing almost linearly from the lower wall to the upper one.

Effect of 𝐵
Recall that  1 represents a measure of the effects of the pressure gradient compared to the viscous effects .The pressure gradient in the x direction is the driving force in this problem.We find that when  1 is more negative, the fluid flows faster everywhere, which is evident in Figure 6a.When  1 is close to zero, there is very limited movement of the fluid.As the absolute value of  1 increases, the velocity of the fluid increases rapidly.Also Figure 6b shows that the gradient of volume fraction is greater in the region near the upper wall corresponding to the greater velocity gradient.When  1 approaches zero, the volume fraction is almost the same everywhere.The gradient of temperature is affected by  1 in a similar manner to the case of variable m and M.  Recall that B 2 is a measure of the effecs of the ratio of the normal stress differences to the force caused by the pressure.Figure 7 indicates that B 2 only has a significant effect on the volume fraction.It is clear from Equation (31) that Since the velocity gradient is not affected too much by  2 , we can observe from Figure 7b that when  2 is relatively large;  ̅ has a larger value in the region close to the upper wall, thus  ̅  ̅ becomes more negative, which agrees with Equation (56).  the velocity and the volume fraction distributions.Recall that B 3 only appears in the energy equation and its effect on velocity and volume fraction is indirect.Figure 8c shows that a larger B 3 makes the temperature distribution more non-linear and helps develop thermal boundary layers adjacent to the walls.When B 3 is close to zero, i.e. negligible viscous dissipation, the temperature variation is linear.B 4 is the ratio of the effect of gravity to the pressure gradient in the direction of gravity.
When B 2 = 0, we can obtain an analytical solution for the volume fraction from Equation (31), that is: When  4 is greater than 1, i.e., the effect of gravity is more evident, the volume fraction varies more non-linearly.When  4 vanishes, volume fraction is almost a constant around N. Figure 9a shows that the velocity profiles are also affected.The slope of the velocity profile near the lower wall becomes milder as  4 increases.The volume fraction is much larger with a larger  4 at the lower wall, and the viscosity is thus larger if we look at the viscous term (1 −  ̅ ) −2.5  − ̅  /2 .( 1 −  0 ) 3   The two dimensionless numbers R 1 and R 5 represent the effects of the coefficients in Equation (15), based on the experimental correlations for the thermal conductivity.We only consider the effects of R 1 and R 5 since they typically denote the highest (third) order effect of temperature.In order to elucidate the effects of R 1 and R 5 more evidently, the values of other parameters (i.e., R 2 ~R4 , R 6 , R 7 ) are selected to be zero so as to minimize or eliminate their effects.From Figure 10 and Figure 11, it is apparent that R 5 influences the velocity, volume fraction and temperature profiles.Figure 10c indicates that the inclusion of R 1 (i.e., when the thermal conductivity is dependent on the shear rate) causes the temperature to vary in a more non-linear fashion.Some differences can also be detected in the velocity profiles (Figure 10a).It can be observed that the value of the maximum velocity increases as R 1 becomes more negative.Figure 11 shows the important role that  5 plays in this problem.When  5 is small, the temperature in the interior of the flow domain is higher everywhere, with all the profiles being non-linear.If we consider the relationship between  5 and the thermal conductivity, we can see that when  5 increases, the thermal conductivity also increases.Thus, with a higher  5 , heat conduction occurs at a higher rate throughout the flow, and the flow cools down from the upper wall temperature to the lower wall temperature.We also find that at  5 = 0, the temperature can vary in a fairly non-linear fashion and a thermal boundary layer develops near the plates.We should notice that when  5 vanishes (provided that the representative values of other   are set to zero), the thermal conductivity becomes a constant.Thus the heat flux in such a case only depends on the temperature gradient.Figure 11a,b show that the velocity and the volume fraction are also dependent on the value of  5 .When  5 increases, the fluid flows more slowly, and correspondingly, the volume fraction varies in a wider range with a larger gradient.

Frictional Effects and Heat Transfer Rate at the Walls
In addition to the above dimensionless numbers and parameters, two other quantities of interest are the skin friction coefficient and the Nusselt number.The skin friction coefficient is related to the shear stress at the walls and is defined by: where   is the wall shear stress: Now,   can be written as: where The Nusselt number is defined by: where  is the convective heat transfer coefficient of the fluid,   is the temperature at the wall.In our problem   is either  0 or  1 at the upper or the lower wall respectively.2-5, to demonstrate the frictional effects and heat transfer at the walls.Table 2 shows that a shear-thickening ( > 0) fluid has a larger velocity gradient at both walls and a shear-thinning ( < 0)fluid has a smaller velocity gradient compared to a Newtonian fluid ( = 0).However, the frictional effect has a different trend.The skin friction coefficient near the lower wall is larger when the fluid is shear-thinning. ̅ ′ (0) increases and  ̅ ′ (1) decreases when  increases, indicating that the temperature variation is more non-linear for a shear-thinning fluid.A similar trend for  ̅ ′ and  ′ can be observed in Table 3 when varying the value of M, while the skin friction coefficients have an opposite trend.Table 4 indicates that when  1 (a measure of the pressure gradient) becomes more negative, the absolute value of the velocity gradient and the skin friction coefficient become larger at the upper and lower walls.For a given value of , , and  3 ,  ̅ ′ (0) increases and  ̅ ′ (1) decreases when  1 becomes more negative.It can be observed in Table 5 that varying  3 does not have a major influence on  ̅ ′ and   .The value of  ′ (0) increases and  ′ (1) decreases when the effect of viscous dissipation is larger.

Concluding Remarks
In this paper, we have suggesteed a new constitutive model for the heat flux vector with the thermal conductivity depending on the shear rate and temperature; we have studied the fully developed flow and heat transfer a of a slurry between two horizontal flat plates, where the upper plate is at a higher temperature.The fluid is assumed to be described by a constitutive relation for a generalized second grade fluid where the shear viscosity is a function of the shear rate, temperature and concentration.The heat flux vector for the slurry is assumed to follow a generalized form of the Fourier's equation where the thermal conductivity k, is based on the experimental results of Lee and Irvine [1], where k depends on the temperature as well as the shear rate.We have performed a parametric study for various dimensionless numbers.The governing equations were solved using Newton's method, and the effects of the dimensionless numbers and various parameters on the velocity, volume fraction and temperature profiles were discussed.The power-law exponent m has a major influence on the velocity.The exponent in the Reynolds viscosity model, M, also affects the velocity, since a larger value of M causes smaller viscosity.We also discussed the effects of the coefficients  1 and  5 in the correlation of thermal conductivity, which significantly affect the variation of the temperature.The frictional effects and the heat transfer rate at the boundaries are considered in terms of the skin friction coefficient and the Nusselt number.Finally we would like to discuss some of the limitations of our model and our formulation of the problem.As mentioned in the introduction, it is possible to model a suspension as a multi-phase fluid where the interaction effects can be studied.Since we have modeled the slurry as a single phase non-linear fluid, even though the concentration of the particles is introduced as a kinematical variable, we really are not studying a two-phase mixture composed of solid particles and a fluid.As a result many interesting phenomena such as particle clustering, particle agglomeration and the stability of the particle motion are not studied here and is beyond the scope of the present study [66,67].The multi-phase approach is necessary at higher volume fractions (see Johnson et al. [68,69] and Massoudi [2,4] or the book by Rajagopal and Tao [5]).As a result one can study the interesting phenomenon related to particle migration [70][71][72][73].Another possible case of modeling flow of densely packed particles is normally referred to as flow of dry granular materials where the effects of interstitial gases or liquids are ignored and the dense suspension is modeled as granular materials.In general, one can use continuum mechanics, statistical mechanics, or numerical simulations to model granular materials [74,75].

6. 5 . 3
Effect of  3 is a measure of viscous dissipation.As we can see from Figure 8a,b B 3 shows minor effects on