Scott W. Sloan

Upper bound limit analysis using discontinuous velocity fields

Abstract A new method for computing rigorous upper bounds under plane strain conditions is described. It is based on a linear three-noded triangular element, which has six unknown nodal velocities and a fixed number of unknown plastic multiplier rates, and uses the kinematic theorem to define a kinematically admissible velocity field as the solution of a linear programming problem. Unlike existing formulations, which permit only a limited number of velocity discontinuities whose directions of shearing must be specified a priori, the new formulation permits velocity discontinuities at all edges shared by adjacent triangles and the directions of shearing are found automatically. The variation of the velocity jump along each discontinuity is described by an additional set of four unknowns. All of the unknowns are subject to the constraints imposed by an associated flow rule and the velocity boundary conditions. The objective function corresponds to the dissipated power, or some related load parameter of interest, and is minimised to yield the desired upper bound. Since plastic deformation may occur not only in the discontinuities, but also throughout the triangular elements as well, the method is capable of modelling complex velocity fields accurately and typically produces tight upper bounds on the true limit load. The formulation is applicable to materials whose strength is cohesive-frictional, purely cohesive and uniform, or purely cohesive and linearly varying, and thus, quite general. The new procedure is very efficient and always requires fewer elements than existing methods to obtain useful upper bound solutions. Moreover, because of the extra degrees of freedom introduced by the discontinuities, the linear elements no longer need to be arranged in a special pattern to model incompressible behaviour accurately.

Refined explicit integration of elastoplastic models with automatic error control

Effective explicit algorithms for integrating complex elastoplastic constitutive models, such as those belonging to the Cam clay family, are described. These automatically divide the applied strain increment into subincrements using an estimate of the local error and attempt to control the global integration error in the stresses. For a given scheme, the number of substeps used is a function of the error tolerance specified, the magnitude of the imposed strain increment, and the non‐linearity of the constitutive relations. The algorithms build on the work of Sloan in 1987 but include a number of important enhancements. The steps required to implement the integration schemes are described in detail and results are presented for a rigid footing resting on a layer of Tresca, Mohr‐Coulomb, modified Cam clay and generalised Cam clay soil. Explicit methods with automatic substepping and error control are shown to be reliable and efficient for these models. Moreover, for a given load path, they are able to control the global integration error in the stresses to lie near a specified tolerance. The methods described can be used for exceedingly complex constitutive laws, including those with a non‐linear elastic response inside the yield surface. This is because most of the code required to program them is independent of the precise form of the stress‐strain relations. In contrast, most of the implicit methods, such as the backward Euler return scheme, are difficult to implement for all but the simplest soil models.

A fast algorithm for constructing Delaunay triangulations in the plane

Abstract This paper describes an algorithm for computing Delaunay triangulations of arbitrary collections of points in the plane. A FORTRAN 77 implementation of the scheme is given. For N points distributed randomly within a square domain, the expected run time for the algorithm is approximately 0(N 5 4 ) . Empirical tests, for N up to 10 000, indicate that the actual run time is substantially less than this prediction and is generally better than 0(N1.1). Excluding the memory required to store the co-ordinates, the algorithm requires slightly greater than 14 N words of integer memory to complete a typical triangulation. The efficiency of the proposed algorithm is verified by comparing its performance with other Delaunay triangulation procedures. Uses of the algorithm include the generation of finite element meshes and the construction of contour plots.

A smooth hyperbolic approximation to the Mohr-Coulomb yield criterion

Abstract The Mohr-Coulomb yield criterion is used widely in elastoplastic geotechnical analysis. There are computational difficulties with this model, however, due to the gradient discontinuities which occur at both the edges and the tip of the hexagonal yield surface pyramid. It is well known that these singularities often cause stress integration schemes to perform inefficiently or fail. This paper describes a simple hyperbolic yield surface that eliminates the singular tip from the Mohr-Coulomb surface. The hyperbolic surface can be generalized to a family of Mohr-Coulomb yield criteria which are also rounded in the octahedral plane, thus eliminating the singularities that occur at the edge intersections as well. This type of yield surface is both continuous and differentiable at all stress states, and can be made to approximate the Mohr-Coulomb yield function as closely as required by adjusting two parameters. The yield surface and its gradients are presented in a form which is suitable for finite element programming with either explicit or implicit stress integration schemes. Two efficient FORTRAN 77 subroutines are given to illustrate how the new yield surface can be implemented in practice.

Adaptive time stepping and error control in a mass conservative numerical solution of the mixed form of Richards equation

Abstract Adaptive time stepping with embedded error control is applied to the mixed form of Richards equation. It is the first mathematically based adaptive scheme applied to this form of Richards equation. The key to the method is the approximation of the local truncation error of the scheme in terms of the pressure head, although, to enforce mass conservation, the principal time approximation is based on the moisture content. The time stepping scheme is closely related to an implicit Thomas–Gladwell approximation and is unconditionally stable and second-order accurate. Numerical trials demonstrate that the new algorithm fully automates stepsize selection and robustly constrains temporal discretisation errors given a user tolerance. The adaptive mechanism is shown to improve the performance of the non-linear solver, providing accurate initial solution estimates for the iterative process. Furthermore, the stepsize variation patterns reflect the adequacy of the spatial discretisation, here accomplished by linear finite elements. When sufficiently dense spatial grids are used, the time step varies smoothly, while excessively coarse grids induce stepsize oscillations.

A fast algorithm for generating constrained delaunay triangulations

Abstract A fast algorithm for generating constrained two-dimensional Delaunay triangulations is described. The scheme permits certain edges to be specified in the final triangulation, such as those that correspond to domain boundaries or natural interfaces, and is suitable for mesh generation and contour plotting applications. Detailed timing statistics indicate that its CPU time requirement is roughly proportional to the number of points in the data set. Subject to the conditions imposed by the edge constraints, the Delaunay scheme automatically avoids the formation of long thin triangles and thus gives high quality grids. A major advantage of the method is that it does not require extra points to be added to the data set in order to ensure that the specified edges are present.

An implementation of Watson's algorithm for computing 2-dimensional Delaunay triangulations

Abstract A FORTRAN 77 implementation of Watsons algorithm for computing two-dimensional Delaunay triangulations is described. The algorithm is shown to have an asymptotic time complexity bound which is better than O(N 1.5 ) by applying it to collections of N points generated randomly within the unit square. The computer code obeys strict FORTRAN 77 syntax. Excluding the memory needed to store the co-ordinates of the points, it requires slightly greater than 9N integer words of memory to assemble and store the Delaunay triangulation.

Undrained stability of a square tunnel in a soil whose strength increases linearly with depth

Abstract This paper examines the undrained stability of a shallow square tunnel in a soil whose strength increases linearly with depth. Rigorous bounds on the loads needed to resist collapse are derived using two numerical techniques which are based on finite element formulations of the classical limit theorems. Both of the numerical procedures assume a linearized perfectly plastic soil model and require the solution of large sparse linear programming problems. For the range of tunnel geometries considered, the numerical results bracket the exact collapse loads closely. The solutions are presented in the form of stability charts which can be used by practising engineers for the purposes of design.

AN AUTOMATIC LOAD STEPPING ALGORITHM WITH ERROR CONTROL

SUMMARY This paper presents an algorithm for controlling the error in non-linear finite element analysis which is caused by the use of finite load steps. In contrast to most recent schemes, the proposed technique is non-iterative and treats the governing load-deflection relations as a system of ordinary differential equations. This permits the governing equations to be integrated adaptively where the step size is controlled by monitoring the local truncation error. The latter is measured by computing the difference between two estimates of the displacement increments for each load step, with the initial estimate being found from the first-order Euler scheme and the improved estimate being found from the second-order modified Euler scheme. If the local truncation error exceeds a specified tolerance, then the load step is abandoned and the integration is repeated with a smaller load step whose size is found by local extrapolation. Local extrapolation is also used to predict the size of the next load step following a successful update. In order to control not only the local load path error, but also the global load path error, the proposed scheme incorporates a correction for the unbalanced forces. Overall, the cost of the automatic error control is modest since it requires only one additional equation solution for each successful load step. Because the solution scheme is non-iterative and founded on successful techniques for integrating systems of ordinary differential equations, it is particularly robust. To illustrate the ability of the scheme to constrain the load path error to lie near a desired tolerance, detailed results are presented for a variety of elastoplastic boundary value problems.

Stability of an undrained plane strain heading revisited

The stability of an idealised heading in undrained soil conditions is investigated in this paper. The heading is rigidly supported along its length, while the face, which may be pressurised, is free to move. The problem approximates any flat wall in an underground excavation. Failure of the heading is initiated by a surface surcharge, acting with the self-weight of the soil. Finite element limit analysis methods, based on classical plasticity theory, are used to derive rigorous bounds on load parameters, for a wide range of heading configurations and ground conditions. Solutions for undrained soils with constant strength, and increasing strength with depth are presented. Recent improvements to finite element limit analysis methods, developed at the University of Newcastle, have allowed close bounds to be drawn in most cases. Previous research in this area has often been presented in terms of a stability ratio, N that combines load and self-weight into a single parameter. The use of a stability ratio for this problem is shown not to be rigorous, a finding that may be applicable to other stability problems in underground geomechanics.