Solid Boundary as Energy Source and Sink in a Dry Granular Dense Flow: A Comparison between Two Turbulent Closure Models

Solid boundary as energy source and sink of the turbulent kinetic energy of the grains, and its influence on the mean and turbulent features of a dry granular dense flow, are investigated by using the proposed zeroand first-order turbulent closure models. The first and second laws of thermodynamics are used to derive the equilibrium closure relations, with the dynamic responses postulated by a quasi-static theory for weak turbulent intensity. Two closure models are applied to analyses of a gravity-driven flow down an inclined moving plane. While the calculated mean porosity and velocity correspond to the experimental outcomes, the influence of the turbulent eddy evolution can be taken into account in the first-order model. Increasing velocity slip on the inclined plane tends to enhance the turbulent dissipation nearby, and the turbulent kinetic energy near the free surface. The turbulent dissipation demonstrates a similarity with that of Newtonian fluids in turbulent boundary layer flows. While two-fold roles of the solid boundary are apparent in the first-order model, its role as an energy sink is more obvious in the zero-order model.


Introduction
Dry granular dense flows are continuous motions of a large amount of discrete solid particles with interstitial space filled by a gas, moving with slow to moderate speed.The grain-grain interaction, in contrast to creeping or rapid flows, results from ong-term enduring frictional contact and sliding, and short-term instantaneous inelastic collision [1]- [4].Two-fold grain-grain interactions induce fluctuations on the field quantities at the macroscopic level, a phenomenon similar to turbulent motion of Newtonian fluids with two distinctions: 1) it emerges from grain-grain interactions, in contrast to those resulted from incoming flow instability, instability in transition region, or flow geometry in Newtonian fluids; and 2) it emerges at slow speed, in contrast to those in Newtonian fluids which are strongly velocity-dependent, characterized by the critical Reynolds' number [5] [6].
Solid boundary has been demonstrated to be an energy source and sink of the turbulent kinetic energy of the grains, and conventional no-slip condition of velocity is not valid [7] [8].Whereas these influences were hardly accounted for in laminar flow formulations, e.g.[9]- [17], they were not appropriately taken into account in the limiting turbulent flow formulations, e.g.[18]- [21].Thus, the goal of the present study is to propose a thermodynamically consistent turbulent closure model to account for these effects, with their influence on the mean and turbulent flow features.Specifically, a zero-and a first-order closure models are proposed, with the focus on the intercomparison of the roles played by the solid boundary, and the influence of velocity slip.
In the following sections, the mean balance equations, state space and entropy inequality are presented for two models.The closure relations are summarized as results from thermodynamic considerations of the first and second laws.Two closure models are applied to analyses of stationary gravitational flows down an inclined moving plane.While solutions of two models demonstrate a qualitative agreement with experimental outcomes in the mean porosity and velocity profiles, the distributions of the turbulent dissipation are similar to those of Newtonian fluids in turbulent boundary layer flows, with their vanishing and finite values obtained on the free surface by the zero-and first-order models, respectively.Increasing velocity slip near the inclined plane tends to enhance the turbulent dissipation nearby, resulting in larger mean porosity and turbulent kinetic energy near the free surface.While boundary as energy source and sink is apparent in the first-order model, its latter role is more obvious in the zero-order model.

Mean Balance Equations and Equilibrium Closure Relations
To account for the distribution of solid volume and its microstructural effect, the (solid) volume fraction ν , defined as the total solid volume divided by the volume of a granular representative volume element (RVE), is introduced, with its time evolution described by the Wilmánski model for dense flows [12] [22].A dense flow is considered a rheological fluid, which must satisfy the basic laws of motion for continuum mechanics.Since in turbulent motion the field quantities experience fluctuations, with solutions random and unpredictable, their statistically averaged values (e.g.Reynolds-averaging) should be defined and investigated.With these, the following mean balance equations must be satisfied [5] The variables and parameters in (1)-( 6) are defined in Table 1.
Equations (1) 1.2 , (2) 1.2 and (3) 1 are respectively the conventional mean balances of mass, linear momentum, angular momentum, internal energy and entropy for a continuum, with the mean density ρ decomposed into  the time evolution of ν , Equation (4) 1 is the phenomenological generalization of the Mohr-Coulomb model for the mean internal friction in a granular material at low energy and high-grain volume fraction [12] [25] [26], while Equations (4) 2 and ( 5) are the balances of turbulent kinetic energy and dissipation, respectively.They are included to denote the influence of the energy cascade from the mean flow scale toward the smallest (dissipation) scale in turbulent flows.In doing so, two turbulent closure models are constructed: Equations ( 1)-( 4) apply for the zero-order model with the turbulent dissipation considered a closure relation, and Equations (1)-( 5) apply for the first-order -k ε model, in which the time evolutions of the turbulent kinetic energy and dissipation are des- cribed independently and separately.
For the application of two models, the quantities are introduced as the primitive mean fields, with the superscripts 0 and 1 denoting the model specification, while the closure relations , , , , , , , , , , ,  ,   =  , , , , , , , , , , ,  , , , are constructed based on the state spaces given by , , , , ,  , , , , ,  ,  , , , with ( ) ( ) and ∇ the Nabla operator.The quantity M ϑ is the material coldness, with T ϑ the granular coldness, a si- milar concept to granular temperature [20] [23]- [25] [27] [28].While it is used to index both variations in turbulent kinetic energy and dissipation in the zero-order model, it is employed only for the variation in turbulent kinetic energy in the first-order model.The quantity 0 ν is the value of ν in the reference configuration, included due to its influence on flowing granular matter.In (9) ∇ denote the influence of turbulent fluctuation, while D and Z are for the viscous and rate- independent effect, respectively.The forms of the closure relations are reduced by the second law of thermodynamics, which is formulated here as the Müller-Liu entropy principle.In its local form, it represents the restriction that the mean entropy production must be non-negative, i.e., T 0 s . A physically admissible process must simultaneously satisfy this Equation, (1)-( 2) and (3) 2 - (5).One can account for all these requirements by using the method of Lagrange multiplier, viz., with the mean balance equations appearing as constraints of the entropy inequality, and Substituting ( 8) and ( 9) into (10) with the assumption of material isotropy and chain rule of differentiation, the corresponding restrictions on forms such as (8) have been defined elsewhere [23] [24].They are expressions for 0 h and 1 h , 0 k and 1 k , as well the dependence of the specific turbulent Helmholtz free energies T0 ψ and T1 ψ , defined by − , for the zero-and first- order models, respectively, viz., with ( ) , , , , and ( ) and ( ) at an thermodynamic equilibrium state denoted by the subscript E | , defined by , , , , , , , in the zero-and first-order models, respectively, are given by . The variables p and β stand for the turbulent thermodynamic and configurational pressures, respectively, viz., for both models.Otherwise, for incompressible grains, p is an independent field and can no longer be determined by Equation (20) 1 ; Equations (11) 1 , ( 18) and ( 19) are simplified to [12] 1 1 C. Fang while (11) 2-3 and ( 12)-( 16) remain unchanged, with (17) reducing to 0 ε = K .

Closure Models
For isothermal flows with incompressible grains, we assume that 0 t , 1 t , 0 K may be decomposed according to with ( ) ; the superscript D indicates dynamic response, which should vanish at thermodynamic equilibrium.Within a quasi-static theory, the dynamic responses are assumed to depend explicitly and linearly on the independent dynamic variables, respectively of the forms, ( ) and the three invariants ( ) , , of D .The functional dependences of ( ) and ( ) λ λ µ µ are the same as above, with ε additionally included.We choose the material vis- cosities ( ) and phenomenological (turbulent) viscosities ( ) with Ξ depending on the three invariants of D ; ( ) , a positive constant; ν ∞ the mean vol- ume fraction corresponding to the possible most dense packing of the grains; and m ν the critical mean volume fraction at which shearing is decoupled from dilatation.Equation (28) asserts that the total stress is larger in turbulent than in laminar flows, and is justified for weak turbulent intensity [29] [30].
For explicit expressions of 0 t , 1 t , 0 R and 1 R , the simplest form of the specific turbulent free energy is proposed following previous works by [12] [17] [26] ( ) ( ) (29) with the plastic contribution confined within T f ψ .Equation ( 29) is justified for weak turbulent flows, and as- serts that smaller granular coldness results in smaller free energy [27]- [30].A hypoplastic form for the mean production of mean internal friction is given by [31] [32] ( ) ( ) to account for rate-independent characteristics, in which ( ) and a a positive constant, relating to the internal friction and frictional angle in the critical state.The scalar functions s f and d f are the stiffness and density factors, denoting strain harding/soft- ing and mean-pressure dependent bulk density, respectively.

Field Equations and Boundary Conditions
Consider a fully-developed, isothermal, two-dimensional stationary turbulent shear flow down an inclined moving plane, as shown in Figure 1.With this, (  [33] [34].Since in the critical state d f is set to be unity, Equation (30) reduces to 2 , ϕ the critical friction angle [34].Since s f does not vanish in general, substituting (37) into (4) 1 yields with s ν the minimum mean volume fraction.With these, the mean field equations reduce to ( ) , where Equations ( 40)-(42) apply for the zero-order model, and (40)-( 41) and ( 43)-( 44) apply for the first-order model.
Due to velocity slip, the grains on the solid plane may not assume the plane velocity.Velocity slip provides extra energy flux toward to, or away from the granular body, which is proportional to the square of the slip velocity and with the same direction of the momentum flux on the boundary [7] [8].On the contrary, due to the experimental setup [37], ν approaches a fixed value on the solid plane.Since experiment is carried out by discharging a constant mass flux on the plane, the flow thickness is fixed, and the shear force on the free surface is negligible due to the significant density difference between the granular body and air.Thus, the boundary conditions are given by

Nondimensionalization and Numerical Method
With the dimensionless parameters defined in Table 2, Equations ( 40)-( 44) are recast respectively by ( ) ( ) with the dimensionless boundary conditions, where The two-point nonlinear BVPs ( 46)-( 50) are solved numerically by means of the method of successive approximation with under-relaxation scheme.For the implementation of numerical simulation, the values of b ν , m ν , ν ∞ and s ν are given by [12] [16] [35]   0.51, 0.
with fixed values of T b ϑ  and b ε , motivated by experimental observations.

Numerical Results
As a parametric study, numerical simulations are carried out for variations in σ , ς and ς ′ , with  is chosen to match the experiments [37].1.0, 0.9, 0.8 σ = indicated by the arrows, with the dashed lines representing laminar flow solutions [12].Decreasing σ tends to increase velocity slip on the solid plane, resulting in more intensive friction near the solid plane with significant turbulent dissipation, as shown in Figure 2 The numerical results approach better to the experimental measurements with smaller values of σ in both models, resulted from the fact that velocity slip is identified near the solid plane in experiments.This finding suggests that the influence of velocity slip needs be taken into account for more accurate numerical prediction.The dimensionless turbulent dissipations, shown in Figure 2(f) and Figure 2(j), decrease from their maximum values on the solid plane toward the minimum values on the free surface with an "exponential-like" tendency in variations in σ in both models.Although such a tendency is similar to that of Newtonian fluids in turbulent boundary layer flows [5] [6], the first-order model is more justified, for finite turbulent dissipation is obtained on the free surface, where the turbulent kinetic energy is maximum (see

Conclusions and Discussions
Boundary as energy source and sink, and the influence of velocity slip near solid boundary on the mean and turbulent features of a dry granular dense flow, were investigated by the proposed zero-and first-order closure models, in which the granular coldness was introduced to index both variations in the turbulent kinetic energy and dissipation in the former model, while they were indexed separately by two independent fields in the latter model.Both models were applied to analyses of isothermal, stationary turbulent shear flows with incompressible grains down an inclined moving plane.Velocity slip near solid boundary tends to enhance turbulent dissipation in both models.The turbulent dissipation profile is similar to that of Newtonian fluids in turbulent boundary layer flows.The first-order model is however more justified, for it asserts that intensive turbulent kinetic energy induces intensive turbulent dissipation, with non-vanishing turbulent dissipation obtained on the free surface, in contrast to vanishing turbulent dissipation identified by the zero-order model.In both models, the mean shearing of the solid plane is less efficiently transferred toward the granular body, and the turbulent dissipation is confined within a thin layer above the solid plane.Outside this thin layer, the grains are dominated by gravity, and collide with one another in a free manner, resulting in significant short-term grain interaction, as reflected by larger mean porosity, velocity and turbulent kinetic energy near the free surface.Two-fold roles played by the solid boundary are more apparent in the first-order model, while boundary as energy source is less apparent in the zero-order model.Comparison with experiments shows qualitative agreement in the ν -and u  -pofiles, and also suggests that velocity slip needs be taken into account for more ac- curate numerical prediction.Although the velocity profiles from the zero-order model approach better to experimental outcomes, the first-order model is better to account for the influence of turbulent eddy evolution by using a - k ε energy cascade.

,
and the symmetry of the mean Cauchy stress is required.Equation (3) 2 is the Wilmánski model for mean velocity component in the - x direction, the mean volume fraction, the granular coldness and the mean internal friction components, respectively; and first-order model.The flow corresponds to the critical state, defined as the state in which of a at the critical state; and c

.Figure 1 .
Figure 1.Gravity-driven stationary flow down an inclined moving plane and the coordinate.

and b ε the boundary values on the inclined plane; 0 V
the velocity of the solid plane; s V the grain velocity on the plane; ϕ and ϕ′ the coefficients relating respectively the contributions of the slip velocity 0 s V V − to the boundary values of T ϑ and ε ; and L the flow thickness.The prescription of ε on the plane applies only for the first-order model.In doing so, solid boundary as an energy source and sink is taken into account by prescribing the values of T ϑ and ε separately in the first-order model, and by a fixed boundary value of T ϑ in the zero-order model.Nonvanishing T b ϑ and b ε correspond to experimental ob- servations at vanishing slip velocity.

Figure 2
Figure 2 illustrates the profiles of 1 ν −  (the mean porosity), u  and the dimensionless turbulent kinetic energy and dissipation, in which

Figure 2 (
Figure 2(c) with Figure 2(g) with Figure 2(a), and Figure 2(d) and Figure 2(h) with Figure 2(b)).The numerical results approach better to the experimental measurements with smaller values of σ in both models, resulted from the fact that velocity slip is identified near the solid plane in experiments.This finding suggests that the influence of velocity slip needs be taken into account for more accurate numerical prediction.The dimensionless turbulent dissipations, shown in Figure2(f) and Figure2(j), decrease from their maximum values on the solid plane toward the minimum values on the free surface with an "exponential-like" tendency in variations in σ in both models.Although such a tendency is similar to that of Newtonian fluids in turbulent boundary layer flows[5] [6], the first-order model is more justified, for finite turbulent dissipation is obtained on the free surface, where the turbulent kinetic energy is maximum (see Figure2(e) and Figure2(f)), in contrast to vanishing turbulent dissipation from the zero-order model (see Figure2(i) and Figure2(j)).This is due to that in the first-order model, the turbulent kinetic energy and dissipation are indexed separately byT

Figure 2 (Figure 2 .
Figure 2(c) with Figure 2(g) with Figure 2(a), and Figure 2(d) and Figure 2(h) with Figure 2(b)).The numerical results approach better to the experimental measurements with smaller values of σ in both models, resulted from the fact that velocity slip is identified near the solid plane in experiments.This finding suggests that the influence of velocity slip needs be taken into account for more accurate numerical prediction.The dimensionless turbulent dissipations, shown in Figure2(f) and Figure2(j), decrease from their maximum values on the solid plane toward the minimum values on the free surface with an "exponential-like" tendency in variations in σ in both models.Although such a tendency is similar to that of Newtonian fluids in turbulent boundary layer flows[5] [6], the first-order model is more justified, for finite turbulent dissipation is obtained on the free surface, where the turbulent kinetic energy is maximum (see Figure 2(e) and Figure 2(f)), in contrast to vanishing turbulent dissipation from the zero-order model (see Figure 2(i) and Figure 2(j)).This is due to that in the first-order model, the turbulent kinetic energy and dissipation are indexed separately by T ϑ

Figure 4
Figure 4 illustrates the profiles of 1 ν −  , u  and the dimensionless turbulent kinetic energy and dissipation from the first-order models, in which [ ] 0.1, 0.5, 0.9 ς ′ = indicated by the arrows.Increasing ς ′ tends to en- hance the energy flux from the granular body toward solid plane, resulting in more intensive turbulent dissipations near the solid plane, with slightly affected turbulent kinetic energy profiles, as shown in Figure 4(d) and Figure 4(c), respectively.Comparing Figure 4(d) with Figure 3(d) illuminates two different mechanisms for the turbulent dissipation: In the latter figure it is induced by the enhanced energy sink on the boundary, yielding enhanced turbulent dissipation near the solid plane.In the former figure it is induced by the intensive turbulent kinetic energy near the free surface, giving rise to the enhanced turbulent dissipation there.In view of these, the profiles of 1 ν −  and u  are only slightly affected by larger values of 2 ς , as displayed in Figure 4(a) and Figure 4(b), respectively.Comparing Figure 4(c) and Figure 4(d) with Figure 3(g) and Figure 3(h) illustrates that the first-order model delivers more justified and accurate estimations on the turbulent kinetic energy and dissipation, and the influence of boundary as energy source can be accomplished.

Figure 4 .
Figure 4. Profiles of ( ) 1 ν −  , u  and the dimensionless turbulent kinetic energy and dissipation from the first-order model,

Table 1 .
Variables and parameters in the mean balance equations.
D symmetric part of ∇v ; α  material derivative of α w.r.t.α , α′ mean and fluctuating values of α , respectively; γ mean density of the solid grains; ε specific turbulent dissipation; η mean specific entropy; ν mean volume fraction; π mean entropy production; φ mean entropy flux; φ′ turbulent entropy flux; Ω any mean orthogonal rotation of a granular RVE;

Table 2 .
Dimensionless parameters in the dimensionless field equations.