Introduction
Recently the authors have been involved with elastoplastic finite element computations of the stress and strain fields around crack tips, with a view to establishing a suitable fracture criterion for crack propagation under these conditions. A particular feature of this class of problem is the large number of elements that are necessary for modelling the singularity at the crack tip. Special crack tip elements, which have been reported for purely elastic solutions, avoid this problem, but they are not suitable for elastoplastic situations because the exact nature of the singularity has not been established. Hence efficient programs are essential for this type of work.
The early work' used an elastoplastic program developed by Nyak,* which was extended to three-dimensions by S a l ~n e n . ~ This program was based on the initial stress method4 and used a banded solution routine (then in common use at Swansea). During this period, Irons developed the 'front' solution ~y s t e r n , ~ which was far more efficient than the banded solution for solving large numbers of equations, both in terms of core requirements and time of execution.
The authors introduced the front solution into the initial stiffness version of Nyak's program.2 This version was preferred to the tangential stiffness method because of its relative efficiency in core utilization (typically 60 per cent of that required by tangential stiffness, as front area was used only once), and the faster execution times for small load increments. Experience with crack geometries6 had shown that small increments were essential in avoiding overestimation of the plastic zone. Hence iterations on the nodal residuals were also small (typically one or two), giving a nett saving in time over the tangential stiffness method.
In the three-dimensional program the stiffness formulation was replaced by an improved technique. Instead of using the full [B]'[D][B]d0 for the integration of the element stiffnesses, it is possible to obtain the correct Ki, terms for that element by rearrangement and multiplication of the [D] matrix. This method, first reported by Gupta, can reduce the computational time for the 20 node parabolic element, when using the full Gaussian quadrature procedure (i.e. 27 integration points), by a factor of 9.