MODELING GRANULAR MATERIALS AS COMPRESSIBLE NONLINEAR FLUIDS : HEAT TRANSFER BOUNDARY VALUE PROBLEMS MEHRDADMASSOUDI

We discuss three boundary value problems in the flow and heat transfer analysis in flowing granular materials: (i) the flow down an inclined plane with radiation effects at the free surface; (ii) the natural convection flow between two heated vertical walls; (iii) the shearing motion between two horizontal flat plates with heat conduction. It is assumed that the material behaves like a continuum, similar to a compressible nonlinear fluid where the effects of density gradients are incorporated in the stress tensor. For a fully developed flow the equations are simplified to a system of three nonlinear ordinary differential equations. The equations are made dimensionless and a parametric study is performed where the effects of various dimensionless numbers representing the effects of heat conduction, viscous dissipation, radiation, and so forth are presented.


Introduction
Multiphase studies present challenging and interesting problems for engineers, physicists, and mathematicians.Different approaches, such as experimental, theoretical, and computational techniques, often combined, have been used.From a mathematical point of view, understanding the behavior and response of each phase, whereby the material is modeled based on sound and physical principles, is the first step towards understanding the more complex multicomponent flows which also include interactions among the different constituents.Some of the early studies were concerned with the engineering and structural design of bins and silos.The inaccuracies of these theories, especially for the dynamic conditions of loading or emptying, occasionally have resulted in failure of the bin or silo [Savage [49]].In many engineering applications such as discharge through bin outlets, flow through hoppers and chutes, flow in mixers, and slurry transports, one needs information on the flow pattern and particle distribution.
A poor understanding of the flow of granular materials can have serious economic consequences in the design and handling of large-scale solids transport plants.Engineers have been forced to resort to costly, cut-and-try methods of design.In order to design equipment such as bins and silos, combustors, hoppers, chutes, hydrocyclones, and so forth in an effective and economical way, a thorough understanding of the various factors governing the flow characteristics of granular materials must be obtained.Heat transfer plays a major role in storage, handling, and processing of bulk solids such as grains, powders, pharmaceutical products, and so forth.Flowing granular streams are being considered for some solar power plants and fusion reactor chambers (see Gudhe et al. [19]).In fossil fuel applications, the self-heating of coal stockpiles has a long history of posing significant problems to coal producers because it lowers the quality of coal and may result in hazardous thermal runaway.Precise prediction of the self-heating process is, therefore, necessary in order to identify and evaluate the control strategies for safe coal mining, storage, and transportation.In a number of applications, granular materials are heated prior to processing, or cooled after processing (see Patton et al. [38]).It can be seen that very little fundamental work, from a mathematical point of view, has been devoted to these types of heat transfer processes in granular materials (see Kaviany [26]).
Granular materials exhibit nonlinear phenomena like yield stress and normal stress differences, the latter is usually referred to as dilatancy (see Reynolds [45]).The normalstress phenomenon is a characteristic of nonlinear fluids and nonlinear elastic solids.For flowing granular materials, where particles are colliding with and sliding over each other, methods in continuum mechanics can be used.In this approach one assumes that the material properties of the ensemble may be represented by continuous functions.One of the early continuum models for flowing granular materials was proposed by Goodman and Cowin (see [17,18]) and Cowin (see [9,10]).Another approach is based on the techniques used in the kinetic theory of gases (see Goldhirsch [16]).In these cases, the effects of the interstitial fluid are neglected and the granular material is modeled as a single phase continuum.In general, however, one should take a multicomponent approach (see Rajagopal and Tao [42]).
The two important constitutive relations needed for the study of flow and heat transfer in granular materials are the heat flux vector and the stress tensor.For most applications the Fourier's law of heat conduction is assumed for the heat flux vector, where an effective thermal conductivity is prescribed.In this paper, we use this approach.For the stress tensor, we use a model developed by Rajagopal and Massoudi [39], where the effects of density gradients are included.This model has been used to study various problems such as flow in a vertical pipe (cf.Gudhe et al. [20]), heat transfer and flow on an inclined plane (see Gudhe et al. [19]), and flow due to natural convection (cf.Massoudi and Phuoc [31]).At the same time this model has been used within the context of mixture theory to study problems of practical interest in multiphase applications (Massoudi et al. [32]).Recent review articles by Savage [49], Hutter and Rajagopal [21], and de Gennes [11], and books by Mehta [34], Duran [12], and Antony et al. [1] address many of the interesting issues in the field of granular materials.
M. Massoudi and T. X. Phuoc 3 In this paper, we review and discuss three important boundary value problems: (i) the flow down an inclined plane with the effects of radiation heat transfer at the free surface; (ii) the natural convection flow between two heated vertical walls; (iii) the shearing motion between two horizontal flat plates, where the effects of viscous dissipation are included.It is assumed that the material behaves like a continuum, similar to a compressible nonlinear fluid.For a fully developed flow the equations are simplified to a system of three nonlinear ordinary differential equations.The equations are made dimensionless and a parametric study is performed, where the effects of various dimensionless numbers representing the effects of heat conduction, viscous dissipation, and so forth are presented.

Governing equations
The balance laws, in the absence of chemical reactions and electromagnetic effects, are the conservation of mass, linear momentum, angular momentum, and energy (see Truesdell and Noll [56]).The conservation of mass in the Eulerian form is given by where ∂/∂t is the partial derivative with respect to time.The balance of linear momentum is where d/dt is the material time derivative, b is the body force, and T is the Cauchy stress tensor.The balance of angular momentum (in the absence of couple stresses) yields the result that the Cauchy stress is symmetric.The energy equation in its general form is where ε denotes the specific internal energy, q is the heat flux vector, r is the radiant heating, and L is the velocity gradient.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 Liu [27, page 130]) where η(x,t) is the specific entropy density, φ(x,t) is the entropy flux, and s is the entropy supply density due to external sources, and the dot denotes the material time derivative.
If it is assumed that where θ is the absolute temperature, then (2.4) reduces to the Clausius-Duhem inequality For a complete study of a thermo-mechanical problem the second law of thermodynamics has to be considered (see Collins and Houlsby [8]).In other words, in addition to other principles in continuum mechanics such as material symmetry, frame indifference, and so forth, the second law also imposes certain restrictions on the type of motion and/or the constitutive parameters.Since, there is no general agreement on the functional form of the constitutive relation and since the Helmholtz free energy is not known, a complete thermodynamical treatment of the present model used in our studies is lacking.Instead, we use the results of Rajagopal et al. (see [40,43]) to obtain restrictions on the material parameters, while using simple and reasonable expressions for the undetermined coefficients.In order to "close" the governing equations, we need constitutive relations for T, q, and r.At this juncture it would be appropriate to point out that in theories for rapid flow of granular materials based on a kinetic theory approach, the fluctuations in the velocity field give rise to the notion of "granular temperature."The convective heat transport, within the context of such theories, is determined by the fluctuations in the velocity field.In the continuum approach that we have taken, the fluctuations in the velocity field are ignored, and the phenomenon of heat transfer is incorporated in the energy equation.To include in addition to the energy equation, the notion of granular temperature would be inconsistent with our approach.The continuum approach is applicable when the packing of the material is reasonably compact, that is, high volume fraction, and the fluctuations from the mean are not significant.

Constitutive equations
3.1.Stress tensor.Reynolds (see [45,46]) observed that for a shearing motion to occur in a bed of closely packed particles, the bed must expand to increase the volume of its voids.He termed this phenomena dilatancy.Reiner [44] proposed a continuum model to describe the mechanics of wet sand.This model does not take into account how the volume fraction affects the stress; however, using his model, Reiner [44] showed that the application of a nonzero shear stress produces a change in volume.McTigue [33] extended the Reiner-Rivlin model to granular materials.Bagnold [3,4] performed experiments on neutrally buoyant, spherical particles suspended in Newtonian fluids undergoing shear in coaxial rotating cylinders.He distinguished three different regimes of flow behavior, which he termed macro-viscous, transitional, and grain-inertia.The interesting phenomenon was the presence of a normal stress proportional to the shear stress, similar to that of the quasistatic behavior of a cohesionless material obeying the Mohr-Coulomb criterion (see Spencer [51], Schaeffer [50]).Many investigators have modeled flow of particles as a fluid phenomenon, even as a compressible fluid (see Tardos [54]).Many researchers have formulated nonlinear fluid-type models to describe the behavior of granular materials (see Reiner [44], McTigue [33], Savage [48], Goddard [15], Astarita and Ocone [2]).For a review see the article by Elaskar and Godoy [13].
In this paper, we use Rajagopal and Massoudi's model [39] where the Cauchy stress tensor T depends on the manner in which the granular material is distributed, that is, the volume fraction ν, its gradient, and the symmetric part of the velocity gradient tensor D. That is, Using standard procedures in mechanics such as frame-indifference, isotropy, and so forth, certain restrictions on (3.1) can be found (cf.Truesdell and Noll [56]).Additional restrictions on the form of the constitutive expression can also be obtained from internal constraints, such as incompressibility and thermodynamics restrictions due to Clausius-Duhem inequality (cf.Müller [36]).A constitutive model that predicts the possibility of normal stress-differences and is properly frame invariant was given by Rajagopal and Massoudi [39]: where where γ is constant.The material parameters β o -β 5 are assumed to have the following forms: (3.4) The above representation can be viewed as Taylor series approximation for the material parameters (see Rajagopal et al. [40]).Such quadratic dependence, at least for the viscosity β 3 , is on the basis of dynamic simulations of particle interactions (see Walton and Braun [57]).Further restrictions on the coefficients have been obtained by using the following argument.Since the stress should vanish as ν → 0, we can conclude that This is a principle of the limiting case.That is, if there are no particles, then ν and grad ν are zero, and the stress should be zero; however, the kinematical terms D, D 2 , and trD, multiplied by β 2 , β 3 , and β 5 do not necessarily go to zero when there are no particles.
Homogenization Motion Therefore, to ensure this we impose the above restrictions.Furthermore, Rajagopal and Massoudi [39] and Rajagopal et al. [43] have shown that as compression should lead to densification of the material.They gave the following rheological interpretation to the material parameters: β 0 (ν) is similar to pressure in a compressible fluid and is to be given by an equation of state, β 2 (ν) is like the second coefficient of viscosity in a compressible fluid, β 1 (ν) and β 4 (ν) are the material parameters connected with the distribution of the granular materials, β 3 (ν) is the viscosity of the granular materials, and β 5 is similar to what is referred to as the "cross-viscosity" in a Reiner-Rivlin fluid.
Looking at (3.2) with (3.4) it can be shown that this model is capable of predicting both normal stress differences in a simple shear flow problem (for details see Massoudi and Mehrabadi [30]).In our analysis in this paper, for simplicity, we assume β 5 is zero and as a result only one of the normal stress differences is nonzero.Furthermore, the volume fraction field ν(x,t) plays a major role in this model.That is, even though we talk of distinct solid particles with a certain diameter, in this theory, the particles through the introduction of the volume fraction field have been homogenized, as shown in Figure 3.1.In other words, it is assumed that the material properties of the ensemble are continuous functions of position.That is, the material may be divided indefinitely without losing any of its defining properties.A distributed volume and a distributed mass M. Massoudi and T. X. Phuoc 7 can be defined, where the function ν is an independent kinematical variable called the volume distribution function and has the property The function ν is represented as a continuous function of position and time; in reality, ν in a granular system is either one or zero at any position and time, depending upon whether one is pointing to a granule or to the void space at that position.That is, the real volume distribution content has been averaged, in some sense, over the neighborhood of any given position.The classical mass density or bulk density, ρ, is related to ρ s and ν through It should be mentioned that in practice ν is never equal to one; its maximum valued, generally designated as the maximum packing fraction, depends on the shape, size, method of packing, and so forth.

Heat flux vector.
The constitutive relation for the heat flux vector is assumed to be given by the Fourier's law of conduction where the heat flux vector q is linearly related to the temperature gradient: where K is an effective or modified form of the thermal conductivity.The dependence of thermal conductivity on material properties and volume fraction is one of the challenging problems in mechanics of microstructure materials.There have been some studies with regard to porous materials but few with regard to flowing granular materials.For example, Kaviany (see [26, page 129]) presents a thorough review of the appropriate correlations for the thermal conductivity for packed beds and the effective thermal conductivity in multiphase flows.Massoudi (see [28,29]) has recently given a brief review of this subject and has also proposed and derived a general constitutive relation for the heat flux vector for a flowing granular media.In this paper, we use the following simple approximations (see Bashir and Goddard [5]): where ψ is the ratio of conductivity of the particle to that of the matrix and K m is the thermal conductivity of the matrix.We use this approximation for two of the boundary value problems considered in this paper, namely, the natural convection problem and the shearing motion.For the inclined flow problem, we use the more nonlinear case due to Jeffrey [24] which includes the second-order effects in the volume fraction: where where where α is the ratio of conductivity of the particle to that of the matrix, κ the effective conductivity of the suspension, κ m the conductivity of the matrix, and ν is the solid volume fraction.

Boundary value problems
Whenever nonlinear constitutive relations are studied, the solution procedures for solving the governing equations, whether analytical or computational, are more complicated.Exact solutions are very rare in heat transfer studies of nonlinear materials or multiphase flows.Next to the exact solutions, finding solutions to simple boundary value problems are extremely desirable.Most of the constitutive relations used in mechanics, whether non-Newtonian models, turbulence models, and so forth, when substituted in the general governing equations, that is, the balance laws, would produce a system of partial differential equations which at times are impossible to solve completely with the numerical techniques currently available.Therefore, from a modeling point of view, it is worthwhile to study problems where due to simplification of the kinematics of the flow or the boundary conditions, one obtains a system of (nonlinear) ordinary differential equations.The solution to these simpler problems would be useful for at least two different reasons: (i) they provide insight into the nature of these nonlinear constitutive relations, and (ii) they provide cases where the accuracy or convergence of solutions to the general multidimensional equations can be tested.Other interesting phenomena such as stability and uniqueness of solutions also sometimes arise.Furthermore, the higher order terms in the constitutive relations require additional boundary conditions, and this itself is an important problem to study.
In most applications the solid particles are transported via a fluid: for example, in coal water slurries, the liquid is either water or oil.In these problems, the driving force is the fluid's pressure gradient, and as a result the motions of the constituents are coupled.In certain cases the driving force is due to the boundary (as in the shearing motion) or due to the gravitational effects (or some other fields such as buoyancy, electric, or magnetic fields), and if the particles are densely packed, the effects of the interstitial fluid can be ignored; in such cases one can use the continuum mechanics of a single constituent to study the flow of the granular media (where u and D are not zero).The three problems selected in this paper represent a few of these special cases where the driving force is due to gravity (Section 4.1), temperature gradient (Section 4.2), and shearing motion (Section 4.3).Free surface flow down an inclined plane.Flow of granular materials down an inclined plane occurs naturally as in the cases of avalanches and mudslides; it is also used for transporting and drying of bulk solids (such as agricultural and pharmaceutical products).It is a viscometric flow (see Truesdell [55]) and one which amends itself to fundamental theoretical and experimental studies.Hutter et al. (see [22,23]) showed that the existence (or nonexistence) of solutions for fully developed flow down an inclined plane depends on the type of boundary conditions that are imposed.Rajagopal et al. [43] also studied existence and uniqueness of solutions to the equations for the flow of granular materials down an inclined plane, where the thermal effects were also included.They used the model developed by Rajagopal and Massoudi [39]; in that study they delineated a range of values for the material parameters, which ensures existence of solutions to the equations under consideration.They also proved rigorously that for a certain range of values of the parameters no solutions exist, while for others there is a multiplicity of solutions.
In this section of the paper, we will look at the flow down a heated inclined plane (see also Gudhe et al. [19]).In many problems, viscous dissipation is ignored.There are many cases in polymer rheology and lubrication, for example, where viscous dissipation cannot be neglected (see Szeri and Rajagopal [53], Szeri [52]).For densely packed granular materials, as particles move and slide over each other, heat is generated due to friction and therefore, in such problems, the viscous dissipation should be included.It is assumed that the flow is fully developed (see Figure 4.1).The free surface is exposed to high ambient temperature and as a result a modified Stefan-Boltzmann correlation for radiation is used at that surface (see Fuchs [14, page 331]).
We make the following assumptions: (i) the motion is steady; (ii) the effects of radiant heating r are imposed at the free surface; (iii) the constitutive equation for the stress tensor is given by (3.2) and the constitutive equation for the heat flux vector is that of Fourier's law, given by (3.11), (3.13)-(3.15); (iv) the density (or volume fraction), velocity, and temperature fields are of the form With the above assumptions, the equation of conservation of mass is automatically satisfied and the momentum and the energy equations in their reduced forms become (i) y-direction momentum: (ii) x-direction momentum: (iii) heat transfer: From (4.2) it is clear that we need two boundary conditions for the volume fraction ν, and (4.3) and (4.4) indicate that two conditions are required for the velocity and the temperature fields, respectively.We can also see that (4.2) is not coupled to the other two equations, and it can therefore be integrated first.Once the volume fraction field is determined, (4.3) can be solved for "u," and finally (4.4) is integrated to find the temperature field "θ."At the surface of the inclined we assume the no-slip condition for the velocity and a constant temperate θ w .Thus, at y = 0 : At the free surface, the no-traction boundary condition is imposed on the stress tensor, and as a result we obtain two expressions for the velocity gradient and the volume fraction (see (4.6), (4.7)), and for the temperature we apply the Stefan-Boltzmann condition (the importance of radiation boundary condition for composite plates is discussed by Miller and Weaver [35], for fluid-particle flow by Chamkha [7], and for packed beds by Wu and Chu [58]; equation (4.8) is really our first approximation, and a more appropriate one for the case of granular materials might be to introduce into the equation a function for the dependence of the volume fraction, for example, q = f (ν)εσ(θ 4 h − θ 4 s )) (see also Saldanha da Gama [47]) when the surrounding temperature is designated as θ ∞ and the temperature at the free surface is θ h .Thus, we have, at y = h, du dy = 0, (4.6) ) We can see that we still need an additional boundary condition for ν.There are at least two options: (i) we can impose a distribution function for ν at the wall, which could have a constant value (this may mean gluing particles to the surface), or (ii) we can give an average value for ν integrated over the cross section (a measure of the amount of particles present in the system).These two conditions can be written as Equations (4.2)-(4.4)are made dimensionless before they are solved numerically.This is a procedure which we will follow throughout the paper.Let us define the dimensionless distance y, the velocity u, and the temperature θ by the following equations (see Na [37]): where u o is a reference velocity, θ w is the wall temperature, and h is the constant height to the free surface.The above equations then become where and the boundary conditions become at y = 0 : dν dy ) In this problem we will use condition (4.14) for the volume fraction.The dimensionless parameters have the following physical interpretations: R 1 could be thought of as the ratio of the pressure force to the gravity force; R 2 is the ratio of forces developed in the material, due to the distribution of the voids, to the force of gravity; R 3 is the ratio of the viscous force to the gravity force (related to the Reynolds number); R 4 is a measure of viscous dissipation, which is the product of the Prandtl number and the Eckert number; and R 5 is a measure of the emissivity of the particles to the thermal conductivity.It follows from Rajagopal and Massoudi [39] that R 1 must always be less than zero for the solution to exist and R 3 and R 4 must be greater than zero, since the viscosity is positive.In addition to these dimensionless numbers, values for N, a, b, and α are to be given a priori.The objective is to conduct a parametric study to see how the various parameters (R 1 , R 2 , R 3 , R 4 , and R 5 ) affect the solution.The volume fraction equation has a nonunique solution as discussed by Rajagopal et al. [43] and Gudhe et al. [19].Such a nonunique solution stems from the boundary condition at y = h by (4.7), where one can see that for a given value of the volume fraction at y = h there exist two possible conditions of the derivative of the volume fraction: one is negative and the other is positive.Thus, there must be two solutions to satisfy these conditions: one solution in which ν must increase monotonically from y = 0 to y = 1 so that the positive condition, (dν/dy) y=1 > 0, is satisfied and one in which ν must decrease monotonically to satisfy the negative condition at the free surface, (dν/dy) y=1 < 0. Typical results of such multiple solutions are shown in Figure 4.2.However, due to gravitational effect particles must settle down toward the surface of the inclined plane rather than floating upward to the free surface.For this reason, the solution in which ν decreases as one approaches the free surface is chosen as the correct solution.
The effects of parameters R 1 , R 2 , and R 3 on ν are shown in Figure 4.3.The distribution of the volume fraction, however, is independent of parameter R 3 .For example, by keeping other parameters constant the volume fraction decreases from 0.355 at the inclined plane to 0.2 at the free surface for all values of R 3 from 0.01 to 1.0.Since R 3 represents the viscous forces, R 1 the pressure force, and R 2 the material forces due to normal stress  effects, this result indicates that the distribution of the volume fraction is dominated by the material forces and the pressure rather than the viscous forces.The effects of R 1 , R 2 , and R 3 on the velocity profiles are shown in Figure 4.4.The velocity profile has a parabolic shape with the maximum value at the free surface.The value of the free surface velocity increases as R 2 increases but it decreases when R 1 increases.Although R 3 does not have a significant effect on ν, it has a strong influence on the velocity distribution.The results indicate that as R 3 increases the flow slows down and the velocity profile approaches a straight line.
For the temperature calculations we chose a = 0.75 and b = 1.The effect of R 3 on the temperature distribution is shown in Figure 4.5 with the surrounding temperature being higher than the surface temperature, θ ∞ = 3.5.For R 3 = 0.01 the maximum temperature is found to be located in a region within the flowfield and the temperature decreases as one approaches the inclined plane and the free stream.This is the effect of viscous heat generation.As R 3 increases the temperature profile becomes linear, increasing monotonically from the inclined surface toward the free stream.This is because the flow slows down as R 3 increases as shown in Figure 4.4, and the viscous dissipation term and the convection term in the energy equation are not as significant as the conduction term.In this case, the flowfield is heated primarily through conduction.
The effects of R 4 on the temperature profile are shown in Figure 4.6, where the temperature profiles and the free surface temperature are calculated as a function of R 4 while other parameters are kept constant.Two different cases are investigated: the cooling case, that is, when the temperature of the environment is lower than the surface temperature and the heating case, that is, when the environment temperature is higher than the surface temperature.An increase in R 4 means that the heat generated by viscous dissipation term is increased.As a result, this generated heat must be transferred to the surface and be removed from the free surface through radiation.Increasing R 4 , therefore increases both the maximum temperature in the flowfield and the free surface temperature.For example, when the environment is hotter than the free surface, as R 4 increases, the free surface      16 Heat transfer studies in granular materials   temperature increases slightly.When the environment temperature is lower than the surface temperature, an increase in R 4 results in a significant increase in the free surface temperature.For all cases, the maximum temperature and the free surface temperature are always higher than both the surface and the surrounding temperatures.The effect of R 5 on the temperature profile for cooling and heating cases is shown in Figure 4.7.For both cases the maximum temperature is found to be inside the flowfield and heat is transferred toward both the inclined surface and the free stream.This is primarily due to the viscous dissipation effect which raises the flowfield temperature to a higher value than both the surface temperature and the environment temperature.The parameter R 5 represents the heat transfer (due to radiation) at the free surface and its magnitude is a measure of how much of the heat generated by the viscous dissipation term can be transferred away.Thus, an increase in R 5 leads to a decrease in the free surface temperature until it reaches the environment temperature.When this condition is reached, the solution becomes independent of R 5 .The free surface temperature could be less than or greater than the inclined surface temperature and the environment temperature depending on the competition between R 4 and R 5 .For a constant R 4 , if R 5 is small, the heat generated due to viscous dissipation cannot be sufficiently removed by the radiation process.In this case the free surface temperature is higher than both the inclined surface and the environment temperature.If R 5 is large, however, the radiative heat removal can sufficiently balance the heat generated.In this case the free surface temperature is low and it decreases as R 5 increases.If R 5 continues to increase, the free surface temperature will reach the surrounding temperature.

Natural convection between two vertical walls.
Flow due to natural (free) convection between two vertical walls has certain applications in the operation of a Clusius-Dickel column for separating isotopes where the combined effects of thermal diffusion and free convection are considered.(See Bird et al. [6, page 300] for the solution of this problem when the fluid is assumed to be a Newtonian one.)Rajagopal and Na [41] studied the same problem for the case of a homogenous incompressible fluid of grade three.The coupling of the momentum and energy equation due to the Boussinesq's approximation presents added complications.In this section, we consider the flow due to natural convection of granular materials.It is assumed that the densely packed material is situated between two infinite vertical walls which are at different temperatures (Figure 4.8).In this arrangement, the particles near the hot wall are heated and the particles near the colder wall are cooled.As a result of such a temperature difference a natural convection flow ensues and the particles move.It is expected that such a flow depends on the competition between various forces such as the pressure force, the viscous and buoyancy forces, the gravity force, the force developed in the materials, and so forth and the transport properties.
We assume: (i) the motion is steady; (ii) radiant heating "r" is ignored; (iii) the constitutive equation for the stress tensor is given by (3.2) and the constitutive equation for the heat flux vector is that of Fourier's law, given by (3.11); (iv) the density, velocity, the temperature fields, and the body force are of the form where θ m is a reference temperature, (e.g., θ m = (θ 1 + θ 2 )/2), g is the acceleration due to gravity, γ is the coefficient of thermal expansion, and j is the unit vector in the ydirection.For thermal conductivity we assume (3.12).Using the following dimensionless parameters: we obtain the following dimensionless forms of the equations, dropping the bars for simplicity: where and the boundary conditions are The dimensionless parameters R 2 , R 3 , and R 4 have the same physical interpretations as those given in the previous section; the additional dimensionless number R 6 is a measure of the buoyancy effects due to thermal expansion.Thus, the solutions to the above equations with respect to R 2 , R 3 , R 4 , R 6 , N, and ξ as parameters will reveal the characteristics of the natural convection flow and heat transfer between the hot and the cold surfaces.We will only focus on the importance of the effects of R 1 and R 3 as typical examples for such an application; a full treatment is given in Massoudi and Phuoc [31].
For a given value of N, (4.17) describes the distribution of the volume fraction in xdirection under the competition between the pressure force and the force developed in the materials due to the distribution of the particles.Such a competition is represented by the dimensionless parameter R 1 which is defined as the ratio of the pressure force to the force developed in the materials due to the volume fraction distribution.Typical results showing the effect of the dimensionless parameter R 1 on the distribution of the volume fraction are plotted in Figure 4.9 for N = 0.9 and R 1 ranging from −1.0 to 1.0.These  values are chosen so that two possible patterns of the distribution of the materials could be demonstrated.The results indicate that the volume fraction ν decreases as the absolute value of R 1 is decreased and it becomes constant (equal to N/2) when R 1 is zero.When R 1 is not zero, however, two distribution patterns are obtained depending on the sign of R 1 .When R 1 is positive (the pressure force has the same sign as the force developed in the materials) ν reaches its maximum value at x = 0 and its minimum values are at the walls.This means that the region near the center has more particles than the regions near the walls.When R 1 is negative (the pressure force has the opposite sign to the force developed in the materials) ν has a minimum value at x = 0 and maximum values are at the walls.More particles are packed near the walls than the region near the center.
The effects of R 1 on the velocity profile are shown in Figure 4.10, when R 2 , R 3 , R 4 , and N are kept constant and R 1 = −1.5, 0.0, and 1.3.Negative R 1 (R 1 = −1.5)represents the situation under which the concentration of the materials is denser near the walls than in the region near the centerline.Positive R 1 (R 1 = 1.3) is for the case where there are more particles in the center region than in the regions near the walls.And R 1 = 0 represents the case that the material is uniformly distributed.The results on the velocity profiles indicate that the particles move with forward velocities in the region near the hotter wall and reversed velocities in the region near the colder wall.
The effects of the dimensionless parameter R 3 on the temperature profiles are shown in Figure 4.11.Parameter R 3 has a strong effect on the heating process.When R 3 is small the heat transfer at the hot wall (x = −1) has negative values.Thus this wall acts as a heat source that heats the particles.When R 3 is large, however, the heat transfer at x = −1 reverses its sign and becomes positive.The hot wall now acts as a heat sink and heat is transferred from the particles to the wall.

Shearing motion and viscous dissipation.
In this section, we will look at the flow and heat transfer to granular materials packed between two long horizontal parallel plates, where the lower plate is fixed and heated and the upper plate is set into motion and is at a lower temperature than the lower plate (see Figure 4.12).As the flow begins the particles roll over each other and slide, and due to the nonlinear effect (dilatancy) the upper plate will have to move outward (upward) to allow for the motion of the particles.This problem can be thought of as a limiting case of flow between two long vertical cylinders, where the inner cylinder is heated and the outer cylinder is rotating with a constant speed.This arrangement could be used to measure the viscosity of the material or to dry the particles, for example.We now make the following assumptions: (i) the motion is steady; (ii) radiant heating r can be ignored; (iii) the constitutive equation for the stress tensor is given by (3.2) and the constitutive equation for the heat flux vector is that of Fourier's law, given by (3.11); (iv) the density (or volume fraction), velocity, and temperature fields are of the form With the above assumptions, the conservation of mass is automatically satisfied and the momentum and the energy equations are given in their reduced dimensionless forms with where u w is the shearing velocity at the top plate and h is the distance between the two plates, as (1 + 3ξν) where and the boundary conditions become at y = 0 : or Case (b) Due to the assumed forms of the volume fraction, velocity, and temperature fields, given by (4.21), the two momentum equations in the x and y directions are not coupled and therefore (4.23) can be solved independently of the velocity and temperature equations.
Once the volume fraction is known, (4.24) is integrated to obtain the velocity field and finally the temperature is obtained by solving (4.25).If the volume fraction, velocity, and temperature are functions of x and y (or z), then we obtain a system of partial differential equations which would have to be integrated simultaneously.The dimensionless parameters have the same physical interpretations as those given by (4.12).
In solving (4.fraction depends strongly on both R 1 and R 2 , the effects of these parameters on the velocity and the temperature distributions seem to be not as significant.For the temperature distribution, however, the effects are due mostly to the parameter ξ which is related to the thermal conductivity and R 4 which is the product of the Prandtl number and the Eckert  number, representing the competition between the viscous dissipation and the heat conduction.Therefore, the shape of the temperature curves and the direction of the heat transfer depend on these parameters.If ξ is large and R 4 is small, the heat generated is conducted away and the temperature decreases from the hotter plate to the colder plate.
On the other hand, if ξ is small and R 4 is large, the heat generated cannot be conducted away.As a result, the flow field is heated significantly and heat is transferred toward both plates.In Figure 4.15, the profiles are obtained for N = 0.4, R 1 = −5, R 2 = 100, and ξ is kept constant at 0.5.For R 4 up to 5.0 the highest temperature is at the bottom plate.For R 4 = 10 and 15.0 the temperature curves have their maximum values at y = 0.2 and 0.3, respectively.
In Figure 4.16 the temperature curves are plotted for N = 0.4, R 1 = −10, and R 2 = 100.It can be seen that for R 4 = 0.5 the temperature decreases from the hot plate to the cold plate.The heating is primarily due to the heat conducted from the hot plate and for higher values of ξ the heating takes place faster.The plot for R 4 = 15 and ξ = 5.0 shows that heat is transferred from the hot plate to the particles.With the same value of R 4 the curve obtained with ξ = 0.25 shows that heat is transferred from the granular materials to the plates and the maximum temperature is located in a region between the plates.

Summary
Three different boundary value problems involving heat transfer in granular materials, namely, (i) the fully developed flow down a heated inclined plane, (ii) flow due to natural convection, and (iii) flow due to shearing motion are studied.The granular material is assumed to behave like a continuum, very similar to a compressible non-Newtonian fluid.The effects of the interstitial fluid are ignored, both in the flow and in the heat transfer analysis.In reality, there will always be some fluid in the pores and in certain cases such as fluidized beds where the fluid interacts with the solid phase; we cannot use an analysis for the single phase dense granular materials as we have presented here.In those cases, multiphase flow theories are more appropriate (see Johnson et al. [25], Rajagopal and Tao [42]).
We can see that the equations of motion and energy, for steady fully developed flows, are of the following compact forms (see Table 5 (5.2) M. Massoudi and T. X. Phuoc 27 For Problem 2,  (5.4) The main reason for doing a parametric study, via nondimensionalizing the equations of motion is that we can gain some insight into a class of problems.Since the material parameters β 0 -β 5 have not been measured experimentally, it is not possible to compare our results to any experiment.However, qualitatively, we can see that since the material parameters are, in general, functions of the volume fraction, there is a stronger nonlinearity in the equations, and therefore, numerically it will be more difficult to obtain solutions.These simple boundary value problems with all the basic assumptions specified should serve as limiting cases for more complicated flow geometries and flow conditions.Obviously, the effects of the interstitial fluid, slip at the wall, particle shape, and so forth, are important issues which need to be studied.

Figure 4 . 4 .
Figure 4.4.Dimensionless velocity as a function of the dimensionless distance: effects of R 1 , R 2 , and R 3 .

Figure 4 . 12
Figure 4.12 22)-(4.26), we use the boundary conditions (4.27) and (4.29), and we guess the values for the unknown conditions of (du/d y) y=0 , (dν/dy) y=0 , and (dθ/dy) y=0 , initially.Typical results for the volume fraction profiles are shown in Figures 4.13 and 4.14.It is clear that both parameters R 1 and R 2 have significant effects on the volume fraction distribution.The range of R 2 for a solution to exist also depends strongly on the value of R 1 .For example, with R 1 = −50 the solution does not exist if R 2 is smaller than 20, but when R 1 = −5 the solution exists for R 2 as low as 2.5.Although the volume 0

Figure 4 . 13 . 1
Figure 4.13.Distribution of the volume fraction as a function of the dimensionless distance: effect of R 1 .

Figure 4 . 14 .
Figure 4.14.Distribution of the volume fraction as a function of the dimensionless distance: effect of R 2 .