%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DEF's %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%your definitions


\documentclass[12pt]{article}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[reqno]{amsmath}
\usepackage{amsxtra}
\usepackage{eucal}
\usepackage{eucalbf}
\usepackage{epsf}
\usepackage[english]{babel}

%TCIDATA{OutputFilter=LATEX.DLL}
%TCIDATA{LastRevised=Tue May 30 10:09:07 2000}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}
%TCIDATA{CSTFile=article.cst}

\textheight=56truepc
\textwidth=150truemm
\hoffset=-10truemm
\voffset=-1.8truecm
\def\centerline#1{\hbox to \textwidth{\hss#1\hss}}
%\input{tcilatex}

\begin{document}

\title{A model for the flow of granular materials and its application to
initial/boundary value problems}
\author{D. Harris}
\date{\textit{${}$Department of Mathematics, UMIST, Sackville Street, Manchester,
M60 1QD}\\
${}$ \smallskip \textrm{david.harris\texttt{@umist.ac.uk }}}
\maketitle

\pagestyle{empty}

\thispagestyle{empty}

\bigskip {\small %your abstract
}

{\small A plasticity type model is here considered for the flow of a dry
granular material such as grain or sand. The physical and kinematic basis
for the model is briefly summarised and the equations governing the model
are presented in terms of the components of the deformation-rate, spin and
stress tensors. The equations comprise a set of six first order partial
differential equations of hyperbolic type for which there are five distinct
characteristic directions. An idealised application to a hopper is
considered for the flow in the vicinity of the upper free surface. A simple
analytic solution is given in which (a) the velocity field is linear in the
space co-ordinates and represents a dilatant or contractant shear, (b) two
possible stress fields are proposed, one linear and one exponential in
space, which satisfy the stress equilibrium equations, the yield condition
and the traction-free condition at the free surface, (c) the density is
homogeneous in space and exponential in time. Finally, a method is proposed
for defining an intrinsic time-scale for the deformation, which enables a
physically realistic density field to be obtained via a sequence of dilatant
and contractant shearing motions. Full advantage is taken of the hyperbolic
nature of the governing equations to allow the solution to have
discontinuities in the field variables, or their derivatives, in crossing
characteristic lines.} \bigskip

\section{Introduction}

The \textit{static} analysis of granular materials such as dry grain or sand
is generally agreed to be well expressed by the Coulomb yield condition
together with the stress equilibrium equations. The text by V.V. Sokolovskii 
\cite{Sokol} gives an excellent treatment both of the theory and application
to the solution of practical problems involving soil mass in a geotechnical
engineering context (for example embankments, soil retained by a wall, wedge
shaped soil mass). By way of contrast there is very little agreement
concerning the equations governing the \textit{flow} of such materials,
indeed it is an open question as to whether or not continuum mechanics is a
reasonable framework for such evidently discrete materials. It is the
contention of the present author that this question may be answered in the
affirmative and that the required equations are summarised in section 1.1
below. Their derivation and associated theory has been presented elsewhere,
see Harris \cite{Harris1}- \cite{Harris5}, and will not be repeated here.

The model presented is one of a class of plasticity models which assumes the
existence of a yield condition limiting the stress states that the material
may sustain but which replaces the standard flow rule obtained from a
plastic potential by a flow rule based upon a kinematic hypothesis that
prescribes the manner in which the material may flow. An early model of this
type was presented by Geniev \cite{Geniev} and was based upon the notion of
a single-slip direction with the material slipping on one of the two stress
characteristic directions. A model based upon simultaneous slipping on both
stress characteristic directions was proposed by de Josselin de Jong \cite
{Jong1}. This model, originally proposed for incompressible materials, also
incorporated the rotation-rate of the slip directions and was extended to
dilatant materials in \cite{Jong2}. In order to overcome an indeterminacy in
the above model, Spencer \cite{Spen1} interpreted the rotation-rate as that
due to the rotation of the in-plane greater principal stress direction. This
introduced the stress-rate into the model. No continuum model has been
universally accepted as describing a good approximation to granular flow,
each model having one or more disadvantages and this is the motivation for
the present work. The model presented here is a single slip model
incorporating a physical rotation of the slip systems and it has been
constructed in a manner designed to avoid the disadvantages suffered by
other models.

For the sake of understanding we shall briefly recapitulate the underlying
physical and kinematic basis of the model. Consider a granular material
confined, by virtue of a combination of containing boundaries and body
force, to occupy a region $\mathcal{R}$ with piecewise smooth boundary $%
\partial \mathcal{R}$. Contact between neighbouring grains is assumed to be
of finite duration and contact forces are assumed non-impulsive. Relative
motion of contacting grains is opposed by solid friction forces and there
may be inter-granular cohesive forces. The grains are treated as rigid
bodies subject to the kinematic constraint of non-overlap. The consequent
reaction forces, invoked to ensure non-violation of this constraint,
dominate the relative (and indeed absolute) motion of the grains.The
existence of relative motion between grains requires the frictional forces
opposing such relative motion to be surmounted and also the re-arrangement
of the grain configuration to prevent grain overlap. This will involve
dilatation for a densely packed configuration as grain re-arrangement makes
room for a relative motion of the grains. In the case of a loosely packed
configuration, for a flow regime that does not entail grain separation,
relative motion of the grains will be accompanied by contraction as the
grain structure collapses and grains move into the voids.

The \textit{local fabric }of the granular material is defined to be the
network, in the vicinity of a point in the granular material, of contact
points between grains, normals to the grain surfaces at contact points and
also the directions of slip tangential to the grain surfaces at contact. The
discrete local fabric of the real granular material is replaced, in the
model, by a \textit{continuum local fabric }at each point $P$ which
comprises a number (possibly infinite) of pairs of unit vector directions $%
\mathbf{t}_{\beta }$ and $\mathbf{t}_{\gamma }$ which we shall term the
macroscopic and microscopic slip directions, respectively. At most one pair
is allowed to be \textit{active }at a given point and time (if no pair is
active the material is rigid in the neighbourhood of the point) and the flow
is assumed to be the resultant of the following three components

\begin{description}
\item[(a)]  a simple shear in the macroscopic slip direction $\mathbf{t}%
_{\beta }$,

\item[(b)]  a dilatation/contraction in the normal direction $\mathbf{n}%
_{\beta }$ to $\mathbf{t}_{\beta }$ such that the resultant of the simple
shear and the dilatation is a dilatant shear in the microscopic slip
direction $\mathbf{t}_{\gamma }$.

\item[(c)]  a local rigid spin of the macroscopic slip direction $\mathbf{t}%
_{\beta }$.
\end{description}

We shall refrain from expressing this kinematic hypothesis mathematically as
this has been done elsewhere, see Harris \cite{Harris1}-\cite{Harris5} but
merely content ourselves with writing down three of the four resulting
second order tensor equations which the kinematic hypothesis gives rise to.

\subsection{The equations governing the model}

Consider a planar deformation in the $Ox_{1}x_{2}$ plane. Let the velocity
vector $\mathbf{v}$ have components $v_{i}$, the stress tensor $\mathbf{%
\sigma }$ have components $\sigma _{ij}$ and let the bulk density be $\rho $%
. Let the deformation-rate tensor (i.e. the symmetric part of the velocity
gradient tensor) be denoted by $\mathbf{d}$ with components $d_{ij},$ $1\leq 
$ $i,j\leq 2.$ Let the spin tensor (i.e. the anti-symmetric part of the
velocity gradient tensor) be denoted by $\mathbf{s}$ with components $s_{ij}$%
. Also, at each point $P$ of $\mathcal{R}$, let the active macroscopic and
microscopic slip directions make angles $\theta =\theta \left(
x_{1},x_{2},t\right) $ and $\chi =\chi \left( x_{1},x_{2},t\right) $,
respectively, with the $x_{1}$-axis, 
\begin{equation*}
\begin{array}{c}
\mathbf{t}_{\beta }=\left( \cos \theta ,\sin \theta \right) , \\ 
\mathbf{t}_{\gamma }=\left( \cos \chi ,\sin \chi \right) .
\end{array}
\end{equation*}
It should be noted that 
\begin{equation*}
\chi =\theta +\nu
\end{equation*}
where $\nu $ denotes the \textit{angle of dilatancy}.

The first equation governing the model states that the dilatation-rate is
proportional to the shear-rate, namely 
\begin{equation}
d_{11}+d_{22}=\sin \nu \left[ -\left( d_{11}-d_{22}\right) \sin \left(
2\theta +\nu \right) +2d_{12}\cos \left( 2\theta +\nu \right) \right]
\label{ke1}
\end{equation}
where the constant of proportionality is $\sin \nu .$ The second equation
relates the angle that the greater in-plane principal value of the
deformation-rate makes with the $x_{1}$-axis to the direction of the
macroscopic slip direction, namely 
\begin{equation}
\left( d_{11}-d_{22}\right) \cos \left( 2\theta +\nu \right) +2d_{12}\sin
\left( 2\theta +\nu \right) =0.  \label{ke2}
\end{equation}
The third equation relates the rotation-rate of $\mathbf{t}_{\beta }$ to the
spin tensor, 
\begin{equation}
2\left( s_{12}+\dot{\theta}\right) =\cos \nu \left[ -\left(
d_{11}-d_{22}\right) \sin \left( 2\theta +\nu \right) +2d_{12}\cos \left(
2\theta +\nu \right) \right] ,  \label{ke3}
\end{equation}
where the superposed dot denotes the material derivative. This completes the
equations governing the flow. The stress field is governed by the Cauchy
equations of motion 
\begin{equation}
\rho \frac{\partial v_{1}}{\partial t}+\rho v_{1}\frac{\partial v_{1}}{%
\partial x_{1}}+\rho v_{2}\frac{\partial v_{1}}{\partial x_{2}}-\frac{%
\partial \sigma _{11}}{\partial x_{1}}-\frac{\partial \sigma _{12}}{\partial
x_{2}}-\rho F_{1}=0,  \label{em1}
\end{equation}
\begin{equation}
\rho \frac{\partial v_{2}}{\partial t}+\rho v_{1}\frac{\partial v_{2}}{%
\partial x_{1}}+\rho v_{2}\frac{\partial v_{2}}{\partial x_{2}}-\frac{%
\partial \sigma _{12}}{\partial x_{1}}-\frac{\partial \sigma _{22}}{\partial
x_{2}}-\rho F_{2}=0,  \label{em2}
\end{equation}
together with the yield condition given by the inequality (\ref{yield1})
below, which is closely related to the Coulomb yield condition. Let $\tau
_{\chi }$, $\sigma _{\chi }$ denote the normal and tangential components of
the traction vector at a point $P$ across the line element along the
microscopic slip direction $\mathbf{t}_{\gamma }$ then we shall assume that
the stress at $P$ satisfies the inequality 
\begin{equation}
\left| \tau _{\chi }\right| \leq k-\left| \sigma _{\chi }\right| \tan \phi
_{\nu },  \label{yield1}
\end{equation}
where $\phi _{\nu }$ denotes the \textit{angle of mobilised internal friction%
} (and takes into account both the resistance to the grain on grain friction
and the resistance to separating the particles during dilatation) and $k$
denotes the \textit{cohesion}.\textit{\ }If $\phi _{0}$ denotes the angle of
internal friction when the state of the material and the stress admits an
isochoric flow then 
\begin{equation*}
\phi _{\nu }=\phi _{0}+\nu .
\end{equation*}
The bulk density of the material $\rho $ is determined by the continuity
equation 
\begin{equation}
\rho \frac{\partial v_{1}}{\partial x_{1}}+\rho \frac{\partial v_{2}}{%
\partial x_{2}}+\frac{\partial \rho }{\partial t}+v_{1}\frac{\partial \rho }{%
\partial x_{1}}+v_{2}\frac{\partial \rho }{\partial x_{2}}=0.  \label{ec}
\end{equation}

\subsection{The characteristic directions}

The system of first order partial differential equations presented in the
previous subsection are hyperbolic in nature, for a proof of this see Harris 
\cite{Harris1}-\cite{Harris3} and \cite{Harris5}. There are five distinct
characteristic directions defined at each point $P$ of the granular material.

\begin{description}
\item[1]  The kinematic equations (\ref{ke1}) and (\ref{ke2}) give rise to
the following pair of velocity characteristic directions, the $\alpha $-
characteristic direction $m_{\alpha }$, 
\begin{equation}
m_{a}=\frac{dx_{2}}{dx_{1}}=\tan \left( \chi +\frac{1}{2}\pi \right) ,
\label{dirb}
\end{equation}
i.e. the direction $\mathbf{n}_{\gamma },$ perpendicular to the microscopic
slip direction, and the $\beta $-characteristic direction $m_{\beta }$, 
\begin{equation}
m_{\beta }=\frac{dx_{2}}{dx_{1}}=\tan \theta ,  \label{dira}
\end{equation}
i.e. the macroscopic slip direction $\mathbf{t}_{\beta }$; the angle between
the $\alpha $- and $\beta $- characteristic directions is $\frac{1}{2}\pi
+\nu .$

\item[2]  The equations of motion (\ref{em1})-(\ref{yield1}) give rise to
the following pair of stress characteristic directions, the $\gamma $%
-characteristic direction $m_{\gamma }$, 
\begin{equation}
m_{\gamma }=\frac{dx_{2}}{dx_{1}}=\tan \chi ,  \label{dirc}
\end{equation}
i.e. the microscopic slip direction $\mathbf{t}_{\gamma }$ and the $\delta $%
-characteristic direction $m_{\delta }$, 
\begin{equation}
m_{\delta }=\frac{dx_{2}}{dx_{1}}=\tan \left( \chi \pm \phi +\frac{1}{2}\pi
\right) ,  \label{dird}
\end{equation}
i.e. an angle $\pm \phi $ to the perpendicular to the microscopic slip
direction; the angle between the $\gamma $ and $\delta $ directions is $%
\frac{1}{2}\pi \pm \phi .$ The ambiguity of sign is removed by the condition
of non-negative work-rate.

\item[3]  The kinematic equation (\ref{ke3}) and continuity equation (\ref
{ec}) both give rise to the same characteristic direction, which is
determined by the pair of ordinary differential equations equations 
\begin{equation}
\frac{dx_{1}}{dt}=v_{1},\qquad \frac{dx_{2}}{dt}=v_{2}  \label{dire1}
\end{equation}
together with the fact that the characteristic is embedded in the flow.
These are parametric equations for the streamline characteristic direction, $%
m_{\varepsilon },$ in the case of time dependent flow and the equation 
\begin{equation}
m_{\varepsilon }=\frac{v_{2}}{v_{1}}=\tan \chi  \label{dire2}
\end{equation}
in the case of time independent flow.
\end{description}

\section{Application: hopper flow near the free surface}

In this section we shall consider an application of the model to a case in
which there is a free surface, for example the emptying of a hopper,
represented schematically in Figure \ref{fig1}. Granular material is
confined between the vertical walls $AC$ and $BD.$ The upper surface of the
granular material is initially the straight line $EFG,$ where $F$ denotes
the mid-point of $EG$. The arrangement whereby the material is let out of
the hopper, somewhere below $CD$, is not shown and we will not attempt to
model the geometry or the flow in the exit region. Instead we shall
concentrate attention on the region in the vicinity of $EFG$ and will
present a simple analytic solution of the governing equations which
represents a simplified flow in the upper portion of the hopper sufficiently
far from the exit. As the hopper empties the material flows downward and the
free surface descends, changing shape as the flow continues until it reaches
a steady state. Let $H$ denote the mid-point of the straight line segment
joining $C$ and $D$. A corner often appears on the free surface on the axis
of symmetry $FH$ of the hopper and this is shown as $E^{\prime }F^{\prime
}G^{\prime }$ in Figure 1. Choose Cartesian co-ordinate axes $Ox_{1}x_{2}$
(where the origin $O$ coincides with the mid-point $F$ in the initial
configuration) with the $x_{1}$-axis horizontal and to the right, the $x_{2}$%
-axis vertically upward. The material properties of the granular material
are assumed to be such that the angle of internal friction $\phi _{0}$ is
constant and that the cohesion $c$ is zero. In this section the angle of
dilatancy $\nu $ will be assumed everywhere constant. We shall assume a
solution in which the active microscopic direction $\mathbf{t}_{\gamma }$ in 
$CEFGDH$ is everywhere parallel to the hopper walls $\overrightarrow{AC}$
and $\overrightarrow{BD},$ i.e the angle $\chi $ is everywhere constant and
equal to $-\frac{1}{2}\pi $ relative to the $Ox_{1}x_{2}$-axes, then 
\begin{equation*}
\theta =-\frac{1}{2}\pi -v.
\end{equation*}
\FRAME{ftbpFU}{1.8031in}{2.6541in}{0pt}{\Qcb{Hopper Flow}}{\Qlb{fig1}}{%
petepic1.eps}{\special{language "Scientific Word";type
"GRAPHIC";maintain-aspect-ratio TRUE;display "USEDEF";valid_file "F";width
1.8031in;height 2.6541in;depth 0pt;original-width 2.7112in;original-height
4.0101in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename
'Petepic1.eps';file-properties "XNPEU";}}The $\alpha $-, $\beta $-,$\gamma $%
- and $\delta $-chracteristic directions are shown in Figure 1 as the
directions$\ \overrightarrow{PL}$, $\overrightarrow{PI}$, $\overrightarrow{PJ%
}$ and $\overrightarrow{PK}$, respectively. If the material in $CEFGDH$ is
nowhere in a state of yield then, as material below $CD$ exits the hopper,
the material in $CEFGD$ will move vertically down the hopper as a rigid
body. If the material in $CEFGD\;$is in a state of yield then it should be
noted that the material placed infinitesimally close either side of the
centre line $FH$ will have no tendency to shear past each other and so the
material on the vertical centre line is not in a state of yield. Thus the
solutions in the two regions in yield $CEFH$, $DGFH$ may be treated
separately and the two solutions patched together along $FH$. It should be
noted that these are not the only two possibilities. Any member of the
family of $\gamma $-characteristics may separate a rigid region from a
yielding region and so the region $CEFGDH$ may be divided up into an
arbitrary number of alternating rigid and deforming regions separated by $%
\gamma $-characteristics. This ambiguity may be removed by consideration of
the exit conditions below $CD$.

If the material interior to the regions $CEFH$, $DGFH$ is everywhere in a
state of yield then we shall model the flow in $DGFH$ by a simple linear
velocity field 
\begin{equation}
v_{1}=0,  \label{vf1}
\end{equation}
\begin{equation}
v_{2}=-k_{\ast }\left( x_{1}-x_{2}\tan \nu \right) +A,  \label{vf2}
\end{equation}
which expresses a dilatant shear of the material vertically downwards
through the hopper. Here $k_{\ast }=k_{\ast }\left( t\right) >0$ is a
constant characterising the magnitude of the shear. The arbitrary constant $%
A=A\left( t\right) $ is determined by the conditions at the exit of the
hopper. This velocity field satisfies the flow rule (\ref{ke1})-(\ref{ke3}).
The flow in $CEFH$ is assumed to be the mirror image of the flow in $DGFH.$
We turn now to the continuum equations of motion in Cartesian form relative
to the $Ox_{1}x_{2}$-axes, then taking into account the above velocity
field, the equations (\ref{em1}) and (\ref{em2}) become 
\begin{equation*}
\frac{\partial \sigma _{11}}{\partial x_{1}}+\frac{\partial \sigma _{12}}{%
\partial x_{2}}=0,
\end{equation*}
\begin{equation*}
\rho \frac{\partial v_{2}}{\partial t}+v_{2}\frac{\partial v_{2}}{\partial
x_{2}}-\frac{\partial \sigma _{12}}{\partial x_{1}}-\frac{\partial \sigma
_{22}}{\partial x_{2}}+\rho g=0.
\end{equation*}
The convection terms are $O\left( k_{\ast }^{2}\right) $ and $O\left(
k_{\ast }A\right) $ and for sufficiently small $k_{\ast }$, $A$ and $\dot{k}%
_{\ast }$, $\dot{A}$ it is usual to neglect the inertia terms in comparison
with the equilibrium stresses. This practice is almost universal in
plasticity theory and we adopt this assumption here. In order that the rate
of working of the stresses be positive for the velocity field (\ref{vf1})
and (\ref{vf2}) it follows that, relative to the $Ox_{1}x_{2}$-axes, $\tau
_{\chi }<0$ and so the yield condition becomes 
\begin{equation*}
\tau _{\chi }=\sigma _{\chi }\tan \left( \phi +\nu \right) .
\end{equation*}
Relative to $Ox_{1}x_{2}$-axes, the yield condition takes the simple form 
\begin{equation*}
\sigma _{12}\cos \left( \phi +\nu \right) -\sigma _{11}\sin \left( \phi +\nu
\right) =0
\end{equation*}
since the $x_{2}$-axis is aligned along the $\gamma $-direction.

Now, the upper surface $EG$ is traction free, 
\begin{equation*}
\mathbf{t}_{EG}=\mathbf{\sigma n}_{EG}=\mathbf{0,}
\end{equation*}
where $\mathbf{t}_{EG}$ denotes the traction vector across $EG$, $\sigma $
denotes the stress tensor and $\mathbf{n}=\left( -\sin \eta ,\cos \eta
\right) ,$ where $\eta $ is the angle the upper surface makes with the $x_{1}
$-axis. Thus, 
\begin{equation*}
\begin{array}{ccc}
\sigma _{12}\cos \eta  & = & \sigma _{11}\sin \eta  \\ 
\sigma _{22}\cos \eta  & = & \sigma _{12}\sin \eta 
\end{array}
\end{equation*}
and, in the case that $\eta \neq \phi +\nu $, this is consistent with the
free surface $EFG$ being in a state of yield only if 
\begin{equation*}
\sigma _{11}=\sigma _{12}=\sigma _{22}=0.
\end{equation*}
In the case that $\eta =\phi +\nu $, there are infinitely many states of
stress that satisfy the traction-free boundary condition and the yield
condition. If $FG$ is a straight line inclined at an angle $\phi +\nu $ to
the $x_{1}$-axis, another slip direction inclined at angle $\nu $ to $FG$
becomes active, and a secondary flow of material from the vicinity of the
hopper walls towards the centre $F^{\prime }$ is activated. We now consider
two forms for the solution of the stress equilibrium equations.

\begin{description}
\item[Case (a)]  Linear stress field.
\end{description}

Suppose that 
\begin{equation*}
\begin{array}{c}
\sigma _{11}=Ax_{1}+Bx_{2}+C, \\ 
\sigma _{12}=Dx_{1}+Ex_{2}+F, \\ 
\sigma _{22}=Gx_{1}+Hx_{2}+I.
\end{array}
\end{equation*}
Substituting into the yield condition gives 
\begin{equation*}
D=A\tan \left( \phi +\nu \right) ,\qquad E=B\tan \left( \phi +\nu \right)
,\qquad F=C\tan \left( \phi +\nu \right)
\end{equation*}
Substituting into the equilibrium equations gives 
\begin{equation*}
A=-B\tan \left( \phi +\nu \right) ,
\end{equation*}
\begin{equation*}
H=\rho g-A\tan \left( \phi +\nu \right) .
\end{equation*}
Substituting into the traction free condition gives $\eta =\phi +\nu $ and 
\begin{equation*}
G=-B\tan ^{3}\left( \phi +\nu \right)
\end{equation*}
\begin{equation*}
I=C\tan ^{2}\left( \phi +\nu \right)
\end{equation*}
Hence, the stress field is 
\begin{equation*}
\sigma _{11}=B\left[ x_{2}-\tan \left( \phi +\nu \right) x_{1}\right] +C,
\end{equation*}
\begin{equation*}
\sigma _{12}=\left\{ B\left[ x_{2}-\tan \left( \phi +\nu \right) x_{1}\right]
+C\right\} \tan \left( \phi +\nu \right) ,
\end{equation*}
\begin{equation*}
\sigma _{22}=\left\{ B\left[ x_{2}-\tan \left( \phi +\nu \right) x_{1}\right]
+C\right\} \tan ^{2}\left( \phi +\nu \right) +\rho g\left[ x_{2}-\tan \left(
\phi +\nu \right) x_{1}\right] .
\end{equation*}
Let the equation of the free surface $F^{\prime }G^{\prime }$, as it
descends, be 
\begin{equation*}
x_{2}=x_{1}\tan \left( \phi +\nu \right) -h
\end{equation*}
then $F^{\prime }G^{\prime }$ is stress-free if $C=Bh$.

\begin{description}
\item[Case (b)]  Exponential stress field.
\end{description}

Consider a stress field of the form 
\begin{equation*}
\begin{array}{c}
\sigma _{11}=A+B\exp \left( k_{1}x_{1}+k_{2}x_{2}\right) , \\ 
\sigma _{12}=C+D\exp \left( k_{1}x_{1}+k_{2}x_{2}\right) , \\ 
\sigma _{22}=E+F\exp \left( k_{1}x_{1}+k_{2}x_{2}\right) +\rho g\left[
x_{2}-\tan \left( \phi +\nu \right) x_{1}\right] ,
\end{array}
\end{equation*}
where $A,B,C,D,E,F,k_{1},k_{2}$ are constants to be determined. Substituting
into the yield condition gives 
\begin{equation*}
C=A\tan \left( \phi +\nu \right) ,\qquad D=B\tan \left( \phi +\nu \right) .
\end{equation*}
Substituting into the equilibrium equations gives 
\begin{equation*}
k_{1}=-k_{2}\tan \left( \phi +\nu \right) ,
\end{equation*}
\begin{equation*}
F=B\tan ^{2}\left( \phi +\nu \right) .
\end{equation*}
Substituting into the traction free condition gives $\eta =\phi +\nu $ and 
\begin{equation*}
E=A\tan ^{2}\left( \phi +\nu \right) .
\end{equation*}
Hence, 
\begin{equation*}
\sigma _{11}=A+B\exp \left\{ k\left[ x_{2}-\tan \left( \phi +\nu \right)
x_{1}\right] \right\} ,
\end{equation*}
\begin{equation*}
\sigma _{12}=\left[ A+B\exp \left\{ k\left[ x_{2}-\tan \left( \phi +\nu
\right) x_{1}\right] \right\} \right] \tan \left( \phi +\nu \right) ,
\end{equation*}
\begin{equation*}
\sigma _{22}=\left[ A+B\exp \left\{ k\left[ x_{2}-\tan \left( \phi +\nu
\right) x_{1}\right] \right\} \right] \tan ^{2}\left( \phi +\nu \right)
+\rho g\left[ x_{2}-\tan \left( \phi +\nu \right) x_{1}\right]
\end{equation*}
where the constants $A$, $B$, $k$ are arbitrary. Note that for the choice 
\begin{equation*}
B=-A
\end{equation*}
the upper free boundary is stress free. The remaining constants $A$ and $k$
are determined by the conditions at the exit to the hopper.

Finally, we consider the density field. Relative to axes $Ox_{1}^{\prime
}x_{2}^{\prime }$ in which the $x_{1}^{\prime }$-axis is directed along the $%
\alpha $-characteristic line the velocity field becomes 
\begin{equation*}
v_{1}^{\prime }=k_{\ast }x_{2}^{\prime },
\end{equation*}
\begin{equation*}
v_{2}^{\prime }=k_{\ast }x_{2}^{\prime }\tan \nu
\end{equation*}
then the continuity equation becomes 
\begin{equation*}
\frac{\partial \rho }{\partial t}+\left( k_{\ast }x_{2}^{\prime }\right) 
\frac{\partial \rho }{\partial x_{1}^{\prime }}+\left( k_{\ast
}x_{2}^{\prime }\tan \nu \right) \frac{\partial \rho }{\partial
x_{2}^{\prime }}+k_{\ast }\rho \tan \nu =0
\end{equation*}
which is equivalent to the set of ordinary differential equations 
\begin{equation*}
\frac{dt}{1}=\frac{dx_{1}^{\prime }}{k_{\ast }x_{2}^{\prime }}=\frac{%
dx_{2}^{\prime }}{k_{\ast }x_{2}^{\prime }\tan \nu }=\frac{d\rho }{-k_{\ast
}\rho \tan \nu }.
\end{equation*}
Hence 
\begin{equation*}
\begin{array}{c}
\rho =\rho _{0}\exp \left( -k_{\ast }t\tan \nu \right)
\end{array}
\end{equation*}
along the characteristic lines 
\begin{equation*}
x_{2}^{\prime }-x_{1}^{\prime }\tan \nu =c,
\end{equation*}
which become 
\begin{equation*}
x_{1}=c
\end{equation*}
relative to the $Ox_{1}x_{2}$-axes. Hence, for a dilating material (i.e. $%
\nu >0$) $\rho \rightarrow 0$ as $t\rightarrow \infty $, whereas for a
compacting material (i.e. $\nu <0$) $\rho \rightarrow \infty $ as $%
t\rightarrow \infty $.

\section{An intrinsic time scale}

The bulk density $\rho $ of the granular material is bounded above by the
density $\rho _{g}$ of the grains which make up the material and below by
the bulk density $\rho _{s}$ at separation packing. In any given flow it is
unlikely that either of these bounds will be attained. The actual maximum
and minimum bulk densities attained during the flow will depend upon several
factors, for example, the initial state of the material, grain size,
confining boundaries and the confining pressure. Let $\rho _{\max }$, $\rho
_{\min },$ respectively, denote these quantities, then 
\begin{equation*}
\rho _{s}\leq \rho _{\min }\leq \rho \leq \rho _{\max }\leq \rho _{g}.
\end{equation*}
At time $t=0$ let $\rho =\rho _{\min }$ then during a deformation governed
by equations (\ref{vf1}) and (\ref{vf2}), the material must initially
densify, i.e. $\nu <0$. Suppose the material continues to densify at a
constant angle of dilatancy $\nu _{c}$ until $\rho =\rho _{\max }$ at time $%
\tau _{c}$, say, then 
\begin{equation*}
\rho _{\max }=\rho _{\min }\exp \left( -k_{\ast }\tau _{c}\tan \nu \right) 
\end{equation*}
i.e 
\begin{equation*}
\tau _{c}=\frac{1}{k_{\ast }\tan \nu _{c}}\ln \left( \frac{\rho _{\max }}{%
\rho _{\min }}\right) .
\end{equation*}
We shall call $\tau _{c}$ the \textit{intrinsic contraction time. }%
Similarly, a dilatant deformation at constant angle of dilatancy $\nu _{d}$ 
\textit{\ }from $\rho =\rho _{\max }$ to $\rho =\rho _{\min }$ gives rise to
an \textit{intrinsic} \textit{dilatation time,} 
\begin{equation*}
\tau _{d}=\frac{1}{k_{\ast }\tan \nu _{d}}\ln \left( \frac{\rho _{\min }}{%
\rho _{\max }}\right) .
\end{equation*}
\textit{\ }

\FRAME{ftbpFU}{4.0491in}{2.5071in}{0pt}{\Qcb{Intrinsic time scale}}{\Qlb{fig2%
}}{petepic2.eps}{\special{language "Scientific Word";type
"GRAPHIC";maintain-aspect-ratio TRUE;display "USEDEF";valid_file "F";width
4.0491in;height 2.5071in;depth 0pt;original-width 3.9998in;original-height
2.4664in;cropleft "0";croptop "1";cropright "1";cropbottom "0";filename
'Petepic2.eps';file-properties "XNPEU";}}We now propose that the velocity
field in the real granular material may be approximated by a sequence of
flows of the type considered in the previous section and alternating between
contraction and dilatation. For definiteness let the initial density $\rho
=\rho _{0}$ be such that the initial flow is contractant with angle of
dilatation $\nu _{c}<0$, followed by a dilatant flow with angle of dilatancy 
$\nu _{d}>0.$ Then the material oscillates between two states, with the
change of state occurring at intervals of time $\tau _{c}$ and $\tau _{d}$.
In Figure \ref{fig2}, $AB$ denotes the value of $\rho _{\min },$ $CD$ the
value of $\rho _{\max }$, $OE$ the value of $\tau _{c}$ and $EF$ the value
of $\tau _{d}$. The deformation is clearly non-smooth, with discontinuities
in both the velocity and stress fields at times when the material flips from
a dilatant to a contractant state and vice versa. Note that, in particular,
the yield condition differs in the two states, with the yield strength
greater in a dilatant state than in a contractant state i.e. 
\begin{equation*}
\left| \sigma _{12}\right| \leq \left| \sigma _{11}\right| \tan \left( \phi
+\nu _{d}\right) ,
\end{equation*}
\begin{equation*}
\left| \sigma _{12}\right| \leq \left| \sigma _{11}\right| \tan \left( \phi
+\nu _{c}\right) .
\end{equation*}

It is a well-known experimental observation that flow in a hopper may suffer
from various discontinuities, for example the pressure at the walls is known
to be discontinuous in time and the flow may either ``pulsate'' or exhibit
``stick-slip''. Such phenomena may be explained in terms of the above
mechanism.

\section{Conclusions}

We have stated the equations of a novel physically based model for the flow
of a granular material and have presented a simple analytic solution based
upon a linear velocity field and a choice of either a linear or exponential
(in space) stress field, together with an exponential (in time) density
field. Despite the simplicity of the solution it is able to satisfy the
yield condition, stress equilibrium equations, and conditions are given when
it may also satisfy a traction-free or stress-free boundary. The solution
has a natural application to hopper flow in the vicinity of the upper free
surface. The non-physical asymptotic behaviour of the density field, i.e.
the fact that dilatation/contraction cannot continue indefinitely suggests
that it is possible to alternate the two possible states and consider a
sequence of initial value problems in which the material successively
dilates and contracts, the final value in a given state giving rise to the
initial conditions of the succeeding state. This idea takes advantage of the
discrete nature of a granular material by relaxing the conditions of
differentiability and continuity usually present and is possible because of
the hyperbolic nature of the governing equations which allow discontinuities
in field variables and their gradients in crossing the characteristics. In
future, it may prove possible to refine the notion of intrinsic time
introduced here and, combining this concept with a typical velocity
magnitude, obtain an intrinsic microscopic length scale associated with the
granular material. It may eventually prove possible for this microscopic
length scale to play the role of a particle length scale. This idea is
speculative but may be worth pursuing. On a more practical level, the next
stage in the development of the model is to overcome the restrictions of a
simple analytic solution and to solve the equations numerically to obtain
numerical approximations to more realistic problems. This work is now in
progress.

\begin{thebibliography}{99}
\bibitem{Geniev}  Geniev, G.A., 1958, Problems of the dynamics of a granular
medium (in Russian), \textit{Akad. Stroit Archit. SSSR, Moscow.}

\bibitem{Harris1}  Harris, D., 1995, A unified formulation for plasticity
models of granular and other materials, \textit{Proc. R. Soc. Lond. A, }%
\textbf{450}, 37-49.

\bibitem{Harris2}  Harris, D., 1997, Modelling mathematically the flow of
granular materials,\textit{\ IUTAM Symposium on Mechanics of Granular and
Porous Materials, }Cambridge July 1996.

\bibitem{Harris3}  Harris, D., 1994, TAM Report No.32, Dept. Mathematics,
UMIST, 1-21.

\bibitem{Harris4}  Harris, D., 1999, Ill- and Well-Posed Models of Granular
Flow, accepted for publication in \textit{Acta Mechanica.}

\bibitem{Harris5}  Harris, D., Characteristic relations for a model for the
flow of granular materials, submitted for publication.

\bibitem{Jong1}  de Josselin de Jong, G., 1959, Statics and Kinematics of
the Failable Zone of a Granular Material\textit{, Thesis, }Uitgeverij, Delft.

\bibitem{Jong2}  de Josselin de Jong, G., 1977, Mathematical elaboration of
the double-sliding, free-rotating, \textit{Archs. Mech., }\textbf{29, }%
561-591.

\bibitem{MandC}  Mehrabadi, M.M. and Cowin, S.C., 1978, Initial Planar
deformation of dilatant granular materials, \textit{J. Mech. Phys. Solids, }%
\textbf{26}, 269-284.

\bibitem{Sokol}  Sokolovskii, V.V., 1965, Statics of Granular Media,
Pergamon Press, Oxford, 3rd. edition.

\bibitem{Spen1}  Spencer, A.J.M., 1964, A theory of the kinematics of ideal
soils under plane strain conditions, \textit{J. Mech. Phys. Solids, }\textbf{%
12}, 337-351.\textbf{\ }

\bibitem{Spen2}  Spencer, A.J.M., 1981, Deformation of ideal granular
materials, \textit{The Rodney Hill 60th. Anniversary Volume, }eds. H.G.
Hopkins and M.J. Sewell, Pergamon Press, Oxford, 607-652.
\end{thebibliography}

\bigskip

\newpage

\section*{Appendix A}

............. 
\begin{equation*}
\begin{array}{c}
\sigma _{11}=Ax_{1}+Bx_{2}+C, \\ 
\sigma _{12}=Dx_{1}+Ex_{2}+F, \\ 
\sigma _{22}=Gx_{1}+Hx_{2}+I+\rho gx_{2}.
\end{array}
\end{equation*}
Substituting into the yield condition gives 
\begin{equation*}
D=A\tan \left( \phi +\nu \right) ,\qquad E=B\tan \left( \phi +\nu \right)
,\qquad F=C\tan \left( \phi +\nu \right)
\end{equation*}
\begin{equation*}
\begin{array}{c}
\sigma _{11}=Ax_{1}+Bx_{2}+C, \\ 
\sigma _{12}=A\tan \left( \phi +\nu \right) x_{1}+B\tan \left( \phi +\nu
\right) x_{2}+C\tan \left( \phi +\nu \right) , \\ 
\sigma _{22}=Gx_{1}+Hx_{2}+I.
\end{array}
\end{equation*}
Substituting into the equilibrium equations gives 
\begin{equation*}
A=-B\tan \left( \phi +\nu \right) ,
\end{equation*}
\begin{equation*}
H=-A\tan \left( \phi +\nu \right) .
\end{equation*}
\begin{equation*}
\begin{array}{c}
\sigma _{11}=-B\tan \left( \phi +\nu \right) x_{1}+Bx_{2}+C=B\left[
x_{2}-\tan \left( \phi +\nu \right) x_{1}\right] +C \\ 
\sigma _{12}=\left\{ B\left[ x_{2}-\tan \left( \phi +\nu \right) x_{1}\right]
+C\right\} \tan \left( \phi +\nu \right) , \\ 
\sigma _{22}=Gx_{1}+Hx_{2}+I.
\end{array}
\end{equation*}
Substituting into the traction free condition gives $\eta =\phi +\nu .$%
\begin{equation*}
-\left[ \left\{ B\left[ x_{2}-\tan \left( \phi +\nu \right) x_{1}\right]
+C\right\} \right] \sin ^{2}\left( \phi +\nu \right) +\left[
Gx_{1}+Hx_{2}+I+\rho gx_{2}\right] \cos ^{2}\left( \phi +\nu \right) =0
\end{equation*}
so 
\begin{equation*}
I=-C\tan ^{2}\left( \phi +\nu \right)
\end{equation*}
\begin{equation*}
-B\left[ x_{2}-\tan \left( \phi +\nu \right) x_{1}\right] \sin ^{2}\left(
\phi +\nu \right) +\left[ Gx_{1}+\left( H+\rho g\right) x_{2}\right] \cos
^{2}\left( \phi +\nu \right) =0
\end{equation*}
\begin{equation*}
\begin{array}{c}
-B\sin ^{2}\left( \phi +\nu \right) +\left( H+\rho g\right) \cos ^{2}\left(
\phi +\nu \right) =0 \\ 
B\tan \left( \phi +\nu \right) \sin ^{2}\left( \phi +\nu \right) +G\cos
^{2}\left( \phi +\nu \right) =0
\end{array}
\end{equation*}
\begin{equation*}
\begin{array}{c}
H=-\rho g+B\tan ^{2}\left( \phi +\nu \right) \\ 
G=-B\tan ^{3}\left( \phi +\nu \right)
\end{array}
\end{equation*}
\begin{equation*}
\begin{array}{c}
\sigma _{11}=-B\tan \left( \phi +\nu \right) x_{1}+Bx_{2}+C=B\left[
x_{2}-\tan \left( \phi +\nu \right) x_{1}\right] +C \\ 
\sigma _{12}=\left\{ B\left[ x_{2}-\tan \left( \phi +\nu \right) x_{1}\right]
+C\right\} \tan \left( \phi +\nu \right) , \\ 
\sigma _{22}=\left\{ B\left[ x_{2}-\tan \left( \phi +\nu \right) x_{1}\right]
-C\right\} \tan ^{2}\left( \phi +\nu \right) +\rho gx_{2}.
\end{array}
\end{equation*}
\begin{figure}[h]
\hbox to \textwidth{\hss \epsffile{newcon.eps}\hss}
\caption{Kelvin's medium}
\label{con}
\end{figure}

\end{document}
