A Granular Fluid as a Limit of a Binary Mixture Theory—Treated as a One-Constituent Goodman–Cowin-Type Material

Abstract

Goodman and Cowin (1972) (Goodman and Cowin (1971). J. Fluid Mech. 45, 321–339. [6]) proposed a continuum theory of a dry cohesionless granular material in which the solid volume fraction $\nu$ is treated as an independent kinematic field for which an additional balance law of equilibrated forces is postulated. They motivated this additional balance law as an equation describing the kinematics of the microstructure and employed a variational formulation for its derivation. By adopting the müller –Liu approach to the exploitation of the entropy inequality we show that in a constitutive model containing $\nu , \dot{\nu }$ and $\mathrm {grad}\,\nu$ as independent variables, results agree with the classical Coleman –Noll approach only, provided the Helmholtz free energy does not depend on $\dot{\nu }$, for which the Goodman–Cowin equations are reproduced. This reduced theory is then applied to analyses of steady fully developed horizontal shearing flows and gravity flows of granular materials down an inclined plane and between parallel plates. It is demonstrated that the equations and numerical results presented by Passman et al. (1980) (Passman, Nunziato, Bailey and Thomas (1980). J. Eng. Mech. Div. ASCE 106, 773–783. [15]) are false, and they are corrected. The results show that the dynamical behavior of these materials is quite different from that of a viscous fluid. In some cases, the dilatant shearing layers exist only in the narrow zones near the boundaries. They motivated this additional balance law as an equation describing the kinematics of the microstructure and employed a variational formulation for its derivation. In an appendix, we present a variational formulation, treating the translational velocity and solid volume fraction as generalized coordinates of a Lagrange an formulation.

This chapter is heavily based on the paper by Wang and Hutter [20].

Appendix 29.A   Variational Principle for a Goodman–Cowin Type Granular Material

29.1.1 29.A.1 Preliminaries

In this appendix, we present a variational principle providing an alternative motivation of the balance law of equilibrated forces, as introduced by Goodman and Cowin [7]. We shall use Cartesian tensor notation; $x_{i}$ $(i=1,2,3)$ and $X_{a}$ $(a=1,2,3)$ will be the coordinates of material points in the present and reference configurations, respectively, and t will denote the time. The volume fraction $\nu ({\varvec{x}})$ $(0\leqslant \nu \leqslant 1)$ and the true density $\gamma ({\varvec{x}})>0$ will denote the averaged volume fraction over an RVE, filling densely the space of the body ${\mathcal B}$ so that the bulk density is expressible as

$$\begin{aligned} \rho = \gamma \nu , \end{aligned}$$

(29.116)

also filling the space ${\mathcal B}$ densely.

Let us consider equilibrium states $x_{i}({\varvec{X}})$, $\nu ({\varvec{X}})$. Neighboring such states will be denoted as $x_{i}({\varvec{X}}, \lambda )$, $\nu ({\varvec{X}}, \lambda )$, parameterized by $\lambda$. The variations of the positions $\delta x_{i}$ and volume fraction $\delta \nu$ may then be defined by

$$\begin{aligned} \delta x_{i}=\frac{\mathrm {d} x_{i}}{\mathrm {d} \lambda }\bigg |_{\begin{array}{l} \scriptstyle \lambda =0 \\ \scriptstyle {\varvec{X}} \;\mathrm {fixed}\end{array}} , \qquad \delta \nu = \frac{\mathrm {d} \nu }{\mathrm {d} \lambda }\bigg |_{\begin{array}{l}\scriptstyle \lambda =0 \\ \scriptstyle {\varvec{X}} \; \mathrm {fixed}\end{array}}. \end{aligned}$$

(29.117)

If the volume fraction is thought of as a function of ${\varvec{x}}$ and t, $\nu =\nu ({\varvec{x}}, t)$ rather than $\nu ({\varvec{X}},t)$ , the variation (29.117)$_{2}$ is represented by

$$\begin{aligned} \delta \nu = \frac{\partial \nu }{\partial \lambda }\bigg |_{\begin{array}{l}\scriptstyle \lambda =0 \\ \scriptstyle {\varvec{x}} \;\mathrm {fixed}\end{array}} + \frac{\partial \nu }{\partial x_{i}} \frac{\partial x_{i}}{\partial \lambda }\bigg |_{\lambda =0}. \end{aligned}$$

(29.118)

A variation of $\nu$ holding the spatial position ${\varvec{x}}$ fixed is denoted by Cowin and Goodman [2] by

$$\begin{aligned} \varDelta \nu = \frac{\partial \nu }{\partial \lambda }\bigg |_{\begin{array}{l}\scriptstyle \lambda =0 \\ \scriptstyle {\varvec{x}} \;\mathrm {fixed}\end{array}} . \end{aligned}$$

(29.119)

Substituting (29.119) into (29.118) and using (29.117)$_{1}$ yields

$$\begin{aligned} \delta \nu = \varDelta \nu + \nu _{,i}\delta x_{i}. \end{aligned}$$

(29.120)

From this representation follows

$$\begin{aligned} \underbrace{(\delta \nu )_{,j} }_{\delta (\nu _{,j})+\nu _{,i}(\delta x_{i})_{,j}} = \underbrace{ (\varDelta \nu )_{,j}}_{\varDelta (\nu _{,j})} + \nu _{,ij}\delta x_{i} + \nu _{,i}(\delta x_{i})_{,j}, \end{aligned}$$

(29.121)

so that

$$\begin{aligned} \delta (\nu _{,j})=\varDelta (\nu _{,j})+\nu _{,ij}\delta x_{i}. \end{aligned}$$

(29.122)

The underbraced term on the left-hand side of (29.121) can be justified by the fact that $(\delta \nu )_{,j}$ differs from $(\delta \nu _{,j})$ by the convective term $\nu _{,i}(\delta x_{i})_{,j}$, since ${\varvec{x}}$ but not ${\varvec{X}}$ is held fixed, whereas in the underbraced term on the right-hand side of (29.121) $\varDelta$ is the variation holding the present configuration fixed.

29.1.2 29.A.2 Variational Principle

We now suppose that in the granular material the interaction of the grains is described by a scalar valued balance law that is added to the force balance by the fact that besides $\delta x_{i}$ also $\delta \nu$ is an independent variation. So, the stored energy function W is for equilibrium states given by a function of the form $W=W(\gamma , \nu , \nu _{,i})$. It follows that the variation of the total energy function of the body can be postulated in the form

$$\begin{aligned} \delta \int _{{\mathcal B}} \gamma \nu W \,\mathrm {d}v= & {} \int _{{\mathcal B}} \left( \gamma \nu b_{i} \,\delta x_{i}+ \gamma \nu \ell \,\delta \nu \right) \,\mathrm {d}v \nonumber \\&+ \int _{\partial {\mathcal B}}\left( t_{i} \,\delta x_{i} + H \,\delta \nu \right) \mathrm {d}a . \end{aligned}$$

(29.123)

The volume term on the right-hand side of this equation consists of the virtual work done by the volumetric body force $\gamma \nu b_{i}$, subjected to the virtual displacement $\delta x_{i}$. In addition the scalar intrinsic equilibrated body force in Goodman -Cowin ’s terminology, performs work when being subjected to volume fraction variations $\delta \nu$. The second term on the right-hand side of (29.123) represents the associated surface work: $t_{i} \,\delta x_{i}$ is this working by the stress tractions when being subjected to the variations of the surface displacement $\delta x_{i}$ and the variation of the equilibrated surface traction H, when being exposed to a variation of the volume fraction at the surface. This terminology follows Goodman and Cowin , however, (29.123) introduces the equilibrated intrinsic body force $\ell$ and the equilibrated traction H (both as scalar quantities). They may equally be introduced in an abstract way without specifying their meaning by a model interpretation. It is given by Goodman in 1969 [5].

We require the mass of a body to be conserved during variations. Thus,

$$\begin{aligned}&\qquad M=\int _{{\mathcal B}} \gamma \nu \,\mathrm {d}v \equiv \mathrm {const.} \nonumber \\&\Longrightarrow \int _{{\mathcal B}}\delta (\gamma \nu )\,\mathrm {d}v +\int _{{\mathcal B}} \gamma \nu \,\delta (\mathrm {d}v) \equiv 0. \end{aligned}$$

(29.124)

Because this identity applies for all ${\mathcal B}$, also infinitesimal ones, we have

$$\begin{aligned} \delta (\gamma \nu ) \,\mathrm {d}v + \gamma \nu \,\delta (\mathrm {d}v) = 0, \end{aligned}$$

from which we deduce

(29.125)

Similarly, (29.124)$_{1}$ can be multiplied with the stored energy function to yield

$$\begin{aligned}&\delta (\gamma \nu )W \,\mathrm {d}v + \gamma \nu W \,\delta (\mathrm {d}v) = 0, \nonumber \\ \\ \mathrm {or}&\int _{{\mathcal B}}\delta (\gamma \nu )W\,\mathrm {d}v + \int _{{\mathcal B}}\gamma \nu W \,\delta (\mathrm {d}v) = 0. \nonumber \end{aligned}$$

(29.126)

These equations imply

$$\begin{aligned} \delta \int _{{\mathcal B}}\gamma \nu W \,\mathrm {d}v= & {} \int _{{\mathcal B}} \gamma \nu W(\delta \,\mathrm {d}v) \nonumber \\= & {} \underbrace{ \int _{{\mathcal B}}\delta (\gamma \nu )W\,\mathrm {d}v + \int _{{\mathcal B}}\gamma \nu W \,\delta (\mathrm {d}v) }_{=0\; \mathrm {because~of~mass~balance}} + \int _{{\mathcal B}}\gamma \nu \,\delta W \,\mathrm {d}v \nonumber \\= & {} \int _{{\mathcal B}}\gamma \nu \,\delta W \,\mathrm {d}v . \end{aligned}$$

(29.127)

The expression on the right-hand side of (29.126) shall now alternatively be written. To this end, we observe that $W=W(\gamma , \nu , \nu _{,i})$ and obtain

Multiplying this relation by ${\gamma \nu }$ yields

$$\begin{aligned} \gamma \nu \,\delta W= & {} \underbrace{ - \gamma ^{2} \nu \frac{\partial W}{\partial \gamma }(\delta x_{i})_{,i} - \gamma \nu \frac{\partial W}{\partial \nu _{,i}}(\delta x_{k})_{,i}\nu _{,k} }_{(i)} \nonumber \\&+ \underbrace{\gamma \nu \frac{\partial W}{\partial \nu _{,i}} (\delta \nu )_{,i} }_{(ii)} \underbrace{-\gamma ^{2} \frac{\partial W}{\partial \gamma } \delta \nu +\gamma \nu \frac{\partial W}{\partial \nu }\delta \nu }_{(iii)}. \end{aligned}$$

(29.128)

Following Cowin and Goodman [2], this relation now suggests the definitions

$$\begin{aligned}&\hat{p}:= \gamma \nu ^{2} \frac{\partial W}{\partial \nu }, \quad \; p := \gamma ^{2} \nu \frac{\partial W}{\partial \gamma }, \quad \; h_{i} := \gamma \nu \frac{\partial W}{\partial \nu _{,i}} \end{aligned}$$

(29.129)

$$\begin{aligned}&P := \frac{1}{\nu }(p - \hat{p}) = \frac{1}{\nu }\left( \gamma ^{2} \nu \frac{\partial W}{\partial \gamma } - \gamma \nu ^{2}\frac{\partial W}{\partial \nu }\right) , \end{aligned}$$

(29.130)

$$\begin{aligned}&T_{ij} := - p\,\delta _{ij} - h_{i}\nu _{,j} = -\gamma ^{2} \nu \frac{\partial W}{\partial \gamma } \delta _{ij} - \gamma \nu \frac{\partial W}{\partial \nu _{,i}}\nu _{,j}. \end{aligned}$$

(29.131)

Thus,

$$\begin{aligned}&T_{ij}(\delta x_{j})_{,i} = - \gamma ^{2} \nu \frac{\partial W}{\partial \gamma } (\delta x_{i})_{,i} - \gamma \nu \frac{\partial W}{\partial \nu _{,i}} \nu _{,j}(\delta x_{j})_{,i} \qquad \qquad (i),\end{aligned}$$

(29.132)

$$\begin{aligned}&h_{i}(\delta \nu )_{,i} = \gamma \nu \frac{\partial W}{\partial \nu _{,i}}(\delta \nu )_{,i} \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \,\, (ii), \end{aligned}$$

(29.133)

$$\begin{aligned}&\quad -P (\delta \nu ) = - \left( \gamma ^{2} \frac{\partial W}{\partial \gamma } - \gamma \nu \frac{\partial W}{\partial \nu }\right) \delta \nu \qquad \qquad \qquad \qquad \,\, (iii). \end{aligned}$$

(29.134)

The symbols (i), (ii), and (iii) correspond to the underscored terms in (29.128), which can now be written as

$$\begin{aligned} \gamma \nu \,\delta W= & {} T_{ij}(\delta x_{j})_{,i} + h_{i}(\delta \nu )_{,i} - P\,\delta \nu \nonumber \\= & {} (T_{ij} \delta x_{j})_{,i }-T_{ij,i}\delta x_{j}+(h_{i}\,\delta \nu )_{,i}-h_{i,i}\,\delta \nu -P \,\delta \nu . \end{aligned}$$

(29.135)

Integrating this expression over the body ${\mathcal B}$ and using the Gauss theorem leads to

$$\begin{aligned} \int _{{\mathcal B}}\gamma \nu \,\delta W \,\mathrm {d}v= & {} - \int _{{\mathcal B}} \left\{ T_{ij,i} \,\delta x_{j}+ (h_{i,i} + P)\delta \nu \right\} \,\mathrm {d}v \nonumber \\&+ \int _{\partial {\mathcal B}}\left\{ T_{ij}n_{i}\,\delta x_{j} + h_{i}n_{i} \,\delta \nu \right\} \mathrm {d}a. \end{aligned}$$

(29.136)

Next, using (29.127) in (29.123) and subsequently employing (29.136) implies

$$\begin{aligned}&\int _{{\mathcal B}}\left\{ (T_{ij,i} + \gamma \nu b_{j})\delta x_{j} + (h_{i,i}+P+\gamma \nu \ell )\delta \nu \right\} \,\mathrm {d}v \nonumber \\&+ \int _{\partial {\mathcal B}} \left\{ (t_{j} - T_{ij}n_{i}) \delta x_{j} + (H - h_{i}n_{i}) \right\} \delta \nu \,\mathrm {d}a \equiv 0, \end{aligned}$$

(29.137)

which must be valid for arbitrary and independent variations of $\delta x_{j}$ and $\delta \nu$, which leads to

$$\begin{aligned} \left. \begin{array}{l} T_{ij,i} + \gamma \nu b_{j} = 0 \\ h_{i,i}+P+\gamma \nu \ell = 0 \end{array}\right\} \;\; \mathrm {in} \; {\mathcal B} , \qquad \left. \begin{array}{l} t_{j} = T_{ij}n_{i} \\ H = h_{i}n_{i} \end{array}\right\} \;\; \mathrm {on} \; \partial {\mathcal B}. \end{aligned}$$

(29.138)

These equations represent the equilibrium equations balancing the stress divergence and the body force and the equation of the balance of the equilibrated forces—divergence of the flux $h_{i}$, the volume fraction- dependent pressure P, and the intrinsic equilibrated body force $\gamma \nu \ell$. The equations on the right-hand side of (29.138) define the “stress boundary conditions”, which would also follow from Cauchy ’s tetrahedral argument.

Cowin and Goodman [2] state at this point that “given the stored energy function W as a function of $\gamma , \nu , \nu _{,i}$ the relations (29.128) and (29.129) together with the equations (29.137) constitute four equations in terms of the two unknowns $\gamma$ and $\nu$”. The system is, therefore, overdetermined as noted by Jenkins [12]. Following this author Cowin and Goodman then propose an equivalent alternative variational principle “which leads directly to a reduction of the equations to form a compatible system”.

Alternative variational formulation (compressible granules) Jenkins [12] did not treat $\delta x_{i}$ and $\delta \nu$ as independent fields, but $\delta x_{i}$ and $\varDelta \nu$, defined in (29.119) and related to $\delta \nu$ via (29.120). The variation of the stored energy function then takes the form

$$\begin{aligned} \gamma \nu \,\delta W= & {} \gamma \nu \frac{\partial _W}{\partial \gamma } \delta \gamma + \gamma \nu \frac{\partial W}{\partial \nu }\delta \nu +\gamma \nu \frac{\partial W}{\partial \nu _{,i}}\delta (\nu _{,i}), \end{aligned}$$

(29.139)

where

$$\begin{aligned}&\delta \gamma = -\gamma (\delta x_{i})_{,i}-\gamma \frac{\delta \nu }{\nu }, \qquad \delta \nu {\mathop {=}\limits ^{(29.120)}} \varDelta \nu +\nu _{,i}\,\delta x_{i}, \\&\delta (\nu _{,i}){\mathop {=}\limits ^{(29.122)}}\varDelta (\nu _{,i})+\nu _{,ij}\,\delta x_{j}, \nonumber \end{aligned}$$

so that (29.139) can be written as

$$\begin{aligned} \gamma \nu \,\delta W= & {} \gamma \nu \frac{\partial W}{\partial \gamma }\left\{ - \gamma (\delta x_{i})_{,i} - \gamma \frac{\delta \nu }{\nu } \right\} + \gamma \nu \frac{\partial W}{\partial \nu } \left\{ \varDelta \nu + \nu _{,i} \,\delta x_{i}\right\} \nonumber \\&+\gamma \nu \frac{\partial W}{\partial \nu _{,i}} \left\{ \varDelta (\nu _{,i})+ \nu _{,ij}\,\delta x_{j}\right\} . \qquad \end{aligned}$$

(29.140)

Next, note that

$$\begin{aligned} W_{,i}=\frac{\partial W}{\partial \gamma } \gamma _{,i} + \frac{\partial W}{\partial \nu }\nu _{,i} + \frac{\partial W}{\partial \nu _{,j}}\nu _{,ji}. \end{aligned}$$

With this and the definitions (29.129) and (29.130) for $p, h_{i}$, and P a somewhat involved but simple identification of individual terms shows that

$$\begin{aligned} \gamma \nu \,\delta W= & {} \gamma \nu \left( p\left( \frac{1}{\gamma \nu }\right) _{,i}+W\right) \delta x_{i} \underbrace{ -p(\delta x_{i})_{,i} }_{-(p\,\delta x_{i})_{,i}+p_{,i}\,\delta x_{i}} -P\varDelta \nu + \underbrace{ h_{i}(\varDelta \nu )_{,i} }_{(h_{i}\varDelta \nu )_{,i} - h_{,i} \varDelta \nu } \nonumber \\= & {} \gamma \nu \Bigg \{ \underbrace{ p\left( \frac{1}{\gamma \nu }\right) _{,i} + \frac{p_{,i}}{\gamma \nu }}_{[p/(\gamma \nu )]_{i}} + W_{,i}\Bigg \}\delta x_{i} \nonumber \\&- (p\,\delta x_{i})_{,i} + (h_{i}\varDelta \nu )_{,i} - P\varDelta \nu - h_{i,i}\varDelta \nu . \end{aligned}$$

Integrating this over the body ${\mathcal B}$, and employing the divergence theorem in the last two terms on the right-hand side yield

$$\begin{aligned} \int _{{\mathcal B}}\gamma \nu \,\delta W \, \mathrm {d}v= & {} \int _{{\mathcal B}} \left\{ \gamma \nu \left( \frac{p}{\gamma \nu }+W\right) _{,i} \,\delta x_{i} - (h_{i,i}+P) \varDelta \nu \right\} \mathrm {d}v \nonumber \\&+ \int _{\partial {\mathcal B}}\left\{ -pn_{i} \,\delta x_{i}+ h_{i,i}n_{i}\right\} \varDelta \nu \,\mathrm {d}a. \end{aligned}$$

(29.141)

The right-hand side of this expression must be identified with the right-hand side of (29.123), viz.,

$$\begin{aligned}&\int _{{\mathcal B}}\left\{ \gamma \nu \left[ \left( \frac{p}{\gamma \nu }+W\right) _{,i} - \ell \nu _{,i} - b_{i}\right] \delta x_{i} - \left[ h_{i,i}+P+\gamma \nu \ell \right] \varDelta \nu \right\} \mathrm {d}v \nonumber \\&+\int _{\partial {\mathcal B}}\left\{ \left( pn_{i} + t_{i} + H\nu _{,i}\right) \delta x_{i} + (h_{i}n_{i} - H)\varDelta \nu \right\} \mathrm {d}a \equiv 0. \end{aligned}$$

(29.142)

For this to hold for all $\delta x_{i}$ and all $\delta \nu$, we obtain

$$\begin{aligned} \left. \begin{array}{l} \displaystyle \frac{\partial }{\partial x_{i}}\left( \frac{p}{\gamma \nu }+W\right) = b_{i}+ \ell \nu _{,i} \\ h_{i,i}+P+\gamma \nu \ell = 0 \end{array} \right\} \; \mathrm {in}\; {\mathcal B}, \quad \left. \begin{array}{l} -p n_{i} = t_{i}+H \nu _{,i} \\ H = h_{i} n_{i} \end{array} \right\} \; \mathrm {in} \; {\partial \mathcal B}.\qquad \quad \end{aligned}$$

(29.143)

These equations possess the advantage that they can directly be integrated, if $b_{i}$ and $\ell$ are expressed in terms of a potential as suggested by Jenkins [12]. If we choose a force potential $\varPhi (x_{i}, \nu )$ such that

$$\begin{aligned} b_{i}=-\frac{\partial \varPhi }{\partial x_{i}}, \quad \ell =-\frac{\partial \varPhi }{\partial \nu }, \end{aligned}$$

(29.144)

then the field equations and boundary conditions take the forms

$$\begin{aligned} \left. \begin{array}{l} \displaystyle \frac{\partial }{\partial x_i}\left( \frac{p}{\gamma \nu }+W+\varPhi \right) = 0 \\ \displaystyle \frac{\partial h_{i}}{\partial x_{i}}+P = \gamma \nu \frac{\partial \varPhi }{\partial \nu } \end{array} \right\} \; \mathrm {in} \;{\mathcal B}, \quad \left. \begin{array}{l} - p n_{i} = t_{i}+H \nu _{,i} \\ H=h_{i}n_{i}, \; t_{j}=T_{ij}n_{i} \end{array} \right\} \; \mathrm {on}\; \partial {\mathcal B}.\qquad \quad \; \end{aligned}$$

(29.145)

The particular form of these field equations and boundary conditions shows that they emerge naturally when one considers variations of $\nu$ by holding ${\varvec{x}}$ fixed rather than ${\varvec{X}}$.

29.1.3 29.A.3 Variational Principle for Density Preserving Granules

For a density preserving material the specific true mass is constant, $\gamma =\mathrm {const.}$, and here, therefore, $\delta \gamma \equiv 0$. From mass balance (29.125), it then follows that

$$\begin{aligned} \delta \nu = -\nu (\delta x_{i})_{,i}. \end{aligned}$$

(29.146)

It follows that the variations of $\delta \nu$ and $\delta {\varvec{x}}$ are no longer independent. Thus, the variational principle (29.123) can now be expressed as

$$\begin{aligned} \delta \int _{{\mathcal B}}\gamma \nu W \,\mathrm {d}v= & {} \int _{{\mathcal B}} \left( \gamma \nu b_{i}\,\delta x_{i} - \underline{ \gamma \nu ^{2} \ell (\delta x_{i})_{,i} } \right) \,\mathrm {d}v \nonumber \\&+\delta \int _{\partial {\mathcal B}}\left( t_{i}\,\delta x_{i} - \nu H(\delta x_{i})_{,i}\right) \mathrm {d}a . \end{aligned}$$

(29.147)

For density preserving granules, $W=W(\nu , \nu _{i})$ and, therefore,

$$\begin{aligned} \delta W = \frac{\partial W}{\partial \nu }\delta \nu + \frac{\partial W}{\partial \nu _{,i}}\delta (\nu _{,i}). \end{aligned}$$

(29.148)

When using

$$\begin{aligned} \delta (\nu _{,j}) = (\delta \nu )_{,j} - \nu _{,i}(\delta x_{i})_{,j},\quad \mathrm {see\; (29.121)} \end{aligned}$$

and (29.146) in the expression (29.147) and employing the definitions (29.129)–(29.131), straightforward computations yield

$$\begin{aligned} -\gamma \nu \,\delta W= & {} \left( - \hat{p}\,\delta _{ij}+\nu h_{k,k}\,\delta _{ij}+ h_{i}\nu _{,j}\right) _{,i}\,\delta x_{j}+ \underline{ [h_{i}\nu (\delta x_{j})_{,j}]_{,i} } \nonumber \\&+ \underline{ [\hat{p}\,\delta _{ij} - \nu h_{k,k}\,\delta _{ij}+h\nu _{,j}]_{,i} } . \end{aligned}$$

(29.149)

This equation is now integrated over the body; hereby, the underlined terms are transformed to surface integrals using the divergence theorem. This step leads to

$$\begin{aligned} \delta \int _{{\mathcal B}}\gamma \nu W\,\mathrm {d}v= & {} \int _{{\mathcal B}} \gamma \nu \,\delta W\,\mathrm {d}v\nonumber \\= & {} - \int _{{\mathcal B}}\left( -\hat{p}\,\delta _{ij}+\nu h_{k,k} \,\delta _{ij} - h_{i} \nu _{,j}\right) _{,i}\,\delta x_j\,\mathrm {d}v - \int _{\partial {\mathcal B}}h_{i}n_{i}\nu (\delta x_{j})_{,j}\,\mathrm {d}v \nonumber \\&+ \int _{\partial {\mathcal B}}n_{i}\left( -\hat{p}\,\delta _{ij}+\nu h_{k,k}\,\delta _{ij} - h_{i}\nu _{,j}\right) _{,i}\,\delta x_{j}\,\mathrm {d}a . \end{aligned}$$

(29.150)

The right-hand side of this variational statement must equal the right-hand side of (29.147), in which we transform the underlined term by writing

$$\begin{aligned} \gamma \nu ^{2}\ell (\delta x_{i})_{,i} = \underline{ (\gamma \nu ^{2} \ell \,\delta x_{i})_{,i} } - (\gamma \nu ^{2} \ell )_{,i}\,\delta x_{i} \end{aligned}$$

and use the divergence theorem for the first term on the right-hand side. These computations lead to

$$\begin{aligned}&\int _{{\mathcal B}}\bigg \{\bigg [\bigg (\underbrace{ -\hat{p}\,\delta _{ij}+\nu h_{k,k}\,\delta _{ij}-h_{j}\nu _{,i}+\gamma \nu ^{2}\ell \,\delta _{ij} }_{:=T_{ij}} \bigg )_{,j}+\gamma \nu b_{i} \bigg ]\delta x_{i}\bigg \} \mathrm {d}v \nonumber \\&+ \int _{\partial {\mathcal B}}\bigg \{\bigg [t_{i}+ \underbrace{ \left( \hat{p}n_{i}-\nu h_{k,k}n_{i} + h_{j} n_{j}\nu _{,i} - \gamma \nu ^{2}\ell n_{i} \right) }_{:=-T_{ij}n_j }\bigg ] \delta x_{i}\bigg \}\mathrm {d}a \nonumber \\&+ \int _{\partial {\mathcal B}} (h_{i}n_{i} - H)\nu (\delta x_{j})_{,j}\,\mathrm {d}a \equiv 0, \end{aligned}$$

(29.151)

which must hold for all $\delta x_{i}$ in ${\mathcal B}$ and on $\partial {\mathcal B}$. However, this requires caution in the evaluation of $(\delta x_{j})_{,j}$ along the free surface $\partial {\mathcal B}$. This quantity is specified by $\delta x_{j}$ everywhere in $\mathcal B\cup \partial {\mathcal B}$ except in the direction normal to the boundary $\partial {\mathcal B}$. Letting $D_{k}$ be the surface gradient and $D_{\perp }$ the directional derivative perpendicular to $\partial {\mathcal B}$, then we have

$$\begin{aligned} (\delta x_{i})_{,i} = (D_{\perp }n_{i} + D_{i})\delta x_{i}, \quad \mathrm {on}\; \partial {\mathcal B}, \end{aligned}$$

(29.152)

and the last term in (29.151) becomes

$$\begin{aligned} \int _{\partial {\mathcal B}}\left[ \left( h_{i}n_{i} - H \right) \nu \left( D_{\perp }n_{j}+D_{j}\right) \delta x_{j}\right] \mathrm {d}a. \end{aligned}$$

(29.153)

Substituting (29.153) into (29.151) yields

$$\begin{aligned} 0\equiv & {} \int _{{\mathcal B}}(T_{ij,j}+\gamma \nu b_{i})\delta x_{i} \,\mathrm {d}v \nonumber \\&+ \int _{\partial {\mathcal B}}\left[ (t_{i} - T_{ij}n_{j})\delta x_{i}+(h_{i}n_{i}-H)\nu D_{\perp }n_{i}\,\delta x_{j}+D_{j}\,\delta x_{j})\right] \mathrm {d}a.\qquad \quad \end{aligned}$$

(29.154)

As the variations $\delta x_{i}$ in ${\mathcal B}\cup \partial {\mathcal B}$ and $D_{\perp }\,\delta x_{j}$ on $\partial {\mathcal B}$ are independent, (29.154) implies

$$\begin{aligned}&T_{ij,j} + \gamma \nu b_{i} = 0, \quad \quad \quad \,\,\, \mathrm {in}\; {\mathcal B}, \end{aligned}$$

(29.155)

$$\begin{aligned}&t_{i} = T_{ij}n_{j}, \quad H = h_{i}n_{i}, \; \mathrm {on}\; {\partial \mathcal B}, \end{aligned}$$

(29.156)

in which $T_{ij}$ is defined in the subbraced terms of (29.151), formally,

$$\begin{aligned} T_{ij}: = -\hat{p}\,\delta _{ij}+\nu h_{k,k}\,\delta _{ij}-h_{j}\nu _{j} +\gamma \nu ^{2}\ell \,\delta _{ij}. \end{aligned}$$

(29.157)

Equation (29.155) is the classical force balance, with the stress tensor defined in (29.157) and the flux boundary conditions in (29.156). They are also given by Goodman and Cowin [7].

Alternative variational formulation (density preserving) For density preserving granules $(\delta \gamma = 0)$ it was shown in (29.146) that $\delta \nu = -\nu (\delta x_{i})_{,i}$. With this, the variational principle (29.123) can be written as

$$\begin{aligned} \delta \int _{{\mathcal B}}\gamma \nu W \,\mathrm {d}v= & {} \int _{{\mathcal B}}\left[ \gamma b_{i}+\gamma \ell \nu _{,i}+(\gamma \nu \ell )_{,i}\right] \nu \,\delta x_{i} \,\mathrm {d}v \nonumber \\&+ \int _{\partial {\mathcal B}}\left[ \left( t_{i}+H\nu _{,i}-\gamma \nu ^{2} \ell n_{i}\right) \delta x_{i} - H(\nu \delta x_{i})_{,i}\right] \mathrm {d}a,\qquad \quad \end{aligned}$$

(29.158)

in which (29.146) and the divergence theorem have been employed. Next, using mass balance (see (29.126)) and the representation (29.148), along with (29.129)–(29.131) we obtain

$$\begin{aligned} \gamma \nu \,\delta W= & {} \left( \frac{\hat{p}}{\nu }+\gamma W-h_{j,j}\right) _{,i} \nu \,\delta x_{i} \nonumber \\&+ \left\{ \left( -\hat{p}+h_{j,j}\nu \right) \delta x_{i} - h_{i} \left( \nu \,\delta x_{j}\right) _{,j}\right\} _{,i}. \end{aligned}$$

(29.159)

Integrating both sides of this equation over the body ${\mathcal B}$ and employing the divergence theorem in the second term in curly brackets yields

$$\begin{aligned} \int _{{\mathcal B}}\gamma \nu \,\delta W \,\mathrm {d}v= & {} \int _{{\mathcal B}} \left\{ \frac{\hat{p}}{\nu }+\gamma W - h_{j,j}\right\} _{,i} \nu \,\delta x_{i}\, \mathrm {d}v \nonumber \\&+ \int _{\partial {\mathcal B}}\left\{ \left( -\hat{p}+h_{j,j}\nu \right) \delta x_{i}\right\} n_{i} \,\mathrm {d}a \nonumber \\&- \int _{\partial {\mathcal B}}h_{i} \nu (\delta x_{j})_{,j}n_{i}\,\mathrm {d}a. \end{aligned}$$

(29.160)

This expression must equal the variation of the stored energy (29.147) expressed in terms of the power of the body forces (classical and equilibrated) plus the corresponding workings of the surface forces; this statement can be written in the form

$$\begin{aligned} 0\equiv & {} \int _{{\mathcal B}} \left\{ \gamma b_{i}+\gamma \ell \nu _{,i}+ (\gamma \nu \ell )_{,i} - \left( \frac{\hat{p}}{\nu }+\gamma W- h_{j,j}\right) _{,i}\right\} \nu \,\delta x_{i} \,\mathrm {d}v \nonumber \\&+\int _{\partial {\mathcal B}}\left\{ t_{i}+H\nu _{,i}+\left( \hat{p}-h_{j,j}\nu - \gamma \nu ^{2}\ell \right) n_{i}\right\} \delta x_{i} \,\mathrm {d}a \nonumber \\&+\int _{\partial {\mathcal B}}\left[ h_{i}n_{i} - H\right] (\nu \,\delta x_{j})_{,j}\,\mathrm {d}a, \end{aligned}$$

(29.161)

an identity, which must hold for arbitrary $\delta x_{i}$ in ${\mathcal B}$ and on $\partial {\mathcal B}$ as well as arbitrary gradients of $\nu \delta x_{j}$ perpendicular to $\partial {\mathcal B}$. Therefore,

$$\begin{aligned}&\frac{\partial }{\partial x_{i}}\left( \frac{\hat{p}}{\nu }+\gamma W - \nu h_{j,j} - \gamma \nu \ell \right) = \gamma b_{i}+\gamma \ell \nu _{,i}, \quad \quad \,\, \mathrm {in} \;{ \mathcal B}, \nonumber \\&t_{i}+H\nu _{,i} = \left( -\hat{p}+\nu h_{j,j} +\gamma \nu ^{2} \ell \right) n_{i}, \quad H=h_{i}n_{i}, \quad \mathrm {on}\; {\partial {\mathcal B}}. \end{aligned}$$

(29.162)

These field equations and boundary conditions are equivalent to (29.156), if the stress tensor (29.157) is substituted.²

29.1.4 29.A.4 Dynamic Case

To apply the principle of virtual work to the dynamic case we follow the Lagrange an method employed in Chap. 26 for liquid crystals with tensorial-order parameter. In this case of a granular assembly the total energy in ${\mathcal B}$ is given in (26.19) as

$$\begin{aligned} {\mathcal F}= & {} \int _{{\mathcal B}}F\,\,\mathrm{d}v, \quad \mathrm{with} \nonumber \\ F= & {} \rho \left( {\varvec{v}}\cdot {\varvec{v}}+\phi +\sigma (\rho )+\kappa (\nu , \dot{\nu }) +\chi (\nu )\right) + W(\nu , \mathrm {grad}\,\nu ), \end{aligned}$$

(29.163)

in which the rank-i tensor $\pmb {\mathbb O}$ from (26.19) has been replaced by the volume fraction $\nu$; analogous to Chap. 26:

• $\rho =\gamma \nu$ is the mass density,

• $\textstyle {\frac{1}{2}}{\varvec{v}}\cdot {\varvec{v}}$, the translational kinetic energy of the grains,

• $\phi$, the potential energy of the body force, ${\varvec{f}}=-\mathrm {grad}\,\phi$,

• $\sigma (\rho )$, the potential energy due to the compressibility of the material,

• $\kappa (\nu , \dot{\kappa })$ the kinetic energy connected with the motion of the volume fraction,

• $\chi (\nu )$, the potential energy of the external actions on $\nu$ (which is physically nonrealistic and will be set to zero),

• W is interpreted as an “elastic” energy of the volume fraction and can be associated with the variation of the volume fraction.

Admittedly, some of these quantities are difficult to realistically identify with specific physical facts.

To evaluate the total time derivative of ${\mathcal F}$, we obtain the following individual expressions

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\int _{{\mathcal B}} {\textstyle \frac{1}{2}} \rho ({\varvec{v}} \cdot {\varvec{v}}) \mathrm {d}v&= \int _{{\mathcal B}}\rho \dot{\varvec{v}} \cdot {\varvec{v}} \,\mathrm {d}v, \nonumber \\ \frac{\mathrm {d}}{\mathrm {d}t} \int _{{\mathcal B}} \rho \dot{\phi }\,\mathrm {d}v&= \int _{{\mathcal B}}\rho \dot{\phi }\,\mathrm {d}v= -\int _{{\mathcal B}}\rho {\varvec{f}} \cdot {\varvec{v}} \,\mathrm {d}v, \nonumber \\ \frac{\mathrm {d}}{\mathrm {d}t}\int _{{\mathcal B}}\rho \sigma (\rho )\,\mathrm {d}v&= \int _{{\mathcal B}}\rho \underbrace{ \frac{\mathrm {d} \sigma }{\mathrm {d} \rho } }_{\sigma ^{\prime }} \dot{\rho }\,\mathrm {d}v = - \int _{{\mathcal B}}\rho ^{2} \sigma ^{\prime }(\rho ) \mathrm {div}\,{\varvec{v}} \,\mathrm {d}v, \nonumber \\ \frac{\mathrm {d}}{\mathrm {d}t} \int _{{\mathcal B}}\rho \chi (\nu )\,\mathrm {d}v&= \int _{{\mathcal B}}\rho \dot{\chi }(\nu )\,\mathrm {d}v = \int _{{\mathcal B}} \rho \frac{\partial \chi }{\partial \nu }\dot{\nu } \,\mathrm {d}v,\\ \frac{\mathrm {d}}{\mathrm {d}t} \int _{{\mathcal B}} W \,\mathrm {d}v \qquad&= \int _{{\mathcal B}}\left\{ \frac{\partial W}{\partial \nu }\dot{\nu } + \frac{\partial W}{\partial \mathrm {grad}\,\nu }(\mathrm {grad}\,\nu )^{\displaystyle \cdot } + W \mathrm {div}\,{\varvec{v}}\right\} \mathrm {d}v \nonumber \\ \frac{\mathrm {d}}{\mathrm {d}t} \int _{{\mathcal B}}\rho \kappa (\nu , \dot{\nu }) \,\mathrm {d}v&=\int _{{\mathcal B}}\rho \left( \kappa (\nu , \dot{\nu }) \right) ^{\displaystyle \cdot } \mathrm {d}v {\mathop {=}\limits ^{*}} \int _{{\mathcal B}} \rho \underbrace{ \left\{ \left( \frac{\partial \kappa }{\partial \dot{\nu }}\right) ^{\displaystyle \cdot } - \frac{\partial \kappa }{\partial \nu }\right\} }_{{\mathfrak m}} \dot{\nu } \,\mathrm {d}v \nonumber \\&=\int _{{\mathcal B}} \rho {\mathfrak m}\dot{\nu } \,\mathrm {d}v, \nonumber \end{aligned}$$

(29.164)

in which the step $({\mathop {=}\limits ^{*}})$ follows from our basic assumption that $\nu$ and $\dot{\nu }$ are interpreted as independent generalized coordinates and velocities in a Lagrange an formulation (see Appendix 26.A in Chap. 26, (26.19) and 26.206), or see [1]).

Evaluating $\dot{\mathcal F}$, with all expressions (29.164) substituted, yields

$$\begin{aligned} \dot{{\mathcal F}}= & {} \int _{{\mathcal B}}\Bigg \{ \rho \left( \dot{{\varvec{v}}}-{\varvec{f}}\right) \cdot {\varvec{v}}+\rho {\mathfrak m}\dot{\nu } + \left( \rho \frac{\partial \chi }{\partial \nu }+\frac{\partial W}{\partial \nu }\right) \dot{\nu }+\frac{\partial W}{\partial \mathrm {grad}\,\nu }(\mathrm {grad}\,\nu )^{\displaystyle \cdot } \nonumber \\&+ \underbrace{ \left( W-\rho ^{2}\sigma ^{\prime }(\rho )\right) \mathrm {div}\,{\varvec{v}} }_{[(W-\rho ^{2}\sigma ^{\prime })v_{i}]_{,i}-(W-\rho ^{2}\sigma ^{\prime })_{,i}v_{i}} + \frac{\partial W}{\partial \mathrm {grad}\,\nu } \underbrace{(\mathrm {grad}\,\nu )^{\displaystyle \cdot } }_{\mathrm {grad}\,\dot{\nu }-\mathrm {grad}\,\nu \cdot \mathrm {grad}\,{\varvec{v}} }\Bigg \}\mathrm {d}v \nonumber \\= & {} \int _{{\mathcal B}}\Bigg \{ \rho (\dot{\varvec{v}}-{\varvec{f}})\cdot {\varvec{v}} + \rho {\mathfrak m} \dot{\nu }+\left( \rho \frac{\partial \chi }{\partial \nu }+\frac{\partial W}{\partial \nu }\right) \dot{\nu } \nonumber \\&+ \underbrace{ \frac{\partial W}{\partial \mathrm {grad}\,{\nu }} \mathrm {grad}\,{\dot{\nu }}}_{[1]} - \underbrace{ \frac{\partial W}{\partial \mathrm {grad}\,\nu }\mathrm {grad}\,\nu \cdot \mathrm {grad}\,{\varvec{v}} }_{[2]} -\left( W-\rho ^{2}\sigma ^{\prime }(\rho )\right) _{,i}v_{i}\Bigg \} \mathrm {d}v\nonumber \\&+\int _{{\partial \mathcal B}}\left[ \left( W-\rho ^{2} \sigma ^{\prime }\right) v_{i}\right] n_{i} \,\mathrm {d}a, \end{aligned}$$

in which

$$\begin{aligned}&[1] = \int _{{\mathcal B}}\frac{\partial W}{\partial \nu _{,i}}\dot{\nu }_{,i} \,\mathrm {d}v \qquad = \int _{{\mathcal B}}\left\{ \left( \frac{\partial W}{\partial \nu _{,i}} \dot{\nu }\right) _{,i} - \left( \frac{\partial W}{\partial \nu _{,i}}\right) _{,i}\dot{\nu } \right\} \mathrm {d}v \\&\qquad \qquad \qquad \qquad \qquad \qquad =\int _{\partial {\mathcal B}}\frac{\partial W}{\partial \mathrm {grad}\,\nu }\cdot {\varvec{n}}\dot{\nu } \,\mathrm {d}a - \int _{{\mathcal B}}\mathrm {div}\,\left( \frac{\partial W}{\partial \mathrm {grad}\,\nu }\right) \dot{\nu }\,\mathrm {d}v \\&{[2]} =\int _{{\mathcal B}}\frac{\partial W}{\partial \nu _{,i}}\nu _{,k}v_{k,i}\,\mathrm {d}v = \int _{{\mathcal B}} \left[ \left( \frac{\partial W}{\partial \nu _{,i}}\nu _{,k}v_{k}\right) _{,i} - \left( \frac{\partial W}{\partial \nu _{,i}} \nu _{,k}\right) _{,i} v_{k}\right] \mathrm {d}v \\&\qquad \qquad \qquad \qquad \qquad \qquad = \int _{\partial {\mathcal B}}\frac{\partial W}{\partial \nu _{,i}} \nu _{,k}v_{k}n_{i}\,\mathrm {d}a \qquad \qquad - \int _{{\mathcal B}}\left( \frac{\partial W}{\partial \nu _{,i}} \nu _{,k}\right) _{,i}v_{k} \,\mathrm {d}v \\&\qquad \qquad \qquad \qquad \qquad \qquad = \int _{\partial {\mathcal B}}\left[ \left( \mathrm {grad}\,\nu \odot \frac{\partial W}{\partial \mathrm {grad}\,\nu }\right) {\varvec{n}} \cdot {\varvec{v}}\right] \mathrm {d}a \\&\qquad \qquad \qquad \qquad \qquad \quad - \int _{{\mathcal B}} \mathrm {div}\,\left\{ \mathrm {grad}\,\nu \odot \frac{\partial W}{\partial \mathrm {grad}\,\nu }\right\} \cdot {\varvec{v}} \,\mathrm {d}v. \end{aligned}$$

Recall that, if ${\varvec{a}}, {\varvec{b}}$ are vectors in ${\mathbb R}^{3}$, then $({\varvec{a}}\odot {\varvec{b}})_{ik} = a_{k}b_{i}$ which equals $({\varvec{a}}\otimes {\varvec{b}})^{T}$. Therefore,

$$\begin{aligned} \dot{\mathcal F}= & {} \int _{{\mathcal B}}\left\{ \rho (\dot{\varvec{v}} - {\varvec{f}}) \cdot {\varvec{v}} + \rho {\mathfrak m}\dot{\nu }+\left( \rho \frac{\partial \chi }{\partial \nu }+\frac{\partial W}{\partial \nu }\right) \dot{\nu } - \mathrm {div}\,\left( \frac{\partial W}{\partial \mathrm {grad}\,\nu } \right) \dot{\nu }\right. \nonumber \\&\left. + \mathrm {div}\,\left( \mathrm {grad}\,\nu \odot \frac{\partial W}{\partial \mathrm {grad}\,\nu }\right) \cdot {\varvec{v}} - \left( W - \rho ^{2} \sigma ^{\prime }(\rho )\right) {\varvec{v}} \right\} \mathrm {d}v \nonumber \\&+ \int _{\partial {\mathcal B}} \left\{ \left[ \left( W- \rho ^{2} \sigma ^{\prime }\right) {\varvec{n}}\cdot {\varvec{v}} + \frac{\partial W}{\partial \mathrm {grad}\,\nu }\cdot {\varvec{n}} \dot{\nu }\right] \right. \nonumber \\&\left. - \left( \mathrm {grad}\,\nu \odot \frac{\partial W}{\partial \mathrm {grad}\,\nu }\right) {\varvec{n}}\cdot {\varvec{v}} \right\} \mathrm {d} a . \end{aligned}$$

(29.165)

In classical mechanics $\dot{{\mathcal F}}$ is the total power input into a mechanical system. Here, this is not so, because (29.165) contains surface terms, which may unduly constrain the potentials W and $\sigma ^{\prime }$ without generalized surface forces. The power ${\mathcal W}^{s}$ of surface forces for a moving boundary ${\partial {\mathcal B}}$ can be written in the form

$$\begin{aligned} {\mathcal W}^{s} = \int _{\partial {\mathcal B}}\left\{ {\varvec{X}}^{s} \cdot {\varvec{v}} + \xi ^{s}\dot{\nu }\right\} \mathrm {d}a, \end{aligned}$$

(29.166)

in which ${\varvec{v}}$ and $\dot{\nu }$ are independent generalized velocities. Moreover, if the surface power can be derived from a potential, i.e., if

$$\begin{aligned} {\mathcal W}^{s} = \frac{\mathrm {d}}{\mathrm {d} t}\int _{\partial {\mathcal B}} W^{s}({\varvec{x}}, \nu ) \,\mathrm {d}a, \end{aligned}$$

(29.167)

in which $W^{s}$ is a scalar potential of position and solid volume fraction, then the relations

$$\begin{aligned} {\varvec{X}}^{s} = \frac{\partial W^{s}}{\partial {\varvec{x}}}, \quad \mathrm {and} \quad \xi ^{s} = \frac{\partial W^{s}}{\partial \nu } \end{aligned}$$

(29.168)

ensue. Adding (29.166) to the expression (29.165) now provides the possibility to relate ${\varvec{X}}^{s}$ and $\xi ^{s}$ to W and $\sigma ^{\prime }$. Explicitly, $\dot{\mathcal F}+{\mathcal W}^{s}$ can now be written as stated already in (26.29), namely as

$$\begin{aligned} \dot{{\mathcal F}}+ {\mathcal W}^{s}= & {} \int _{{\mathcal B}}\left\{ {\varvec{X}}\cdot {\varvec{v}} + \xi ^{s} \dot{\nu }\right\} \mathrm {d}v \nonumber \\&+ \int _{\partial {\mathcal B}}\left\{ \left( {\varvec{X}}^{b}+{\varvec{X}}^{s}\right) \cdot {\varvec{v}}+ \left( \xi ^{b}+\xi ^{s}\right) \dot{\nu }\right\} \mathrm {d}a, \end{aligned}$$

(29.169)

in which

$$\begin{aligned}&\left. \begin{array}{ll} \displaystyle {\varvec{X}} = \rho ({\varvec{v}} - {\varvec{f}}) - \mathrm {grad}\,(W-\rho ^{2}\sigma ^{\prime }) + \mathrm {div}\,\left( \mathrm {grad}\,\nu \odot \frac{\partial W}{\partial \mathrm {grad}\,\nu }\right) \\ \displaystyle \xi = \rho \left( {\mathfrak m}+\frac{\partial \chi }{\partial \nu }\right) + \frac{\partial W}{\partial \nu } - \mathrm {div}\,\left( \frac{\partial W}{\partial \mathrm {grad}\,\nu }\right) \end{array} \right\} \; \mathrm {in} \;{\mathcal B},\nonumber \\ \\&\left. \begin{array}{ll} \displaystyle {\varvec{X}}^{b} = \left( W-\rho ^{2}\sigma ^{\prime }\right) {\varvec{n}} - \left( \mathrm {grad}\,\nu \odot \frac{\partial W}{\partial \mathrm {grad}\,\nu }{\varvec{n}}\right) \\ \displaystyle \xi ^{b} = \frac{\partial W}{\partial \mathrm {grad}\,\nu }\cdot {\varvec{n}} \end{array}\right\} \; \mathrm {on}\; \partial {\mathcal B} . \nonumber \end{aligned}$$

(29.170)

We repeat, as we already pointed out in Chap. 26, p. 288: It is physically significant to note that the surface integral in (29.169) with generalized forces ${\varvec{X}}^{s}$ and $\xi ^{s}$ was introduced in (29.168). In the absence of dissipative terms we will prove that ${\varvec{X}}^{s}+{\varvec{X}}^{b} = {\varvec{0}}$, $\xi ^{s}+\xi ^{b} = 0$. So, when $W^{s} \equiv 0$ then ${\varvec{X}}^{b}$ and $\xi ^{b}$ would have to separately vanish, which would constrain the functions W and $\sigma$.

To avoid this situation, we now introduce, following Lord Rayleigh [18] or Sonnet and Virga [17], see Chap. 26, (26.31), the dissipation potential as a frame indifferent functional of the stretching tensor ${\varvec{D}}$, the volume fraction $\nu$ and its rate $\dot{\nu }$; explicitly

$$\begin{aligned} {\mathcal R} = \int _{{\mathcal B}}R\,\mathrm {d}v, \quad R = R(\nu , \dot{\nu }, {\varvec{D}}). \end{aligned}$$

(29.171)

We suppose R to be a bilinear function of $\dot{\nu }$ and ${\varvec{D}}$; writing it in terms of $\dot{\nu }$ and $\mathrm {grad}\,\nu$, we have

$$\begin{aligned} \delta {\mathcal R} = \int _{{\mathcal B}}\Bigg \{ \frac{\partial R}{\partial \dot{\nu }} \delta \dot{\nu }+\underbrace{ \frac{\partial R}{\partial \mathrm {grad}\,\nu }\cdot \underbrace{ \delta (\mathrm {grad}\,{\varvec{v}}) }_{\mathrm {grad}\,(\delta {\varvec{v}}) }}_{ [1]}\Bigg \} \,\mathrm {d}v, \end{aligned}$$

where

$$\begin{aligned}{}[1] = \frac{\partial R}{\partial v_{i,j}}(\delta v_{i})_{,j} = \left( \frac{\partial R}{\partial v_{i,j}} \delta v_{i}\right) _{,j} - \left( \frac{\partial R}{\partial v_{i,j}}\right) _{,j} \delta v_{i}, \end{aligned}$$

implying

$$\begin{aligned} \delta {\mathcal R}= & {} \int _{\partial {\mathcal B}}\left( \frac{\partial R}{\partial \mathrm {grad}\,{\varvec{v}}} {\varvec{n}}\right) \cdot \delta {\varvec{v}} \,\mathrm {d}a \nonumber \\&+\int _{{\mathcal B}}\left\{ \frac{\partial R}{\partial \dot{\nu }} \,\delta \dot{\nu } - \mathrm {div}\,\left( \frac{\partial R}{\partial \mathrm {grad}\,{\varvec{v}}}\right) \cdot \delta {\varvec{v}} \right\} \mathrm {d}v, \end{aligned}$$

(29.172)

or, since $\partial R/\partial \mathrm {grad}\,{\varvec{v}} = \partial R/\partial {\varvec{D}}$,

$$\begin{aligned} \delta {\mathcal R}= & {} \int _{\partial {\mathcal B}}\left( \frac{\partial R}{\partial {\varvec{D}}}{\varvec{n}}\right) \cdot \delta {\varvec{v}} \,\mathrm {d}a \nonumber \\+ & {} \int _{{\mathcal B}}\left\{ \frac{\partial R}{\partial \dot{\nu }}\,\delta \dot{\nu } - \mathrm {div}\,\left( \frac{\partial R}{\partial {\varvec{D}}}\right) \cdot \delta {\varvec{v}} \right\} \mathrm {d}v. \end{aligned}$$

(29.173)

Here, ${\varvec{n}}$ is the outward unit normal vector on $\partial {\mathcal B}$, and the divergence theorem has been used to obtain (29.173).

With (29.169), (29.170) and (29.173), all elements are now at disposal to apply the principle of virtual power, which requires

$$\begin{aligned} \delta (\dot{\mathcal F}+{\mathcal W}^{s}) + \delta {\mathcal R} = 0, \quad \forall \{\delta {\varvec{v}}, \delta \dot{\nu }\}, \end{aligned}$$

(29.174)

with the understanding that the generalized forces and their power of working remain constant during the variation, see also Chap. 26, Eqs. (26.8) and (26.11). Consequently, with (29.169), (29.170), and (29.173), we have

$$\begin{aligned} \left. \begin{array}{l} \displaystyle {\varvec{X}} - \mathrm {div}\,\left( \frac{\partial R}{\partial {\varvec{D}}}\right) = {\varvec{0}} \\ \displaystyle \xi + \frac{\partial R}{\partial \dot{\nu }} = 0 \end{array} \right\} \quad \mathrm {in}\; {\mathcal B}, \end{aligned}$$

(29.175)

$$\begin{aligned} \left. \begin{array}{l} \displaystyle {\varvec{X}}^{b} + {\varvec{X}}^{s} + \frac{\partial R}{\partial {\varvec{D}}} \cdot {\varvec{n}} = {\varvec{0}} \\ \xi ^{b}+ \xi ^{s} = 0\end{array} \right\} \quad \mathrm {on}\; {\partial {\mathcal B}}. \end{aligned}$$

(29.176)

These equations, with the interpretations (29.170), are the equations of motion in ${\mathcal B}$ and the boundary conditions on $\partial {\mathcal B}$. Note that ${\varvec{X}}^{b}+{\varvec{X}}^{s} = {\varvec{0}}$ and $xi^{b}+\xi ^{s} = 0$ if effects of dissipation are ignored. In this limited case, ${\varvec{X}}^{s}$ and $\xi ^{s}$ are nonzero, since ${\varvec{X}}^{b} \ne {\varvec{0}}$ and $\xi ^{b} \ne 0$. According to (29.170)$_{3,4}$, as stated earlier, their vanishing would severely constrain the potentials W and $\sigma$ on the free surfaces.

Explicitly, the momentum equations (29.175)$_{1}$ take the forms

$$\begin{aligned} \rho \dot{\varvec{v}}= & {} \rho {\varvec{f}} + \mathrm {div}\,\left( (W-\rho ^{2}\sigma ^{\prime }) {\varvec{I}}\right) \nonumber \\&- \mathrm {div}\,\left( \mathrm {grad}\,\nu \odot \frac{\partial W}{\partial \mathrm {grad}\,\nu }\right) + \mathrm {div}\,\left( \frac{\partial R}{\partial {\varvec{D}}}\right) , \end{aligned}$$

(29.177)

from which the stress tensor can directly be read off:

$$\begin{aligned} {\varvec{T}} = (W-\rho ^{2}\sigma ^{\prime }){\varvec{I}}- \mathrm {grad}\,\nu \odot \frac{\partial W}{\partial \mathrm {grad}\,\nu } + \frac{\partial R}{\partial {\varvec{D}}}. \end{aligned}$$

(29.178)

Since W must be objective, it can be written as $W= \hat{W}(\nu , g)$, $g:=\textstyle {\frac{1}{2}}\mathrm {grad}\,\nu \cdot \mathrm {grad}\,\nu$; all three terms on the right-hand side of (29.178) are now symmetric, so that ${\varvec{T}} = {\varvec{T}}^{T}$. In particular,

$$\begin{aligned}&\mathrm {grad}\,\nu \odot \frac{\partial W}{\partial \mathrm {grad}\,\nu } = \frac{\partial \hat{W}}{\partial g} \mathrm {grad}\,\nu \otimes \mathrm {grad}\,\nu \nonumber \\&\qquad \qquad \qquad {\mathop {=}\limits ^{\hat{W}=f(\nu )g}} f(\nu )\mathrm {grad}\,\nu \otimes \mathrm {grad}\,\nu . \end{aligned}$$

(29.179)

Moreover, the dissipative term $\partial R/\partial {\varvec{D}}$ does not need to be restricted to Newtonian behavior, viz.,

$$\begin{aligned} \frac{\partial R}{\partial {\varvec{D}}} = \kappa (\mathrm {div}\,{\varvec{v}}){\varvec{I}} + 2 \mu \left( {\varvec{D}}-\textstyle {\frac{1}{3}}(\mathrm {div}\,{\varvec{v}})\right) {\varvec{I}} \end{aligned}$$

with bulk, $\kappa$, and shear, $\mu$, viscosity. For isotropic behavior the function R in (29.171) takes the form

$$\begin{aligned} R = R(\nu , \dot{\nu }, I_{{\varvec{D}}}, I\!I_{{\varvec{D}}}, I\!I\!I_{{\varvec{D}}}), \end{aligned}$$

(29.180)

where $I_{{\varvec{D}}}, I\!I_{{\varvec{D}}}, I\!I\!I_{{\varvec{D}}}$ are the invariants of the stretching tensor ${\varvec{D}}$. It is convenient to separate the $(\nu , \dot{\nu })$-dependence of R from the dependence on ${\varvec{D}}$. Instead of (29.180), we shall, therefore, write

$$\begin{aligned} R = r(\nu , \dot{\nu })\tilde{R}(I_{\varvec{D}}, I\!I_{\varvec{D}}, I\!I\!I_{\varvec{D}}) , \end{aligned}$$

(29.181)

and then may write

$$\begin{aligned} {\varvec{T}}_{\mathrm {diss}} = r(\nu , \dot{\nu })\left\{ \frac{\partial \tilde{R}}{\partial I_{\varvec{D}}} \frac{\partial I_{\varvec{D}}}{\partial {\varvec{D}}} + \frac{\partial \tilde{R} }{\partial I\!I_{\varvec{D}}} \frac{\partial I\!I_{\varvec{D}}}{\partial {\varvec{D}}} +\frac{\partial \tilde{R}}{\partial I\!I\!I_{\varvec{D}}} \frac{\partial I\!I\!I_{\varvec{D}}}{\partial {\varvec{D}}}\right\} . \end{aligned}$$

(29.182)

For the principal invariants one may deduce (see Hutter and jöhnk [8], pp. 46–49)

$$\begin{aligned}&\frac{\partial I_{\varvec{D}}}{\partial {\varvec{D}}} = {\varvec{I}}, \qquad \frac{\partial I\!I_{\varvec{D}}}{\partial {\varvec{D}}} = I_{{\varvec{D}}}{\varvec{I}} - {\varvec{D}}, \nonumber \\ \\&\frac{\partial I\!I\!I_{\varvec{D}}}{\partial {\varvec{D}}}= {\varvec{D}}^{-1}I\!I\!I_{{\varvec{D}}} = \left\{ {\varvec{D}}^{2} - I_{{\varvec{D}}}{\varvec{D}}+ I\!I_{{\varvec{D}}}{\varvec{I}}\right\} , \nonumber \end{aligned}$$

(29.183)

in which ${\varvec{D}}$ must be symmetric and the Caley -Hamilton theorem has also been used. Substituting (29.183) into (29.182) leads to the following expression for the dissipative Cauchy stress

$$\begin{aligned} {\varvec{T}}_{\mathrm {diss}}= & {} r(\nu , \dot{\nu }) \left[ \phi _{0}{\varvec{I}}+\phi _{1}{\varvec{D}} +\phi _{2}{\varvec{D}}^{2}\right] , \end{aligned}$$

(29.184)

where

$$\begin{aligned}&\phi _{0} = \frac{\partial \tilde{R}}{\partial I_{{\varvec{D}}}}+\frac{\partial \tilde{R}}{\partial I\!I_{{\varvec{D}}}}I_{{\varvec{D}}}+\frac{\partial \tilde{R}}{\partial I\!I\!I_{{\varvec{D}}}} I\!I_{{\varvec{D}}}, \nonumber \\&\phi _{1} = -\left( \frac{\partial \tilde{R}}{\partial I\!I_{{\varvec{D}}}}-\partial \frac{ \tilde{R}}{\partial I\!I\!I_{{\varvec{D}}}}I_{\varvec{D}}\right) , \\&\phi _{2} = \frac{\partial \tilde{R}}{\partial I\!I\!I_{{\varvec{D}}}}. \nonumber \end{aligned}$$

(29.185)

This stress representation possesses the structure of a Reiner –Riwlin fluid. Notice that in case $\tilde{R}$ does not depend upon the third invariant $I\!I\!I_{{\varvec{D}}}$, the dependence on ${\varvec{D}}^{2}$ in (29.184) drops out and the $\phi$-dependences simplify to

$$\begin{aligned} \phi _{0}= \frac{\partial \tilde{R}}{\partial I_{{\varvec{D}}}}+\frac{\partial \tilde{R}}{\partial I\!I_{{\varvec{D}}}}I_{{\varvec{D}}}, \;\; \phi _{1} = \frac{\partial \tilde{R}}{\partial I\!I_{{\varvec{D}}}}, \;\; \phi _{2} = 0,\;\; \mathrm {if}\; \tilde{R}\ne \tilde{R}(\cdot , I\!I\!I_{{\varvec{D}}}).\qquad \quad \end{aligned}$$

(29.186)

Note, moreover, a density preserving granular material ($\gamma =\mathrm {const.}$) is not necessarily volume preserving. Its bulk viscosity is $\mu _{\mathrm {bulk}}r(\nu , \dot{\nu })\phi _{0}$, while its shear viscosity is given by $\mu _{\mathrm {shear}}= r(\nu , \dot{\nu })\phi _{1}$. However, if we ignore both $I_{{\varvec{D}}}$- and $I\!I\!I_{{\varvec{D}}}$-dependences of $\tilde{R}$, then the bulk viscosity vanishes. In this case, the dissipative stress reduces simply to

$$\begin{aligned} {\varvec{T}}_{\mathrm {diss}} = r(\nu ,\dot{\nu })\phi _{1}{\varvec{D}}, \quad \phi _{1}=\frac{\partial \tilde{R}}{\partial I\!I_{{\varvec{D}}}^{\prime }} \quad I\!I_{{\varvec{D}}}^{\prime }= \textstyle {\frac{1}{2}} I_{{\varvec{D}}^{2}}. \end{aligned}$$

(29.187)

Except for the $r(\nu , \dot{\nu })$-dependence this law corresponds to the classical viscous shearing behavior, known in fluid mechanics as dilatant and pseudoplastic behavior (see e.g., [11], pp. 366–371). Popular examples are so-called power laws.

The second equation of motion (for the evolution of the volume fraction) is obtained by combining (29.175)$_{2}$ with (29.170)$_{2}$ and using the definition of ${\mathfrak m}$ given in (29.164). This yields

$$\begin{aligned} \gamma \nu \left\{ \left( \frac{\partial \kappa }{\partial \dot{\nu }}\right) ^{\displaystyle \cdot } - \frac{\partial \kappa }{\partial \nu }\right\} + \frac{\partial W}{\partial \nu } - \mathrm {div}\,\left( \frac{\partial W}{\partial \mathrm {grad}\,\nu }\right) + \frac{\partial R}{\partial \dot{\nu }} = 0. \end{aligned}$$

(29.188)

With the stored energy function $W=\hat{W}(\nu , g)$, $g:=\textstyle {\frac{1}{2}}\mathrm {grad}\,\nu \cdot \mathrm {grad}\,\nu$ and the volume fraction kinetic energy $\kappa$ given by

$$\begin{aligned} \hat{W} = {\mathfrak f}(\nu ) g \quad \text{ and } \quad \kappa (\nu , \dot{\nu })=\textstyle {\frac{1}{2}} k(\nu )\dot{\nu }^{2}, \end{aligned}$$

(29.189)

(both are quadratic in $\mathrm {grad}\,\nu$ and $\dot{\nu }$) one easily deduces

$$\begin{aligned} \gamma \nu \left\{ k(\nu )\ddot{\nu } + \textstyle {\frac{1}{2}}k^{\prime }(\nu ) \dot{\nu }^{2}\right\} - {\mathfrak f}^{\prime }(\nu )g - {\mathfrak f} (\nu )\varDelta \nu + \frac{\partial R}{\partial \dot{\nu }} = 0 . \end{aligned}$$

(29.190)

The Eqs. (29.188) and (29.190) look more general than that introduced in [7],

$$\begin{aligned} \gamma \nu k \ddot{\nu }+\mathrm {div}\,{\varvec{h}} - \gamma \nu f = 0. \end{aligned}$$

(29.191)

It is evident that the coefficient of equilibrated inertia must be a material constant as a condition to reach agreement between the field equations derived by the balance law and variational approaches. This may well be the reason why Eringen in his writings of micro-morphic continua postulates explicitly a conservation law of equilibrated inertia $k_{{\varvec{\mathfrak a}}}= 0$ $({\mathfrak a}=1,\ldots )$ in mixture formulations (see also Chap. 30, (30.5)). On the other hand, it is easy to verify that

$$\begin{aligned} {\varvec{h}}= & {} - \frac{\partial \hat{W}}{\partial g} \mathrm {grad}\,\nu = - {\mathfrak f}(\nu ) \mathrm {grad}\,\nu , \nonumber \\ \\ f= & {} \frac{1}{\gamma \nu } \frac{\partial \hat{W}}{\partial \nu }g = {\textstyle \frac{1}{2}} {\mathfrak f}'(\nu ) \mathrm {grad}\,\nu \cdot \mathrm {grad}\,\nu . \nonumber \end{aligned}$$

(29.192)

Thus, there is a one-to-one correspondence of the differential equations for the fluid velocity ${\varvec{v}}$ and the volume fraction $\nu$ between the variational and balance law derivation of the field equations. Incorporating the dissipation potential ${\mathcal R}$ the variational method has gone beyond the balance law approach, since it delivered explicit formulae for the dissipative stress given in (29.184) and (29.185). Moreover, with the introduction of the dissipative surface potential $W^{s}$, it also delivers in (29.170)$_{3,4}$ boundary conditions for the surface forces ${\varvec{X}}^{b}$ and $\xi ^{b}$.