## Abstract

This paper develops analytical solutions for shakedown limits of a cohesive-frictional half-space under a three-dimensional moving surface load. Melan's lower-bound shakedown theorem has been adopted as the theoretical basis for deriving shakedown limits. Rigorous lower-bound solutions are obtained for shakedown limits by establishing a self-equilibrated residual stress field that, together with the applied elastic stress fields, lies within the Mohr–Coulomb yield criterion throughout the half-space. By searching through the half-space, this study shows that the most critical location for satisfying the yield condition lies on the central plane. The analytical solutions derived in the paper can be used to benchmark numerical shakedown results, as well as to serve as a theoretical basis for the development of an analytical design method for pavements under moving traffic loads.

## 1. Introduction

Over the last two decades, there has been growing interest in the application of the fundamental shakedown theorems of Melan (1938) and Koiter (1960) to study the behaviour of elastic–plastic structures subjected to repeated or cyclic loading conditions (e.g. Johnson 1985, 1992; Ponter *et al*. 1985; Konig 1987). In the field of geotechnical engineering, the theory of shakedown is particularly useful for solving problems in the design of roads and pavements under moving traffic loads. This has been demonstrated by many researchers including Sharp (1983), Sharp & Booker (1984), Collins & Cliffe (1987), Raad *et al*. (1988), Collins *et al*. (1993), Yu & Hossain (1998), Shiau & Yu (2000*a*,*b*) and Collins & Boulbibane (2000).

When an elastic–plastic structure is subjected to a cyclic load, three distinctive situations may occur. First, if the load magnitude applied is so low (i.e. lower than the elastic limit of the structure) that nowhere in the structure is deforming plastically, then the behaviour will be entirely elastic. Second, if the load is larger than the first yield limit so that some part of the structure is deforming plastically but is less than a critical limit, then the plastic deformation will cease to occur after a number of load cycles. In other words, after a number of cycles with plastic deformation, the whole structure will respond purely elastically to the remaining load cycles. If this happens, then we consider the structure to have ‘shaken down’, and the critical limit (below which shakedown can occur) is termed a ‘shakedown (or elastic shakedown) limit’. Third, if the applied load is greater than the shakedown limit, then the structure will continue to exhibit plastic strains for however long the load cycles are applied. If this occurs, then the structure would eventually fail owing to excessive plastic deformation.

The essential calculation of applying the theory of shakedown is to determine theoretically the shakedown limit for a given structure with known material properties. Lower and upper bounds on the true shakedown limit can be estimated for elastic perfectly plastic materials using the bound theorems derived by Melan (1938) and Koiter (1960) assuming associated flow rules. Many numerical solutions have been obtained for the shakedown limit of pavements under moving surface loads. However, almost all existing lower-bound shakedown solutions are for pavements under two-dimensional plane-strain moving surface loads (e.g. Sharp & Booker 1984; Raad *et al*. 1988; Radovsky & Murashina 1996; Yu & Hossain 1998). As a key contribution in the field, Ponter *et al*. (1985) developed an upper-bound approach for shakedown analysis of three-dimensional rolling/sliding contact problems in cohesive materials. As a useful extension, Collins & Cliffe (1987) and Collins & Boulbibane (2000) later presented limited three-dimensional upper-bound solutions for shakedown of frictional pavements under moving surface loads. In the absence of good lower-bound solutions for three-dimensional loading, it is not possible to judge how effective the existing two-dimensional lower-bound and three-dimensional upper-bound solutions really are. The objective of this paper is to fill this gap by developing rigorous analytical lower-bound solutions for shakedown limits of cohesive-frictional soils under three-dimensional moving surface loads. Details of how shakedown limits can be used in pavement analysis and design will not be repeated here, but can be found in Sharp (1983), Collins *et al*. (1993) and Shiau & Yu (2000*a*,*b*).

## 2. Problem definition

We consider a cohesive-frictional half-space that is subjected to a surface contact loading limited to a circle of radius, *a* (i.e. *x*^{2}+*y*^{2}≤*a*^{2}), as shown in figure 1. If tensile stresses are treated as positive, then the boundary stresses on the surface are given as follows:(2.1)(2.2)(2.3)where *P* is the total normal load applied in the *z*-direction (i.e. the vertical direction) and *Q* is the total shear force applied in the *x*-direction (i.e. the moving load direction). All of the boundary stresses outside the circle of contact are assumed to be zero. It is also common to link the normal and shear loads by a frictional coefficient *μ* so that(2.4)

The stress distribution defined by equations (2.1)–(2.4), often referred to as the three-dimensional Hertz load distribution (see Johnson 1985; Collins & Cliffe 1987), is found to model the loaded region between a type and a pavement reasonably well (Sharp & Booker 1984). If the circular load moves along the *x*-direction, then we wish to determine the load condition (i.e. the value of *P* for a given *μ*) at which shakedown of cohesive-frictional materials can occur.

## 3. The lower-bound shakedown theorem

Melan's lower-bound shakedown theorem states that if any self-equilibrated residual stress field can be found, which, when combined with the elastic stress field produced by the applied loads, does not violate the yield condition anywhere, then shakedown will occur. If the applied load is denoted by *λp*_{0} (where *p*_{0} may be conveniently set as the unit pressure in the actual calculation), and *λ* is a dimensionless scale parameter, then all the induced elastic stress components are also proportional to *λ*. Melan's shakedown theorem hence demands(3.1)where is the elastic stress field owing to the applied pressure *p*_{0} and *f*(*σ*_{ij})=0 is the yield condition for the soil. The largest value of *λ* obtained by searching *all possible* self-equilibrated residual stress fields, , will give the actual shakedown limit .

## 4. The elastic stress field due to the three-dimensional Hertz stress distribution

The analytical solutions for elastic stress fields at any point (*x*, *y*, *z*) in the half-space, due to the three-dimensional Hertz stress distribution, defined in equations (2.1)–(2.4), were given by Hamilton (1983). The stress expressions that are relevant to this study are given below.

### (a) The stresses due to the normal load

The elastic stresses due to the normal load *P* are given as follows:(4.1)

(4.2)

### (b) The stresses due to the tangential (or shear) load

The elastic stresses due to the tangential shear force *Q* are given below:(4.3)(4.4)where the quantities *A*, *S*, *M*, *N*, *G*, *H*, *r* and *Φ* in equations (4.1)–(4.4) are defined below:

## 5. The residual stress field and shakedown condition

To determine a lower-bound shakedown limit, we need to consider a residual stress field. Any residual stress fields must satisfy the equations of equilibrium and stress boundary conditions upon removal of applied surface loads. A general material point may have all the six components of residual stresses, but symmetry and other considerations impose some constraints. This has been discussed in detail by Sharp (1983), Hills & Sackfield (1984), Johnson (1985) and Kapoor & Johnson (1992). As noted by Hills & Sackfield (1984), it is also useful to remember that, in some cases, one can conceive of alternative sets of residual stresses that result in the same shakedown limit, although some will be more physically reasonable than others. The closeness of a lower-bound solution to the true shakedown limit depends on our ability to propose a residual stress field as close as possible to the true one.

For the problem considered here, in which the material is assumed to be isotropic and homogenous, the resulting permanent deformation and therefore the residual stress field will be independent of the travel (*x*)-direction. It is reasonable to assume that under a moving three-dimensional Hertz pressure distribution, the most critical plane is one of the *xz* planes defined by *y*=const. On these planes, the only non-zero residual stress that may increase the shakedown limit would be . This is because for the case of normal loading, the shear stress is antisymmetric in *x*, as it has peaks that are equal in magnitude but opposite in sign on either side of *x*=0. It is therefore not possible to increase the shakedown limit by introducing a residual shear stress distribution , which must be independent of *x* (Johnson 1985; Kapoor & Johnson 1992; Collins & Boulbibane 2000). Of course, in the *y*-direction, a non-zero normal residual stress may well exist. The validity of using such a residual stress field containing both and in the half-space has recently been confirmed numerically by a three-dimensional numerical lower-bound shakedown study supervised by the author (Shiau 2001). In this numerical study, finite elements and mathematical programming were used to search for an optimum residual stress field that would lead to a best possible lower-bound solution for shakedown limits. The numerical results suggest that the optimum residual stress field contains non-zero and . Given that the residual stress field must be independent of the travel (*x*) direction, the equilibrium equations would require that be a function of *y* and *z* and be a function of *z*. Other than this requirement, there are no other constraints on these two residual stress components.

Therefore, the total stress field owing to a moving Hertz load distribution for an element on any *xz* plane (i.e. a *y*=const. plane) at any given moment can be defined as the sum of elastic stresses and residual stresses:(5.1)(5.2)(5.3)

According to Melan's lower-bound shakedown theorem, the material shakedowns if the total stress field does not violate the yield condition anywhere. If we assume that the yield condition for the half-space can be described by the Mohr–Coulomb criterion, and the normal stress in the *y*-direction, then , is the intermediate stress (which is possible as the residual stress is an arbitrary function of *z*), then the necessary condition for shakedown to occur is that the stress field defined by equations (5.1)–(5.3) satisfies the following inequality:(5.4)where *c* and *ϕ* are material cohesion and angle of friction. The yield condition (5.4) may be rewritten as follows:(5.5)withIt then follows that in order to satisfy the shakedown condition (5.5), the following must be met:(5.6)

We need to find the maximum value of by searching through every location in the half-space in order to determine the shakedown limit *λ*_{sd}. The shakedown condition (5.5) cannot be satisfied if the condition (5.6) is violated, but can simply be satisfied with *Z*=0 if the residual stress were chosen to make the first term of (5.5) vanish (i.e. *X*+*Y*=0; which should be possible as the residual stress is an arbitrary function of *y* and *z*). The limiting condition for shakedown occurs at the point in the half-space where the magnitude of is a maximum. By using the elastic stress solutions (4.1)–(4.4), it is a simple matter to search for the most critical plane (defined by the *y* coordinate) on which is a maximum at a critical point (defined by the *x* and *z* coordinates). Because the elastic stresses are symmetric along the central plane (*y*=0), it is sufficient only to search for either *y*≤0 or *y*≥0 for the most critical plane (i.e. the critical *y* value).

The above procedure satisfies all the requirements of Melan's shakedown theorem and therefore is a rigorous lower-bound analysis for cohesive-frictional materials. By making an *a priori* assumption that shakedown is controlled by the stress condition on the central plane (*y*=0), some shakedown solutions were obtained previously by Johnson & Jefferis (1963) and Johnson (1985) for three-dimensional rolling/sliding contact problems in a purely cohesive material. It is stressed that no such assumption is needed in the analysis presented in this paper, which is valid for a more general cohesive-frictional material.

Radovsky & Murashina (1996) applied a similar analysis to study the condition of shakedown for a cohesive-frictional half-space under a two-dimensional moving surface load. It is interesting to note that an alternative analysis (called the method of conics) used by Sharp & Booker (1984) gives similar shakedown solutions for pavements under two-dimensional moving loads. In an important contribution, Collins & Cliffe (1987) extended the upper-bound approach of Ponter *et al*. (1985) to frictional materials and obtained the same condition obtained from lower-bound shakedown theorem for two-dimensional surface loading. Therefore, it can be concluded that the shakedown solutions determined by the condition (5.6) for the case of a two-dimensional surface loading are true shakedown limit solutions. For the special case of purely cohesive materials (*ϕ*=0), this true shakedown limit under a rolling cylinder (i.e. two-dimensional plane-strain loading) was first established by Belyakov (see Gokhfeld & Cherniavski 1980) who obtained an upper-bound solution that was identical to the lower-bound solution derived by Johnson (1962) (see also Johnson 1985, 1992).

## 6. Numerical results and discussion

For a Hertz stress distribution, the maximum compressive pressure occurs at the centre of the loading area (*x*=*y*=*z*=0), and is obtained from equation (2.1) as . The shakedown limit determined by the condition (5.6) may be represented by a dimensionless shakedown limit parameter . Alternatively, the shakedown limits may also be defined in terms of the mean or average pressure acting on the circular area, . In the latter case, the dimensionless shakedown limit parameter is . These dimensionless shakedown limit parameters are dependent on material properties including the angle of internal soil friction *ϕ* and tangential frictional coefficient *μ*.

Table 1 and figure 2 present the numerical results demonstrating the dependence of the shakedown limit on both the angle of friction and tangential frictional coefficient. It is evident that the shakedown limit of a cohesive-frictional half-space increases with increasing friction angle. In addition, the presence of surface shear stress tends to reduce the vertical shakedown load significantly. This effect is important as a real traffic load generally has a non-zero tangential component.

A most interesting and useful conclusion obtained from this study is that the numerical search for the most critical point throughout the half-space that controls the yield condition suggests that it always lies on the central plane *y*=0. In fact, this condition was assumed previously by Johnson & Jefferis (1963) and Johnson (1985), without rigorous proof, in deriving lower-bound solutions for rolling/sliding contact problems in purely cohesive materials. Some justification was later provided by Ponter *et al*. (1985) who showed, using a mechanism-based upper-bound analysis, that the central plane is indeed the most critical plane for determining elastic shakedown limits in a cohesive material. However, the present paper provides a complete theoretical justification for assuming that the stress condition on the central plane controls shakedown analysis of rolling/sliding contact problems in a general cohesive-frictional material.

To illustrate the need to ensure that the yield condition is not violated anywhere in the half-space in shakedown analysis, table 2 presents shakedown limit results for two different cases:

The yield condition is satisfied everywhere in the half-space ; and

View this table:The yield condition is satisfied in a large part of the half-space defined by and ,

where *a* is the radius of the surface contact loading area. As expected, table 2 shows that the shakedown limits for case (1) is generally smaller than the results for case (2) in which the yield criterion is violated in the region . Obviously, the difference will increase if the region where the yield condition is violated becomes larger. For situations where the critical stress location lies on the loading surface *z*=0 (i.e. for cases when *μ* is close to 1.0), the difference becomes very small.

For the special case of zero friction angle (i.e. purely cohesive materials with the Tresca yield criterion), the solutions obtained in this paper are identical to the well-known solutions presented by Johnson & Jefferis (1963) and Johnson (1985). In particular, the shakedown limit parameter for the case of normal loading only is found to be *k*_{max}=4.68. This analytical shakedown limit is also very close to a numerical lower-bound solution of *k*_{max}=4.6 obtained by Shiau (2001) using finite elements and nonlinear programming techniques. The analytical lower-bound solution of 4.68 is almost identical to an upper-bound solution of *k*_{max}=4.7 for this special case derived by Ponter *et al*. (1985). This shakedown solution, which may be regarded as the exact solution, can also be compared with the upper-bound solution *k*_{max}=7.0 obtained by Collins & Cliffe (1987) using a simpler sliding mechanism with ‘V’ cross-section. When the tangential friction coefficient *μ* is greater than 0.3, the present lower-bound solutions are practically identical to the upper-bound solutions of Ponter *et al*. (1985). The difference between the present lower-bounds and Ponter *et al*.'s upper-bounds is largest (by some 20%) when *μ* is about 0.2, and then decreases quickly when *μ* is either reduced or increased.

For cohesive-frictional materials, the analytical lower-bound shakedown solutions obtained in this paper for a three-dimensional Hertz loading are close to the numerical lower-bounds obtained by Shiau (2001) using finite elements and nonlinear programming techniques. For *ϕ*=10°, the analytical lower-bound for *k*_{max} is 6.53 (from table 1) and the numerical lower-bound of Shiau is about 6.7. When the friction angle increases to 20°, the analytical lower-bound of *k*_{max}=9.25 from this paper is comparable to 9.7 from Shiau's numerical method. Of course, the numerical lower-bound solutions presented in Shiau (2001) are approximate only because the elastic stress field derived from the displacement finite-element method is not rigorous. Given that shakedown limits are sensitive to the accuracy of elastic stress fields, some errors would be introduced into the results of shakedown limits in any numerical lower-bound approach. Most recently, using a newly developed numerical linear matching technique presented in Boulbibane & Ponter (2003), limited three-dimensional upper-bound solutions have been obtained (A. R. S. Ponter 2004, personal communication) for a cohesive-frictional material, which agree well with the lower-bound solutions derived in this paper.

It is instructive to compare the solutions of shakedown limits for a cohesive-frictional half-space under a three-dimensional Hertz moving surface load with those for a two-dimensional plane-strain moving surface load (by a rolling cylinder), as shown in figure 3. The shakedown limit solutions for a cohesive-frictional half-space under a two-dimensional moving surface load can also be obtained using the condition (5.6). The elastic stress field owing to a two-dimensional plane-strain moving Hertz surface load can be found in Johnson (1985). In this case, the Hertz stress distribution on the surface with the normal load *P* per unit length in the *y*-direction and the tangential frictional coefficient *μ* is given by the following equations:(6.1)(6.2)The maximum vertical stress occurs at *x*=*z*=0, which is obtained from equation (6.1) as . On the other hand, the mean or average stress is . Therefore, the dimensional shakedown limit parameters in terms of both the maximum pressure and mean pressure are related by

The shakedown limits obtained in this way for a two-dimensional plane-strain moving load are presented in table 3 and figure 4. Some solutions for this two-dimensional shakedown problem were also given by Radovsky & Murashina (1996). A direct comparison between table 1 and three-dimensional suggests that the shakedown limits for a three-dimensional moving surface load is generally larger than those for a two-dimensional, plane-strain moving surface load. The difference increases with the angle of friction and decreases with the tangential frictional coefficient. In particular, when the tangential frictional coefficient becomes very large (with a large shear stress component), the difference between two- and three-dimensional solutions becomes very small as in these cases is a maximum on the surface. For the case of normal load only, the difference varies from 17 to 27% when the angle of friction increases from 0° to 45°. For a realistic case (e.g. *ϕ*=25°–45° and *μ*=0.2–0.3), the assumption of a two-dimensional traffic load distribution for the contact area between a tyre and a pavement would lead to a lower shakedown limit and therefore a conservative design.

A similar difference between two- and three-dimensional solutions has also been observed numerically by Shiau (2001) when finite elements and nonlinear programming techniques were used to determine shakedown limits. It is important to note that since the shakedown limit is sensitive to the accuracy of elastic stresses, a major limitation of using finite elements to determine shakedown limits is that a very large number of displacement finite elements are needed to ensure that the elastic stresses obtained are reasonably accurate. Therefore, this could make the resulting mathematical programming problem too large to solve.

## 7. Concluding remarks

In this paper, rigorous analytical solutions are developed for the necessary condition of shakedown of a cohesive-frictional half-space under a three-dimensional Hertz moving surface load. The solutions are based on Melan's lower-bound shakedown theorem and, strictly speaking, can only be regarded as lower bounds. However, the same lower-bound approach was shown to give exact shakedown limits for a two-dimensional loading. It is also shown that for the case of normal loading, the lower-bound solution obtained is almost identical to the upper-bound solution of three-dimensional rolling contact problems in a cohesive material. It is thus believed that the lower-bound solutions obtained in this paper are very good estimates of the exact solutions for three-dimensional loading, and this suggestion is supported by comparisons with a number of recent numerical lower- (Shiau 2001) and upper-bound solutions (Ponter *et al*. 1985; Boulbibane & Ponter 2003; A. R. S. Ponter 2004, personal communication) for three-dimensional loading problems.

Although the analytical solutions presented in this paper are for a single layered medium, the procedure developed is equally applicable to a multilayered system. For the case of a multilayered medium, however, analytical treatments are no longer possible, as the elastic stress fields caused by a surface loading will have to be determined numerically.

Finally, it should be pointed out that this study is concerned only with load limits where elastic shakedown will occur for which Melan's lower-bound and Koiter's upper-bound theorems are valid. Repeated applications of a load in excess of the elastic shakedown limit will cause repeated plastic deformation. As explained by Ponter *et al*. (1985) and Johnson (1992), it may take one of two forms: (i) a closed cycle of plastic strain (known as reverse or alternating plasticity), which may be expected to lead to failure by fatigue; or (ii) a steady accumulation of plastic strain (known as ratchetting or incremental collapse), which may be expected to lead to failure by excessive deformation. If the first case occurs, the structure is said to have undergone a plastic shakedown. As pointed out by Ponter *et al*. (1985), the load limit for plastic shakedown is not governed by Melan's lower-bound and Koiter upper-bound shakedown theorems. Using a heuristic argument presented in Ponter & Karadeniz (1985), Ponter *et al*. (1985) were able to derive some load limits for plastic shakedown (termed as ratchetting limits), which were found to be significantly higher than those for elastic shakedown with low surface friction coefficients. However, it must be stressed that the upper-bound solutions obtained by Ponter *et al*. (1985) for elastic shakedown agree well with the lower-bound elastic shakedown limits obtained in this paper. Furthermore, as rightly pointed out by Johnson (1992), the elastic shakedown limit provides a more rational design criterion for cyclically loaded structures, since, eventually, a structure that is subjected to a steady reverse plasticity (i.e. plastic shakedown) will probably fail by fatigue.

## Acknowledgments

The research described in this paper forms part of an EPSRC-funded research programme on stability and shakedown analysis in geotechnical engineering at the University of Nottingham. The author wishes to thank Professor Alan Ponter of the University of Leicester and Professor Ian Collins of Auckland University for useful discussions on shakedown theory. The author is very grateful to the referees for their constructive comments and helpful queries, which have considerably improved the rigour and clarity of this paper.

## Footnotes

- Received June 28, 2004.
- Accepted December 9, 2004.

- © 2005 The Royal Society