Second-order work in barodesy

Second-order work analyses, based on elasto-plastic models, have been frequently carried out leading to the result that failure may occur before the limit yield condition is encountered. In this article, second-order work investigations are carried out with barodesy regarding standard element tests and finite element applications. In barodesy, it was shown—like in hypoplasticity and elasto-plasticity—that second-order work may vanish at stress states inside the critical limit surface. For boundary value problems, an end-to-end shear band of vanishing second-order work marks situations, where failure is imminent.


Introduction
Investigations, based on elasto-plastic models, have been carried out by several authors [2-4, 10, 17, 25, 35] leading to the result that failure may occur before the limit yield condition is encountered. Second-order work investigations with hypoplasticity showed similar results [2,7,14,28], among other things that for loose soil second-order work vanishes at stress states inside the critical limit surface.
The literature on uniqueness, stability, bifurcation and failure is vast. Stability refers to systems, as characterized by their boundary conditions, and not only to materials. In particular, the nature of tractions on the boundary plays an important role, e.g. the question whether they are dead or follower loads or not. In this article, we consider the second-order work expressed as trð TDÞ. For symbols and notation, see Sect. 2. Many relevant citations can be found in Hill [9]. As they refer mainly to elastic solids, we cite some publications from soil mechanics.
Lade [15], based on experimental results, concludes that the violation of the stability criteria of Hill, trð TDÞ ! 0, and Drucker, trð TD pl Þ ! 0, does not necessarily evoke an observable collapse of the sample.
Nova [29] investigated the controllability of element tests of soils obeying the constitutive law D ¼ C T, with C being the fourth-order compliance matrix. As loading program, he denotes the prescription of n components of D (or linear combinations thereof) and the remaining 6 À n components of T (or linear combinations thereof), with 0 n 6. The response to this loading program consists of the corresponding 6 À n components of D and the remaining n components of T. He concludes that the constitutive relation determines a unique response, if the symmetric matrix C s is positive definite: or, equivalently, if the 'second-order work' W 2 is positive: The unique invertibility (in the sense of unique response to any loading program) was called by Nova 'controllability'. This term is, however, somehow misleading, as it can be easily conceived as the condition for obtaining a unique deformation of a sample by application of boundary tractions and displacements (in the sense that the deformation of a body is controlled by the boundary tractions and displacements). Such deformation is in many cases expected to be homogeneous, the corresponding tests are then called 'element tests'. The unique solution of a boundary value problem with prescribed boundary displacements does, however, not follow from W 2 [ 0 but from another condition, firstly derived by Hill [9], as shown here in a simplified form: we consider the velocity field v as the solution of a boundary value problem and investigate whether this solution is unique. Assume that there exists also another solution v 6 ¼ v. Denoting differences by a prime, e.g. v 0 ¼ v À v, we observe that v 0 vanishes at the boundary. The equilibrium equation reads r Á T ¼ 0, and continued equilibrium reads r Á _ T ¼ 0. The same equations hold also for the stress difference ÞdV and apply the theorem of Gauss. We thus obtain that this integral vanishes: because v 0 ¼ 0 on the surface S. Further, The second integral on the right-hand side vanishes due to continued equilibrium. Thus, for non-uniqueness must hold: implies uniqueness. For the special case v ¼ 0 we have: Thus we have and However, the statement that trð _ TDÞ [ 0 implies uniqueness cannot be generally inferred, e.g. uniqueness is implied if the material is incrementally linear.
Negative second-order work denotes softening, i.e. the tangential stiffness has at least one negative eigenvalue. For some static boundary conditions of dead loads, this implies increasing in kinetic energy and thus collapse. In fact, Nicot et al. [26] correlate vanishing second-order work with increase of kinetic energy and corresponding failure.

Symbols and notation
We use the symbolic notation for Cauchy effective stress T and stretching D. In the figures, the more familiar symbol r i instead of T i is used for the principal stresses. Normal stresses are defined negative for compression. Tensors are written in bold capital letters (e.g. X). jXj :¼ ffiffiffiffiffiffiffiffi ffi trX 2 p is the Euclidean norm of X, trX is the sum of the diagonal components of X. The superscript 0 marks a normalized tensor, i.e. X 0 ¼ X=jXj. 1 denotes the second-order unit tensor. Stresses are considered as effective ones, the normally used dash is omitted. The stretching tensor D is the symmetric part of the velocity gradient. Stretching D is only approximately equivalent to the strain rate _ e. For rectilinear extensions however, D equals _ e, considering the logarithmic strain e. p :¼ À 1 3 trT is the mean effective stress, e vol ¼ tre is the volumetric strain. For compressive strain, e i is defined negative.
The deviatoric stress is written as q ¼ Àðr 1 À r 3 Þ and the deviatoric strain reads e q ¼ À2=3 Á ðe 1 À e 3 Þ for axisymmetric conditions. Note that several definitions for dilatancy can be found. In barodesy is used d : A common definition of the angle of dilatancy w for axisymmetric triaxial compression, i.e. _ e 1 \0, is, see e.g. [23]: This ratio of volumetric to deviatoric strain rate can be expressed in dependence of d, (for d [ 0) as: From Eqs. 10 and 11 follows: 3 Barodesy Similar to hypoplasticity, barodesy 1 is a constitutive relation of the form T ¼ hðT; D; eÞ, which is based on the asymptotic behaviour of soil [12,13,20]. The general form of the barodetic constitutive relation is: R 0 and T 0 are the normalized tensors of proportional stress paths R and stress T, respectively. The tensor R depends on D, and the relation RðDÞ ¼ À expða D 0 Þ includes a stress-dilatancy relation [21], with a as a scalar quantity. The scalar quantities h, f and g depend on the actual stress level, the actual void ratio e and a stress-dependent critical void ratio e c and, thus describe barotropy and pyknotropy. 2 Appendix 1 summarizes all equations of barodesy.
Regarding second-order work, we have to check whether the expression can become negative or equal to zero. In barodesy, the second-order work is zero when Y (Eq. 15) is zero: The variables f, g and R 0 depend on D 0 . We therefore may rewrite equation 15 as: To determine the boundary of the region with W 2 [ 0 in stress space, we consider a stress state T and ask whether there is a D such that Eq. (16) be just fulfilled. 'Just fulfilled' means that there is only one solutionD. This implies Equation (16) introduces a relation between D 0 and T 0 , that is not necessarily unique. However, the additional requirement (17) imposes uniqueness and, thus, establishes a functionD As in Eq. 18 T 0 is a function ofD 0 , it can be written as only valid for vanishing second-order work. It therefore follows for W 2 ¼ 0, that (for isotropic materials) T 0 andD 0 are coaxial. As in Sect. 4 only rectilinear extensions are examined, the co-rotational, objective stress rate T coincides with _ T. This follows from T ¼ _ T þ TW À WT and W ¼ 0. In general cases (e.g. the finite element applications in Sect. 5), we do consider the Zaremba-Jaumann rate T.
Several rectilinear deformations represented by the principal stresses T 1 , T 2 , T 3 will be numerically examined as to whether there can be found D 0 -tensors such that trð _ TD 0 Þ ¼ 0. The boundary of the region in stress space with trð _ TD 0 Þ [ 0 is the surface of vanishing second-order work.

Element tests
We investigate the second-order work for specific loading paths in standard element tests (Sects. 4.1-4.2) and give a more general perspective in Sect. 4.3.

Undrained triaxial test
The undrained triaxial test is an illustrative example to explain vanishing second-order work inside the critical stress surface. As it is a rectilinear extension, we set D ¼ _ e. For an undrained test applies trD ¼ 0 From Eq. 20 follows trð _  Table 1 for parameters. The solid line marks the points of vanishing second-order work, the dotdashed line is the critical state line. The maximum of q ( _ q ¼ 0 and trð _ TDÞ ¼ 0) is clearly visible. Note that the mobilized friction angle at the maximum of q is lower than at critical state (marked with the crosses þ in Fig. 1).
In Appendix 2, we add the drained triaxial test as an illustrative example to investigate second-order work in barodesy.

Non-conventional drained triaxial tests
We consider drained triaxial tests with reduction of p at q ¼ const. For normally consolidated Weald clay (for parameters see Table 1) we set p ini ¼ 50 kPa, p ini ¼ 100 kPa and p ini ¼ 200 kPa. The tests start as conventional drained triaxial tests and at r 1 ¼ 1=K 0 Á r 2 the mean effective stress p is decreased by increasing the pore pressure, cf. similar experiments by Lade [15] and simulations by Wan and Pinheiro [34]. A reduction in the mean stress is obtained e.g. in the case of an excavation [8]. For the tests in Fig. 2a, the deviatoric stress remains constant (q ¼ 75 kPa for test A, q ¼ 150 kPa for test B, q ¼ 300 kPa for test C), hence _ q ¼ 0. The second-order work according to equation 37 simplifies thus to trð _ TDÞ ¼ _ p Á _ e vol . With decreasing p, i.e. _ p 6 ¼ 0, the secondorder work vanishes for this specific loading path at _ e vol ¼ 0. Simulations with barodesy show that in the nonconventional drained triaxial tests of Fig. 2, the secondorder work vanishes inside the critical limit surface. Figure 2a shows the stress paths of the non-conventional triaxial tests. The stress states with trð _ TDÞ ¼ 0 are marked with circles () and connected with a line. In Fig. 2b, c the solid line shows the volumetric strain and the dashed line the second-order work of test A. It is visible that secondorder work is zero at the local maximum of e vol , i.e. _ e vol ¼ 0, which has experimentally been confirmed: a sudden collapse is reported to occur [3,17] at the local maximum of volumetric strain. The mobilized friction angle at trð _ TDÞ ¼ 0 is smaller than the mobilized friction angle at critical state, cf. Fig. 2a.

Investigations in the deviatoric plane
The following analysis has been carried out numerically. We consider the deviatoric plane trT ¼ À500 kPa ¼ const. in the principal stress space spanned by T 1 ; T 2 ; T 3 and we search for the boundary of the region where trð _ TDÞ [ 0. We examine stress rays starting from the hydrostatic axis On each ray, we step forward with small increments of deviatoric stress. At each step, we check whether there are D tensors such that the condition trð _ TDÞ ¼ 0 is fulfilled. To this end, it is sufficient to check tensors D 0 , coaxial to T, distributed in all directions of the space D 1 ; D 2 ; D 3 , with the magnitude jDj ¼ 1, the polar angle 0\h\p, and azimuth angle 0\/\2p. The principal values of D are: We vary h and / independently on each stress ray for every step and search for minimum values of second-order work.
As soon as trð _ TDÞ ¼ 0 is encountered, the stress state T belongs to the searched boundary.  The Mohr-Coulomb yield locus in the deviatoric plane corresponds to a hexagon and the mobilized friction angle u m reads 3 : In barodesy, the cone of critical stress states practically coincides with the locus according to Matsuoka-Nakai [5,21] and is characterized through Eq. 25: In Fig. 3, vanishing second-order work is investigated with barodesy for Weald clay (u c ¼ 24 ). The following results are obtained: • For normally consolidated soil 4 (OCR ¼ 1 and e [ e c ) this cone lies inside the cone of critical stress states (Eq. 25). • For e ¼ e c (p e =p ¼ 2 respectively), the cone of critical states and the cone of trð _ TDÞ ¼ 0 differ only slightly: However, the cone of vanishing second-order work lies inside the cone of critical stress states.
• For highly overconsolidated soil e\e c (OCR = 6), the cone with trð _ TDÞ ¼ 0 lies outside the cone of critical stress states.
In Fig. 5a-c, stress locations in the deviatoric plane (trT ¼ const), characterized by the Lode angle a r (Eq. 26), are further investigated. The Lode angle a r is defined as follows: a r ¼ 30 holds for triaxial compression and a r ¼ À30 holds for triaxial extension. The angle b between the directions of stress T 0 and stretching D 0 reads, cf. [7]: For the undrained triaxial tests (with a prescribed D, i.e. Fig. 1b b is approximately 63:6 and the mobilized friction angle is u W 2 % 21:3 , when W 2 ¼ 0. A variation of D according to Eqs. 21-23 in order to find vanishing values of W 2 results in a slightly lower mobilized friction angle of u W 2 ¼ 21:1 . This variation was carried out at axisymmetric stress states. The critical friction angle of London clay is 22:6 . For the non-conventional triaxial tests in Fig. 2, b is approximately 86:4 and the mobilized friction angle u W 2 is approximately 22:340 , which is also lower than the critical friction angle of Weald clay with 24 . A variation of D according to Eqs. 21-23 in order to find vanishing values of W 2 results in a slightly lower mobilized friction angle of u W 2 ¼ 22:336 . In this article, the mobilized friction angle when second-order work vanishes is denoted as u W2 . 4 It is common to define the overconsolidation ratio (OCR = p e =p) by means of the so-called Hvorslev's equivalent consolidation pressure p e ¼ exp NÀlnð1þeÞ k Ã , divided by the actual mean stress p.
p e is the value of mean stress on the isotropic normal consolidation line which refers to the current specific volume ð1 þ eÞ. If p=p e ¼ const in a deviatoric plane (p ¼ const), then the void ratio e is constant. Figures 4 and 5 give a more general view of the vanishing second-order work locus: • In Fig. 4 a 3D representation of surfaces formed by trð _ TDÞ ¼ 0 for three different overconsolidation ratios (OCR ¼ 1, OCR ¼ 2 and OCR ¼ 6 according to Fig. 3) is shown. The cross section of the critical stress surface of barodesy (Eq. 25) with the deviatoric plane trT ¼ À500 kPa is added.
• Figure 5a shows the mobilized friction angles u W 2 (obtained with sin u m ¼ T min ÀT max T min þT max ) along the trð _ TDÞ ¼ 0 locus versus a r . For normally consolidated samples (OCR ¼ 1), the minimum mobilized friction angle is u W 2 % 18 , which is only 3 higher than the mobilized friction angle under oedometric conditions estimated with Jáky's relation. 5 Similar results have been obtained with hypoplasticity [7]. For the OCR ¼ 2, the cone of vanishing second-order work lies slightly inside the cone of critical states, cf. Fig. 5a. For highly overconsolidated soil, u W 2 is higher than u W 2 of the critical stress surface, cf. Fig. 5a. • Furthermore for OCR ¼ 1 the angle b between normalized stress T 0 and stretching D 0 according to Eq. 27 is 63 \b\69 , cf. Fig. 5b, the lower the void ratio (the higher the OCR), the higher the angle b. For highly overconsolidated soil, the angle b (77 \b\82 ) in Fig. 5b is higher than for slightly overconsolidated or normally consolidated soil. • Figure 5c shows the dilatancy d ¼ trD 0 in dependence of the Lode angle a r . For normally consolidated clay (OCR ¼ 1), the behaviour is slightly contractant (trD 0 % À0:2). Note that trD 0 ¼ 0 describes isochoric deformation and trD 0 ¼ À1 applies for oedometric compression. In addition, the angle of dilatancy w is also shown in Fig. 5c. 6 For an overconsolidation ratio of 2, clay is slightly dilatant (trD 0 % 0:1), for overconsolidated samples (OCR ¼ 6), trD 0 % 0:4, cf. Fig. 5c. Arthur et al. [1] (cited in [33]) report that the angles of dilatancy w in the shear plane in dense biaxial tests with sand were about 9 w 30 . Simulations of overconsolidated samples (OCR ¼ 2. . .6) with barodesy result in angles of dilatancy in the range of 3 \w\14 , see  Fig. 4 In barodesy, trð _ TDÞ ¼ 0 is described by a cone. A 3D representation of surfaces formed by trð _ TDÞ ¼ 0 for OCR ¼ 1, OCR ¼ 2 and OCR ¼ 6 according to Fig. 3 is shown. The critical stress surface of barodesy (Eq. 25) is added for trT ¼ À500 kPa. In this plot, Weald clay is simulated with barodesy [20] 5 From the earth pressure coefficient at rest K 0 ¼ 1 À sin u c ¼ 1 À sin u m 1 þ sin u m we obtain with u c ¼ 24 the mobilized friction angle under oedometric conditions u W2 ¼ arcsin sin u c 2 À sin u c % 15 . 6 Note that the angle of dilatancy w in Fig. 5c strictly applies for axisymmetric compression (i.e. a r ¼ 30 ) only. Experimental studies by Nakai [24] showed that 0\a r \15 for plane strain conditions. These findings correspond to results obtained with barodesy [21]. Attention should be paid to normally consolidated/slightly overconsolidated soils, where second-order work may vanish inside the critical stress surface. The results provide a basis for finite element applications, as shown below.

Finite element calculations
State of the art in geotechnical engineering are calculations of stress and strain fields with finite element approaches. Commercial finite element programs often allow assessment of stability by means of so-called strength reduction analyses (u-c reduction). A frequently used approach to assess stability is to reduce the shear parameters until loss of convergence in the numeric calculation. This approach is ambiguous, as convergence depends not only on the stability but also on numerical issues as, e.g. incrementation. Investigations on the occurrence of trð _ TDÞ 0 in finite element calculations could give a more clear identification of instability [17,22,25,27]. The here presented finite element calculations have been performed using Abaqus. For Abaqus, a user subroutine UMAT for the material model barodesy is available [31]. As barodesy is not formulated in the framework of elasto-plasticity, a special strength reduction approach has been developed [32]. Second-order work has been made available as additional output variable in the user material subroutine. Thus, second-order work can be visualized easily in the abaqus framework for barodesy.
Note that in the finite element applications addressed here, W 2 is evaluated with the actual _ T and D tensors. In other words, a search for the D-tensor that minimizes W 2 has not been carried out, as it would render the calculations extremely lengthy. Consequently, the condition W 2 ¼ 0 could be encountered even earlier. An analytical solution for the W 2 surfaces in stress space as developed by Niemunis [28] for hypoplasticity and applied by Meier et al. [22] is not yet archived for barodesy. However, the so obtained instability is still useful as indicator of failure.

Biaxial tests
The capability of modelling shear bands is an important property of material models. A first approach of visualizing shear bands in finite element calculations can be done on fine-meshed biaxial tests [14]. Finite element calculations of biaxial tests with barodesy have already been performed by Schneider-Muntau et al. [31], and the appearance of shear bands has been discussed. The same example is used in this article for shear band visualization with the secondorder work criterion. For a biaxial test with a homogeneous void ratio distribution over 200 elements (e ini ¼ 0:55), all elements have the same deformation, see Fig. 6 for the stress-strain relationship and Fig. 7 for the second-order work distribution. Second-order work vanishes for all elements at the same calculation step, which corresponds to the peak at an axial strain of e 1 ¼ 6:7%. Note that in laboratory tests inhomogeneous deformation exists from the very beginning of the test, be it small or pronounced.
As to be expected, biaxial tests with an imperfection (a single looser element with e ini ¼ 0:57) show an inhomogeneous deformation from the very beginning of the test. In Fig. 6b, the stress-strain relationship of every single element is displayed. A similar behaviour for all elements is obtained until the peak. With continued deformation, the stress-strain relationships of each element get chaotic. At the peak, second-order work becomes negative for some elements. As can be seen in Fig. 8a, those elements form a band of vanishing/negative second-order work. Note that dead load controlled deformation would not be possible at this stage any more. The first occurrence of a continuous shearband is at 6.7% axial strain, which corresponds to the peak of the homogeneous void ratio distribution. At this stage, the shear band is not visible in the void ratio distribution, Fig. 8b, and can only be guessed in the distribution of deviatoric strains, see Fig. 8c. Visualizing shear bands in terms of deviatoric strain or void ratio distribution only indicate strain accumulation and are dependent on the visualization scale.

Slope stability
Slope stability analyses in the framework of elasto-plasticity have been carried out, in order to identify unstable situations [4,16,27,30]. Meier et al. [22] determine unstable areas of shallow slopes on the basis of a hypoplastic second-order work stability criterion according to Niemunis [28]. Slope stability finite element calculations with a strength reduction method for barodesy have been presented by Muntau et al. [32]. The material parameters u c (critical friction angle) and N (ordinate intercept of the NCL) have been reduced gradually until the appearance of shear bands. The initial void ratio was set to 0.5 in all elements. Evaluating the shear band in terms of secondorder work 0 yields the hatched areas in Fig. 9a. The band of vanishing/negative second-order work from bottom to top of the slope is clearly visible, contrary to the shear band in terms of void ratio distribution in Fig. 9b.

Conclusions
Vanishing second-order work appears to be an suitable criterion for a situation where failure may occur. To evaluate trð _ TDÞ, all possible D-tensors and the pertinent _ Ttensors (resulting from a particular constitutive relation) should be investigated. This search is time-consuming and has been applied for the element tests in Figs. 3, 4 and 5: we varied the stretching tensor in the deviatoric plane in order to search for minimum values of second-order work. As soon as second-order work vanishes, the investigated stress state belongs to the searched boundary. For a constant overconsolidation ratio, vanishing second-order work is described by a cone. For normally consolidated to slightly overconsolidated soil, these cones lie inside the cone of critical stress states.
In barodesy, the second-order work approach showedlike in hypoplasticity and elasto-plasticity-that secondorder may vanish at stress states inside the critical limit surface. It is obtained that for rather loose soils, secondorder work vanishes inside the critical stress surface. Evaluating W 2 ¼ 0 in the non-conventional triaxial tests in Fig. 2 with the actual D-tensor and with a variation of D, almost led to the same results for mobilized friction angles.
For finite element calculations, the variation of D is not feasible, so second-order work has been evaluated with the D-tensor which is obtained from the equations of motion. For boundary value problems, an end-to-end band of vanishing second-order work marks the state where failure is imminent. In this article, second-order work has been investigated for the following boundary value problems: (1) Strain-controlled drained biaxial test with an initial imperfection (a slightly higher void ratio in one element). The stress-strain behaviour is similar for all elements until the peak, but then gets chaotic. At a certain strain-corresponding to the strain at the peak of a homogeneous sample-an area of vanishing or negative second-order work appears forming an end-to-end shear band. (2) Slope stability has been investigated by means of second-order work. The first appearance of an endto-end band of W 2 0 can define system failure.  Fig. 10 A drained triaxial test (e ini ¼ 0:55 and p ini ¼ 100 kPa) of Weald clay is simulated with barodesy. W 2 ¼ 0 when _ q ¼ _ p ¼ 0. In a mean stress and deviatoric stress are plotted versus deviatoric strain, in b second-order work is plotted versus the deviatoric strain e q