## Abstract

A theoretical framework of rigid plasticity is presented that is based on optimization and includes frictional dissipation. It has been used in this paper as a foundation for existing and modified models of granular materials consisting of rigid granules. It has the major advantage that it enables yield criteria to be created numerically, which is particularly useful when analytical expressions cannot be found. This framework is constructed by first postulating: (i) a dissipation function that can depend on the current components of stress, but is always homogeneous of degree one and positive definite, (ii) a volume constraint function that is also homogeneous of degree one, and (iii) a balance of the rate of doing work and the rate at which energy needs to be dissipated. A mathematical process similar to the construction of a dual norm in convex analysis then leads to: a flow rule; a single natural representation of the yield surface; and a useful constitutive inequality involving the components of stress and strain rate.

Dissipation functions appropriate for the Drucker–Prager and Matsuoka–Nakai yield surfaces are investigated when a simple dilation rule is the volume constraint. These include a case where an explicit expression for the yield function is not found and, instead, the yield function is found numerically. Such numerical yield functions have been checked graphically against carefully constructed envelopes and found to be consistent with them.

## 1. Introduction

The deformation of a compact assembly of hard, free-flowing particles suffering compressive principal stresses is largely independent of the rate of deformation and exhibits what can usefully be described as a yield strength (Bishop & Henkel 1962). At sufficiently low levels of the principal stresses, the particles can be thought of as rigid, and they therefore have no characteristic stress. Under those conditions, yield can be characterized by a yield surface that takes the form of a cone centred at the origin in principal stress space.

The issue of the direction of plastic flow, however, presents more difficulty. Although for metals, the direction of plastic strain is in the direction normal to the yield surface, this normality-based flow rule is not always appropriate for granular materials, as has been demonstrated by, for example, Lade (1988). In particular, for effectively rigid particles, the actual dilation is much less than would be predicted by assuming that the flow was normal to a conical yield surface. This feature of the behaviour of granular materials is known as non-associated flow—as the association with the yield surface in describing the direction of flow is lost. Instead, it is common for the direction of flow to be modelled as normal to a scalar ‘plastic potential’ function of the principal stresses that can be chosen appropriately to provide a good fit to the available data. A typical example is provided by Gajo & Muir Wood (1999).

As an alternative to working with a yield surface and plastic potential, Houlsby (1982, 1992) started with two thermodynamic potentials and a kinematic constraint (or dilation rule) and used Ziegler's principle of orthogonality (Ziegler 1983) to justify the application of non-associated flow to the elasto-plastic deformation of granular materials. Collins & Houlsby (1997) used the singular case of the Legendre transformation to extend the method. Chandler (1985), also starting with a dissipation function and kinematic constraint, used the theory of envelopes (Courant 1934; Bruce & Giblin 1984), rather than Legendre transforms, to achieve the same results for rigid plasticity and to produce some new yield functions.

These approaches have enabled physically based understanding of the material behaviour to be encapsulated within the mathematical formulation of material models. Such investigations (Chandler 1985; Collins & Houlsby 1997) have shown that the flow is non-associated because frictional dissipation (i.e. dissipation that is dependent on the *current* state of stress carried by the assembly) is present. Models based on these approaches have been used to simulate the behaviour of granular materials and have been able to capture many of the features of granular material behaviour that have been demonstrated experimentally (Chandler 1990; Houlsby 1992; Puzrin *et al*. 2001; Chandler & Sands 2007).

There remain some limitations to this approach, however. Reasonable choices of dissipation functions and dilation rules do not always result in analytical expressions for the yield surface, and this is limiting progress in this area. For example, there is a dissipation function that, with the incompressibility condition, leads to the Matsuoka–Nakai yield criterion (Collins & Houlsby 1997). If, however, the volume strain rate is proportional to the square root of the second invariant of the deviatoric strain rate tensor (a dilation rule commonly used in conjunction with the Drucker–Prager yield criterion), an appropriate expression for the yield surface cannot, as yet, be found. The paper presented here derives the formulations described above using an optimization approach for rigid plasticity. This not only provides a more straightforward method than used previously, but also enables numerical solutions to be obtained when analytical solutions are not available.

Optimization has an important history in the theory of plasticity (see Johnson & Mellor (1973) for applications to metalworking). With the upper and lower bound theorems of associated plasticity as the basis, optimization methods have long been useful in producing approximations to boundary value problems, and progress in this area is still being made. For example, Makrodimopoulos & Martin (2006) have recently made use of second-order cone programming for cohesive soils and Hjiaj *et al*. (2005) have introduced the dual norm to the literature of soil plasticity in the context of the Cam Clay model.

Much less use has been made of other possible optimization methods. Some mixed objective functions have been constructed where a maximum is found with respect to both the stress and strain rate tensors simultaneously. One has been used to find approximate solutions to boundary value problems (Chandler & Song 1990) with a model of powder compaction using associated plasticity. These solutions involved both strain hardening in the material and frictional sliding on the boundaries.

Similar mixed methods, however, have been discovered for non-associated plasticity. In the development of a particular critical state soil model (Chandler 1990) making use of the theory of envelopes, an interesting inequality was found. This could be written in terms of the rate of doing work, the dissipation function and the yield function. This inequality, which reduced to the Cauchy–Schwarz inequality in terms of the stress and strain rate tensors, was augmented by exact penalty functions (Fletcher 1987) in order to provide an objective function. It was then used, in conjunction with the evolution of internal parameters, to find the values of the non-specified components of stress and strain rate in the simulation of triaxial tests.

In this paper, an optimization approach to the formulation of constitutive equations not only confirms earlier results found using the theory of envelopes, but also gives new ones. These show the pivotal role of the theory of optimization in providing an insightful structure for both associated and non-associated plasticity. It also generates, for each combination of dissipation function and dilation rule, a preferential form of the yield function and a potentially useful inequality. The paper concludes with some examples of the method, including when numerical forms of the yield function are produced as well as algebraic ones.

## 2. Underlying optimization problem

Consider a unit volume of material in Cartesian space (*x*_{1}, *x*_{2}, *x*_{3}) subject to tractions on the boundary giving rise to a stress tensor with components *σ*_{ij}. Each point moves with velocity *v*_{i} and the components of the strain rate tensor are given by(2.1)It is useful to partition the strain rate into volumetric and deviatoric parts(2.2)and(2.3)where *δ*_{ij} is the Kronecker delta. Similarly, the components of the stress tensor can be partitioned into the mean stress (*σ*) and the components of the deviatoric stress tensor (*s*_{ij})(2.4)and(2.5)The rate of doing work can then be written as .

To construct a constitutive relationship for rigid plastic flow using the method proposed by this paper, the first step is to postulate a dissipation function, , that is positive definite and homogeneous of degree one in the components of the strain rate tensor. This restriction ensures that energy is dissipated and the resistance to deformation is independent of the strain rate. In this paper, we consider only smooth dissipation functions, treating the important case of non-smooth dissipation functions (Frossard 1983) elsewhere. It is possible to use dissipation functions that are independent of the rate of volume change. In these circumstances, it is necessary to apply a volume, or kinematic, constraint, which can take the form of a dilatancy rule. In order to maintain rate independence, this kinematic constraint, , must be homogeneous and, without loss of generality, is constructed to be of degree one in the components of the strain rate tensor. While a kinematic constraint on the rate of volume strain is always necessary for the dissipation functions used in this paper, that need not always be the case, depending on the choice of dissipation function (Chandler 1985).

For a rigid plastic material, a power balance takes the form of the rate of doing work minus the rate of energy dissipation and can be written as(2.6)So when *ϕ*<0, the rate of doing work is insufficient to supply the demands of energy dissipation. But if *ϕ*=0, the boundary tractions just provide enough power to drive the plastic deformation. In this work, we make the assumptions that: (i) the stress state is initially such that plastic flow cannot occur and (ii) as soon as the stress changes sufficiently, plastic flow will take place.

The question ‘what combinations of *σ* and *s*_{ij} can induce plastic deformation?’ is therefore answered by specifying *σ* and *s*_{ij} and then determining whether there are values of the components of strain rate that make *ϕ*≥0. Clearly, to do this, the values of and that maximize *ϕ*, subject to any constraints, need to be found.

These constraints are important. There may be a kinematic constraint, as mentioned earlier, but the constraint , where is a positive constant, is always needed. Otherwise, if *ϕ* is positive, it can always be increased by a factor, *λ*, simply by multiplying all the components of strain rate by *λ*. Hence, takes the place of what is commonly called ‘the plastic multiplier’ (usually denoted by in the flow rule), and its value is found as part of a specific boundary value problem. This value will take account of any hardening. It will be seen that, like the plastic multiplier, does not appear in the yield function. As will be seen in the examples, for some combinations of dissipation function and kinematic constraint, the constraint is necessary so that and can be treated as independent variables. If this constraint is not applied explicitly, the optimization can sometimes predict a solution that has a deviatoric strain rate with a non-zero trace.

The problem can therefore be summarized as: maximize(2.7)with respect to and subject to(2.8)(2.9)and(2.10)To facilitate this operation, a Lagrangian function can be constructed with *ψ*, *ρ* and *κ* as Lagrangian multipliers(2.11)At a constrained maximum or minimum, where and take the values and , respectively, and the Lagrangian multipliers take the values *ψ*^{*}, *ρ*^{*} and *κ*^{*}, the optimality conditions require that(2.12)(2.13)(2.14)(2.15)and(2.16)As , and *ν* are all homogeneous of degree one, Euler's theorem can be used to help find the values of the Lagrangian multipliers in terms of the basic variables (Lasserre & Hiriart-Urruty 2002). First, multiply equation (2.12) by and then contract the terms in equation (2.13) with and add the results to give(2.17)At the constrained maximum, the constraints are satisfied and this reduces to(2.18)Furthermore, the contraction of equation (2.13) with the Kronecker delta (equivalent of taking the trace) gives(2.19)Equations (2.18), (2.12) and (2.19) reduce to(2.20)(2.21)and(2.22)respectively.

It should be noted that and are functions of *σ* and *s*_{ij}, so the values of the Lagrangian multipliers are also functions of *σ* and *s*_{ij}. In particular, *ψ*^{*} takes the role of a yield function defining a yield surface. It is negative inside the yield surface, where and positive outside the yield surface where . Only when , does(2.23)This expression could be manipulated to give a family of other functions that share the zero-level surface described by equation (2.23), but have a different functional form. An example is(2.24)In the literature of plasticity, it has not been considered important to specify which of these functions is specifically chosen. As represents one of the Lagrangian multipliers of the maximization problem, it takes on a particular significance among all those other functions sharing the zero-level surface. In this paper, will be known as the *natural* yield function and, as will be seen later, it is important in the construction of constitutive inequalities.

The conditions for a natural yield surface to exist are that: (i) for all values of *σ* and *s*_{ij}, has a unique maximum (subject to the constraints) with respect to and , and (ii) this maximum must be zero for the set of contiguous points in *σ*, *s*_{ij} space that form the yield surface. For some choices of dissipation function and kinematic constraint, this may not hold, and the precise restrictions required for it to do so have yet to be established. It is however clear that many dissipation functions in the form of a matrix norm (e.g. Boyd & Vandenberghe 2004) appear to be satisfactory. It is also worthwhile noting that, in cases where kinematic constraints are not required, the procedure presented here parallels closely the process of finding the dual norm from the norm in convex analysis (Boyd & Vandenberghe 2004; Hjiaj *et al*. 2005). Specifically, when the norm is the dissipation function, the dual norm is 1+*ψ*^{*}.

Despite these uncertainties, it is useful to proceed on the basis that a unique maximum can be found such that *ψ*^{*} is a single-valued scalar function of the components of the stress tensor. For some choices of and *ν*, the natural yield function, *ψ*^{*}(*σ*, *s*_{ij}), can be written down explicitly and, in §3, we assume that that is the case, leaving the discussion of an exclusively numerical procedure to §4.

## 3. Mixed maximum principle

Normally, in using constitutive equations, only six independent components are specified from the six independent components of stress and the six independent components of strain rate. The remaining unknown six need to be found. We desire a convenient method for finding the unknown values. A mixed maximum principle with respect to the known and unknown values of both the components of stress and strain rate is developed here with this purpose in mind. A suitable principle can be formulated by first recalling that is a constrained maximum for any values of *σ* and *s*_{ij}. Secondly, non-optimal values of the components of the strain rate, which nonetheless obey the kinematic constraints, are written as (where *α* is the corresponding value of the dissipation function). As is a maximum, it follows that the inequality(3.1)holds for the values of the components of strain rate represented by . This inequality can be rewritten using the definition of the natural yield function,(3.2)and the property of homogeneous functions of degree one to become(3.3)

The requirement that has a specific value is now relaxed. Let , without square brackets, be the components of the strain rate that obey the kinematic constraint but are now allowed to give an arbitrary positive value to . As both and , equation (3.3) can be replaced by(3.4)Furthermore, upon again recalling that *ϕ* is homogeneous of degree one, this becomes(3.5)

In the interests of making connections between plasticity and convex analysis, it is interesting to compare this to the inequality of norm and dual norm (e.g. Boyd & Vandenberghe (2004) for a general treatment and Hjiaj *et al*. (2005) for an application to soil mechanics).

As will be seen in the examples, finding values of stress and strain rates that give the maximum value of this does not ensure that the components of stress lie on the yield surface or that is non-zero. To ensure the first of these, we can impose the yield criterion by subtracting the square of the yield surface as a penalty function to give(3.6)and satisfy the second by imposing that one component of the strain rate tensor has a non-zero value.

The kinematic constraint also needs to be satisfied. When a kinematic constraint for the volume strain rate is written as(3.7)all that is required is to replace by throughout, giving(3.8)In order for to be found, a quadratic penalty term is included to give(3.9)This is an inequality containing the unknown and known values of the components of the stress and strain rate tensors. This has been found previously by another route for a small number of specific cases (Chandler 1990).

We are now in a position to deliver the desired algorithm promised at the beginning of this section. Specifically, if any 6 of the 12 independent components are known, then the remaining unknown 6 can be found by a straightforward optimization algorithm. (Often in problems involving granular material, at least one diagonal term in the stress tensor must be one of the values specified, and for rigid perfectly plastic granules, at least one non-zero component of strain rate must also be prescribed.) Importantly, the penalty terms are exact penalty functions (Fletcher 1987), which means that neither *ρ*_{1} nor *ρ*_{2} need to be large in the context of the problem.

## 4. Numerical alternative

As mentioned earlier, it is not always possible to obtain an explicit form for the yield function and, in that situation, the optimization method described in §3 would require a numerical approximation to the natural yield function. In this section, an alternative approach is investigated, and that is minimizing the sum of the squares of the deviations from the flow rule. To see how this might work, firstly rearrange equations (2.12) and (2.13),(4.1)(4.2)

As we wish to find values of the components of stress and strain rate consistent with yield, the Lagrangian multipliers *ρ*^{*} and *κ*^{*} can be calculated from equations (2.21) and (2.22) with the other Lagrangian multiplier (*ψ*^{*}) set to zero. All the terms on the left-hand side of equations (4.1) and (4.2) can then be evaluated if the components of the stress and strain rate tensor are either prescribed or guessed. With *ρ*^{*} determined by equation (2.21), the left-hand side of equation (4.1) can be seen to be identically zero.

Both sides of equation (4.2) can be contracted with themselves and added to the squares of both sides of equation (4.1) to give(4.3)The right-hand side is only zero when *ψ*^{*} is zero as(4.4)is positive unless and are both zero and, as is positive definite and homogeneous of degree one, this can only occur when and are *all* zero.

As *ψ*^{*} is only zero when the components of stress are on the yield surface, the left-hand side of equation (4.3) is also only zero when the components of stress lie on the yield surface. Additionally, a minimum value of zero for the left-hand side implies that the deviatoric part of the flow rule must also be obeyed. The volumetric kinematic constraint can be imposed using a quadratic penalty function as before, completing the solution.

## 5. Examples

For illustrative purposes, two examples of dissipation functions are examined. One leads to the Drucker–Prager yield function (Drucker & Prager 1952), if used with the simple dilation rule specified below. The other leads to the Matsuoka–Nakai yield function (Matsuoka & Nakai 1974), if used with an incompressibility constraint, but results in a different yield surface (obtained numerically and checked graphically), if used with the same simple dilation rule.

### (a) Drucker–Prager

In order to reproduce the Drucker–Prager yield surface, we start with the dissipation function(5.1)and the dilation rule (kinematic constraint) of , where *ν* is a constant. This gives a Lagrangian function(5.2)leading to the optimality conditions,(5.3)(5.4)(5.5)(5.6)and(5.7)Using the methods presented in §2, the Lagrangian multipliers can be found explicitly as(5.8)(5.9)and(5.10)with representing the rate of doing work at the optimum. For this dissipation function and dilation rule, the deviatoric rates of strain can, noting that the mean stress is negative and using , be written explicitly as(5.11)providing a flow rule.

Using the above values for the Lagrangian multipliers and components of deviatoric strain rate allows the rate of doing work to be written in terms of the components of stress and the dissipation function,(5.12)Substitution of this in equation (5.8) gives(5.13)which, for the case of compressive stresses, leads directly to an explicit form for the natural yield function,(5.14)The mixed extremum principle (equation (3.9)) reduces to the Cauchy–Schwarz inequality augmented by quadratic penalty functions(5.15)

Effectively, the Cauchy–Schwarz inequality becomes an equality when the deviatoric stress and deviatoric strain rate tensors are proportional with a positive constant of proportionality; and the penalty functions are zero if both the yield condition and the dilatancy rule are obeyed.

### (b) Matsuoka–Nakai with no volume change

The restricted case where we assume *a priori* that the principal axes of the stress and strain rate tensors coincide is considered in this section. In this case, in order to reproduce the Matsuoka–Nakai yield surface in terms of principal stresses, we start with the dissipation function(5.16)

This is a more compact version, used by Nixon (1999), of the dissipation function given by Houlsby (1992) and is positive definite when the principal values of stress all have the same sign. While the case, where there is zero volume change, will be shown to correspond to the Matsuoka–Nakai yield function (see the outer function in figure 1), the option of dilation is included in the following for later use. The seemingly most straightforward choice for a dilation rule was that used in the Drucker–Prager case, but with *Δ* replacing , as shorthand, to give(5.17)This choice is consistent with the experimental work of Lanier & Zitouni (1989) and the numerical results of Thornton (2000). With dilation included (and *I*_{1} replacing (*σ*_{1}+*σ*_{2}+*σ*_{3})), the Lagrangian function can be written as(5.18)leading to the optimality conditions,(5.19)(5.20)(5.21)(5.22)and(5.23)with the subscript *I* taking value 1, 2 or 3 and no sum implied on repetition. For this dissipation function and dilation rule, and using the methods presented in §2, the Lagrangian multipliers can be found explicitly as(5.24)(5.25)and(5.26)and the principal deviatoric rates of strain can be written as(5.27)We require that , and for the special case when *ν*=0, this reduces to(5.28)which can be solved for *κ*^{*} in terms of the principal stresses to give(5.29)This allows equation (5.27) to become a flow rule where principal deviatoric strain rates are expressed in terms of , *ψ*^{*} and the values of the principal stresses. This, in turn, allows the rate of doing work to be written as(5.30)which leads directly to an explicit form for the natural yield function when combined with equation (5.24) to give(5.31)The inequality(5.32)follows from equation (3.9).

This seemingly unlikely inequality is clearly true if(5.33)Noting that , by the inequality of the arithmetic and harmonic means, and recalling *D*_{MN}>0, both sides of equation (5.33) can be squared. Subsequent rearrangement gives(5.34)The right-hand side can be expanded and shown to be(5.35)which is clearly never positive if the principal stresses all have the same sign.

### (c) Matsuoka–Nakai with dilation

The case in §5(*b*), where *ν*≠0, provides an opportunity to exploit the least squares method described in §4. The yield surface for a specific value of the pressure is shown in figure 1 as a line, where *ψ*^{*}(*σ*_{1}, *σ*_{2}, *σ*_{3})=0 augmented by vectors representing the direction of the deviatoric strain rate, for two combinations of *ω* and *ν*.

Each point on the curves was found by: (i) specifying the direction of a vector in principal stress space lying in the constant pressure plane out from the point where *σ*_{1}=*σ*_{2}=*σ*_{3}, (ii) choosing a magnitude for this vector and then finding the least squares minimum of the left-hand side of equation (4.3) with respect to the principal strain rates, and (iii) finding the magnitude of the vector that gives *ϕ*=0. The directions of the strain rate vectors then follow from the flow rule.

These results were confirmed numerically. When *ν*=0, this process gives identical curves to those obtained from equation (5.31). For non-zero dilations, the yield surfaces can be compared with those produced by constructing envelopes (Chandler 1985), and an example of this is shown in figure 2, which can be compared with the inner function shown in figure 1. Figure 3 shows the effect of negative and positive dilations.

## 6. Concluding remarks

For models representing granular materials using smooth dissipation functions and dilation rules, the optimization method (which has much in common with the norm : dual-norm transformation in convex analysis) can be used to establish flow rules and (what we have called) natural yield functions. This approach can provide benefits for the mathematical formulation of rigid plastic constitutive equations, as well as providing a flexible numerical implementation allowing a diverse range of dissipation functions and kinematic constraints to be used. Specifically:

The optimization approach presented here allows for the development of yield surfaces and flow rules that are guaranteed to be thermodynamically consistent even when no analytical form of the yield surface can be found. It follows that physical insight can be incorporated into the dissipation function and dilation rule without regard to the mathematical difficulty of formulating an explicit expression for the yield function. The ability to formulate numerical representation of natural yield functions, in this way, provides a workable route to developing a very wide range of yield functions, some of which might better reflect the behaviour of real materials.

The optimization approach generates specific yield functions which we have called the natural yield functions. This is the only form that guarantees the ability to develop constitutive inequalities as shown in §3. Methods using either envelopes or the singular case of the Legendre transform only produce yield surfaces. These are represented by the zero-level surface of scalar functions whose values remote from the yield surface have a degree of arbitrariness. These scalar functions may not have a form that is suitable for the development of useful inequalities.

The constitutive inequalities that can be constructed within this framework can be used for the approximate solution of boundary value problems. For example, Chandler & Song (1990) used such a method to simulate the compaction of granular material in a parallel-sided bin under the action of gravity.

Although the Lagrangian multiplier,

*κ*^{*}, associated with the constraint is zero in the Drucker–Prager model, it is non-zero in the Matsuoka–Nakai case. The inclusion of this constraint directly at the start avoids any need to insert it later in an ad hoc manner (Collins 2003).The models used as examples in this paper use only two adjustable parameters to encapsulate the features of an admittedly limited range of material behaviour. The formulation, however, could be extended to include terms that reflect particle deformation and these have been incorporated into the dissipation function by Chandler (1985) for the Drucker–Prager case and by Collins (2003) for the Matsuoka–Nakai case.

## Footnotes

- Received February 12, 2007.
- Accepted May 4, 2007.

- © 2007 The Royal Society