Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

54, 4, pp. 1245-1256, Warsaw 2016
DOI: 10.15632/jtam-pl.54.4.1245

ON THE TURBULENT BOUNDARY LAYER OF A DRY GRANULAR

AVALANCHE DOWN AN INCLINE. II – CLOSURE MODEL AND

NUMERICAL SIMULATIONS

Chung Fang

National Cheng Kung University, Department of Civil Engineering, Tainan City, Taiwan

e-mail: cfang@mail.ncku.edu.tw

Dynamic responses of the closure relations, specific turbulent Helmholtz free energy and
turbulent viscosity are postulated followed by experimental calibrations. The established
closure model is applied to analyses of a gravity-driven stationary avalanche with incom-
pressible grains down an incline.While themean velocity and volume fraction increase from
their minimum values on the plane toward maximum values on the free surface exponen-
tially, two-fold turbulent kinetic energies and dissipations evolve in a reversemanner. Most
two-fold turbulent kinetic energies and dissipations are confined within the thin turbulent
boundary layer immediately above the plane, with nearly vanishing two-fold turbulent kine-
tic energies andfinite two-fold turbulent dissipations in the passive layer.The two layers are
similar to those of Newtonian fluids in turbulent boundary layer flows, and are preferable
recognized by the distributions of turbulent kinetic energies and dissipations.

Keywords: closure model, gravity-driven flow, passive layer, turbulent boundary layer

1. Introduction

This paper continues Fang (2016b), hereafter referred to as Part I. The balance equations of the
mean fields for isothermal flows with incompressible grains are summarized in the following

0= ˙̄ν + ν̄∇· v̄ 0= γ̄ν̄ ˙̄v− div(̄t+R)− γ̄ν̄b̄
0= γ̄ν̄ℓ¨̄ν −∇· (h̄+H)− γ̄ν̄f̄ 0=˚̄Z− Φ̄ (̊Z̄≡ ˙̄Z− [Ω̄,Z̄])
0= γ̄ν̄k̇−R ·D̄−∇·K+ γ̄ν̄ε 0= γ̄ν̄ṡ− ℓH ·∇ ˙̄ν −∇·L+ γ̄ν̄H

(1.1)

for which

P = {p̄, ν̄, v̄,Z̄,ϑM,ϑT ,ϑG} C = {t̄,R, h̄,H, f̄,Φ̄,k,s,K,L,ε,H} (1.2)

are introduced respectively as the primitive mean fields and closure relations based on the
turbulent state space given by

Q= {ν0, ν̄, ˙̄ν,g1, γ̄ = c1,g2,ϑM = c2,g3,ϑT ,g4,ϑG,g5,D̄,Z̄} C= Ĉ(Q) (1.3)

with c1 and c2 are constants, and g2 = g3 = 0. Quantities in (1.1)-(1.3) have been defined
in Part I. Müller-Liu entropy principle has been investigated to derive the equilibrium closure
relations, with the results summarized in Table 1, in which the subscript E denotes that the
indexed quantity is evaluated at an equilibrium state, defined viz.,

Q
∣

∣

E
≡ (ν0, ν̄,0,g1,c1,0,c2,0,ϑT ,0,ϑG,0,0,Z̄) QD ≡ ( ˙̄ν,g3,g4,g5,D̄) (1.4)

withQD the dynamic sub-state space, uponwhich the dynamic closure relations should depend.



1246 C. Fang

Table 1.Thermodynamically consistent equilibrium closure relations (Fang, 2016b)

ψT = ψ̂T(ν0, ν̄,∇ν̄, γ̄ = c1,ϑM = c2,ϑT ,ϑG,Z̄) β̄ = γ̄ν̄ψTν̄
γ̄ν̄k = γ̄ν̄ϑMψT

,ϑT
γ̄ν̄s = γ̄ν̄ϑMψT

,ϑG

γ̄ν̄ε
∣

∣

E
=0 γ̄ν̄H

∣

∣

E
=0 Φ̄

∣

∣

E
=0

K
∣

∣

E
=(ϑM −ϑT)γ̄ν̄ε,g4

∣

∣

E
+(ϑM −ϑG)γ̄ν̄H,g4

∣

∣

E
− γ̄ν̄ϑMψT

,Z̄
· Φ̄,g4

∣

∣

E

L
∣

∣

E
=(ϑM −ϑT)γ̄ν̄ε,g5

∣

∣

E
+(ϑM −ϑG)γ̄ν̄H,g5

∣

∣

E
− γ̄ν̄ϑMψT

,Z̄
· Φ̄,g5

∣

∣

E

ℓ(ϑMh̄+ϑGH)= γ̄ν̄ϑMψT,g1 H= ℓRg1

f̄
∣

∣

E
=(ℓ)−1

{

(p̄− β̄)/(γ̄ν̄)+(1−ϑT/ϑM)ε, ˙̄ν
∣

∣

E
+(1−ϑG/ϑM)H, ˙̄ν

∣

∣

E
−ψT
,Z̄
· Φ̄, ˙̄ν
∣

∣

E

}

t̄
∣

∣

E
=−ν̄p̄I− γ̄ν̄ψT,g1 ⊗g1+ γ̄ν̄ψ

T
,Z̄
· Φ̄,D̄

∣

∣

E

R
∣

∣

E
=−(ϑM/ϑT −1)γ̄ν̄ε,D

∣

∣

E
− (ϑM/ϑT −ϑG/ϑT)γ̄ν̄H,D

∣

∣

E

In Section 2, the dynamic responses of the closure relations are postulated by a quasi-static
theory, followed by the specific postulates of the turbulent Helmholtz free energy, viscosities
and the hypoplastic model for rate-independent characteristics. The established closure model
is applied to analyses of a gravity-driven stationary avalanche down an incline in Section 3.
Numerical simulations are compared with laminar flow solutions. The study is concluded in
Section 4.

2. Zero-order closure model

2.1. Dynamic response

It is assumed that the closure relations consist of the equilibrium and dynamic parts viz.

C = C
∣

∣

E
+CD C ∈ {t̄,R,Φ̄,K,L, f̄, γ̄ν̄ε, γ̄ν̄H} (2.1)

Specifically, t̄D,RD, f̄D, γ̄ν̄εD, γ̄ν̄HD,KD andLD are assumed to be the quasi-static expres-
sions of QD given by

0= t̄D− ǫM ˙̄νI−λM(trD̄)I−2µMD̄ 0=RD− ǫT ˙̄νI−λT(trD̄)I−2µTD̄
0= fD+ζ ˙̄ν + δ(trD̄) 0= γ̄ν̄εD−f1 ˙̄ν −f2(trD̄)−f3(g4 ·g4)
0= γ̄ν̄HD−f4 ˙̄ν −f5(trD̄)−f6(g5 ·g5) 0=KD+f7g4 0=LD+f8g5

(2.2)

with ǫM, ǫT , ζ,f1, f3, f4, f6−8 being scalar functions of (ν0, ν̄, γ̄, ϑ
M, ϑT , ϑG); and λM, λT , δ,

µM, µT , f2, f5 scalar functions depending additionally on the three invariants (I
1
D̄
,I2
D̄
,I3
D̄
) of D̄

given by I1
D̄
≡ trD̄, I2

D̄
≡ 0.5(tr2D̄− trD̄2) and I3

D̄
≡detD̄.

In equation (2.2), Truesdell’s equi-presence principle is used, by which (ϑM t̄D+ϑTRD) and
f̄D depend explicitly and linearly on ˙̄ν and D̄; γ̄ν̄εD and γ̄ν̄HD depend explicitly and linearly on
˙̄ν, D̄,g4 andg5;K

D andLD depend explicitly and linearly ong4 andg5, respectively,motivated
by the Fourier law. Thus, a dry granular avalanche is considered a Stokes or Reiner-Rivlin fluid.
Scalar functions µM and µT are respectively the material viscosity and phenomenological (tur-
bulent) viscosity induced by turbulent fluctuation. Equation (2.2) has been applied for creeping,
dense and rapid laminar flows, and for weak turbulent dense flows as the first approximation
(Fang, 2008, 2009, 2016a; Fang and Wu, 2014; Kirchner and Teufel, 2002; Wang and Hutter,
1999).



On the turbulent boundary layer of a dry granular avalanche... 1247

2.2. Turbulent Helmholtz free energy, material and turbulent viscosities

It is assumed that ψT consists of an elastic part, ψTe , and a rate-independent part, ψ
T
f , viz.

ψT = ψTe (ν0, ν̄,g1, γ̄,ϑ
M,ϑT ,ϑG)+ψTf (I

1
Z̄
,I2
Z̄
,I3
Z̄
) (2.3)

with the irreversible effect confined within ψTf (Fang, 2009; Kirchner and Teufel, 2002; Wang

and Hutter, 1999). Following the previous works, ψTe is assumed to be expanded in a Taylor
series about ν̄ = ν̄m and |g1|=0, with ν̄m the critical mean volume fraction at which shearing
is decoupled from dilatation, see Fang (2009), Savage (1993), Wang and Hutter (2001)

ψTe =
{

α0(ν̄ − ν̄m)2+β0
( ν̄m
ν̄∞− ν̄

)2
g1 ·g1

}

Fc

Fc =
n=2
∑

n=0

1

n!

{(ϑT

ϑM

)n

+
(ϑG

ϑM

)n}

−1
(2.4)

with ν̄∞ the value of ν̄ corresponding to the denst possible packing of the grains, and {α0,β0}
dependingon {ν̄m, γ̄,ϑM}. Equation (2.4) is an extension of its laminar flow counterpartwithFc
accounting for the influence of turbulent fluctuation, motivated by the nonlinear characteristics
of rapid flows (Pudasaini and Hutter, 2007; Rao and Nott, 2008; Wang and Hutter, 2001). It
asserts that smaller two-fold granular coldnesses result in smaller free energy.

The specific forms of the material viscosity µM and turbulent viscosity µT are given by

µM = µ0γ̄
2
( ν̄m
ν̄∞− ν̄

)8
Ξ µT = µ0γ̄

2(Fc−1)
( ν̄m
ν̄∞− ν̄

)8
Ξ (2.5)

with Ξ= Ξ̂(I1
D̄
,I2
D̄
,I3
D̄
), and µ0 = µ̂0(ν0,ϑ

M), a positive constant. They are postulated followed
by the previous works (Fang, 2009; Kirchner andTeufel, 2002), with the power 8 a curve-fitting
(Savage, 1993), and the dependency of µT on ϑT and ϑG motivated by Newtonian fluids in
turbulent flow. Both µM and µT assert that total stress is larger in turbulent flows than in
laminar flows. For laminar flows, both ϑT and ϑG vanish, yielding the vanishing µT .

2.3. Hypoplasticity

A hypoplastic model of Φ̄ is given by (Fellin, 2013; Fuentes et al., 2012; Niemunis et al.,
2009)

Φ̄= ˆ̄Φ(ν̄,D̄,Z̄)= fs(ν̄,I
1
Z̄
)
{

a2D̄+ ˇ̄Ztr(ˇ̄ZD̄)+fd(ν̄)a(
ˇ̄Z+ ˇ̄Z

∗

)‖D̄‖
}

(2.6)

for rate-independent characteristics, with ˇ̄Z = Z̄/tr(Z̄), the versor of Z̄; ˇ̄Z
∗

= ˇ̄Z− I/3, the
deviator of ˇ̄Z; ‖D̄‖ =

√
trD̄2; and a a positive constant. The scalar functions fs and fd are

respectively the stiffness and density factors. The constant a is related to the stress state Z̄c
and frictional angle ϕc in the critical state, and can be determined empirically (Bauer and
Herle, 2000; Buscamera, 2014; Marcher et al., 2000). Equation (2.6) is used to account for the
rate-independent features of dry granular systems,with the benefits that (1) distinction between
loading andunloading is automatically accomplished, and (2) elastic/inelastic deformations need
not a priori be distinguished; information about yield surface and plastic potential is no longer
necessary.



1248 C. Fang

2.4. Closure relations

With these, the closure relations of an isochoric flow are given by

0= γ̄ν̄k− γ̄ν̄
[

α0(ν̄ − ν̄m)2+β0
( ν̄m
ν̄∞− ν̄

)2
(g1 ·g1)

](

1+
ϑT

ϑM

)

0= γ̄ν̄ε−f1 ˙̄ν −f2(trD̄)−f3(g4 ·g4)
0= γ̄ν̄H −f4 ˙̄ν −f5(trD̄)−f6(g5 ·g5)

0= γ̄ν̄s− γ̄ν̄
[

α0(ν̄ − ν̄m)2+β0
( ν̄m
ν̄∞− ν̄

)2
(g1 ·g1)

](

1+
ϑG

ϑM

)

0= ℓ(ϑMh̄+ϑGH)−2β0γ̄ν̄ϑMFc
( ν̄m
ν̄∞− ν̄

)2
g1

0=K+f7g4 0=L+f8g5

0= f̄ −
p̄

γ̄ν̄ℓ
+
2

ℓ

[

α0(ν̄ − ν̄m)+
β0ν̄
2
m

(ν̄∞− ν̄)3
(g1 ·g1)

]

Fc−
(

1−
ϑT

ϑM

) f1
γ̄ν̄ℓ

−
(

1−
ϑG

ϑM

) f4
γ̄ν̄ℓ
+ ζ ˙̄ν + δ(trD̄)

0= t̄− (−ν̄p+ ǫM ˙̄ν +λM trD̄)I−fs(ζ1I+ ζ2Z̄+ ζ3Z̄2)

+2β0γ̄ν̄Fc
( ν̄m
ν̄∞− ν̄

)2
g1⊗g1−2µ0γ̄2

( ν̄m
ν̄∞− ν̄

)8√

|I2
D
|D̄

0=R−
[

−
(ϑM

ϑT
−1
)

f2−
(ϑM

ϑT
− ϑ

G

ϑT

)

f5+ ǫ
T ˙̄ν +λT trD̄

]

I

−2µ0γ̄2(Fc−1)
( ν̄m
ν̄∞− ν̄

)8√

|I2
D̄
|D̄

(2.7)

where the Cayley-Hamilton theorem and the notations

c1 =ψ
T
f,I1
Z̄

c2 = ψ
T
f,I2
Z̄

c3 = ψ
T
f,I3
Z̄

(I1
Z̄
),Z̄ = I (I

2
Z̄
),Z̄ = I

1
Z̄
I− Z̄ (I3

Z̄
),Z̄ = I

3
Z̄
Z̄
−1

ζ1 = a
2(c1+ c2I

1
Z̄
)− c3a2I2Z̄I

3
Z̄

ζ3 = c3a
2(I3
Z̄
)2

ζ2 =(c1+ c2I
1
Z̄
)(I1
Z̄
)−1− c2(a2+(I1Z̄)

−2trZ̄2)+ c3I
3
Z̄
(3(I1

Z̄
)−2+a2I1

Z̄
)

(2.8)

have been used. For laminar flows, (2.7) reduce to those in previouswork (Fang, 2009). Two-fold
turbulent kinetic energies in (2.7)1 and (2.7)4 are determined once ν̄ is known.
The field equations of {p̄, v̄, ν̄,ϑT ,ϑG,Z̄} are obtained by substituting (2.7) into (1.1). Since

the system ismathematically likely well-posed, one has the chance to obtain the primitivemean
fields. In applying (2.7), the phenomenological parameters α0, β0, f1−8, fs, fd, a, ǫ

M, λM, ǫT ,
λT , µ0, ζ1−3, ν̄m and ν̄∞ need be prescribed. Since detailed information of them is insufficient,
numerical simulation is restricted to a parametric study.

3. Gravity-driven flow

3.1. Field equations and boundary conditions

Consider a fully-developed, isochoric, two-dimensional stationary avalanche down an incline,
as shown in Fig. 1. It is assumed that

v̄= [ū(y), v̄(y),0] ν̄ = ν̄(y) p̄ = p̄(y)

ϑT = ϑT(y) ϑG = ϑG(y) Z̄ij = Z̄ij(y)
(3.1)



On the turbulent boundary layer of a dry granular avalanche... 1249

with v̄/ū ∼ 0; u′ 6= 0, v′ 6= 0; {i,j} = (x,y), motivated by the assumptions that α,x ≪ α,y for
any quantity α in simple turbulent shear flows of Newtonian fluids (Batchelor, 1993). The flow

corresponds to the critical state defined as the state in which ˙̄ρ =0 and ˚̄Z= 0 (Ai et al., 2014;
Kirchner andTeufel, 2002). Since in the critical state fd is set to be unity, equations (1.1)4 and
(2.6) reduce to

0= fs
{

a2cD̄+
ˇ̄Ztr(ˇ̄ZD̄)+ac(

ˇ̄Z+ ˇ̄Z
∗

)‖D̄‖
}

(3.2)

with ac =
√

8/27sinϕc, and ϕc the critical friction angle (Kirchner and Teufel, 2002). Since fs
does not vanish generally, substituting (3.1) into (3.2) yields

0= ˇ̄Zxx
ˇ̄Zxym+ac

(

2 ˇ̄Zxx−
1

3

)

0= ˇ̄Zyy
ˇ̄Zxym+ac

(

2 ˇ̄Zyy −
1

3

)

0= a2cm+
ˇ̄Z2xym+2ac

ˇ̄Zxy

(3.3)

with m ≡ D̄xy/‖D̄‖ = D̄xy/|D̄xy|. A non-trivial solution to (3.3) is only that Z̄xx = Z̄yy
and Z̄xy = −m

√

8/3sinϕcZ̄yy. Thus, equation (1.1)4 is decoupled from other mean balance
equations. For further identification, a specific form of fs is given by (Bauer and Herle, 2000)

fs =
(1− ν̄s
1− ν̄

)m

m =1 (3.4)

with ν̄s theminimummean volume fraction; and unity power justified formost cases (Herle and
Gudehus, 2000; Marcher et al., 2000; Niemunis et al., 2009).

Fig. 1. Gravity-driven stationary avalanche down an incline and the coordinate

With these, the field equations are obtained

0=
d

dy

{

1− ν̄s
1− ν̄

(ζ2Z̄xy + ζ3Z̄
2
xy)+µ0γ̄

2Fc
( ν̄m
ν̄∞− ν̄

)8(dū

dy

)2
}

+ γ̄ν̄bsinθ

0=
dt̄yy
dy
− γ̄ν̄ cosθ

0=
d

dy

{

2β0γ̄ν̄Fc
ℓ

( ν̄m
ν̄∞− ν̄

)2dν̄

dy

}

+
1

ν̄ℓ

{

− t̄yy+
1− ν̄s
1− ν̄

l(ζ1+ ζ2Z̄yy + ζ3Z̄
2
yy)

−2α0γ̄ν̄2(ν̄ − ν̄m)Fc−2β0γ̄ν̄
( ν̄m
ν̄∞− ν̄

)2(dν̄

dy

)2 ν̄∞Fc
ν̄∞− ν̄

}

0= µ0γ̄
2(Fc−1)

( ν̄m
ν̄∞− ν̄

)8(dū

dy

)3
−f7

d2ϑT

dy2
−f3
(dϑT

dy

)2

0=−f8
d2ϑG

dy2
−f6
(dϑG

dy

)2

(3.5)

for ū(y), ν̄(y), t̄yy(y), ϑ
T(y) and ϑG(y), in which p̄(y) is replaced by t̄yy(y) for simplicity.



1250 C. Fang

Although solid boundaries have been demonstrated to act as energy sources and sinks of
the turbulent kinetic energy of the grains (e.g. Pudasaini and Hutter, 2007; Richman andMar-
ciniec, 1990), and the conventional no-slip condition for ū does not hold (Savage, 1993). The
field observations suggest that the no-slip condition can still be used as the first approximation
(Pudasaini andHutter, 2007; Rao andNott, 2008), with which a fixed value of ν̄ is assumed.As
motivated by the finite turbulent kinetic energies and dissipations on solid boundaries of New-
tonian fluids in turbulent boundary layer flows (Batchelor, 1993; Tsinober, 2009), finite two-fold
turbulent kinetic energies and dissipations on the plane are assumed through the prescriptions
of ϑT and ϑG, respectively. On the other hand, field observations suggest that the grains on the
free surface interlock with one another to form a kind of inelastic network. Due to a significant
density difference between the grains and the air, the shear stress is negligible, yielding vanishing
normal gradients of ū and ν̄. However, the air entrainment on the free surface provides a finite
normal stress, giving a rise to nonvanishing t̄yy and normal gradients of ϑ

T and ϑG (Pudasa-
ini and Hutter, 2007; Rao and Nott, 2008; Wang and Hutter, 2001). Thus, BCs for (3.5) are
postulated by:
— for y =0

ū =0 ν̄ = ν̄b ϑ
T = ϑTb ϑ

G = ϑGb (3.6)

— for y = L

dū

dy
=0

dν̄

dy
=0

dϑT

dy
= aT

dϑG

dy
= aG t̄yy = t̄b (3.7)

with the superscript b denoting the boundary values.

3.2. Nondimensionalisation

With the dimensionless parameters defined in Table 2, in which V0 is the characteristic
velocity of the flow, equations (3.5) are recast in dimensionless forms

0=
d

dỹ

{

1− ν̄mν̃s
1− ν̄mν̃

Ξ1+
F̃c

(ν̃∞− ν̃)8
(dũ

dỹ

)2
}

+S1ν̃ sinθ

0=
dπ̃

dỹ
+S2ν̃ cosθ

0=
d

dỹ

{

2ν̃F̃c
(ν̃∞− ν̃)2

dν̃

dỹ

}

+
1

ν̃

{

π̃+
1− ν̄mν̃s
1− ν̄mν̃

Ξ2−2ν̃2(ν̃ −1)−
2ν̃ν̃∞F̃c
(ν̃∞− ν̃)3

(dν̃

dỹ

)2
}

0=
χ1(F̃c−1)
(ν̃∞− ν̃)8

(dũ

dỹ

)3
− d
2ϑ̃T

dỹ2
−χ2

(dϑ̃T

dỹ

)2

0=−
d2ϑ̃G

dỹ2
−χ3
(dϑ̃G

dỹ

)2

(3.8)

for ũ(ỹ), ν̃(ỹ), π̃(ỹ), ϑ̃T(ỹ) and ϑ̃G(ỹ),with F̃c =1+ϑ̃T+ϑ̃G+(ϑ̃T)2+(ϑ̃G)2, andthedimensionless
BCs:
— for ỹ =0

ũ =0 ν̃ = ν̃b ϑ̃
T = ϑ̃Tb ϑ̃

G = ϑ̃Gb (3.9)

— for ỹ = L̃

dũ

dỹ
=0

dν̃

dỹ
=0

dϑ̃T

dỹ
= ãT

dϑ̃G

dỹ
= ãG π̃ = π̃b (3.10)



On the turbulent boundary layer of a dry granular avalanche... 1251

Table 2.Dimensionless parameters

ξ2 =
α0
β0
=
1

ℓ2
ỹ = ξy L̃ = ξL ν̃ =

ν̄

ν̄m
ũ =

ū

V0

ν̃∞ =
ν̄∞
ν̄m

ñub =
ν̄b
ν̄m

ϑ̃T = ϑ
T

ϑM
ϑ̃G =

ϑG

ϑM

Ξ1 =
ζ2Z̄xy+ ζ3Z̄

2
xy

µ0γ̄2V
2
0 ξ
2

Ξ2 =
ζ1+ ζ2Z̄yy + ζ3Z̄

2
yy

α0γ̄ν̄
3
m

S1 =
ν̄mb

µ0γ̄V
2
0 ξ
3

S2 =
b

α0ν̄2mξ
χ1 =

µ0γ̄V
3
0 ξ

f7ϑM
χ2 =

f3ϑ
M

f7
χ3 =

f6ϑ
M

f8

π̃ =
−t̄yy
α0γ̄ν̄3m

ν̃s =
ν̄s
ν̄m

π̃b =
−t̄b

α0γ̄ν̄3m
ϑ̃Tb =

ϑTb
ϑM

ϑ̃Gb =
ϑGb
ϑM

ãT =
aT
ϑMξ

ãG =
aG
ϑMξ

Equations (3.8)-(3.10) define a nonlinear BVP, with L̃ being the effect of flow thickness;
S2 the effect of gravity; S1 the influence of viscosity for fixed S2; Ξ1 and Ξ2 the hypoplastic
effect; χ1 the relative contribution between viscosity and turbulent kinetic energy flux; χ2 and
χ3 the relative significancesbetween two-fold turbulentkinetic energyfluxesanddissipations.For
implementation of numerical simulation, the values of ν̄m, ν̄b, ν̄∞ and ν̄s are given by ν̄b =0.51,
ν̄m =0.555, ν̄∞ =0.644, ν̄s =0.25 (thus, ν̃b =0.919, ν̃∞ =1.16, ν̃s =0.451), with fixed values
of ϑTb , ϑ

G
b and t̄b assumed as a first approximation (Bauer andHerle, 2000; Fang andWu, 2014;

Savage, 1993; Wang and Hutter, 1999).

Two-point nonlinear BVP (3.8)-(3.10) is solved numerically by using the iterative methods
with a successive under-relaxation scheme (Fang andWu, 2014;Wang andHutter, 1999;Wendt,
2009). So, sequences of the primitivemean fields are calculated at each iteration step, which are
incorporated into the next step, until demanded convergence is reached. Moreover, integrating
last Eq. (3.8) yields an analytical solution of ϑ̃G(ỹ) under a fixed value of χ3, viz., ϑ̃

G = ϑ̃Gb +
χ−13 ln

(

(1− ãGχ3(L̃− ỹ))/(1− ãGχ3L̃)
)

, indicating that ϑ̃G increases logarithmically from the
plane toward the free surface, and corresponding to the previousworks (Fang, 2009, 2016b; Fang
andWu, 2014).

3.3. Numerical results

Numerical tests show that only the relative magnitudes of the ν̃-, ũ-, ϑ̃T- and ϑG-profiles
are influenced by the variations in S1 and χ1−3, but the tendencies remain unchanged. Thus,
S1 = 0.02 and χ1 = χ2 = χ3 = 0.01 are used, with ϑ̃

T
b = ϑ̃

G
b = 0.1 for finite turbulent kinetic

energies and dissipations on the boundary (Pudasaini and Hutter, 2007; Rao and Nott, 2008).
Since Ξ1 and Ξ2 are of equal importance (Fang, 2009; Kirchner and Teufel, 2002), they are set
equal. In all figures, the normalized calculated values are displayed for comparison.

Figure2 illustrates theprofilesof ν̃, ũ, γ̄ν̃s, γ̄ν̃k, γ̄ν̃H and γ̄ν̃ε forvariations in L̃ = [10,15,20]
indicated by the arrows, with ãT = ãG = 0.1, Ξ1 = Ξ2 = 0.01, π̃b = 0.01 and S2 = 0.02. The
solid lines are the simulated results; the dashed lines are the laminar flow solutions from Fang
(2009); the dotted line are Newtonian fluid characteristics in a laminar flow. Increasing L̃ tends
to enlarge the difference in ν̃ between the free surface andplane, as shown inFig. 2a.This results
from the weight of the granular body: when the flow is thicker, larger compressive stress applies
on the grains in the thin layer immediately above the plane (the turbulent boundary layer, TBL),
with maximum shearing there, causing the grains to collide intensively with one another and
resulting in smaller values of ν̃. Above this thin layer, there exists a relatively thick layer (the
passive layer, PL), in which the grains form a kind of inelastic network and behave as a lump



1252 C. Fang

Fig. 2. Normalized profiles of ν̃, ũ, γ̄ν̃s, γ̄ν̃k, γ̄ν̃H, γ̄ν̃ε, with L̃ = [10,15,20] indicated by the arrows.
Dashed lines: laminar flow solutions; dotted line: laminar Newtonian flow

solid with nearly uniform ν̃ and ũ, as displayed in Figs. 2a and 2b. As L̃ increases, the TBL
becomes thinner with larger ũ-gradients at the interface between two layers. When compared
with laminar flow solutions, the ν̃- and ũ-profiles aremore convex, with larger amplitudes in the
PL. These are due to the influence of turbulent kinetic energy and dissipation, to be discussed
later.

Two-fold turbulent kinetic energies in Figs. 2c and 2d decrease from their maximum values
in the plane toward nearly vanishing values on the free surface exponentially. Similar tendencies
appear for the profiles of γ̄ν̃H and γ̄ν̃ε in Figs. 2e and 2f, except for finite values on the free
surface. As L̃ increases, γ̄ν̃s, γ̄ν̃k, γ̄ν̃H and γ̄ν̃ε decreasemore obviously. These correspond not
only to those of Newtonian fluids in turbulent boundary layer flows, but also are justified, for
turbulent kinetic energy and dissipation should assume maximum values in the regions where
shearing is maximum, and a larger turbulent kinetic energy induces larger turbulent dissipation
(Batchelor, 1993;Tsinober, 2009).Although inNewtonianfluidsanddrygranular avalanches the
turbulent kinetic energies and dissipations evolve in a similar manner, their vanishing values on
the free surface are identified forNewtonian fluids, while it is not so for dry granular avalanches.
These reflect the discrete nature of dry granular systems.

Although theTBLandPLcanbe identifiedby theprofiles of ν andu in laminar formulations
(e.g. Fang, 2009; Wang and Hutter, 1999), they are preferably recognized by the distributions
of γ̄ν̃s, γ̄ν̃k, γ̄ν̃H and γ̄ν̃ε. In the PL, the dominant grain-grain interaction is the long-term
one, causing the grains to form a kind of inelastic network to yield nearly vanishing γ̄ν̃s and
γ̄ν̃k, and finite γ̄ν̃H and γ̄ν̃ε. On the other hand, the grains in the TBL are dominated by the
short-term interaction, giving a rise to intensive turbulent fluctuation with significant γ̄ν̃s, γ̄ν̃k,
γ̄ν̃H and γ̄ν̃ε, resulting in larger ν̃ and ũ in thePL,when comparedwith laminar flow solutions.

Figure 3 illustrates the profiles of ν̃, ũ, γ̄ν̃s, γ̄ν̃k, γ̄ν̃H and γ̄ν̃ε for variations in
S2 = [0.01,0.035,0.07] indicated by the arrows, with ãT = ãG = 0.1, L̄ = 15, Ξ1 = Ξ2 = 0.01
and π̃b =0.01. Increasing S2 tends to enhance the gravitational effect, resulting inmore convex
ν̃- and ũ-profiles in Figs. 3a and 3b. This goes back to the influence of a larger grain weight. As



On the turbulent boundary layer of a dry granular avalanche... 1253

Fig. 3. Normalized profiles of ν̃, ũ, γ̄ν̃s, γ̄ν̃k, γ̄ν̃H, γ̄ν̃ε, with S2 = [0.01,0.035,0.07] indicated by the
arrows. Dashed lines: laminar flow solutions; dotted line: laminar Newtonian flow

S2 increases, most γ̄ν̃s, γ̄ν̃k, γ̄ν̃H and γ̄ν̃ε are confined within even thinner TBLs, resulting in
more energetic grain collisions there, as shown inFigs. 3c-3f. In addition, γ̄ν̃s, γ̄ν̃k, γ̄ν̃H and γ̄ν̃ε
evolve with similar tendencies described in Fig. 2. Due to the distributions of two-fold turbulent
kinetic energies, the ν̃- and ũ-profiles are more convex than their laminar-flow counterparts.

The influence of π̃b is summarized in Fig. 4, with ãT = ãG = 0.1, Ξ1 =Ξ2 = 0.01, L̃ = 15,
S2 = 0.02 and π̃b = [0.01,0.1,0.25] indicated by the arrows. Increasing π̃b is to apply larger
normal traction on the free surface, exciting the grains in the TBL to collide with one another
more vigorously. This reduces base friction, resulting inmore convex ν̃- and ū-profiles in Figs. 4a
and 4b, and thicker PLs. Figures 4c-4f show that as π̃b increases, γ̄ν̃s, γ̄ν̃k, γ̄ν̃H and γ̄ν̃ε are
confinedmostly in the TBL and decrease exponentially from the plane toward the free surface.
Theprofiles of ν̃ and ũaremoreconvex than their laminarflowcounterparts due to the influences
of γ̄ν̃s, γ̄ν̃k, γ̄ν̃H and γ̄ν̃ε.

Simulations for variations in Ξ1 and Ξ2 are summarized in Fig. 5, with Ξ1 = Ξ2 =
[0.01,0.05,0.1] indicated by the arrows; ãT = ãG = 0.1, π̃b = 0.01, L̃ = 15 and S2 = 0.02.
When Ξ1 and Ξ2 increase, the hypoplastic effect inside a granular RVE is enhanced, which we-
akens the frictional contact between the grains in the TBL, giving a rise to reduced γ̄ν̃s, γ̄ν̃k,
γ̄ν̃H and γ̄ν̃ε with thinner TBLs. With an enhanced hypoplastic effect, most γ̄ν̃s, γ̄ν̃k, γ̄ν̃H
and γ̄ν̃ε are confined within even thinner TBLs, which are equally recognized in the ν̃- and
ũ-profiles.

Calculations for variations in ãT and ãG are given inFig. 6,with ãT = ãG = [0.01,0.025,0.05]
indicated by the arrows; Ξ1 =Ξ2 =0.01, π̃b =0.01, S2 =0.02 and L̃ =15. Equal values of ãT
and ãG are used for simplicity. Increasing ãT and ãG allows more fluxes of γ̄ν̃s and γ̄ν̃k enter
into the granular body from the free surface, inducing more γ̄ν̃H and γ̄ν̃ε for counter-balance.
This yieldsmore thicker PLs, illustrated bymore convex profiles of ν̃, ũ, γ̄ν̃s, γ̄ν̃k, γ̄ν̃H and γ̄ν̃ε.
These correspond to field observations, for in the PL, the grains are interlocked and behave as
a lump solid, causing γ̄ν̃H and γ̄ν̃ε to overcome γ̄ν̃s and γ̄ν̃k.



1254 C. Fang

Fig. 4. Normalized profiles of ν̃, ũ, γ̄ν̃s, γ̄ν̃k, γ̄ν̃H, γ̄ν̃ε, with π̃b = [0.01,0.1,0.25] indicated by the
arrows. Dashed lines: laminar flow solutions; dotted line: laminar Newtonian flow

Fig. 5. Normalized profiles of ν̃, ũ, γ̄ν̃s, γ̄ν̃k, γ̄ν̃H, γ̄ν̃ε, with Ξ1 =Ξ2 = [0.01,0.5,0.1] indicated by the
arrows. Dashed lines: laminar flow solutions; dotted line: laminar Newtonian flow



On the turbulent boundary layer of a dry granular avalanche... 1255

Fig. 6. Normalized profiles of ν̃, ũ, γ̄ν̃s, γ̄ν̃k, γ̄ν̃H, γ̄ν̃ε, with ãT = ãG = [0.01,0.025,0.05] indicated by
the arrows. Dashed lines: laminar flow solutions; dotted line: laminar Newtonian flow

4. Conclusions and discussions

The derived equilibrium closure relations in Part I (Fang, 2016b) have been implemented to
obtain a zero-order closuremodel, which has beem applied to analyses of a stationary avalanche
down an incline; numerical simulations have been compared with laminar flow solutions.

While ν̃ and ũ evolve from their minimum values on the plane toward maximum values on
the free surface, γ̄ν̃s, γ̄ν̃k, γ̄ν̃H and γ̄ν̃ε distribute in a reverse manner, with most of them
confined within the TBL immediately above the plane. Above this, there exists a PL in which
the grains behave as a lump solid with nearly uniform ν̃ and ũ. In the TBL, the grains are
dominated by the short-term interaction, giving a rise to intensive turbulent fluctuation with
significant turbulent kinetic energy and dissipation, while those in the PL are dominated by the
long-term interaction to form a kind of inelastic network. Two layers are preferable recognized
from the turbulent kinetic energy and dissipation profiles.

The TBL and PL of a dry granular avalanche are similar to those of Newtonian fluids in
turbulent boundary layer flows. Although the turbulent kinetic energies and dissipations evolve
in a similar manner, their vanishing values on the free surface are found for Newtonian fluids,
while nearly vanishing turbulent kinetic energies and finite turbulent dissipations are obtained
for granular avalanches, resulted from their discrete nature and different dominant grain-grain
interactions in the TBL and PL. Discrepancies in the estimated ν̃- and ũ-profiles from the
laminar flow solutions suggest that the energy cascade induced by turbulent fluctuation needs
to be considered for better estimations on the characteristics of dry granular avalanches.

Acknowledgements

The author is indebted to theMinistry of Science and Technology, Taiwan, for the financial support

through the projectMOST103-2221-E-006-116-.The author also thanks the editor andProfessor Janusz

Badur for the detailed comments and suggestions which led to improvements.



1256 C. Fang

References

1. Ai J., Langston P.A., Yu H.-S., 2014, Discrete elementmodeling of material non-coaxiality in
simple shear flows, International Journal for Numerical and Analytical Methods in Geomechanics,
38, 6, 615-635

2. Batchelor G.K., 1993, The Theory of Homogeneous Turbulence, Cambridge University Press,
Cambridge NewYork

3. Bauer E., Herle I., 2000, Stationary states in hypoplasticity, [In:] Constitutive Modeling of
Granular Materials, Kolymbas D. (ed.), Springer Verlag, Berlin, 167-192

4. Buscamera G., 2014, Uniqueness and existence in plastic models for unsaturated soils, Acta
Geotechnica, 9, 313-327

5. Fang C., 2008, Modeling dry granular mass flows as elasto-visco-hypoplastic continua with mi-
crostructural effects. II. Numerical simulations of benchmark flowproblems,ActaMechanica, 197,
191-209

6. Fang C., 2009, Gravity-driven dry granular slow flows down an inclined moving plane: a compa-
rative study between two concepts of the evolution of porosity,Rheologica Acta, 48, 971-992

7. Fang C., 2016a, A k-ε turbulent closuremodel of an isothermal dry granular dense matter,Con-
tinuum Mechanics and Thermodynamics, 28, 4, 1048-1069

8. Fang C., 2016b, On the turbulent boundary layer of a dry granular avalanche down an incline.
I. Thermodynamic analysis, Journal of Theoretical and Applied Mechanics, 54, 3, 1051-1062

9. Fang C., Wu W., 2014, On the weak turbulent motions of an isothermal dry granular dense
flowwith incompressible grains. Part II. Complete closuremodels and numerical simulations,Acta
Geotechnica, 9, 739-752

10. Fellin W., 2013, Extension to barodesy to model void ratio and stress dependency of the K0
value,Acta Geotechnica, 8, 561-565

11. FuentesW.,Triantaftllidis T., LizcanoA., 2012,Hypoplasticmodel for sandswith loading
surface,Acta Geotechnica, 7, 177-192

12. Herle I., Gudehus G., 1999, Determination of parameters of a hypoplastic constitutive model
from properties of grain assemblies,Mechanics of Cohesive and Frictional Materials, 4, 461-485

13. Kirchner N., Teufel A., 2002, Thermodynamically consistent modeling of abrasive granular
materials. II: Thermodynamic equilibrium and applications to steady shear flows, Proceeding of
Royal Society London A, 458, 3053-3077

14. MarcherTh., VermeerP.A.,WolffersdorfP.-A., 2000,Hypoplastic and elastoplasticmo-
deling – a comparisonwith test data, [In:]Constitutive Modeling of GranularMaterials, Kolymbas
D. (ed.), Springer Verlag, Berlin, 353-371

15. Niemunis A., Grandas-Tavera C.E., Prada-Sarmiento L.F., 2009, Anisotropic visco-
hypoplasticity,Acta Geotechnica, 4, 293-314

16. Pudasaini S., Hutter K., 2007,Avalanche Dynamics, Springer Verlag, Berlin Heidelberg

17. Rao K.K., Nott P.R., 2008, Introduction to Granular Flows, CambridgeUniversity Press, Lon-
don NewYork

18. RichmanM.W.,Marciniec R.P., 1990,Gravity-driven granular flows of smooth, inelastic sphe-
res down bumpy inclines, Journal of Applied Mechanics, 57, 1036-1043

19. Savage S.B., 1993, Mechanics of granular flows, [In:] Continuum Mechanics in Environmental
Sciences and Geophysics, Hutter K. (ed.), Heidelberg, Springer, 467-522

20. Tsinober A., 2009,An Informal Conceptual Introduction to Turbulence, Springer, Heidelberg

21. Wang Y., Hutter K., 1999, A constitutive theory of fluid-saturated granular materials and its
application in gravitational flows,Rheologica Acta, 38, 214-223

22. Wang Y., Hutter K., 2001, Granular material theories revisited, [In:] Geomorphological Fluid
Mechanics, Balmforth N.J., ProvenzaleA. (eds.), Heidelberg, Springer, 79-107

23. Wendt J.F., 2009,Computational Fluid Dynamics: An Introduction, Springer, Berlin Heidelberg

Manuscript received January 29, 2015; accepted for print February 23, 2016