Particle mechanics approach to continuum constitutive modelling

A new type of constitutive model is proposed for soils, based on three foundations. The first foundation involves a new approach to Cauchy’s concept of stress, based on a way of organizing particle mechanical information. This leads to a description of effective stress that is consistent with Terzaghi’s principle but contains more information than a tensor. The second foundation involves the load-spreading and displacement-spreading effects of particles in an aggregate. The third foundation takes account of recent proposals in the thermomechanics of granular media. A simple model is then proposed based on expectations inferred from published data of asymptotic proportional straining. Validation calculations are presented showing realistic predictions for monotonic and cyclic plasticity, critical states, dilation, induced anisotropy and non-coaxiality as a result of principal axis rotation. The model suggests that a yield envelope, plastic flow and critical states are not fundamental properties but emergent ones. The paper concludes with a brief discussion of possibilities for further developments.


Particle mechanics approach to continuum constitutive modelling
x, y, z relating to the corresponding coordinate directions y associated with interparticle contact forces in direction y 1, 2, 3 relating to the corresponding coordinate directions Introduction Background A constitutive model for soils for engineering purposes is a relation among stress, strain, their increments and their histories. It may also include strength and bifurcation. Its aim is to provide engineers with tools to understand better how to classify soils, characterise and measure their properties and make predictive or forensic calculations for the responses of soil bodies to various scenarios of applied loads, displacements or other perturbations. Terzaghi's (1923Terzaghi's ( , 1936 Principle of Effective Stress forms the basis of modern models of fully saturated soils (Clayton et al., 1995). In some contexts, soil can be modelled as isotropic and linear elastic, and useful practical solutions are provided by Davis and Selvadurai (1996), Poulos and Davis (1974), and others. However, soils are not isotropic, not linear and not elastic (Atkinson, 2000;Burland, 2012;Graham et al., 1988;Jardine et al., 2004;Mitchell and Soga, 2005;Oda et al., 1985). At the very least, a modern model would be expected to incorporate some aspects of elastoplasticity, in which a soil element will not fully return to its original shape when a stress or load increment is applied and then removed.
The first credible elastoplastic constitutive models for soils were proposed about 50 years ago. Among those, Drucker and Prager's (1952) model, Drucker et al.'s (1957) approach to work-hardening and the Cam-Clay models (Roscoe and Burland, 1968;Schofield and Wroth, 1968) are in routine use in finite-element codes today. Their concepts of plasticity, hardening, dilation and critical states underlie modern codes of geotechnical practice. These early models omitted effects of cyclic loading, anisotropy and principal axis rotation, which are now regarded with increased importance (Casagrande and Carillo, 1944;Scott, 1985Scott, , 1987. Data of asymptotic behaviour in compressive proportional straining, including but not limited to one-dimensional (1D) compression, also provide challenges (e.g. Chu and Lo, 1994;Gudehus, 2011;Ibraim et al., 2010;Topolnicki et al., 1990).

Aims of this paper
The present paper proposes an approach that is consistent with this new focus, while avoiding much of its potential complexity.
Ideally, to develop a model based on particle mechanics, one might hope to start with ideas about particle properties and interactions, then deduce consequences for continuum mechanical behaviour and then test predictions against data. The present paper adopts a different, iterative approach. Ideas of particle mechanics are used to interpret data, and the interpretations are then used to develop the ideas further, the data are used again and so on.
The result is significantly different from previous concepts. The differences can be illustrated by considering the following three different approaches to the constitutive modelling of granular media.
(a) Tensor-based approaches are exemplified by the Cam-Clay models and their implementation in the CRISP finite-element program (Britto and Gunn, 1987). Constitutive equations relate components of the stress tensor, stress increment tensor, strain increment tensor and/or invariants thereof, and a description of state such as specific volume. The maths involves a small number of tensor quantities and invariants. (b) Numerical (discrete-element) modelling of particle assemblies is described by Cundall and Strack (1979), Nguyen et al. (2015), Tordesillas et al. (2009), and others. Stress is related to interparticle forces. Equations involve these discrete forces, their increments and incremental changes of position and orientation and shapes of discrete particles. (c) In the approach proposed in this paper, the relation between stress and interparticle force is used to develop a continuum representation of the forces. The representation is a scalar function of directions in physical space; the directions may be represented, for example, by latitude and longitude. Other scalar functions are also developed, and the maths is the maths of functions and of operators on functions.
Foundation 1: organising particle mechanical information

Preamble
In their Cam-Clay model, Schofield and Wroth (1968) considered a 'random aggregate of irregular solid [soil] particles of diverse sizes which tear, rub, scratch, chip, and even bounce against each other during the process of continuous deformation'. The idea of a 'random' aggregate suggests an absence of organisation at the microscopic scale.
However, experiments by Al-Tabbaa (1984), Arthur and Menzies (1972), Arthur et al. (1977) and others showed that even initially isotropic soils could develop induced anisotropy. This had been defined by Casagrande and Carillo (1944) as anisotropy associated with a history of shearing strain. Baker and Desai (1984) used the phrase 'isotropic state', and considered how a soil state could change from isotropic to anisotropic as a result of strains. Calladine (1971) showed that it was possible to understand macroscopic plasticity in terms of microstructural properties. Concepts of 'fabric' were developed by Oda (1972Oda ( , 1978 and others, and computationally intensive approaches to soil modelling are now available in the discrete-element method (O'Sullivan, 2013;Radjai and Dubois, 2011). That method seeks to replace continuum mechanics with particle mechanics. The present paper seeks instead to replace particle mechanics with a more appropriate continuum mechanics.

Previous approaches to stress
Cauchy's (1822/3) general theory of stress underpins Terzaghi's (1923Terzaghi's ( , 1936 principle and is the basis of modern continuum mechanics (Maugin, 2013;Truesdell, 1968Truesdell, , 1977. His theory preceded Rankine's (1857) work and followed Coulomb's (1773) Memoir, which is generally regarded as having laid the foundations of modern soil mechanics (Heyman, 1972(Heyman, , 1997. In Figure 1, a closed orientable surface surrounds a volume of continuum material. Tractions act on the surface. For 'uniform' conditions or for a small enough volume (strictly an infinitesimal one if conditions are non-uniform), the traction vector t can be obtained by applying a linear operator r to the outwards unit vector n normal to the surface 1.

= t n σ
The operator r is the total stress tensor, expressed as a matrix so that t is the product of a matrix and a vector. It is a 'Cauchy' tensor in the sense that it relates directly to tractions on a unit volume in a current configuration, independently of any reference configuration.
A special case is a uniaxial stress tensor Sy of magnitude Sy and with the major principal axis in the direction of a unit vector ( ) , , x y z n n n ψ ψ ψ ψ = T n . Appendix A shows that Sy can be expressed as follows: x z y z n n n n n n n n n n n n n n n n n N Ny is proportional to the dyadic product of a vector with itself, and is a tensor. Since ny is a unit vector, the sum of the diagonal elements of Ny is 3. The trace tr(Sy) (the sum of its diagonal terms) is Sy. Terzaghi (1923) recognised that pore stress plays no role in the stress-strength behaviour of many fully saturated granular soils. The tensor of effective stress takes account of this by subtracting the pore stress, giving 4.
where t¢ represents compressive tractions that equilibrate interparticle contact forces, excluding forces due to pore fluid pressures, and r ¢ is the effective stress tensor obtained by subtracting pore stress from the total stress tensor. Figure 2 shows the sign convention for positive directions of tensor stress components used in this paper.

A measurement frame for directions
This paper seeks a way to organise information about interparticle forces in terms of their directions. For that purpose, it is helpful to start by developing a way to describe direction. Figure 3 shows a sphere of surface area a, centred on some point O. The arrow labelled y represents a direction. A curvilinear grid could be marked on the sphere, for example, lines of latitude and longitude, and the direction could be specified by coordinates on the grid. The direction can also be specified as a unit vector ny in relation to a coordinate frame {x,y,z} or in other ways. A direction is independent of the way it is specified, so a discussion of y should be possible without saying which way is used.
Like the coordinate frame {x,y,z}, the sphere is a concept. It can act as a fixed reference frame for directions just as well as a set of fixed axes {x,y,z} can be used for position. Let fy be some quantity that depends on direction y, and let á fyñ denote its average taken over all directions. The average can be calculated by integrating over the surface of the sphere and dividing by the surface area.
This process of integration can be considered as follows. Suppose the spherical surface of area a is divided into many small convex patches, such as patch P with area day. The ratio day /a is 4p times the solid angle subtended by the patch P. The average is obtained by multiplying the value of fy by the area day of a patch, adding the results for all patches, dividing by the total surface area of the patches that have been included in the integral, and then taking the limit for finer and finer patchworks. This gives In practice, the average can be approximated by carrying out calculations for a sufficiently fine patchwork, and an example is given in Pande and Sharma's (1983) description of numerical procedures in their multilaminate model.

Effective stress and interparticle forces
Cauchy's original definition of stress 'at a point' can be said to refer to an infinitesimal volume of material. Lambe and Whitman (1979) commented that a 'large' point must be considered when this concept is applied to soils -large enough to contain many soil particles. The width of a rupture surface is sometimes considered to be ten or so median particle diameters (e.g. Roscoe, 1970;Sadrekarimi and Olson, 2010), so a cube with side dimensions of 30-40 diameters, containing several tens of thousands of particles, might need to be examined in order to check that a rupture surface was not passing through a ʻlarge pointʼ. Figure 4 shows such a ʻlarge pointʼ. A wavy surface is constructed around the point, passing through voids and interparticle contacts but not particles. The surface as sketched is nowhere more than one particle size distant from the surface of a perfect sphere. It would be the equivalent of the wavy surface and Euler cut concept of previous explanations of effective stress (de Boer and Ehlers, 1990).
In Figure 3, a direction y ¢ might be said to be ʻclose toʼ y if the arrow for y ¢ intersects the same patch P as the arrow for y. In Figure 4, consider the subset of interparticle forces that act on the wavy surface from outside in directions y ¢ that are close to a given direction y. Such a subset is indicated by arrows. An approximately uniaxial stress tensor dr ¢ can be defined by forming a crosssectional cut that is normal to direction y, summing the forces of the subset, and dividing by the cross-sectional area of the cut. Other forces act in other directions, and the whole soil body is assumed to be in equilibrium.
The major principal value of this uniaxial stress will be proportional to day /a for small values of this ratio, so a magnitude ψ σ ′ can be defined as  In principle, all quantities except the stress ψ σ ′ can be known independently, so this serves as a definition of ψ σ ′ . To be consistent with the idea of a continuum, the limit should perhaps be taken as the area tends to zero, but such a concept is obviously constrained by the same considerations as those mentioned by Lambe and Whitman (1979).
The stress ψ σ ′ thus represents a simple, continuum mechanical way to represent interparticle forces. In general, it will be a function of direction y for any given soil state, and the function will change during a stress-strain process. Summing over all patches, then using the familiar concept of taking limits, in this case as day /a tends to zero, gives a tensor r ¢ 7.
The process of taking limits represents a conceptual replacement of the particulate material with an equivalent continuum. It would be expected to reduce to zero any effects of non-uniaxiality due to small differences between y ¢ and y.
All the effective forces acting on the wavy surface have been included in the summation and integration, and Appendix B confirms that r ¢ is indeed a tensor. The present paper interprets it as the familiar Cauchy effective stress tensor. The quantity ψ σ ′ is here called 'outer stress'. It is a function of direction y for a given soil state.
The above equation shows that the six independent components of the tensor are differently weighted averages of this function. The function contains more information than the tensor, and many different functions can produce the same values for the six weighted averages. Taking the trace of both sides of Equation 7, and using Equation 3, gives Hence the average outer stress always equals the familiar mean normal effective stress p ¢.

Continuum descriptions of work and strain
Work calculations in continuum mechanics are described by Chadwick (1999), Spencer (1980) and others in terms of a deformation gradient tensor F, and the calculation for soils has been explored by Houlsby (1979Houlsby ( , 2005 and others. The gradient describes an overall deformation from a reference configuration to a current configuration. It is exact for large as well as small strains applied to a continuum. Let dF be the change of deformation gradient associated with an increment of stress-strain behaviour in which the current configuration changes to a new current configuration. It is well known that the product dFF -1 can be separated into a symmetric part representing incremental strain and an antisymmetric part representing incremental rotation. If these parts are divided by an appropriate time increment, the results are the deformation (strain) rate tensor and the spin (or 'vorticity') tensor, respectively (Spencer, 1980).
In the absence of fluid flow, the effective work dW¢ done on the soil per unit particle volume is given by 9.
T 1 d tr( ) The Cauchy stress tensor is symmetric, but the incremental rotation is antisymmetric, implying that no work is associated with incremental rotation. Using Equation 7 to substitute for the   ( ) The first result shows that dW ¢ is the average of increments dW ψ ′ associated with directions y. The second shows that dW ψ ′ can be calculated as the product of the outer stress with the negative of an increment dVy. This is comparable to the expression -p¢dV for work associated with an increment dV of specific volume under hydrostatic effective stress p¢.
By expanding the right-hand side of Equation 12 and expressing dFF -1 in terms of incremental strains and incremental rotations, dVy is found to be independent of incremental rotation. It is also independent of reference configuration, because dFF -1 is. However, this does not imply that a reference-independent quantity Vy exists. Instead, a different reference-independent quantity will be developed later, denoted with a prime as V ψ ′. By expanding the matrix product in the trace in Equation 12, the average ádVyñ is found to be equal to the increment dV of overall specific volume V.
The name 'outer work (per unit particle volume)' is proposed for dW ψ ′, and 'outer (specific) volume' for Vy. Both quantities are functions of y, just as outer stress is a function of y. They can depend on the stress-strain history of the soil. In general, Vy depends on reference configuration, but dVy does not. For simplicity, the words in brackets are usually omitted in the remainder of this paper.

Discussion
The outer stresses and outer incremental volumes are continuum variables that can be used to construct new models of the constitutive behaviours of soils. They are directly associated with particle mechanics. Outer stresses contain a two-dimensional (2D) infinity of information, which is far more information than in a stress tensor. This suggests that the maths for outer stresses is the maths of functions, and so of real Hilbert spaces, described by Akhiezer and Glaxman (1993), Stone (1932) and others.
Consider any description of micromechanics that can be expressed in terms of a function fy of the directions of interparticle forces.  ( ) fq is the root mean square deviation from the average. An example would be an outer deviatoric stress q σ ′ calculated by replacing the symbol f with s′. There is no general relation between this and the deviator of the stress tensor r ¢, because Equation 7 implies that no element of the stress tensor involves a squared stress. The information represented by q σ ′ is absent from the tensor.
In classical continuum mechanics, Cauchy stresses are often associated with small-strain processes, and Piola-Kirchoff stresses are used to model large strain effects (Spencer, 1980). Appendix C of the present paper shows that outer quantities are perfectly adequate for large-strain processes, provided the following is used for the incremental form of Equation 7 15.
On the right, d ψ σ ¢ is an increment due to constitutive behaviour, while d ψ σ ¢ G is due to incremental changes of geometry. A similar result applies for increments of tensor f. The equation is independent of any arbitrary assumptions about a reference configuration.
It remains to demonstrate that these proposals can be usefully incorporated into a constitutive model. The next two sections develop some general concepts for that task, and a simple constitutive model is then proposed and validated.

Preamble
Many existing models use composite stress and strain quantities, formed from dimensionally consistent combinations of measurable quantities. An example is the triaxial deviator stress, which equals measurable axial less radial stress. The deviator is convenient, relates easily to shear strength and can be related to fundamental aspects of the elastic responses of isotropic soils (Schofield and Wroth, 1968).
Composites can also be formed from outer quantities, and the present paper calls them 'inner' quantities. They are interpreted as descriptions of what happens inside a body of soil (Figure 4), in contrast to the outer quantities that were defined in terms of interparticle forces acting on the surface of a body. The treatment presented below relates to reversible aspects of particle interactions. Inelastic processes are considered later.

Inner specific volumes
Figure 5 sketches some of those interactions. The force F1 is equilibrated by forces F2 and F3, which are in different physical directions. In effect, the particle spreads the load F1. Of course the particle can also be interpreted as providing a way to spread F2 and F3. This 'load-spreading' role is matched by a displacementspreading role, because a movement at contact 1 necessarily entails movements at other contacts, at least if the particle is relatively rigid.
The displacement-spreading phenomenon suggests that a useful measure of geometry change inside the sphere can involve combinations of increments from different directions y. The present paper explores combinations * dV ψ that can be calculated as follows 16.
where fy is some function of direction and soil state. * dV ψ is herein called an 'inner (specific) volume increment', with the word 'specific' omitted for simplicity. The symbol y inside the averaging operator á¼ñ is a dummy variable over which an integration is done. The average represents an influence, on direction y, of other directions.
A meaningful function fy might be a material constant or a description of fabric that depends on direction. Or it might depend on invariants like relative density, or on overconsolidation ratio (OCR), or it may be more complicated. In some increments and some directions, Equation 16 can imply that the inner volume increment may be dilative ( * 0 dV ψ > ) while the outer increment is contractive (dVy < 0), or vice versa. However, taking the average of both sides gives 17.
The á fy dVyñ terms cancel, and the right-hand side reduces to ádVyñ.
Hence for inner volumes defined as in Equation 16, * dV ψ will always equal the increment of overall specific volume.

Work conjugacy and inner stress
The idea of work conjugacy is well established in continuum mechanics and in soil mechanics (Peric et al., 1990;Schofield and Wroth, 1968;Spencer, 1980). The following proposals seek a way to apply work conjugacy to outer and inner quantities. Let * ψ σ ′ represent inner effective stresses that are work-conjugate with * dV ψ . More precisely, let * ψ σ ′ be a function that is work-conjugate with * dV ψ in the sense that work equations for inner quantities have the same form as Equations 10 and 11 18.
This does not imply that the inner increment * dW ψ ′ equals the outer increment dW ψ ′, only that their averages are the same. Using Equation 16 to substitute for * dV ψ , and then comparing the result with Equations 10 and 11, gives 20.
The average of the last term on the right equals * ψ σ ′ multiplied by á fy dVyñ, and so is the same as the average of * f dV ψ ψ ψ σ ′ . Hence the above equation simplifies as follows Particle at end of increment Particle in its position at start of increment  Now the quantities that are being averaged on the left-and righthand sides can only differ by a quantity whose average is always zero. Hence

22.
( possible geometric increments. A stress that cannot do work may in principle influence behaviour (Houlsby, 2005;Peric et al., 1990). However, the present paper assumes that ** ψ σ ′ is always zero.
A feature of the above result is that the average outer stress ψ σ ′ is not in general equal to the average inner stress * ψ σ ′ . However, there are two special cases. If fy represents a material constant or an invariant (independent of y, although not necessarily of stressstrain history), then The second case is if the soil happens to be in an isotropic state with all directions equivalent and none special. In this state fy must be the same for all directions y ; otherwise the state would not be isotropic. For both special cases, if the workless stress is zero, , and both equal the mean normal effective stress p¢.

The no-tension condition
The no-tension condition for uncemented soils is part of Schofield's (2005) soil map. So, if the principal values of the tensor on the lefthand side of Equation 7 cannot be tensile, then something on the right-hand side of that equation must enforce this.
One possibility is that the same no-tension condition applies to the outer stresses in all directions. This would be consistent with the concept of Inspection shows that one way to ensure that this applies for all possible inner stresses would be to ensure that fy is always between 0 and 1, and that the inner stresses also satisfy a no-tension condition.
The no-tension condition is also related to invertibility, as follows. Putting ** 0 ψ σ ′ = in Equation 22, dividing both sides by (1 -fy), and then taking averages gives 24.

Discussion
The method of work conjugacy does not consider loss of work or energy due to irreversibility. For this reason, the function fy is interpreted as a property that relates to reversible interactions in the soil aggregate, and the motions sketched in Figure 5 are restricted to reversible changes in particle arrangements. A different treatment is needed for irreversible processes such as frictional slip at interparticle contacts or particle breakage. This will be considered later in the paper.
The treatment above is not general but is adequate for present purposes. For example, the treatment started with the geometric Equation 16 and ended with the stress Equation 22, but it is equally possible to start with a stress equation and end with a geometric one. This gives a slightly different result for workless quantities.
For simplicity in calculations later in this paper, fy is assumed to be a material constant (denoted as a), and the no-tension condition then implies that a is between 0 and 1.

Foundation 3: thermodynamics Preamble
The importance of thermodynamics has long been recognised in theories of plasticity and elastoplasticity (Hill, 1950;Schofield and Wroth, 1968). Lemons (2013) describes thermodynamics as 'the science of macroscopic objects composed of many parts' -which clearly relates to soils. Specific assumptions for the 'thermomechanics' of solids were proposed by Ziegler (1983) and others.
Thermodynamics can also be described as the science of what is possible and what is not. For this reason, this part of the paper focuses on looking for general opportunities. It is not always possible to know what the physical interpretation of an opportunity may be. Collins and Houlsby (1997) showed that, for an elastoplastic granular soil under isothermal conditions, the rate of work done on an element is the sum of a reversible rate of change of Helmholtz free energy and an irreversible, strictly dissipative part. In incremental terms, this can be expressed as

An existing work equation
where H¢ is the Helmholtz energy of the soil per unit particle volume, and d dW ′ is an increment of work dissipated (or energy irreversibly lost), again per unit particle volume. Irreversibility implies that d dW ′ is non-negative for the sign convention used in this paper.
The increment of Helmholtz free energy was taken as the sum of the increment of internal energy and the increment of reversible heat flow out of the soil. Elasticity is typically associated only with changes of internal energy (Spencer, 1980), and unrecovered strain can cause the material energy to change even if the overall change of stress is zero. For these reasons, dH¢ and d dW ′ are not in general the same as elastic and plastic increments.

An inner work balance equation
Proposals earlier in this paper considered reversible aspects of interactions between different physical directions. To provide an opportunity to account for irreversible aspects, the present paper proposes to adapt the above equation as follows. For a given direction y, the work on the left-hand side is replaced by the sum of external work and energy that has been transferred internally from other directions, as follows where all quantities are for direction y and are per unit particle volume, dE ψ ′ is an increment of energy received from all other directions, H ψ ′ is a Helmholtz free energy and d dW ψ ′ is a nonnegative incremental dissipated work. This equation might be interpreted physically as a balance between inputs to direction y, on the left-hand side, and uses of those inputs, on the right-hand side.
The receipt dE ψ ′ is considered to be the result of energy transfers from other directions. A particle mechanical explanation might be as follows. Suppose that slip occurs at a particle contact where the force is in a given direction. The slip will alter the equilibrium of the particles involved, changing the other forces on those particles, which in general will be in different directions. These changes will then propagate to neighbours of these particles, and to neighbours of neighbours, affecting more interparticle forces in more directions. If interparticle force relates to stress, then the energies in different directions will be changed compared to what would have happened in the absence of slip.
This explanation suggests that the magnitudes of energy transfers that cause the receipts might be related to dissipated work, although the transfers will not themselves be dissipative. Since a transfer of energy from one direction to another is the negative of the transfer from the second to the first, the average receipt is zero Taking averages of both sides of Equation 26 then shows that Collins and Houlsby's (1997)

Testing the proposed work balance
Consider a soil sample that is prepared in an isotropic state and then subjected to purely hydrostatic effective loading and unloading. No direction is special, so all inner quantities equal their averages (assuming that there is no internal directional instability). All inner volumes * V ψ would evolve as the overall specific volume V, and all outer and inner stresses would evolve as the mean normal effective stress p¢.
Changes of energy could then be measured by measuring work done and dissipative losses and using Equation 25. These would be changes of an isotropic energy iso H ′ , which would be some function H of p¢, V and perhaps of variables describing isotropic history It could then be consistent to assume that in more complex stressstrain processes, the energy H ψ ′ associated with direction ψ can be obtained using the same function, but using variables associated with direction ψ 29.
where V ψ ′ is a specific volume for direction y. In the model developed later, V ψ ′ is not the same as the outer volume Vy or the inner volume * V ψ . Instead it is associated with hardening and is herein called the hardening volume in direction y. Figure 6 shows results for isotropic compression and swelling by Roscoe and Burland (1968) of isotropically reconstituted kaolin clay. The large dots represent data points that have been read and replotted from Roscoe and Burland's paper. The curves represent calculations using the model presented later, but for the present discussion, it is convenient to consider these as representing the data too.
The kaolin responses were bounded by curve ABG representing the asymptotic isotropic compression curve. On unloading from B, then reloading, a hysteresis loop occurs. At the low-stress end of the loop, the unloading segment CD is steeper than the reloading segment Particle mechanics approach to continuum constitutive modelling Dean DE, indicating that the modulus in reloading from D is higher than in unloading to D. On reloading, the (p¢,V) state of the soil returns asymptotically to the bounding curve. Many soils show similar features of cyclic hysteresis and of asymptotic compressive responses Been, 2000, 2006;Pestana and Whittle, 1995).
It is convenient to analyse the data by first assuming that the function H depends only on p¢ and V. If this is the case, the increment of energy in an incremental process is given by The first partial differential is carried out at constant stress, and the second at constant specific volume. Hence using Equation 25, noting that incremental work is -p¢dV for these isotropic increments, and rearranging The left-hand side represents a modulus under isotropic conditions. On the right-hand side, energy is expected to increase with increasing stress, so the denominator will be positive. Also, the incremental dissipated work is non-negative. So this result says that the modulus in compression (dV negative) can never be greater than the modulus in swelling (dV positive) for the same soil state. This directly contradicts the data. In Figure 6, the modulus in compression from D to E is greater than the modulus in swelling from C to D.
It follows that there must be at least one more variable involved, apart from p¢ and V. The extra variable cannot be preconsolidation stress, since that is a function of p¢ and V if the bounding curve ABG is a soil property. In similar results for a sand, the variable could not be relative density or relative dilatancy index, since those too are functions of one or both of p¢ and V (Bolton, 1986). Nor could it be the 'state parameter' described by Been (2000, 2006), since that is a function of stress and void ratio e = V -1 for a given soil.
Several authors discuss the idea of 'back-stress'; see for example Collins (2003), Gudehus (2011) and others. The present paper proposes something similar, as follows.

A way forward
The present paper assumes that the missing variable is a stress that depends on stress-strain history, although it would not be a preconsolidation stress. By subtracting this missing stress from p¢, another new stress is obtained. The two new stresses sum to p¢ by definition.
Given the relations among tensor, outer and inner quantities, a consistent assumption would be that inner stress is also the sum of two stresses. Denoting these with subscripts P and Q gives 32.
Although this equation has been developed from isotropic data, it is proposed as a general equation valid for all stress-strain processes.
Both types of inner stress may depend differently on direction y and depend differently on stress-strain history.
The subscripts P and Q are not intended to match concepts of mean and deviatoric stresses, although a limited analogy is considered below. Both types of stress can be involved in volumetric behaviours in isotropic processes, and in volumetric and shear behaviours for more general stress-strain processes.

A generalised no-tension condition
A natural question to ask is, what might distinguish a P stress from a Q stress? The answer proposed herein relates to the no-tension condition for uncemented granular soils. Schofield's (2005) 'no-tension condition' implies that mean normal effective stress is never tensile, and a simple extrapolation would be that the inner stresses are never tensile. If the P-stress is restricted to being positive, then the following is sufficient to ensure * 0 This suggests that, like average and deviatoric effective stresses, the P-stress cannot be negative, whereas the Q-stress can be negative. This gives a limited analogy with average and deviatoric quantities.
This also suggests that the P-and Q-stresses play different roles in a particle aggregate and will need to play different roles in a constitutive model. A slightly more general limit will be assumed in the model developed later in this paper, as follows 34.
where m is a positive material constant. If * 0 P ψ σ ′ ≥ and m is between 0 and 1, then Schofield's condition is achieved.

Implications for energy and work
If there are two stresses, then it seems feasible that there may be two energies. The present paper assumes that the energy H ψ ′ is the sum of two energies, denoted as P ψ ′ and Q ψ ′

13
The existence of two types of stress and two types of energy suggests that two different mechanical processes are involved, and these are called the P-and Q-processes in this paper.
With two processes, it is now necessary to extend some of the previous proposals. One way to extend the calculation of incremental work in Equation 19 would be to define incremental works for each process as 36.
Positive work for the P-process occurs only for inner compression ( * 0 dV ψ < ) in direction ψ. Positive work for the Q-process can occur for inner compression if the Q-stress is positive, and for inner dilation ( * 0 dV ψ > ) if the Q-stress is negative.
Each process will need a balance equation consistent with Equation 26. The following proposal is explored herein 38.
where QP dW ψ ′ represents an increment of work done by the Q-process for direction ψ on the P-process for that direction. It is included for generality and will be used in the model developed later. Incremental receipts P dE ψ ′ and Q dE ψ ′ would sum to the net receipt dE ψ ′ for direction ψ, and their individual averages would be zero. Incremental dissipated works dP dW ψ ′ and dQ dW ψ ′ would each be non-negative and would sum to d dW ψ ′ . Figure 7 illustrates the proposed equations schematically as a flow diagram for a material direction y. Incoming incremental work, * dW ψ ′ , is split into components for the P-process (on the left) and the Q-process (on the right). Each process also receives energy from other material directions, and there is a transfer QP dW ψ ′ from the Q-process to the P-process (a negative value would represent a transfer in the opposite direction). Each process dissipates some of the net work and energy it receives, and stores the remainder as Helmholtz energy.
Negative changes of Helmholtz energy are possible, and energy can reduce in some processes through dissipated work. The schematic diagram confirms that no work or energy is unaccounted for, and that all opportunities for work or energy flows have been included.

Exploratory expressions for energies
The macroscopic data do not provide a particle mechanical explanation for the P-and Q-processes. Even so, it is useful to develop exploratory expressions for the energies, so that calculations can be done and predictions can be assessed. Table 1 gives some proposals. In group A, the expression for P ψ ′ is adapted from the expression for isotropic energy in the original Cam-Clay model (Schofield and Wroth, 1968). A different dimensionless constant m is used here. The expression for Q ψ ′ involves a quadratic because the Q-stress can be negative while the energy cannot be. The constant k here is also dimensionless and not the same as in Cam-Clay. The denominator ensures that the equation is dimensionally consistent.
The equations in group B are adapted from expressions that relate to the double-log form of isotropic compression curve discussed by Hashiguchi (1995) The familiar formulae for the asymptotic compression line for the two groups are given under 'Further inferences'. The present paper proposes that this line is evidence of a more fundamental curve that applies in terms of inner stress and hardening volumes.
The formula for this curve can be inferred directly, as listed. By using these equations to eliminate stresses from the equations for , asy H ψ ′ , an expression is obtained for asymptotic energy in terms of the hardening volumes. The symbol A ψ ′ is used to make it clear that the model in the present paper assumes that this energy is formally a function solely of hardening volume and relevant material constants.
By replacing the hardening volumes with specific volume, an isotropic energy R¢ is obtained. This will serve as a reference energy in the model developed later. The isotropic compression curves are obtained from the condition A R ψ ′ ′ = in circumstances when all stresses equal the mean normal effective stress and all hardening volumes equal the overall specific volume. More generally, it is possible to deduce the asymptotic value for the energy ratio / ψ ψ ¢ ¢ P A , as listed.
The significance of these inferences is that they provide guidance on the structure of a constitutive model. Data by Topolnicki et al. (1990) and others show that, in compressive proportional straining, the relation between effective stress and specific volume   The same data show that effective stress ratios tend to become constant in compressive proportional straining. This is well established for 1D compression, where the constant is represented by the coefficient of lateral earth pressure at rest. Different types of proportional straining give different asymptotic stress ratios, suggesting that different types produce different steady relationships between asymptotic energy and direction y. This idea will now be discussed from a different starting point.

Asymptotic energy surfaces
If energy is related to stress, then the familiar idea of yield surfaces might be usefully extended to surfaces in energy space or in terms of hardening volumes.
The present paper explores the concept of an 'asymptotic energy surface'. This is a relation among average asymptotic energy, reference energy and possibly other statistical features of the relation between asymptotic energy and direction y. If asymptotic energy varies solely with hardening volume, then the surface can also be interpreted as a relation between specific volume and statistical aspects of the hardening volumes. It would not necessarily be equivalent to a yield envelope, limiting envelope or other previous concept.
The present paper explores a surface described by the following relation among the deviator q A′ of the asymptotic energy, the average asymptotic energy and the reference energy R¢ 40.
where b is a positive material constant. By using the properties of a deviator (Equation 14), the equation can be rearranged as follows 41.
The surface is plotted in Figure 8 as an ellipse in a 2D space with axes of the average and the deviator. There is an obvious analogy with the yield envelope for Roscoe and Burland's (1968) modified Cam-Clay. However, as indicated above, the surface does not function as a yield envelope or a bounding surface in the model developed later.
In order to be able to relate the surface to directions of interparticle forces, it proves useful to define a dimensionless quantity ψ ξ′ as follows 43. It is proposed to call the average ψ ξ′ a 'surface label', since it identifies a surface in the asymptotic energy space. The label is invariant in the sense of being independent of direction, but variable in the sense that its value may change during a stress-strain process. ψ ξ′ is the contribution to the label provided by particle mechanical features associated with direction y.
where dT ψ ′ depends on the incremental contributions d ψ ξ′ and is given in Appendix D. The first term on the right-hand side is interpreted as volumetric hardening, because reference energy depends only on specific volume (Table 1). The remainder dT ψ ′ is referred to as 'transient incremental hardening' in the present paper. In general, it will have different algebraic forms for different equations for the asymptotic surface.
As noted earlier, data by Topolnicki et al. (1990) and others suggest that / / ψ ψ = ¢ ¢ ¢ ¢ dA A dR R along asymptotic compression curves for proportional straining. In terms of the proposed concepts, this suggests that asymptotic behaviour involves only volumetric hardening. Hence non-zero transient hardening would be interpreted as associated with non-proportional straining and with pre-asymptotic proportional straining.

Preamble
This part of the paper develops a simple mathematical model based on interpretations about material behaviour developed above. The aim is primarily to confirm that a model can indeed be developed that has some of the known features of the engineering constitutive behaviours of soils.
The behaviours targeted herein are the asymptotic ones. The surface label ψ ξ′ is assumed to move asymptotically to 1 during compressive proportional straining, and the increments d ψ ξ′ should move asymptotically to zero. The energy ratios / ψ ψ ¢ ¢ P A and / ψ ψ ¢ ¢ Q P are expected to tend to asymptotes P L′ and Q L′ , respectively, as given in Table 1.
The primary task is therefore to find incremental evolution laws that achieve these asymptotic behaviours. The following algebraic identities prove useful 45.
The identities show how the energy increments can be interpreted in terms of incremental energy ratios and volumetric and transient incremental hardening.
To achieve the expected behaviours, the present paper introduces a new incremental work quantity, a dW ψ ′ , which it is proposed to call 'active' work, hence the subscript 'a'. This has two key features. The first relates to Noll's (1958) principle of determinism and is based on the fact that the incremental active work for a (P or Q) process can be known in an incremental calculation without using the constitutive laws for the process, provided the strain increments are the givens. The second is that, in compressive proportional straining, incremental active work tends to become asymptotically equal to incremental dissipated work. Appendices E and F describe this further. The latter feature means that the difference between active and dissipated work can be used in a model to control behaviour, and this is illustrated below in three ways. ( )

Transient incremental hardening
where h is a proposed material constant, and Bb is an invariant derived in Appendix D.
The proposal satisfies the expectation, inferred earlier from data by Topolnicki et al. (1990)

Interprocess work
The interprocess work in Equations 38 and 39 is incremental work done by the Q-process on the P-process. For the purposes of a simple model, this was assumed to be related to the transient incremental hardening in the following way. First, hardening was split as follows 49.
dT ψ ′ is independent of the incremental active and dissipated works for the Q-process, and Q dT ψ ′ is independent of increments for the P-process. The following assumption was made for the interprocess work 50.
The assumption means that the balance equation (Equation 38) for the P-process simplifies, using Equations 44 and 45, as follows 51.
( ) The first equation has the form of an evolution law for the incremental energy ratio, and it does not contain any increments from the Q-process. Some assumptions that simplify the right-hand side are discussed later. The increment sP dW ψ ′ in the second equation is the 'source work' for the P-process. It represents the part of the external work increment that remains after the effect of volumetric hardening has been accounted for.

The Q-process
Using Equation 46 to simplify the work balance equation (Equation 39) for the Q-process, a result is obtained that can be further simplified by defining incremental source work sQ dW ψ ′ for the Q-process, as follows

53.
( ) This is the external work for the Q-process less a factor that depends on the P-process and the transient incremental hardening. Using this with Equations 39, 46 and 50 gives 54.
( ) This has the same form as Equation 51 for the P-process.
Since the aim is to find an evolution law for / ψ ψ ¢ ¢ Q P , it seems appropriate to find assumptions that simplify the right-hand side. Details are given in Appendix E. The result is 55. To achieve the expected asymptotic behaviour for the energy ratio, the present paper assumes that incremental dissipated work is related to the absolute value of incremental active work through a variable wyQ wyQ is expected to be an increasing function of / ψ ψ ¢ ¢ Q P and to be 1 when the energy ratio reaches its asymptote. The following is explored in the present paper 57.
where DQ is a proposed positive material constant.
If h < 1, and if the incremental source work is positive in direction y, these proposals make the increment of / ψ ψ ¢ ¢ Q P positive if / ψ ψ ¢ ¢ Q P is less than its asymptote, but zero if / ψ ψ ¢ ¢ Q P is at the asymptote. Hence the ratio will increase towards the asymptote in this case but will not go higher than the asymptote. If the incremental source work is negative, Equation 55 makes the incremental ratio negative. Hence, there is no possibility of / ψ ψ ¢ ¢ Q P moving from a value below its asymptote to a value above.

The P-process
The balance equation for the P-process was reduced previously to the evolution law in Equation 51. A formula is needed for the A is a mathematical operator whose form is determined by the equations for the asymptotic surface and incremental hardening. It operates on one function of direction y to produce another function of the direction. In this case, it operates on the difference between incremental active work aP dW ψ ′ for the P-process and incremental dissipated work dP dW ψ ′ for the P-process. Appendix F gives details.
The energy ratio is expected to tend to its asymptote in compressive proportional straining, so the left-hand side of the above equation will tend to zero in this kind of straining. To achieve this, the present paper assumes that, under all conditions, incremental dissipated work for the P-process is related to the absolute value of incremental active work as follows 59.
where wyP is a positive scalar multiplier. This is like Equation 56 for the Q-process, but the multiplier need not be the same. The present paper assumes that 60.
where DP is a proposed positive material constant.
The asymptotic behaviours implied by these proposals are complicated by the involvement of the operator A, but may be expected to be similar in general nature to the behaviour for the Q-process.

Preamble
The simple model was programmed in Microsoft Excel VBA. Appendix G outlines the calculation procedure. The aims were (a) to confirm that the proposals are feasible and internally self-consistent, (b) to show that behaviours relevant to engineering practice are achieved and (c) to compare results with a selection of published experimental data.  Figure 16) were selected by fitting the calculated responses to the limited data available, and so are not necessarily definitive. Particle mechanics approach to continuum constitutive modelling Dean the units for stress because it is the value of the specific volume on the asymptotic isotropic compression line (ICL) at a stress of 1 unit.
For the calculations below, the initial soil state in the calculations was assumed to be isotropic. The initial value of the surface label was taken as 1, and the initial values of the stress ratios * were all m. The initial soil state was on the asymptotic compression curve at a specific volume of 4 and a mean normal effective stress of 1·43 kPa. Figure 6 showed the results for calculations of isotropic loading and cycling for the model. Also shown are points representing the data for kaolin clay by Roscoe and Burland (1968) that was discussed earlier.

Isotropic loading and cycling
The results confirm that the model can reproduce the data well. Figure 9 shows calculated developments of work and energy plotted against the increment number. The points A, B, … correspond to the points in Figure 6. Those calculations started considerably before point A. The results confirm that cumulative dissipated work never reduced, as required in thermodynamics. Changes of energy for the Q-process were relatively small. The difference between the upper two curves represents change of energy for the P-process since the start of the calculations. These changes are approximately m/l of the overall cumulative work. This is about 0·4 in the calculations, representing a much higher fraction of recoverable work compared to an interpretation of the data in terms of Cam-Clay. Figure 10 shows some further details. Figure 10( against the increment number. During asymptotic compression, the active and dissipated works were about 85% of the work done on the material. During unloading, the work and the active work were both negative, so the active work rate remained positive while the dissipated work rate was negative. ′ ′ ′ for the Q-process. This is due to the small value of m 2 k/m for these calculations. Figure 10(c) confirms that the normalised inner stress ratio can pass through zero and change sign smoothly. have a noticeable effect, which will normally be larger for less stiff materials that deform more in shear.

One-dimensional loading and cycling
The calculations confirm that the model results have the same general characteristics as data for clay (e.g. Al-Tabbaa, 1984;Hight et al., 1985). The asymptotic curves for oedometer conditions in the upper diagram are of the same shape as the asymptotic curve for isotropic conditions, and would plot as a straight line if the horizontal axis was a log scale. There is a hysteresis loop on large unloading and reloading. There is no sharp yield point on reaching the previous maximum stress. This is consistent with data for many soils (Wesley, 2010).
The stress paths in the lower diagram are consistent in general shape with published data of many soils. The loading curves are consistent with an asymptotically constant value of the coefficient of lateral earth pressure. The unloading curves are broadly consistent with the familiar empirical relation between Ko and OCR by Mayne and Kulhawy (1982); see also Wroth and Houlsby (1985).
Points B, J and K have been marked. B is at the end of the first loading. J is when the deviator stress reaches zero during unloading. K is at the end of unloading. Figure 12 shows rose diagrams for these states (Meier et al., 2008). In each diagram, the magnitude of a quantity is plotted radially, in a direction corresponding to the direction y of interparticle forces in physical space. Vertical represents the axial direction of triaxial compression.
Figure 12(a) shows the variation of outer stress. The stress is a measure of the product of contact forces and numbers of contacts with forces in the relevant direction. The plot shows that this product is relatively larger for the axial direction of compression compared to the radial directions. This may be broadly consistent with expectations. Figure 12(b) shows the variation of hardening volumes V ψ ′. These are expected to be related to volumetric aspects of the fabric of particle packing. At the initial isotropic state, the diagram would have shown a circle of radius equal to the specific volume. The absence of circularity at the end of 1D loading shows that the particle structure  Particle mechanics approach to continuum constitutive modelling Dean 20 has become anisotropic. On unloading to state K, the shape of the rose curve does not change much, indicating that the anisotropy that was induced during loading persists during unloading.
Figure 12(c) shows the variation of contributions ψ ξ′ with direction.
The initial isotropic state with 1 ψ ξ′ = would be represented by a circle of radius 1. After loading, the contributions are highest in the direction of compression, indicating more hardening in that direction compared to the lateral direction of zero strain.

Undrained triaxial compression: isotropically prepared samples
Model calculations were carried out to represent samples prepared by isotropic compression to various specific volumes, followed by isotropic unloading to a common specific volume V = 2·42. These and all subsequent calculations include strain-related effects of changes of orientations of interparticle force. The stress-volume paths are shown in Figure 13. Figure 14 shows calculated responses for undrained triaxial compression applied to samples prepared in this way. The full curves represent calculation results. They are broadly consistent with expectations for a remoulded clay. The results at low OCR indicate contractive response (reduction of mean normal effective stress with shearing), while those at high OCR show dilative responses. The results also confirm the approximate existence of critical states in the model. Since critical states were not assumed in the model's equations, the results suggest that these states are emergent properties, not fundamental ones.
Also shown on the left are scaled data points for reconstituted London Clay from Bishop et al. (1965). The stresses have been obtained by dividing the stresses in the original data by 3·3. This was done so that the shape of the experimental stress path could be compared with the calculation results (the magnitudes of stresses predicted by the model can be readily adjusted by changing the model constant Giso). These data were also used by Schofield and Wroth (1968) to support their proposals for original Cam-Clay. While the comparison in Figure 14 is far from a full validation of the present model, it does indicate that the ideas underlying the model are on the right lines. Figure 15 shows rose plots for sample B of Figure 14. In each diagram, the four curves correspond to axial strains of 0%, 5%, 20% and 60%, respectively. In almost all diagrams, there is a clear difference between compressive responses for interparticle forces with orientations above 30° to the radial, and extensional behaviour for orientations below 30° to the radial. This is broadly consistent  Particle mechanics approach to continuum constitutive modelling Dean with the pattern of linear strains in this kind of stress-strain process, with extension in directions close to radial, and compression close to vertical.
Figure 15(a) shows that outer stresses reduce close to the radial direction but increase close to the vertical direction. This seems broadly consistent with expectations for interparticle forces. Normalised energy ratios for the P-process collapse below 30° (see Figure 15(b)), and inner stress ratios in these directions collapse to their minimum values (Figure 15(c)). Normalised energy ratios for the Q-process are close to 1 except near the changeover angle (see Figure 15(d)). A consistent response occurs for the contributions (see Figure 15(e)). Figure 16 shows the results for undrained triaxial shearing after a history of 1D compression followed by partial unloading.

Undrained triaxial responses after 1D preparation
The full curves are model calculations, while the dotted curves are the data for Lower Cromer Till from Hight et al. (1985). Eight samples were tested. All were initially one-dimensionally consolidated to point K. They were then variously one-dimensionally unloaded   Particle mechanics approach to continuum constitutive modelling Dean Start points for shearing phase Scaled data for London Clay (Bishop et al., 1965;also used by Schofield and Wroth, 1968)  There are some differences of detail between the measured paths and the calculated stress paths. However, the calculation results show the same general features as the data. This confirms that the ideas underlying the proposed model are on the right track. Figure 17 shows drained triaxial responses for the model. The plots in Figures 17(a) and 17(b) show the relation among the deviatoric stress ratio q/p ¢, axial strain and volumetric strain during shearing. The same initial states A, B, C and D were used as in Figure 14. The calculated stress paths are shown in Figure 17(c), and the stressvolume behaviour is shown in Figure 17(d).

Drained triaxial compression
The calculation results confirm that the model can reproduce the familiar drained behaviours of typical soils. The normally consolidated sample A and the lightly overconsolidated samples B and C experienced compressive volume strains, with stable plastic Particle mechanics approach to continuum constitutive modelling Dean yielding and the deviator stress increasing steadily to a critical state. The heavily overconsolidated sample D experienced dilative strain, with a peak deviator stress at relative small axial strain, followed by a calculated drop in deviator stress with further axial strain.
An estimate from the model's critical state line is also shown in Figures 17(c) and 17(d). The ratio of the pressure on the ICL to that on the critical state line at the same specific volume is about 2·3 for these calculations. This is broadly consistent with much data for clays; see for example data of London Clay by Parry (1958) discussed by Schofield and Wroth (1968).
Additional calculations were also carried out to determine the relation between peak deviator stress and mean normal stress for highly overconsolidated samples. These resulted in a curved 'failure' envelope shown dotted in Figure 17(c). The existence of a curved envelope is consistent with Schofield's (2005) 'soil map' and with general experience of soil behaviour (e.g. Lambe and Whitman, 1979). Figure 18 shows model calculation results for undrained simple shear applied to one-dimensionally prepared samples. All samples were initially compressed from a specific volume of 4, then (except for A) partially unloaded, so that all reached a specific volume of 2·42 at the start of shearing. Figure 18(a) shows the calculated stress paths in terms of the vertical effective stress and the shear stress on a horizontal plane (Figure 18(c)). The stress paths show the expected variations with OCR. The normally consolidated sample A shows a reduction of vertical effective stress during shearing, which is consistent with expected contractive behaviour. The highly overconsolidated sample D experienced an increase in vertical effective stress during shearing, consistent with dilative response.

Undrained simple shear: monotonic loading
For these calculations, it seems that an engineering shear strain of 30° was not sufficient to reach a critical state. Additional calculations confirmed that the contour representing a shear strain of 30° is curved in the model, as indicated in Figure 18(a).
Principal directions of incremental strain in this kind of shearing are at 45° to the vertical at all stages of shearing. Hence the principal directions of the in-plane effective stresses are expected to move towards a 45° orientation at large strain. However, this would be expected to be modified by anisotropy induced during the history of 1D compression. The calculation results for the major principal stress direction, shown in Figure 18 Cyclic simple shear data are important in the design of many offshore structures and for offshore and onshore design against seismic loading (Anderson, 1991;Ishihara, 1996;Luong, 1980;Muir Wood, 1991;Youd et al., 2001). Typically many cyclic tests may be performed and the results integrated into design in complex ways. A constitutive model capable of accurately predicting cyclic response can therefore have considerable utility. Figure 19 shows some calculation results for the simple model of the present paper. On the left is a record of the stress-strain response for undrained two-way symmetric cycles of simple shear applied after a history of 1D normal consolidation. The first cycle has the shape of a familiar hysteresis loop, of the type sometimes modelled using Masing's rule (Chiang, 1999;Kramer, 1996). As cycling progresses, larger cyclic strains develop due to a form of collapse of the stress-strain response in the middle of each half-cycle.
The diagrams on the right show developments with the cycle number. The vertical effective stress gradually reduces as cycling continues and will eventually reach near-zero. This is part of what causes the softening or 'degradation' in the load cycles. The magnitudes of the cyclic shear strains increase, as do the values of the peak angles of friction mobilised in each cycle on the horizontal plane.
The calculated responses are broadly consistent with much data of both sands and clays, as presented for example by Andersen (1991), Ishihara (1996), Jefferies and Been (2006), Seed and Lee (1966), Vaid and Chern (1985) and others.

Practical relevance of this research
As noted in the Introduction, the focus of much modern research in constitutive modelling is shifting from continuum mechanics to particle mechanics (Bolton, 2000;Cundall, 2001;Mesarovic et al., 2012;Seyedi Hosseininia, 2013). The research reported herein offers a mathematical basis that is consistent with this focus but which can be viewed as simpler.
The research started with a mathematical implementation of the familiar interpretation of Terzaghi's (1923Terzaghi's ( , 1936 Principle of Effective Stress. This resulted in Equation 7, in which ψ σ ′ is a continuous function of directions y of interparticle forces. The equation shows that the six independent components of the effective stress tensor consist of six differently weighted averages of the function. The domain of y is a 2D continuum that may be represented by the surface of a measurement sphere (Figure 3). The function therefore contains much more information than the tensor. For example, the tensor does not contain any information about the function's deviator (Equation 14).
If soil happened to be an elastic material, Kirchoff's (1850) uniqueness theorem might imply that the six averages would be adequate descriptions of stress for purposes of constitutive modelling. However, irreversible behaviour including frictional sliding at interparticle contacts, and particle crushing, means that Kirchoff's theorem does not apply to soils. Equation 7 therefore indicates that there is no reason to imagine that the effective stress tensor is an adequate description of stress for purposes of constitutive modelling. Particle mechanics approach to continuum constitutive modelling Dean Implications were developed for descriptions of work (Equation 10) and strain (Equation 12), and a preliminary approach was made in modelling the load-spreading and displacement-spreading roles of particles (Equations 16 and 22). Thermodynamic approaches were explored based on Collins and Houlsby's (1997) results for isothermal conditions. An assessment of data of kaolin clay published by Roscoe and Burland (1968) led to a concept of two particle mechanical processes, the P-and Q-processes.
New work balances were developed (Equations 38 and 39), which embody the first and second laws of thermodynamics. Data of proportional straining allowed behaviours under asymptotic conditions to be inferred, for the energy ratios / P A Constitutive equations were developed primarily to satisfy the requirement to achieve the expected asymptotic behaviours.
Validation calculations confirmed that the model performs well, including in comparisons with small datasets of kaolin clay ( Figure  6), London Clay ( Figure 14) and Lower Cromer Till ( Figure 16). The familiar concepts of a yield locus, plastic potential and critical states appeared in the results, even though they were not part of the model's assumptions. Just as in Calladine's (1971) model, these concepts are emergent properties of the model, not fundamental ones. There was no particular problem in modelling induced anisotropy, cyclic loading or effects of principal axis rotation, and the stress-strain responses in multiple cycling reproduce aspects of the known cyclic degradation of stiffness ( Figure 19).

Potential applications
An understanding of the constitutive behaviour of soils underlies all geotechnical design and forensic work, and shapes standards, codes of practice, methods of soil sampling and testing and data interpretation (Baligh et al., 1987;Lunne et al., 2006) and methods of teaching. The potential applications of this research are therefore wide-ranging, but there is also considerable work to be done to realise the opportunities that the work provides.
New particle mechanical explanations are needed for the new phenomena identified in this paper: the two energies, transient incremental hardening, receipts and active and source work. It is expected that approaches such as those pioneered by Rowe (1962Rowe ( , 1969) may be helpful, and much may be expected from the numerical modelling of granular assemblies, which also has a long history (Cundall et al., 1982;Tordesillas et al., 2009). It is expected that rose plots from the model, such as in Figures 12 This is just the familiar tensor transformation law for effective stress.

Appendix C. Effects of deformations on directions
The Cauchy stress tensor in Equation 7 can be interpreted as describing effective stress on the sides of a unit cube of soil. During an incremental deformation, a unit cube of material will change shape and orientation. Procedures for correcting the stresses as a result of this are well established and can be inferred by considering the Polar-Kirchoff stress tensors or by other means (Hashiguchi and Yamakawa, 2013;Spencer, 1980). However, the proposals in the present paper involve directions y. Material directions will change in any incremental rotation or shearing. Also, if shear strains occur, different material directions will change by different amounts. A need therefore arises for a new procedure to correct for these changes.
The procedure used in this paper may be described as follows.
One imagines placing a measurement sphere (Figure 3) around a spherical region of soil (Figure 4). For any patch P and direction y, one imagines marking the patch and the direction in the soil. As a result of an incremental deformation, the soil sphere will be an ellipsoid and may also have rotated. The marks in the soil will no longer line up with marks on the fixed measurement sphere. Now the outer stress at the start of the increment is a function of y and so can be plotted on a graph. Figure 20 shows a plot in which, for simplicity, y is represented on the 1D horizontal axis. Point A represents the outer stress at the start of the increment at a given y. Point B represents the stress at the end of the increment, in the direction as it is marked in the soil. The change of outer stress between A and B is due to constitutive behaviour associated with the material direction that was initially in direction y, but this material direction has now changed. For the purpose of the next incremental calculation, the value of stress is needed at point C in the original direction y.
Let dz be the difference between the directions as shown. Then for the simplified graph, the stress at C can be calculated from the stress at B and the gradient of stress with respect to direction 67.  ( ) The right-hand side is the scalar product of two vectors. The gradient has become the vector ψ σ Ñ ¢ , representing the conventional gradient of a function on the surface of a sphere. The change of direction has become the vector dFF -1 ny, which represents the incremental change of ny.
Similar calculations were used for the inner stresses and for the hardening volumes. Consequences for energies can then be calculated directly (Table 1).
The procedure is material-independent, since constitutive assumptions are not involved. dFF -1 is independent of reference configuration, so too is d G ψ σ ¢ . The increment dFF -1 ny already includes the familiar effects of strain, so there is no need to use the Piola-Kirchoff methods. The procedure does not require decomposition of F or dFF -1 , so questions do not arise about which decomposition is appropriate. Also, the calculation does not change the area under the graph of the quantity that is corrected. This is vital to preserve thermodynamic validity, because the Helmholtz energy in Collins and Houlsby's (1997) thermodynamics is represented by the average H ψ ′ in the present model.

Appendix D. Maths for the asymptotic surface
Three different possible surfaces were investigated as part of the present study. Table 3 lists the surfaces and the corresponding formulae for transient incremental hardening. Option C was used for the calculations presented in the main text, and its equations were derived as follows.   Table 3.
The second term in the brackets comes from Equation 73. The first term is then obtained by applying the condition that the average receipt is zero. A further simplification may then be achieved if incremental source work for the Q-process is assumed to be related to incremental active work as follows

76.
( ) where S ψ ′ and U ψ ′ are given in Table 3. Substituting these proposals into Equation 54 gives Equation 55. Detailed calculations confirm that this expression satisfies the two considerations mentioned at the start of this appendix.

Appendix F. Maths for the P-process
The evolution law for / P A where S ψ ′ and U ψ ′ are given in Table 3. A may be interpreted as a mathematical operator that, in the above equation, operates on one function of direction ( dP dW ψ ′ ) to produce another (defined by the right-hand side of the equation).
The operator can in principle be applied to different increments. The present paper assumes that incremental source work for the P-process is related to incremental active work as follows 81.
( ) The reason for making this assumption in the present paper is that it simplifies the maths. Using it to substitute for the source work in Equation 78 gives Equation 58.
where the quadratic nature of Q ψ ′ ensures that the first term on the right-hand side remains finite as * Q ψ σ ′ passes through zero, and Equations 53, 55 and 56 ensure that the last term on the right-hand side remains finite as * Q ψ σ ′ passes through zero.
6. Incremental changes of P-stress, Q-stress and hardening volumes can be determined using the balance equations and the energy equations (Table 1). These can then be re-referenced to the directions y using the procedure of Appendix C. 7. The new values of outer stresses at the end of the increment can then be determined (Equation 22). The new outer stresses are used to determine the new value of the effective stress tensor (Equation 7).