Numerical integration of the Glasgow Coupled Model (GCM)

. This paper presents the derivation of the incremental constitutive relations of the Glasgow Coupled Model (GCM), an established elasto-plastic model for unsaturated soils intended to represent their mechanical and water retention behaviour. The set of incremental equations derived defines an initial value problem (IVP) that needs to be solved at each integration step once the initial (or previous) state and the increments of suction and strain are known. To solve this IVP effectively, a small reformulation of the GCM (not involving any modification of the model, but simply presenting it differently) has been adopted with the aim to facilitate the development of a robust algorithm capable to identify unambiguously the correct model response activated by a given stress path (i.e. elastic, mechanical yielding, water retention yielding or simultaneous mechanical and water retention yielding). Subsequent to the description of this algorithm, the paper presents an explicit substepping integration scheme with automatic error control able to integrate the incremental constitutive equations of the model for each possible response, including unsaturated and saturated conditions. The numerical capabilities of the substepping scheme in terms of accuracy and efficiency are then demonstrated in the context of the performance maps proposed originally for saturated soil models.


Introduction
The Glasgow Coupled Model (GCM) is a coupled mechanical-water retention constitutive model for unsaturated soils developed by Wheeler et al. (2003) [1]. The model was originally proposed for isotropic stress conditions and was then extended to include the effects of deviatoric stresses. Firstly, for the simplified triaxial stress conditions [2] and later for the general stress state conditions [3]. In both cases the formulation proposed assumed the Modified Cam Clay (MCC) as underlying model for saturated conditions. Extensive work on the use of the GCM has been carried out since its development in 2003, demonstrating its capabilities in representing relevant features of unsaturated soil behaviour such as the existence of two unique normal compression planar surfaces for specific volume and degree of saturation when soil states are at the intersection of mechanical (M) and wetting retention (WR) yield surfaces [4], the unification of the representation of the occurrence of plastic compression in unsaturated soils during loading, wetting (collapse compression) and drying (irreversible shrinkage) as yielding on a single M yield curve [5] and the consistent representation of transitions between saturated and unsaturated conditions via the existence of two unique expressions to predict the occurrence of saturation and de-saturation in the soil [6].
While all these findings have facilitated better understanding of the behaviour of unsaturated soils and, more in particular, of the range of applications of the GCM, a robust algorithm capable to numerically integrate the incremental constitutive relations of the GCM in an efficient manner has been a challenging task because of the coupled nature of its mathematical formulation. In this context, [7] presented a complete explicit substepping integration scheme with automatic error control which main features and capabilities are briefed and further investigated in this conference paper, including supplementary numerical examples.

Reformulation of the GCM
With the aim of simplifying the numerical integration of the GCM, certain aspects of the model are reformulated in this section. This reformulation does not involve any modification of the model, simply a change in how it is presented. As in [4], the version of the GCM considered here predicts no elastic changes of (the gradient of elastic scanning curves in the water retention plane is zero i.e. κ s = 0 in the original model of [1]). This consideration is to ensure consistent behaviour across transitions between unsaturated and saturated states [4].

Stress variables
To represent the "average soil skeleton stress" [8] of a soil under unsaturated conditions, the GCM uses the "Bishop's" stress tensor defined by where is the total stress tensor, = (1,1,1,0,0,0) an auxiliary vector, the degree of saturation, the pore air pressure, the pore water pressure.
Similarly, the stress variable for the water retention behaviour is the modified suction defined by: where is the porosity.

Mechanical yield curve
For the simplified stress conditions of the triaxial test, the expression of the mechanical yield curve (M) of the GCM f M is: where p * is the mean Bishop's stress, q is the deviatoric stress, 0 * is the mechanical yield stress (or preconsolidation stress) is the slope of the critical state line in the q:p * plane. Under unsaturated conditions, the mechanical yield stress 0 * increases with decreasing degree of saturation according to the following expression: where ′ 0 is the saturated mehcanical yield stress. 1 and are soil constants as defined in the original [1]. Equation 3 indicates that the mechanical yield curve is elliptical in shape of aspect ratio M when plotted in the q:p * plane (see Figure 1a). The size of this ellipse is given by the current value of the mechanical yield stress 0 * which varies with the degree of saturation according to Equation 4 as illustrated in Figure 1b. In particular, when the soil becomes saturated, the mechanical yield curve reverts to the conventional ellipse of the MCC (as indicated in Figure 1 for S r = 1). Equation 4 resembles the expression proposed by Jommi and Di Prisco [9], with the difference here that the variation of S r is represented within the GCM constitutive framework. Representing the mechanical yield condition in the GCM in terms of Bishop's stresses and degree of saturation means that there is no movement of f M until the soil state reaches the surface (see Figure 1) which differs with the original presentation of the GCM in [1], where coupled movements of the mechanical yield surface (expressed there in terms of Bishop's stresses and modified suction s * ) occur during yielding on water retention yield surfaces. The new formulation presented here has advantages in with regards to the identification of which is the correct model response activated by a stress path. In addition to this, Equation 4 also provides a very simple representation of the transitions between saturated and unsaturated conditions as discussed in Lloret-Cabot & Wheeler [10].

Water retention yield curves
The water retention behaviour is described by two yield functions, the wetting f WR and drying f DR retention. Variations of s * occurring inside f WR and f DR result in no changes of S r (i.e. dS r = d = 0). Yielding on f WR causes plastic increases of S r (i.e. dS r = d > 0) and yielding on f DR produces plastic decreases of S r (i.e. dS r = d < 0). The expression of the wetting retention yield curve is: where s 1 * is the wetting yield stress controlling the occurrence of yielding on f WR (equivalent to p 0 * for mechanical yielding) and varies with the occurrence of mechanical yielding according to: where k 2 is a coupling parameter, p′ 0 is the mechanical hardening parameter and ∆ indicates plastic decreases of specific volume from a reference state. s 10 * and p′ 00 are, respectively, the values of s 1 * and p′ 0 at the reference states when = 0. The expression of the drying retention yield curve is: where s 2 * is the drying yield stress for f DR which varies according to: where s 20 * and p′ 00 are, respectively, the values of s 2 * and p′ 0 when ∆ = 0.

Equation 5 (and Equation 7
) indicate that f WR forms a straight line when plotted in the lns * :lnp 0 ′ plane which is parallel to f DR (see Figure 2). The positions of these straight lines and their gradient with respect to lnp 0 ′ are given by Equations 6 and 8. The current values of the parameters s 10 * and s 20 * (which correspond, respectively, to the values of s 1 * and s 2 * at a reference state with p′ 0 = p′ 00 ) fix the position of f WR and f DR respectively, and the gradient is given by k 2 . Therefore, the parameters s 10 * and s 20 * are referred hereafter as the hardening parameters of the water retention response. Equations 5-8 are still active under fully saturated conditions, because they track the influence of mechanical yielding on the potential occurrence of desaturation on drying (i.e. air-entry point) and re-saturation on wetting or loading (i.e. air-exclusion point).
The spacing between f WR and f DR is assumed constant when plotted in terms of lns * (i.e. s 2 * =R‧s 1 * , where R is a soil constant [4]). This spacing defines the elastic domain of the water retention behaviour (see Figure 3). Yielding on the drying retention yield curve reduces the values of S r and causes a coupled movement of the wetting retention yield curve [1]. Equivalent comments apply when yielding on f WR .

Model responses
To represent the mechanical and water retention behaviour in unsaturated soils, the GCM considers the following six possible responses: Desaturation occurs whilst on f DR and saturation occurs whilst on f WR (see [4], [6], [11] for more details).

Formulation of the problem
The numerical integration of the incremental relationships of a constitutive model defines an initial value problem (IVP) that can be solved knowing the initial (or current) state, the model parameters and, for a strain-driven formulations, a given pair of ∆ɛ and ∆s ( denotes a finite variation). In the context of the GCM, this IVP is written in terms of two mechanical equations (Equations 9 and 10 representing Bishop's stress -strain relations) and two water retention equations (Equations 11 and 12 including relations between modified suction and degree of saturation): where the subscript M indicates mechanical response and the subscript R indicates retention response (f WR or f DR ), Δλ M and Δλ R are the respective plastic multipliers, p 0 ' and R0 * are the respective hardening parameters, aM is the gradient of f M , a R is the gradient of f WR (or f DR ), B M and B R are a scalar functions.

Elastic and elasto-plastic behaviour
Elastic behaviour under saturated or unsaturated conditions is a particular case of the IVP defined by Equations 9-12. Elastic behaviour is represented in the GCM in terms of the secant bulk ̄ and shear ̄ moduli as further described [7].
Yielding on one water retention curve alone involve no mechanical yielding and consequently p 0 ' remains unchanged (i.e. Δλ M = 0). Similarly, yielding on the mechanical yield curve alone means that Δλ R = 0 (and, hence, ΔS r = 0) which simplifies the expression for Δλ M (see [7]). The two remaining cases (simultaneous yielding on f M and f R ) require the derivation of an expression for Δλ M and Δλ R which can be found by imposing consistency condition on both yield curves as detailed in [7].

Identification of the model response activated by the trial stress path
Once the model response is known, all variables are updated using the appropriate set of incremental relations. In such update, the algorithm automatically checks if the stress path intersects a yield curve and finds the corresponding intersection if necessary, using the Pegasus algorithm [12] extensively tested for saturated soil models ([13]- [14]). Figure 3 illustrates the various steps carried out by the algorithm to decide how to integrate the given increments of ∆ɛ and ∆s correctly. The case illustrated corresponds to simultaneous yielding on f M and f DR (a similar algorithm applies for the other cases). A maximum of three different trials is needed to handle correctly this problem, meaning that the algorithm needs to break ∆ɛ and ∆s in three parts. All other cases (i.e. initial stress point on one or two yield curves) are a simplified version of this one and follow a similar logic. Figure 3 is in two parts. Part a shows the full sequence of steps in the lnp 0 ′:lns * plane and Part b illustrates their counterparts in the S r :lnp * plane. The current stress point is indicated by i and is found inside the three yield curves of the model ( i f WR is not plotted for clarity, but its location is to the left of point i). Trial 1 (indicated by t 1 ) is purely elastic (ΔS r = 0 and Δp 0 ' = 0) and ends up outside both i f DR (see Figure 3a) and i f M (see Figure 3b). It is then required to know which of these two yield curves is hit first by trial 1. This problem involves finding two scalars (α 1 for f DR and α 2 for f M ), both between 0 and 1. The scalars indicate the portion of ∆ɛ and ∆s required to move, elastically, the stress point i to the corresponding intersection point (indicated as i R1 for f DR and i M1 for f M ). The lower value between α 1 and α 2 corresponds to the yield curve hit first by trial 1. In the example here, f M is the yield curve hit first (i.e. α 2 < α 1 ). This means that a purely elastic update of the stress point from i to the intersection point i M1 is carried out using the corresponding portion of the given increments (i.e. α 2 ∆ɛ and α 2 ∆s). The next step is to compute Trial 2 (indicated as t 2 ) starting from i M1 (also indicated as i M in Figure 3) and now assuming yielding on only f M . Note that Trial 2 uses the not yet integrated part of the increments of strains and suction i.e. (1-α 2 )∆ɛ and (1α 2 )∆s. Given that yielding on only f M is the model response assumed in computing t 2 , the mechanical hardening parameter p 0 ' is not constant (see Figure 3a) but S r is (see Figure 3b). A second intersection problem arises, now with i f DR (Figure 3). This second intersection problem involves finding a scalar β (also between 0 and 1) that defines the portion of (1-α 2 )∆ɛ and (1-α 2 )∆s required to move, under yielding on only f M , the stress point from i M to i R (also indicated as i Y in Figure 3 as the stress point lies on both yield curves). Once β has been found, the stress point is updated from i M to i R assuming yielding on only f M and using the relevant portion of strain and suction increments i.e. β(1-α 2 )∆ɛ and β(1α 2 )∆s.

Second order accurate explicit substepping integration scheme
This section presents a second order accurate explicit substepping integration scheme for the numerical integration of the GCM (typically referred to as the second order accurate modified Euler with substepping. Given a pseudo-time step/substep i (∆T) with 0 < i (∆T) ≤ 1, the forward Euler and modified Euler updates for σ * , p 0 ′, S r and s R0 * are described in the following by adopting the Butcher tableau (Dormand and Prince, [15]). The coefficients for the two methods are summarised in Table 1. The subscripts i and i+1 denote quantities evaluated at pseudo-times i T and i+1 T= i T + i (∆T) respectively: where the coefficients k b are summarised in Table 1 and the coefficients kj a are summarised in Table 1.

Verification
A numerical test (Test 1) is carried out to study how the error in σ * and S r propagates during the numerical integration. Adopting the same terminology used in [7], the error incurred by the numerical scheme in a single substep (or step in case of no substepping) is indicated by e whereas E indicates the cumulative error over a number of substeps. The numerical test considered here, assumes axisymmetric conditions and an initial unsaturated stress state lying on both f M and f WR , at zero deviatoric stress. The soil constants and initial state considered are summarised in Tables 2 and 3, respectively. The tolerance is assumed equal to 10 -12 .  In particular, the decrement of suction analysed varies from -∆s = 10 -04 to 10.0 kPa (with ∆ε v = 0). Accuracy is assessed by plotting the error in σ * and S r against the size of the input suction variations using logarithmic scales. This representation provides a form of verification of the integration scheme, because the gradient obtained for the best-fitted straight line through a particular set of error results should be in correspondence with the order of accuracy of the numerical integration method (e.g. [16], [17]). Figure 4 illustrates the behaviour of e in terms of Bishop's stress (Figure 4a) and degree of saturation (Figure 4b). Similar results were obtained for the mechanical hardening parameter and the water retention hardening parameter. Symbols indicate the computed relative error and the dashed lines indicate the best-fitted straight line through the computed relative error for the corresponding integration scheme. The respective gradients of each best-fitted straight line plotted in both plots match the expected order of accuracy of the method, suggesting that both substepping schemes work correctly at a single step/substep level. Once a substepping integration scheme has been verified at a single step level, the verification process should study the numerical performance over several substeps in terms of the cumulative error E, as defined in [16] in the context of the performance maps. Various performance maps can be plotted to study specific computational aspects of the substepping integration process. For example, Figure 5 illustrates and demonstrates that, in fact, the variation of cumulative error with the number of substeps follows approximately a straight line of gradient -2 for the modified Euler substepping integration scheme (see [16] for further details).

Conclusions
The formulation of the modified Euler substepping scheme to integrate the incremental constitute relations of the GCM has been presented and its numerical performance has been studied for a simple numerical test involving wetting at constant volume. The results show that both, the error over an individual step/substep and the cumulative error over multiple substeps, behave as expected. This applies to the relative error of both the mechanical (stresses and mechanical hardening parameter) and water retention (degree of saturation and water retention hardening parameter) components of the problem.