The authors have very effectively presented upper and lower bound plasticity solutions for the classical problem of the bearing capacity of two-layered clays.
In referring to alternative methods of solving this problem however, the authors state the ‘experience has indicated that results from the displacement finite element method tend to overestimate the true limit load and, in some instances, fail to provide a clear indication of collapse altogether’.
This is not the experience of others working in the field, who know that the locking phenomenon to which the authors refer is easily avoided. The elasto-plastic displacement finite element method has been shown repeatedly (e.g. Zienkiewicz et al., 1975; Griffiths, 1982; Griffiths & Lane, 1999) to give reliable solutions to a range of collapse problems in geomechanics, especially plane strain analyses involving ϕu = 0 soil. In the 1982 paper mentioned above, the two-layered clay bearing problem was solved quite successfully using a relatively crude finite element mesh.
To avoid any further misunderstanding on this issue, a more thorough parametric study of the two-layer problem has been performed using a version of Program 6·0 (8-node quadrilaterals, reduced integration) from the published software of Smith & Griffiths (1988).
A typical mesh for this study is shown in Fig. 19. In all the analyses, a vertical displacement, δv, was applied incrementally to a smooth footing, with nodal reactions beneath the footing back-figured after each increment from the converged stress field. Convergence after each increment was defined as having occurred when the nodal displacements from one iteration to the next were changing by less than 0·1%.
Bearing capacity failure was deemed to have occurred when the nodal reactions reached a maximum and levelled out to within a tolerance of 0·1%. A typical plot of the development of bearing resistance with vertical displacement is shown in Fig. 20.
The shear strength of each layer of undrained clay was governed by Tresca's failure criterion, defined by the dimensionless function
Positive values of F generated within the mesh were considered ‘illegal’ and redistributed to neighbouring regions that still had reserves of strength.
Using the authors' definition of the bearing capacity factor
computed values of this quantity by elasto-plastic displacement finite elements are shown in Tables 4 and 5 and Fig. 21, for a range of cases, together with Fmax, the maximum value of F observed within the mesh at convergence. It should be noted that the limit solutions for H/B = 0·2 were obtained by linear interpolation of results from the authors' paper.
The upper and lower bounds bracket the displacement finite element results, as might be expected, the one exception being when H/B = 1 and cu1/cu2 = 5, where the upper bound solution appears to drift slightly below the displacement finite element result.
The reason for this discrepancy is unclear; however, the tables indicate consistently small values of Fmax, indicating a high level of convergence and accuracy in the displacement finite element results.
Figure 22 shows typical nodal displacement patterns from the displacement finite element analyses at failure for the cases of ‘weak on strong’ and ‘strong on weak’. The contrasting nature of the failure mechanisms in each case is clearly indicated.
In summary, elasto-plastic displacement finite element methods, in a single analysis, can be relied upon to give robust and accurate solutions to a wide range of geotechnical ‘failure’ problems.
Authors' reply
We should like to thank Goss and Griffiths for their interest in our paper. Their discussion does not focus on the results presented, but instead is chiefly concerned with the following statement: ‘In practice, great care must be exercised when finite element analysis is employed to predict limit loads. Even for quite simple problems, experience has indicated that results from displacement finite element method tend to overestimate the true collapse load and, in some instances, fail to provide a clear indication of collapse altogether’.
It is a well-known fact that displacement finite element analysis of undrained geotechnical problems can encounter severe numerical difficulties. In particular, the accuracy of the stresses computed from conventional finite elements is often reduced dramatically as the compressibility approaches zero. This phenomenon, which leads to an erroneous stiffening of the load deformation response, is widely known as ‘locking’, and has been reported in the literature by many researchers (e.g. Herrmann, 1965; Christian, 1968; Zienkiewicz et al., 1971; Naylor, 1974; Nagtegaal et al., 1974; Sloan, 1981; Sloan & Randolph, 1982; de Borst & Vermeer, 1984; Burd & Houlsby, 1990; Yu et al., 1993; Yu & Netherton, 2000).
In 1974, Nagtegaal et al. published a landmark paper on the difficulties associated with finite element calculations in the fully plastic range involving incompressible behaviour. By considering the limiting case of a very fine mesh, they proved that most displacement finite elements that employ low-order polynomials to model the displacement field are not suitable for incompressibility computations, particularly for axisymmetric problems. This is because the incremental incompressibility condition imposes a large number of constraints on the nodal velocities, which effectively reduces the available number of degrees of freedom. Since these constraints may multiply at a faster rate than the new degrees of freedom as the mesh is refined, it may not be possible to ensure that there are sufficient degrees of freedom available to accommodate the constant volume condition, regardless of how many elements are used in the grid.
One of the earlier approaches used to overcome this problem is the reduced integration rule suggested by Zienkiewicz et al. (1971). The element most commonly used in this method is the 8-noded rectangle with 4-point integration. As discussed by Naylor (1974) and Sloan & Randolph (1982), reduced integration has the beneficial effect of decreasing the total number of incompressibility constraints on the nodal degrees of freedom. This is clearly seen by noting that the maximum number of constraints per element must be less than, or equal to, the total number of integration points used in the calculation of the element stiffness matrices. A quasi-theoretical justification for using reduced integration in analysing incompressible materials has been given by Malkus & Hughes (1978). They proved that displacement formulations with reduced integration are, in certain cases, equivalent to mixed formulations where both stresses and displacements are treated as variables. Although conceptually appealing, this equivalence does not guarantee that spurious deformation modes will not occur.
Although it was once widely used in the finite element community, the reduced integration method can produce spurious stress and displacement oscillations. To illustrate the limitations of the reduced integration method, Sloan & Randolph (1983) presented examples of footings and vertical cuts in which the reduced integration approach leads to incorrect or unacceptable deformation predictions. More recently, Naylor (1994) demonstrated that even a high-order element (cubic triangles), when used with six integration points (reduced integration), produces a zero-energy mechanism. These shortcomings are well known in the area of computational geomechanics and a number of other important cases have been discussed by Sloan (1981) and de Borst & Vermeer (1984) among others. At a more fundamental level, the major limitation of using the reduced integration technique is that the incompressibility condition is satisfied only at a limited number of integration points, rather than everywhere within the element.
The authors did not mean to suggest that displacement finite elements should be avoided for collapse analysis, merely that they should be used with extreme care.
Finally, a major advantage of the methods used in the paper is that the upper and lower bounds automatically provide an error estimate for the limit loads. This feature is invaluable for collapse analysis where the exact solution is unknown.
