An Efficient Model for Numerical Simulation of the Mechanical Behavior of Soils. Part 1: Theory and Numerical Algorithm

. A numerical model to simulate the mechanical behavior of soils is presented in this first part of the present work. This paper reports some advances obtained in the first stage of an ongoing research, which aims to develop a robust and efficient algorithm for consolidation analyses of saturated/unsaturated soils. Hence, the present work is mainly concerned with basic issues, such as efficiency of the finite element formulation and accuracy of the constitutive formulation. Since geotechnical materials exhibit elastoplastic characteristics, the theory of plasticity is applied here by means of the critical state concept using the Cam-Clay formulation. The constitutive equation is integrated using an explicit algorithm where the strain increment is divided into a number of sub-steps defined automatically by the numerical scheme. Eight-node hexahedral finite elements with one-point quadrature are employed in the spatial discretization of the geometrical domain. In order to avoid excitation of spurious modes, an efficient hourglass control is utilized in conjunction with a corotational formulation, which also contributes to the treatment of geometrical and physical nonlinearities.


Introduction
The use of numerical models to simulate geotechnical problems has strongly increased in the last decades. This growth is mainly due to expressive advances obtained in constitutive modeling of soils and other granular materials and the great reliability achieved by present numerical methods, such as the Finite Element Method (FEM). Since the pioneering works of Coulomb (1776) and Rankine (1857), the theory of plasticity has been applied to describe analytically the mechanical behavior of soils, which is entirely justified considering experimental observations that suggest irreversible behavior for strains, yield phenomena and shear-induced dilatancy. In this sense, computational efforts to analyze such problems are usually high, owing to nonlinear characteristics referring to the material and kinematical descriptions of the problem. Therefore, robust low-order finite element formulations using the one-point quadrature technique are welcome.
One of the most popular constitutive formulations based on critical state soil mechanics (CSSM) is the Cam-Clay model. The original Cam-Clay formulation was presented by Roscoe & Schofield (1963) and Schofield & Wroth (1968) using an elastic-plastic constitutive framework. Later, Roscoe & Burland (1968) proposed the Modified Cam-Clay model, where some drawbacks of the original formulation were eliminated. The original Cam-Clay yield surface presents a discontinuity for isotropic compression stress states, which leads to numerical difficulties to determine strain increments when an associated flow rule is employed.
The Modified Cam-Clay model has been widely used as a constitutive model to describe essential mechanical properties of clays. The main advantages of the model are referred to its relative simplicity and ability to represent strength and deformation characteristics of clays realistically with a limited number of parameters to be defined. On the other hand, despite its great popularity, Cam-Clay formulation seems to be inadequate to reproduce very complex soil behavior. Shortcomings associated to the model are: (a) soils are assumed to be isotropic. It is well known that natural soils may be anisotropic due to the action of deposition; (b) time effect on soil deformation is not taken into account (viscoelastic/plastic behavior is neglected); (c) failure stress states to the left of the critical state line are overestimated due to the form adopted for the yield surface; (d) shear strains are not well predicted within the yield surface because either the shear modulus or the Poisson's ratio is assumed to be constant; (e) the behavior of sands cannot be precisely predicted because sands do not follow exactly the principle of normality and experimental data show that the critical state point lies to the left of the peak of the yield locus; (f) soils under cyclic loading is another problematic subject because Cam-Clay models estimate large plastic strains for primary loading but for the remaining unloadreload cycles within the yield surface only purely elastic strains are obtained. Experimental evidences show that all unload-reload cycles result in hysteretic behavior.
Nevertheless, most of the basic features referring to mechanical behavior of soils can be covered using the Cam-Clay model, such as increasing stiffness while the material undergoes compression, hardening/softening and compaction/dilatancy behaviors, and a tendency to reach the critical state eventually. Moreover, the Basic Barcelona Model (BBM) by Alonso et al. (1990) utilizes the Cam-Clay formulation to develop a mathematical model to analyze unsaturated soil behavior, which is the final goal of the present research. Some other important critical state models may be found in Zienkiewicz & Naylor (1973), Carter et al. (1982), Naylor (1985), Gens & Potts (1988), Pastor & Zienkiewicz (1990), Crouch et al. (1994), Collins & Houlsby (1997), Yu (1998), Sheng et al. (2000) and Zhao et al. (2005).
Numerical procedures to simulate three-dimensional nonlinear problems lead usually to very time-consuming algorithms because an iterative process is required in order to satisfy the mechanical equilibrium. The stiffness matrix as well as the internal and external force vectors must be updated regularly and a numerical scheme is needed to integrate the incremental constitutive equation. As models based on the FEM utilize quadrature techniques to evaluate element matrices and force vectors numerically, the computational work for each iterative step is multiplied by the number of quadrature points to perform the numerical integration. Therefore, element formulations with full integration have been frequently replaced by one-point quadrature elements in order to obtain more efficient schemes for quadrature computations. In this case, hourglass control techniques are required to avoid excitation of spurious modes (hourglass modes). A hexahedral element formulation using one-point quadrature and a corotational coordinate system with hourglass control was presented by Hu & Nagy (1997). More recently, Duarte Filho & Awruch (2004) extended that formulation to geometrically nonlinear analysis using the model proposed by Liu et al. (1998) and uniform reduced integration. Numerical investigations based on the numerical model introduced by Duarte Filho & Awruch (2004) may be found in Andrade et al. (2007), Braun & Awruch (2008) and Braun & Awruch (2009). An application of three-dimensional low-order finite elements on soil mechanics was previously performed by de Borst & Groen (1999), but without dealing with reduced integration techniques.
In a finite element analysis where nonlinear materials are considered, stress updates are obtained from numerical integration of the constitutive equation, which is performed using implicit or explicit algorithms. The main difference between implicit and explicit integration schemes lies on how variables are evaluated: at known stress states for explicit integration methods or at unknown stress states for implicit ones. From previous investigations (see, for instance, Gens & Potts, 1988;Potts & Ganendra, 1994;Sloan et al., 2001;Zhao et al., 2005), one verifies that explicit al-gorithms are more indicated to reproduce complex stressstrain relations, such as those observed for critical state models.
In this first part of the present work, the numerical model presented by Duarte Filho & Awruch (2004) is extended to cover problems with elastic-plastic constitutive models, particularly the Modified Cam-Clay model. The critical state theory for soil mechanics is briefly explained using the elastoplastic framework to introduce the constitutive relations and a corotational reference system is presented in order to describe kinematically the motion of the continuum. A finite element formulation based on reduced integration for the eight-node hexahedral element is described, where stabilization is performed using hourglass control techniques to avoid numerical instabilities such as volumetric locking and shear locking. In addition, an explicit integration scheme, similar to that presented by Sloan et al. (2001), is adopted in order to update stress states at the center of the finite elements.

Analytical Model for Critical State Soil Mechanics
Soil materials present a very complex nature where a multiphase system is usually considered to describe the mechanical behavior of the soil mass (see Fredlund & Rahardjo, 1993). The phases observed within the soil volume are the soil skeleton, or the solid phase, and the pores or void spaces, which may be filled with gaseous (air and/or other chemical species) and/or liquid (water and/or other chemical species) matter. In order to describe the soil behavior using the theory of continuum mechanics, all phases are assumed to be continuous and linear such that all elements may be considered together by superposition. Momentum, mass and energy balance equations must be satisfied and a constitutive equation that relates energetically conjugated stress and strain measures must be also defined. In soil mechanics, the momentum balance for each infinitesimal element is described in terms of effective stresses and in the absence of temperature changes, the energy conserving equation may be disregarded. A measure for the stress tensor must be chosen such that the stress-strain relation maintains its objectivity in the nonlinear range. In this work, the Cauchy stress tensor and the small strain tensor are utilized, which will be defined in a corotational reference system.

Stress-strain relations
In this work, a nonlinear hypoelastic constitutive equation is considered to relate strain and stress measures in the elastic regime using a rate form as follows: with: where s ij and e ij are components of the Cauchy stress tensor s and small strain tensor e, D ijkl e represents components of the elastic constitutive matrix D e and K and G are the bulk and shear moduli, respectively. For critical state models, the tangential form of these moduli is usually assumed to be dependent on the effective mean normal stress p, which are usually expressed as: where e v e and e d e are the volumetric and deviatoric parts of the elastic strain tensor e e , q is the deviatoric stress, e 0 is the initial voids ratio, v s is the specific volume, n is the Poisson's ratio and k is the slope of the unloading-reloading lines (URL) on a ln p-v plane (see Schofield & Wroth, 1968 for additional information). It is worth to notice that the stress dependence on K and G leads to a nonlinear elastic constitutive matrix D e . The secant forms of Eqs. 3 and 4 are obtained by integrating Eq. 3, which yields: and, consequently: where p 0 is the effective normal stress at the beginning of the volumetric strain increment De v e .
In the elastoplastic range, a yield function f, a plastic potential g, a flow rule and a hardening law must be introduced in order to define the elastoplastic formulation (see Owen & Hinton, 1980). The yield function for the Modified Cam-Clay model may be written as: where M is the slope of the critical state line and p c is the preconsolidation pressure, which is related to the maximum effective pressure experienced by the soil mass. Hence, if a stress state leads to f < 0, it is assumed to be elastic and if a stress state leads to f ³ 0, it is assumed to be plastic (see also Fig. 1). Critical state models usually describe rate relations of volumetric plastic strains and the preconsolidation pressure according to: where l is the slope of the normal compression lines (NCL) on a ln p-v plane (see Schofield & Wroth, 1968 for additional information) and dl is the plastic multiplier of the flow rule, which is defined according to the principle of normality.
The incremental stress-strain equation in the elastoplatic range may be written as follows: and: where A is a hardening parameter and a f and a g are the flow vectors based on the yield function f and plastic potential g, respectively. The flow vectors are evaluated using the generalized formulation presented by Nayak & Zienkiewicz (1972). By assuming an associative flow rule, where f º g, a g and a f lead to the same result.

Motion description in the corotational reference system
For elastoplastic formulations, strain increments are imposed taking into account the last incremental displacement field obtained from the solution of the equilibrium equation, considering a generic iterative step. At each iterative step, the respective stress updates are obtained according to the material behavior, using either Eq. 1 in the elastic range or Eq. 9 in the plastic range. Since any motion of a continuous medium may be decomposed into rigid body and deformation motions, assuming that the finite element discretization is fine enough, the pure deformation portion of the motion is a small quantity compared with the element dimensions if the motion decomposition is performed in a corotational coordinate system. Consequently, the small strain hypothesis can be considered appropriately.
In order to update strains in the corotational system, the following expression is adopted: where n + 1 and n indicate initial and final positions of the time interval [t n , t n+1 ] and $ e is the strain tensor calculated in the corotational system. The strain increment D$ e is obtained using the mid-point integration of the strain rate tensor as follows: where $ x n +1 /2 is the geometric configuration of the finite element in the corotational system at the mid-point of the time interval [t n , t n+1 ] and D $ u def is referred to the deformation part of the total displacement increment D $ u in the corotational system. The total displacement increment D $ u is decomposed into a part owing to pure deformation D $ u def and a part owing to pure rotation D $ u rot , such that The increment of deformation displacements in the corotational system is obtained using: where $ x n and $ x n 1 + are geometric configurations of the finite element in the corotational system at t n and t n+1 , respectively (see Fig. 2). The geometric configurations $ x n , $ x n 1/2 + and $ x n 1 + are obtained from the following transformations: where R n , R n+1/2 and R n+1 are orthogonal transformation matrices, performing rotations from the global coordinate system to the corotational coordinate system, and x n , x n+1/2 and x n+1 are geometric configurations defined in the global coordinate system. The subscripts n, n + 1/2 and n + 1 denote positions in the time interval [t n , t n+1 ]. The components of the transformation matrix R are given by: ; where x, h and x j are vectors containing local nodal coordinates and global nodal coordinates associated to the eight-node hexahedral element (see Fig. 3), respectively. Since the local coordinate system rotate attached to each element in a corotational coordinate system, stress measures are not affected by rigid motions (objectivity). Therefore, the Cauchy stress tensor is employed in this work to calculate stress values in the corotational system. The Cauchy stress tensor in the global system is obtained from an objective tensor transformation as follows: where s and $ s are the Cauchy stress tensor evaluated in the global and corotational coordinate systems, respectively.
On the other hand, stress rate measures are performed in this work using the Truesdell rate tensor in order to maintain objectivity of the stress updates in the corotational system: where I is the unit tensor and where & $ e and & $ w are the strain rate tensor and the spin tensor, respectively, which are evaluated in the corotational system.
where t is the traction vector, W E is the volume referred to element E and G E represents its boundary surfaces, both considered in the FEM context. Spatial coordinates and displacements are approximated at element level using the eight-node hexahedral finite element formulation, which may be expressed as: where x, du and u are the coordinate, virtual displacement and displacement vectors evaluated within the element domain and x E , du E and u E are their respective nodal values. The column matrix N contains the shape functions of the eight-node hexahedral element, which is usually presented as (see also Fig. 3): where x, h and z represent the axis directions of the local coordinate system (corotational system) and j denotes the local node numbers of the hexahedral element.
The equilibrium equation (Eq. 21) must be iteratively satisfied using the incremental approach, where the stiff-ness matrix and the internal force vector are considered as functions of the current element configuration as well as the current stress state. The Newton-Raphson method is employed in order to linearize all nonlinearities of the model. Consequently, the internal force vector is submitted to a Taylor series expansion within the time interval [t n , t n+1 ], which gives: where i and i-1 indicate current and previous iterative steps within the time interval [t n , t n+1 ] and K n +1, i -1 tan represents the tangent stiffness matrix.
The tangent stiffness matrix and the internal force vector are evaluated in the corotational coordinate system at time instant t and iteration i as follows: where the gradient matrix $ B and $ W E are referred to the current element configuration in the corotational system, D and T are fourth order tensors related to the elastic constitutive equation and Truesdell rate terms, repectively, and s i is the corotational Cauchy stress tensor evaluated at the iterative step i (see Braun & Awruch, 2008;Duarte Filho & Awruch, 2004 for further details). It is important to notice that T vanishes for geometrically linear problems as well as D becomes D e for elastic materials and D ep for elastoplastic materials.
In order to solve the equilibrium equation, the tangent stiffness matrix and the internal force vector are brought back to the global coordinate system using the following objective transformations: where R is the transformation matrix defined in Eqs. 16 and 17. The final matrix format of Eq. 21 in the global system, considering the incremental approach, is given as follows: where subscripts n+1 denote current position in the time marching with i and i-1 indicating current and previous iterative steps in the Newton-Raphson method. The vector of external forces f ext represents the right-hand side terms of Eq. 21.

Reduced integration and element stabilization
Element formulations based on reduced integration must be stabilized using hourglass control techniques in order to avoid numerical instabilities. Volumetric locking, for instance, is remedied using reduced selective integration, where the gradient matrix B is decomposed as follows: where B( ) 0 corresponds to volumetric terms of the strain tensor, which is evaluated at the center of the element (x = 0, h = 0, z = 0), and~( , , ) B x h z is referred to deviatoric terms of the strain tensor.
In addition,~( , , ) B x h z must be expanded using Taylor series at the center of the element up to bilinear terms. Consequently, Eq. 28 can be re-written as: where B(0) is a sum of the volumetric and deviatoric parts of the gradient matrix, both obtained from one-point quadrature.
The vector representation of the stress tensor s is also submitted to a Taylor series expansion over the center of the element, which leads to: where all the terms above have similar interpretations to those described in Eq. 29. By substituting Eqs. 29 and 30 into the left-hand side term of Eq. 21, considering the relation between strain and displacement components and the decomposition performed in Eq. 28, the internal force vector must be re-written as follows: where f int,0 and f int,hg are the one-point quadrature and hourglass control parts of the internal force vector, respectively, which may be expressed as: where K stab is the element stabilization stiffness matrix and W E is the volume of the element. Derivatives of the stress vector are obtained employing the stabilization matrix E proposed by Hu & Nagy (1997), which leads to a constitutive-like relation between stress and strain derivatives (see Duarte Filho & Awruch, 2004;Braun & Awruch, 2008). For elastoplastic behavior, Reese (2005) proposed a modified stabilization matrix constituted by an optimized stabilization parameter m. In this work, the following procedure is adopted: where A 0 is the hardening parameter given by Eq. 11 and evaluated only at the first yielding of the respective element. The stabilization matrix E is obtained using the stabilization parameter m as follows: Shear locking is removed describing the shear components of the strain tensor in an orthogonal corotational coordinate system. In addition, all shear components are linearly interpolated in a single direction of the reference system as follows:  The internal force vector is not well evaluated for distorted elements if one-point quadrature is used. In order to correct this deficiency the gradient matrix obtained with reduced integration B(0) must be replaced by uniform gradi-ent submatrices ¢ B a ( ) 0 defined by Flanagan & Belytschko (1981): where the subscript a corresponds to the local node number of the element. It is worth to notice that all equations presented in this section are implicitly referred to the corotational coordinate system, although the symbol^adopted previously has been omitted in the respective expressions.

Integration of the stress-strain relation
In elastoplastic analysis, strain increments are imposed at element level based on the incremental solution obtained from the equilibrium equation. On the other hand, stress increments are dependent on the material behavior. A trial elastic stress increment is usually considered as an initial estimative for the new stress state. If this new stress state does not lead to plastic yielding, the corresponding stress increment is taken as true. However, if plastic yielding is observed, the following system of ordinary differential equations must be solved: where T is the artificial time, which is defined in the range 0 £ T £ 1, t 0 is the time at the starting point of the load increment (Dt), t is the time in a position within the load increment (Dt) and B is a hardening parameter.
In the remainder of the present section, numerical procedures to solve the elastoplastic problem are briefly described. A step-by-step algorithm of the present numerical model may be found in Appendix A.
where D e is the secant constitutive matrix, which is defined using K and G (see Eqs. 5 and 6) owing to the nonlinear elastic behavior observed in critical state models. It is important to notice that K and G are evaluated using the stress state at the starting point of the current strain increment (s 0 ) and the total volumetric strain increment De v . For conventional constitutive models, such as the Mohr-Coulomb yield criterion, these parameters are linear elastic constants evaluated using classical elastic relations, which leads to the tangent form of the constitutive matrix D e . In this case, D e must be replaced by D e in Eq. 41 and hereafter.
(b) Determination of the elastic fraction of the stress increment The elastic trial stress increment Ds e is now utilized to check the yield function (Eq. 7) with the corresponding trial stress state ( ) s s 0 + D e . If no plastic yielding occurs ( ( , ) ) f s s , the stress state is updated according to the trial stress state. On the other hand, if the trial stress state causes plastic yielding, the following alternatives must be considered (see also Fig. 4): : the stress state has changed from elastic to plastic. The elastic fraction of the total strain increment De must be determined in order to advance the initial stress state s 0 to the stress state on the yield surface s yld ; 2) if f ( , ) : the stress state is initially lying on the yield surface and the trial stress state causes plastic yielding. In this case, elastoplastic unloading may occur, which is observed if the angle q between the flow vector a 0 and the tangential elastic increment Ds e (K,G) is larger than 90º. The angle q may be calculated using: where a 0 is evaluated at the initial stress state s 0 and Ds e is evaluated using the incremental form of Eq. 1. If no elastoplastic unloading occurs, the strain increment De is assumed to be totally plastic. Nevertheless, if the elastoplastic unloading conditions are satisfied, the elastic fraction of the total strain increment De must be determined.
The stress state on the yield surface s yld may be determined considering the following nonlinear equation: where a is a scalar to be determined in order to satisfy Eq. 43. The secant constitutive matrix D e is evaluated using the initial stress s 0 and the strain increment aDe. If a = 0, the strain increment De is assumed to be totally plastic and if a = 1, the strain increment De is assumed to be totally elastic. Elastic to plastic transition is characterized by 0 < a < 1, where the elastic fraction of the total strain increment De is given by aDe. In order to solve Eq. 43, the Pegasus algorithm introduced by Dowell & Jarratt (1972) is adopted in the present work (see Appendix B for a detailed description).
(c) Integration of the stress-strain equation: Once the elastic fraction of the stress increment is obtained from the Pegasus algorithm, the plastic fraction of the strain increment D ¢ e is determined using: In order to define initial conditions for the system of equations described by Eq. 39, the stress state s 0 is updated to the stress state at the onset of the plastic yielding s yld and the hardening parameter k 0 is considered at the start of the strain increment D ¢ e where T = 0 and t = t 0 . By using the explicit Euler method, a substepping technique, as proposed by Sloan et al. (2001), is performed over the strain increment D ¢ e to find the stress state and the hardening parameter at the end of D ¢ e , where T = 1. Subincrements of D ¢ e are defined according to the artificial time increment DT, which is calculated taking into account a local error measure obtained by the difference between a second order accurate modified Euler solution and a first order accurate Euler solution. Consequently, the size of the subincrements is automatically modified along the integration process.
In the explicit Euler method, stress state s n and hardening parameter k n at the end of an artificial time increment DT n are obtained as follows: where n-1 and n are referred to previous and current artificial time positions T n-1 and T n , respectively, such that T n = T n-1 + DT n .
In the modified Euler method, stress state s n and hardening parameter k n at the end of an artificial time increment DT n are obtained as follows: where the Euclidian norm s n and the absolute value k n are calculated using Eq. 49. The current subincrement D ¢ e is accepted if both the errors given by Eq. 52 are less than a prescribed tolerance value. The next artificial time increment DT n+1 is obtained taking into account the relative errors of the current artificial time T n , which permits the size of the subincrements to vary during the integration process. At the end of each subincrement, corrections in the stress state and hardening parameter may be needed in order to satisfy the yield condi- tion accurately and avoid yield surface drift. The integration process is carried on until the plastic increment D ¢ e is completely applied, i.e. T n = 1. Suitable tolerance values for all the iterative procedures utilized in the present scheme are found in Sloan et al. (2001), where a sensitive analysis is presented in order to investigate the influence of the tolerance values over the numerical predictions.
For a complete description of the numerical model adopted in this work, readers are addressed to the algorithm presented in Appendix A.

Conclusions
A numerical model based on one-point quadrature and critical state formulation was proposed in this work to simulate the mechanical behavior of soils. The critical state formulation was briefly described using a classical elastoplastic approach where the Modified Cam-Clay model was emphasized. The finite element implementation of the analytical model was performed using the eight-node hexahedral element formulation and reduced integration techniques. A corotational reference system was utilized to stabilize the element formulation as well as to describe stress-strain relations and motion, especially for geometrically nonlinear analysis. The constitutive equation was integrated in this paper using an explicit algorithm with an automatic procedure to split the strain increment into a number of subincrements. Numerical algorithms to determine the yield intersection point, to handle elastoplastic unloading and to restore uncorrected stresses to the yield surface were also presented. One can observe that the formulation presented here leads to a highly efficient numerical model, especially for nonlinear analysis, where the computational effort is substantially reduced when compared to other formulations employing standard approaches. In future works, following this on-going research, the present formulation will be extended to analyze the mechanical behavior of saturated/unsaturated soils.

Appendix A
All numerical procedures adopted in this work may be summarized according to the following algorithm:

Initialization procedures
-Read input data (algorithmic parameters, physical constants and mesh configuration) -Evaluate the transformation matrix R and the gradient matrix B at element level considering the initial geometric configuration in the corotational system. For geometrically linear problems, those matrices will be maintained throughout the analysis. -Evaluate the internal force vector f int considering the imposed initial stress state. Iterative loop: i = i + 1 1) For geometrically nonlinear problems, update the gradient matrix B and the transformation matrix R at element level considering the current geometric configuration in the corotational system and evaluate the Truesdell rate tensor T, which must be added to the constitutive matrix D. 2) Evaluate the tangent constitutive matrix D: in the elastic regime, consider D = D e ; in the plastic regime, consider D = D ep . 3) Evaluate the tangent stiffness matrix K tan in the corotational system including the element stabilization terms. 4) Solve the equilibrium equation in the global system: K if f ( , ) s n n k ¹ 0 then correct the stress state and the hardening parameter to satisfy the yield condition precisely.

End of substepping loop (T n = 1)
11) Update the stress state and the hardening parameter at the end of the strain increment De. 12) For geometrically nonlinear problems, determine the final geometric configuration x n+1 based on the last incremental solution Du n+1,i and update the gradient matrix B and the transformation matrix R at element level considering x n+1 in the corotational system. 13) Compute the internal force vector considering the last updated stress state. 14) Verify convergence of the current load step.

Appendix B
The Pegasus algorithm utilized in this work may be described as follows: 1) Retain the stress state s 0 and the hardening parameter k 0 evaluated at the start of the strain increment De. 2) Set a 0 = 0 and a 1 = 1. £ then leave the iterative loop retaining the last evaluation of a, else go to the next step of the iterative loop. TOL is referred to a tolerance criterion adopted to define the convergence of the numerical procedure. e) If F N presents an opposite sign to F 0 , set a 1 = a and F 1 = F N , otherwise update F 1 using F F F F F N 1 1 0 0 = and set a 0 = a with F 0 = F N .

End of iterative loop
Exit with a.
For elastoplastic unloading, different values must be considered for a 0 and a 1 to initiate the iterative procedure in order to ensure that the Pegasus algorithm finds the intersection point on the yield surface precisely. Therefore, a 0 and a 1 should satisfy the following conditions: