Evaluate the performance of an accelerated initial stiffness method in three dimensional finite element analysis

Newton’s method is a commonly used algorithm for elasto-plastic finite element analysis and has three common variations: the full Newton–Raphson method, the modified Newton–Raphson method and the initial stiffness method. The Newton–Raphson methods can converge to the solution in a small number of iterations when the system is stable; however, the methods can be quite computationally expensive in some types of problems, for example where the tangent stiffness matrix is unsymmetric or the plasticity is highly localized. The initial stiffness method is robust in those cases but requires a larger number of iterations. This prompted the formulation of many acceleration techniques in literature. In this paper, those techniques will be briefly discussed. This will be followed by the development of a modified acceleration technique for the initial stiffness method. The performance of the modified accelerated initial stiffness method will be examined in elasto-plastic analyses, using both direct and iterative matrix solvers. The results will be compared – in terms of the required number of iterations and the computation time – with an existing accelerated initial stiffness method, the non-accelerated initial stiffness method and the Newton–Raphson tangent stiffness method.     

Implicit numerical integration for a kinematic hardening soil plasticity model

Soil models based on kinematic hardening together with elements of bounding surface plasticity, provide a means of introducing some memory of recent history and stiffness variation in the predicted response of soils. Such models provide an improvement on simple elasto‐plastic models in describing soil behaviour under non‐monotonic loading. Routine use of such models requires robust numerical integration schemes. Explicit integration of highly non‐linear models requires extremely small steps in order to guarantee convergence. Here, a fully implicit scheme is presented for a simple kinematic hardening extension of the Cam clay soil model. The algorithm is based on the operator split methodology and the implicit Euler backward integration scheme is proposed to integrate the rate form of the constitutive relations. This algorithm maintains a quadratic rate of asymptotic convergence when used with a Newton–Raphson iterative procedure. Various strain‐driven axisymmetric triaxial paths are simulated in order to demonstrate the efficiency and good performance of the proposed algorithm. Copyright © 2001 John Wiley & Sons, Ltd.     

Finite difference techniques for modelling geothermal reservoirs

Two finite difference algorithms suitable for long-time simulation of the exploitation of a two-phase geothermal reservoir are presented. One is based on the hopscotch method proposed by Saul'yev¹ and analysed further by Gordon² and Gourlay.³ The other is based on the well-known ADI method. Both methods use a Newton–Raphson iterative technique in order to obtain accurate solutions of the non-linear difference equations involved. Rapid convergence of the iterative schemes occurs both for single-phase and two-phase reservoir problems. One- and two-dimensional model problems are presented.     

Numerical study on finite element implementation of hypoplastic models

A comprehensive numerical study on finite element implementation of hypoplastic models is presented. Two crucial aspects, local integration of the constitutive equations (the local problem) and forming tangent operators for Newton–Raphson iteration (the global problem), are investigated. For solving the local problem, different integration algorithms, including explicit and implicit methods, are examined using tri-axial compression tests and incremental stress response envelopes, as well as typical boundary value problems. For solving global problems, three different ways of generating the tangent operator are compared. The numerical evidences indicate that, in terms of accuracy, efficiency and robustness, explicit methods with substepping and error control are the best choices for constitutive integration of hypoplastic models while the so-called continuum tangent operators have certain advantages over two other types of numerically-generated consistent tangent operators.     

Isogeometric analysis for unsaturated flow problems

Unsaturated flow problems in porous media often described by Richards’ equation are of great importance in many engineering applications. In this contribution, we propose a new numerical flow approach based on isogeometric analysis (IGA) for modeling the unsaturated flow problems. The non-uniform rational B-spline (NURBS) basis is utilized for spatial discretization whereas the stable implicit backward Euler method for time discretization. The nonlinear Richards’ equation is iteratively solved with the aid of the Newton–Raphson scheme. Owing to some desirable features of an efficient numerical flow approach, major advantages of the present formulation involve: (a) numerical oscillation at the wetting front can be avoided or facilitated, simply by using either an h-refinement or a lumped mass matrix technique; (b) higher-order exactness can be obtained due to the nature of the IGA features; (c) the approach is straightforward to implement and it does not need any transformation, e.g., Kirchhoff transformation or filter algorithm; and (d) in contrast to the Picard iteration scheme, which forms linear convergences, the proposed approach can however yield quadratic convergences by using the Newton–Raphson method for solving resultant nonlinear equations. Numerical model validation is analyzed by solving a three-dimensional unsaturated flow problem in soil, and its derived results are verified against analytical solutions. Numerical applications are then studied by considering three extensive examples with simple and complex configurations to further show the accuracy and applicability of the present IGA.     

Implicit integration of a mixed isotropic–kinematic hardening plasticity model for structured clays

In recent years, a number of constitutive models have been proposed to describe mathematically the mechanical response of natural clays. Some of these models are characterized by complex formulations, often leading to non‐trivial problems in their numerical integration in finite elements codes. The paper describes a fully implicit stress‐point algorithm for the numerical integration of a single‐surface mixed isotropic–kinematic hardening plasticity model for structured clays. The formulation of the model stems from a compromise between its capability of reproducing the larger number of features characterizing the behaviour of structured clays and the possibility of developing a robust integration algorithm for its implementation in a finite elements code. The model is characterized by an ellipsoid‐shaped yield function, inside which a stress‐dependent reversible stiffness is accounted for by a non‐linear hyperelastic formulation. The isotropic part of the hardening law extends the standard Cam‐Clay one to include plastic strain‐driven softening due to bond degradation, while the kinematic hardening part controls the evolution of the position of the yield surface in the stress space. The proposed algorithm allows the consistent linearization of the constitutive equations guaranteeing the quadratic rate of asymptotic convergence in the global‐level Newton–Raphson iterative procedure. The accuracy and the convergence properties of the proposed algorithm are evaluated with reference to the numerical simulations of single element tests and the analysis of a typical geotechnical boundary value problem. Copyright © 2007 John Wiley & Sons, Ltd.     

An improved accelerated initial stress procedure for elasto-plastic finite element analysis

The most flexible and generally applicable methods for elasto-plastic analysis are those based on an incremental-iterative form of the initial stress approach, but such methods often exhibit slow convergence. The acceleration procedure known as the alpha-constant stiffness method is reconsidered and some modifications are proposed. The principal difference in the present approach lies in the use of a single acceleration parameter, rather than a diagonal matrix of acceleration coefficients. The new scheme shows a significant improvement in numerical stability and converges three times faster than the standard initial stress method. Some practical aspects associated with the method are discussed and a number of applications are presented.     

Explicit stress integration of complex soil models

In this paper, two complex critical‐state models are implemented in a displacement finite element code. The two models are used for structured clays and sands, and are characterized by multiple yield surfaces, plastic yielding within the yield surface, and complex kinematic and isotropic hardening laws. The consistent tangent operators—which lead to a quadratic convergence when used in a fully implicit algorithm—are difficult to derive or may even not exist. The stress integration scheme used in this paper is based on the explicit Euler method with automatic substepping and error control. This scheme employs the classical elastoplastic stiffness matrix and requires only the first derivatives of the yield function and plastic potential. This explicit scheme is used to integrate the two complex critical‐state models—the sub/super‐loading surfaces model (SSLSM) and the kinematic hardening structure model (KHSM). Various boundary‐value problems are then analysed. The results for the two models are compared with each other, as well with those from standard Cam‐clay models. Accuracy and efficiency of the scheme used for the complex models are also investigated. Copyright © 2005 John Wiley & Sons, Ltd.     

Improved numerical algorithms for frictional contact in pile penetration analysis

This paper presents a numerical formulation for frictional contact problems associated with pile penetration. The frictional contact at the soil–pile interface is formulated using the theory of hardening/softening plasticity, so that advanced models for the interface can be dealt with. A smooth discretisation of the pile surface is proposed using Bézier polynomials. An automatic load stepping scheme is proposed, which features an error control algorithm and automatic subincrementation of the load increments. The numerical algorithms are then used to analyse the installation process of pushed-in axial piles. It is shown that the smooth discretisation of the pile surface is effective in reducing the oscillation in the predicted pile resistances and the automatic load stepping scheme outperforms the classical Newton–Raphson scheme for this type of problem.     

Stress‐driven integration strategies and m‐AGC tangent operator for Perzyna viscoplasticity and viscoplastic relaxation: application to geomechanical interfaces

The paper proposes a stress‐driven integration strategy for Perzyna‐type viscoplastic constitutive models, which leads also to a convenient algorithm for viscoplastic relaxation schemes. A generalized trapezoidal rule for the strain increment, combined with a linearized form of the yield function and flow rules, leads to a plasticity‐like compliance operator that can be explicitly inverted to give an algorithmic tangent stiffness tensor also denoted as the m‐AGC tangent operator. This operator is combined with the stress‐prescribed integration scheme, to obtain a natural error indicator that can be used as a convergence criterion of the intra‐step iterations (in physical viscoplasticity), or to a variable time‐step size in viscoplastic relaxation schemes based on a single linear calculation per time step. The proposed schemes have been implemented for an existing zero‐thickness interface constitutive model. Some numerical application examples are presented to illustrate the advantages of the new schemes proposed. Copyright © 2016 John Wiley & Sons, Ltd.     

Analysis of single rock blocks for general failure modes under conservative and non‐conservative forces

After describing the kinematics of a generic rigid block subjected to large rotations and displacements, the Udwadia's General Principle of Mechanics is applied to the dynamics of a rigid block with frictional constraints to show that the reaction forces and moments are indeterminate. Thus, the paper presents an incremental‐iterative algorithm for analysing general failure modes of rock blocks subject to generic forces, including non‐conservative forces such as water forces. Consistent stiffness matrices have been developed that fully exploit the quadratic convergence of the adopted Newton–Raphson iterative scheme. The algorithm takes into account large block displacements and rotations, which together with non‐conservative forces make the stiffness matrix non‐symmetric. Also included in the algorithm are in situ stress and fracture dilatancy, which introduces non‐symmetric rank‐one modifications to the stiffness matrix. Progressive failure is captured by the algorithm, which has proven capable of detecting numerically challenging failure modes, such as rotations about only one point. Failure modes may originate from a limit point or from dynamic instability (divergence or flutter); equilibrium paths emanating from bifurcation points are followed by the algorithm. The algorithm identifies both static and dynamic failure modes. The calculation of the factor of safety comes with no overhead. Examples show the equilibrium path of a rock block that undergoes slumping failure must first pass through a bifurcation point, unless the block is laterally constrained. Rock blocks subjected to water forces (or other non‐conservative forces) may undergo flutter failure before reaching a limit point. Copyright © 2007 John Wiley & Sons, Ltd.     

A consistent formulation of a dilatant interface element

The paper considers a plane joint or interface element suitable for implementation into a standard non-linear finite element code. The element is intended to model discontinuities with rough contact surfaces, such as rock joints, where dilatant behaviour is present. Of particular concern is the formulation of a constitutive model which fully caters for all possible histories of opening, closing and sliding (accompained by dilation or contraction) in any direction. The non-linear incremental constitutive equations are formulated in a manner appropriate for a back-ward difference discretization in time along the path of loading. The advantage of such an approach is that no essential distinction need be drawn between opening, closing and sliding. Further, a convenient formulation of the constitutive equations is facilitated by representing the different contact conditions in relative displacement space. The state diagram in relative displacement space, however, changes from one time step to the next, and evolution equations for the updating must be formulated. These concepts are illustrated for two rock-joint models: a sawtooth asperity model and a limited dilation model. The models are based on a penalty formulation to enforce the contact constraints, and explicit equations for the tangent stiffness matrix and for the corrector step of the standard Newton–Raphson iterative algorithm are derived. These equations have been implemented as an user element into the finite element code ABAQUS⁷. Three examples are presented to illustrate the predictions of the formulation.     

A residual force finite element approach to soil–structure interaction analysis

An iterative process based upon a hybrid ‘residual force’ method is presented for solving elasto–plastic soil–structure interaction problems. In this approach the soil and the structure are treated as separate bodies and related only by compatibility of displacements and equilibrium of forces at the soil–structure interface. This scheme enables a significant improvement in numerical stability and rate of convergence over the conventional initial stress method. It is also shown that various interface conditions such as shear failure, slip and breakaway, and frictional and dilatant behaviour can be readily accounted for. Some practical aspects associated with the proposed scheme are emphasized for a number of numerical examples.     

Zero thickness interface elements—numerical stability and application

Many methods have been proposed to model joints in rocks or the interface between soil and a structure. Many analysts have reported numerical problems when using zero thickness interface elements while others have presented satisfactory results without comment of such difficulties. The numerical behaviour of zero thickness interface elements is further investigated in this paper. Some simple examples illustrate the application of interface elements to practical situations and highlight the numerical difficulties that may be encountered. Both ill-conditioning of the stiffness matrix and high stress gradients were found to cause numerical instability. Ill-conditioning can be reduced by careful selection of the size of the 2D elements adjacent to the interface. The problem of steep stress gradients is entirely one of inadequate mesh design. Contrary to other reports, this paper shows that the Newton–Cotes integration scheme has no benefit over Gaussian integration. Analyses of a retaining wall using interface elements confirm the analytical values of active and passive earth pressure coefficients which are commonly used in analysis and design of retaining walls.     

Consistent tangent matrices for density‐dependent finite plasticity models

The consistent tangent matrix for density‐dependent plastic models within the theory of isotropic multiplicative hyperelastoplasticity is presented here. Plastic equations expressed as general functions of the Kirchhoff stresses and density are considered. They include the Cauchy‐based plastic models as a particular case. The standard exponential return‐mapping algorithm is applied, with the density playing the role of a fixed parameter during the nonlinear plastic corrector problem. The consistent tangent matrix has the same structure as in the usual density‐independent plastic models. A simple additional term takes into account the influence of the density on the plastic corrector problem. Quadratic convergence results are shown for several representative examples involving geomaterial and powder constitutive models. Copyright © 2001 John Wiley & Sons, Ltd.     

A study on iterative methods for solving Richards’ equation

This work concerns linearization methods for efficiently solving the Richards equation, a degenerate elliptic-parabolic equation which models flow in saturated/unsaturated porous media. The discretization of Richards’ equation is based on backward Euler in time and Galerkin finite elements in space. The most valuable linearization schemes for Richards’ equation, i.e. the Newton method, the Picard method, the Picard/Newton method and the L-scheme are presented and their performance is comparatively studied. The convergence, the computational time and the condition numbers for the underlying linear systems are recorded. The convergence of the L-scheme is theoretically proved and the convergence of the other methods is discussed. A new scheme is proposed, the L-scheme/Newton method which is more robust and quadratically convergent. The linearization methods are tested on illustrative numerical examples.     

Mixed finite element discretization and Newton iteration for a reactive contaminant transport model with nonequilibrium sorption: convergence analysis and error estimates

We present a numerical scheme for reactive contaminant transport with nonequilibrium sorption in porous media. The mass conservative scheme is based on Euler implicit, mixed finite elements, and Newton method. We consider the case of a Freundlich-type sorption. In this case, the sorption isotherm is not Lipschitz but just Hölder continuous. To deal with this, we perform a regularization step. The convergence of the scheme is analyzed. An explicit order of convergence depending only on the regularization parameter, the time step, and the mesh size is derived. We give also a sufficient condition for the quadratic convergence of the Newton method. Finally, relevant numerical results are presented.     

Estimating Intrinsic Formation Constants of Mineral Surface Species Using a Genetic Algorithm

The application of a powerful evolutionary optimization technique for the estimation of intrinsic formation constants describing geologically relevant adsorption reactions at mineral surfaces is introduced. We illustrate the optimization power of a simple Genetic Algorithm (GA) for forward (aqueous chemical speciation calculations) and inverse (calibration of Surface Complexation Models, SCMs) modeling problems of varying degrees of complexity, including problems where conventional deterministic derivative-based root-finding techniques such as Newton–Raphson, implemented in popular programs such as FITEQL, fail to converge or yield poor data fits upon convergence. Subject to sound a priori physical–chemical constraints, adequate solution encoding schemes, and simple GA operators, the GA conducts an exhaustive probabilistic search in a broad solution space and finds a suitable solution regardless of the input values and without requiring sophisticated GA implementations (e.g., advanced GA operators, parallel genetic programming). The drawback of the GA approach is the large number of iterations that must be performed to obtain a satisfactory solution. Nevertheless, for computationally demanding problems, the efficiency of the optimization can be greatly improved by combining heuristic GA optimization with the Newton–Raphson approach to exploit the power of deterministic techniques after the evolutionary-driven set of potential solutions has reached a suitable level of numerical viability. Despite the computational requirements of the GA, its robustness, flexibility, and simplicity make it a very powerful, alternative tool for the calibration of SCMs, a critical step in the generation of a reliable thermodynamic database describing adsorption equilibria. The latter is fundamental to the forward modeling of the adsorption behavior of minerals and geologically based adsorbents in hydro-geological settings (e.g., aquifers, pore waters, water basins) and/or in engineered reactors (e.g., mining, hazardous waste disposal industries).     

Finite element formulation for poroelastic problem with zero effective stress boundary condition

Zero effective stress boundary condition along with constant fluid flux is commonly encountered in geotechnical applications such as uncased borehole stability, fluid injection and production at an uncased borehole, hydraulic fracturing and sand production. This complex boundary condition introduces high nonlinearity in the numerical simulation. Conventional iterative methods such as Newton–Raphson method are required to solve this nonlinear problem iteratively, which involve huge computing time and also pose numerical difficulties on the convergence. To overcome this numerical difficulty and hence reduce the computing time, a novel numerical technique is proposed in this paper. Its performance is evaluated using a numerical example simulating fluid injection around an uncased borehole. Copyright © 2009 John Wiley & Sons, Ltd.     

Development of the distinct lattice spring model for large deformation analyses

This study develops the distinct lattice spring model (DLSM) for geometrically nonlinear large deformation problems. The formulation of a spring bond deformation under a large deformation is derived under the Lagrange framework using polar decomposition. The results reveal that the DLSM's stiffness matrix under small deformations is the tangent stiffness matrix of the DLSM under large deformations. The formulation of the spring bond internal force under a given configuration is also presented and can be used to calculate the unbalanced force. Using these formulations, three nonlinear solving methods (the Euler method, modified Euler method, and Newton method) are developed for the DLSM with which to tackle large deformation problems. To investigate the performance of the developed model, three numerical examples involving large deformations are presented, the results of which are also in good agreement with the analytical and finite element method solutions. Copyright © 2013 John Wiley & Sons, Ltd.