Precomputation of Critical State Soil Plastic Models

Abstract

In this paper, a simple precomputing procedure is proposed to improve the numerical performance of the technological application of critical state soil models. In these models, if associated plasticity is assumed, the normalization of the stress space allows both the yield surface and the plastic components of the elastoplastic matrix to be defined as a function of a single variable. This approach facilitates their parameterization and precomputation, preventing the repetition of calculations when the boundary value problems appear at the yield surface with the calculation of plastic strain. To illustrate the scope of the procedure, its application on a modified Cam Clay model is analysed, which shows that the method allows a significant reduction of about 50% (as compared with the conventional explicit integration algorithm) in the computational time without reducing the precision. Although it is intended for critical state models in soils, the approach can be applied to other materials and types of constitutive models provided that parameterization is possible. It is therefore a methodology of practical interest, especially when a large volume of calculations is required, for example when studying large-scale engineering systems, performing sensitivity analysis, or solving optimization problems.

1. Introduction

Simulation of soil yield processes is often computationally expensive, especially when analysing active clays such as MX-80 bentonite [1], FEBEX bentonite [2], or GMZ bentonite [3]. These materials, which usually have a high smectite content, when hydrated unconfined from partially saturated conditions, can undergo deformations of 500%. Under confined conditions, they can lead to swelling pressures approaching 10 MPa. This makes the simulation of their hydromechanical behaviour always computationally demanding, particularly when simulating complex domains or when simulations have to be repeated a large number of times. Even considering the development experienced in recent years by microelectronic technology, which has allowed the multi-core and many-core hybrid heterogeneous parallel computing platform to facilitate a very important advance in computing power, the efficiency of the calculation algorithm continues to be a key issue in the application of massive calculation processes. This is shown, for example, in [4] where the application of the Generalized Likelihood Uncertainty Estimation method in the probabilistic estimation of parameters in hydrology is analysed. In the field of soil mechanics, the improvement of computational performance is especially important in design or parameter identification processes, where the computational time can determine the viability of the study [5].

In some constitutive models, such as the critical state soil models, both the stress space [6] and the plastic terms D^p of the elastoplastic matrix D^ep [7,8,9] can be normalized. Thus, different stress states are associated with the same point of the normalized yield surface and with the same D^p. The result could be that when solving boundary value problems, for a nonnegligible number of times, the calculations associated with the plastic behaviour of the same point (or in the vicinity of the same point) of the normalized stress space are repeated. It is therefore interesting to precompute D^p for a grid of normalised stresses and to interpolate any value within the points of that grid when solving boundary value problems.

This approach is the resolution strategy proposed in the present work. In the following sections, the normalisation of the yield surface and D^p is first described. Then, the CPU time savings that the precomputation entails are evaluated. Later, the criteria are defined to ensure that the proposed method does not compromise the quality of the simulation, and finally, an inspection exercise is performed to evaluate its performance in solving boundary problems.

2. Theoretical Background

In addition to saturated conditions, this paper also assumes axisymmetric conditions. These hypotheses do not reduce the scope of the proposed formulation and yet allow it to be explained more clearly. Furthermore, in line with this quest for clarity, the modified Cam Clay model (MCCm) [6] is adopted as the reference critical state soil model. In this way, the complexity of the constitutive formulation is reduced to a minimum, and the present work is focused on the description of the proposed calculation method. However, the widespread use of the MCCm in computational geotechnics [10,11,12] means that the methodology presented here is not trivial.

Usually, the yield surface f of the MCCm is formulated as

$$f\equiv {\left(\frac{q}{M}\right)}^2-p\left(p_{\mathrm{O}}-p\right)=0,$$

(1)

where q is the von Mises deviatoric stress, M is the slope of the critical state line, and the isotropic preconsolidation pressure p_O is the hardening parameter. This parameter is a good normalisation parameter since by introducing the dimensionless stresses P = p/p_O and Q = q/p_O, f projects onto F (Figure 1), as defined by the expression

$$F\equiv {\left(\frac{Q}{M}\right)}^2-P\left(1-P\right)=0.$$

(2)

All the points (p,q)_i located on the same axis η, although they are associated with a different hardening parameters p_Oi, define the same dimensionless stress (P,Q) (Figure 1). Moreover, in all of them, the vector n, normal to f, given by the expression

$$n=\left(\frac{\partial f}{\partial p},\frac{\partial f}{\partial q}\right)=\left(2p-p_{\mathrm{O}},\frac{2q}{M^2}\right)=p_{\mathrm{O}}\hspace{0.17em}\left(2P-1,\frac{2Q}{M^2}\right)=p_{\mathrm{O}}\hspace{0.17em}N$$

(3)

has the same direction of N, with its modulus varying according to the p_O factor. Therefore, assuming an associated plasticity model, the plastic terms of the elastoplastic matrix will be the same for all yield stresses which project onto the same dimensionless stress. Indeed, if the elastoplastic matrix D^ep defines the increase in the effective stress dσ associated with an increase in the strain dε (the Voigt notation is used for both stress and strain tensors), then we have

$$d\mathsf{\sigma }=D^{\mathrm{ep}}\cdot d\mathsf{\epsilon }$$

(4)

If the consistency equation is applied to the definition of the yield surface f and both the flow rule and the hardening law are considered, then D^ep can be expressed as (see, for example, [13,14])

$$D^{\mathrm{ep}}=D^{\mathrm{e}} · \left(I-\frac{{\left(\frac{\partial f}{\partial \mathsf{\sigma }}\right)}^{\mathrm{t}}\cdot \frac{\partial f}{\partial \mathsf{\sigma }} · D^{\mathrm{e}}}{\frac{\partial f}{\partial \mathsf{\sigma }} · D^{\mathrm{e}}· {\left(\frac{\partial f}{\partial \mathsf{\sigma }}\right)}^{\mathrm{t}}-\frac{\partial f}{\partial x} \cdot \frac{\partial x}{\partial {\mathsf{\epsilon }}^{\mathrm{p}}} \cdot {\left(\frac{\partial f}{\partial \mathsf{\sigma }}\right)}^{\mathrm{t}} } \right)=D^{\mathrm{e}} · \left(I-D^{\mathrm{p}}\right)$$

(5)

where D^e is the elastic strain matrix, I is the identity matrix with the dimension of D^e, the symbol (·)^t indicates the transpose operator, and x defines the vector of hardening parameters. As noted, in the MCCm, x is equal to the scalar p_O. Its value is usually defined in critical state models by the hardening law

$$d{p}_{\mathrm{O}}=\frac{1+e}{\lambda -\kappa } p_{\mathrm{O}} d{\epsilon }_{\mathrm{V}}^{\mathrm{p}}$$

(6)

where λ is the slope of the virgin compression line of the soil, and dε_V^p is the volumetric plastic strain. For the axisymmetric conditions adopted, where dε = (dε_V, dε _S) in which dε _S is the deviatoric strain, D^p can be calculated as

$$D^{\mathrm{p}}=\frac{ \left(\begin{array}{cc}{\left(2p-p_{\mathrm{O}}\right)}^2& \left(2p-p_{\mathrm{O}}\right)\frac{2q}{M^2}{f}_{\nu }\\ \left(2p-p_{\mathrm{O}}\right)\frac{2q}{M^2}& {\left(\frac{2q}{M^2}\right)}^2{f}_{\nu }\end{array}\right)}{ {\left(2p-p_{\mathrm{O}}\right)}^2+f_{\nu }{\left(\frac{2q}{M^2}\right)}^2+\frac{\kappa }{\lambda -\kappa }{p}_{\mathrm{O}}\left(p-p_{\mathrm{O}}\right)}$$

(7)

where f_ν is a function of Poisson’s ratio given by

$$f_{\nu }=\frac{9\left(1-2\nu \right)}{2\left(1+\nu \right)}$$

(8)

Dividing by (p_O)² results in the dimensionless expression

$$D^{\mathrm{p}}=\frac{ \left(\begin{array}{cc}{\left(2P-1\right)}^2& \left(2P-1\right)\frac{2Q}{M^2}{f}_{\nu }\\ \left(2P-1\right)\frac{2Q}{M^2}& {\left(\frac{2Q}{M^2}\right)}^2{f}_{\nu }\end{array}\right)}{ {\left(2P-1\right)}^2+f_{\nu }{\left(\frac{2Q}{M^2}\right)}^2+\frac{\kappa }{\lambda -\kappa }\left(2P-1\right)}$$

(9)

Additionally, D^e is defined as

$$D^{\mathrm{e}}=\left(\begin{array}{cc}K& 0\\ 0& 3G\end{array}\right)$$

(10)

where K is the bulk modulus

$$K=\frac{1+e}{\kappa }p$$

(11)

and G is the shear modulus

$$G=\frac{3\left(1-2\nu \right)}{2\left(1+\nu \right)}K$$

(12)

These moduli depend on the void ratio and are updated as the soil is loaded and deformed. Therefore, the elastic stiffness matrix depends on the strain path followed, as well as the preconsolidation stress, which must be updated according to Equation (6) from the volumetric plastic strain. However, Equation (9) shows that D^p depends on P and Q only. Actually, it depends on only one parameter, since both P and Q are on the surface F = 0 (Figure 1), which, like D^p, can be parametrised by the angle θ shown in Figure 2. Furthermore, in this figure, it is proposed to select N equispaced points (with an Δθ angle between them) at F = 0, calculating D^p_i for each of them. Subsequently, when solving a boundary value problem, the value of D^p is approximated by the interpolation D^p* for angle values between two previously calculated values

$$D^{\mathrm{p}*}\left(\theta \right)=\frac{\left({\theta }_{i+1}-\theta \right)}{\left({\theta }_{\mathrm{i}+1}-{\theta }_{\mathrm{i}}\right)}{D}^{\mathrm{p}}{}_{\mathrm{i}}\hspace{0.17em}+\frac{\left(\theta -{\theta }_{\mathrm{i}}\right)}{\left({\theta }_{\mathrm{i}+1}-{\theta }_{\mathrm{i}}\right)}{D}^{\mathrm{p}}{}_{\mathrm{i}+1}\hspace{0.17em},\hspace{0.17em}{\theta }_{\mathrm{i}}\le \theta <{\theta }_{\mathrm{i}+1}$$

(13)

obtaining D^ep* by

$$D^{\mathrm{ep}*}\left(\theta \right)=D^{\mathrm{e}} · \left(I-D^{\mathrm{p}*}\left(\theta \right)\right)$$

(14)

The calculation procedure proposed in this work is based on precomputing the vector of matrices {D^p₁, D^p_2, …, D^p_N}, thus avoiding repeating the calculation of D^p during the resolution of the boundary value problems.

3. Evaluation of the Precomputation Efficiency

The efficiency of the proposed procedure depends on the cost of calculating D^p. If this cost is not an important part of the total cost of the calculation process, the performance associated with using D^p* might not be of interest. Figure 3a summarises the process followed to solve a computational step when integrating the constitutive model.

Table 1. Equivalent sum costs of the main stages that conform to a conventional explicit integration algorithm and an algorithm with the proposed procedure.

Stage	Sums	Multiplications	Divisions	Computational Cost
De matrix, CDe	3	6	2	17
Dp matrix, CDp	8	25	9	69
Dp* matrix, CDp*	4	16	2	29

lists the number of operations associated with the calculation of D^e and D^p. Although these numbers can vary depending on the way in which the algorithm is programmed, the changes are not significant in terms of the relative computational cost of each magnitude. Using the results of the benchmark defined in Appendix A (where it is indicated that the estimated cost of each multiplication is equivalent to 1.101 sums, with the cost of division being 3.698 sums), the costs shown in the last column of Table 1 were obtained. While it is advisable to be cautious about the reliability in absolute terms of the representativeness of these numbers, they do characterise the weight that each calculation carries in relative terms. Figure 3b shows the strategy followed when applying the precomputational procedure. As noted, D^e must be calculated as in the conventional method. In addition, interpolation of the D^p_i matrices must be made to obtain D^p* by performing the operations indicated in Table 1. Using again the values of Appendix A, the computational cost of Table 1 is obtained, which is approximately 42% of the cost of calculating D^p. The computational time savings are remarkable, and the application of the approach proposed here is therefore of interest.

4. Precomputation Density

To define the precomputation strategy, the minimum number of points N to be taken at F = 0 must be determined in such a way that the accuracy of the calculation when using D^p* is comparable to the accuracy obtained when calculating D^p. To accomplish this goal, a set of constant strain rate paths with a final strain of 20% and with different ratios between volumetric and shear strains, as in Figure 4, were inspected. In all of these inspection exercises, initial spherical conditions were assumed with an effective mean stress of p = p_O/2 = 200 kPa. The analysis was performed for the three soils, namely Weald clay, Klein Belt Ton, and kaolin, which are characterised in

Table 2. Parameters used to characterise the behaviour of the clays assumed in the models (adapted from Schofield and Wroth [15]).

Parameter	Symbol	Weald Clay	Klein Belt Ton	Kaolin
Slope of the critical state line in the p − q effective stress space	M	0.95	0.845	1.02
Slope of the virgin compression line in the ln p − e space	λ	0.093	0.356	0.260
Slope of the unloading–reloading line in the ln p − e space	κ	0.035	0.184	0.050
Void ratio at p = 1 kPa	e1	1.088	3.677	2.764
Poisson’s ratio	ν	0.3	0.3	0.3
Liquid limit (%)	LL	43	127	74
Plastic limit (%)	PL	18	36	42
Plasticity index (%)	PI	25	91	32

. The differences between their properties make the conclusions reached not limited to a particular material.

As a reference, all paths were simulated using an explicit Euler integration method. An intensive substepping process was applied, i.e., the number of calculation steps was divided into smaller substeps successively until the solution differed from the previous solution by only the sixth significant digit, which is a strategy similar to that used by Sołowski and Gallipoli [16]. These results were compared with the precomputed–interpolated solution using the same number of substeps during the explicit integration.

For each path, for each value of α in Figure 4, the number of points N to be computed on the yield surface was determined in such a way that the normalised root-mean-square deviation between the stress solution obtained by explicit Euler integration and the solution reached by calculation with the precomputation–interpolation procedure was less than 10⁰, 10⁻¹, 10⁻², 10⁻³, and 10⁻⁴. Figure 5 shows the values of N obtained for each of these tolerances and for each value of α.

It is important to note that, as illustrated in Figure 6, when imposing a total strain of 20%, in all cases, a significant yield was attained. In other words, trajectories that were highly sensitive to the quality of the precomputation results were analysed.

In all paths, as the tolerance decreased, N increased. This trend can be clearly seen in Figure 7, where for each tolerance value, the maximum value of N obtained for any α was plotted. The behaviour of N was well defined and not erratic; additionally, Figure 7 can be used to define the number of points N to be taken to avoid exceeding a certain value of tolerance.

5. An Inspection on Solving Boundary Value Problems

Although the comparison in Table 1 highlights the interest of the precomputational approach, it does not illustrate its scope clearly enough, because it is associated with the simulation of individual stress–strain paths. Considering a finite element model in which conventionally the stresses are a state function, those paths would illustrate the behaviour of a single Gaussian point. It is reasonable to question what will happen when, in a real problem, the number of Gaussian points is in the hundreds or thousands.

It is not easy to answer this question. Savings will depend on the structure of the finite element code used and the calculation and information management strategies implemented. Each code will have a different response. However, in general, it is reasonable to assume that if N_GP is the number of Gaussian points considered in the mesh, the cost of a conventional explicit integration, χ_C, can be estimated as

$${\chi }_{\mathrm{C}}=N_{\mathrm{GP}}\hspace{0.17em}\left({\eta }_{\mathsf{\Delta }\mathsf{\epsilon }\mathrm{e}}\hspace{0.17em}{C}_{\mathrm{De}}+{\eta }_{\mathsf{\Delta }\mathsf{\epsilon }\mathrm{ep}}\hspace{0.17em}\left(C_{\mathrm{De}}+C_{\mathrm{Dp}}\right)\right)$$

(15)

where the computational costs of D^e and D^p, and C_De and C_Dp have been defined (as well as the computational cost of D^p* and C_Dp*) in Table 1, and for the total number of computational time steps needed to solve the boundary problem, η_Δεe defines the weighted value of the elastic increments, and η_Δεep defines one of the elastoplastic increments. In turn, the cost of the precomputation–interpolation strategy, χ_PI, is estimated as

$${\chi }_{\mathrm{PI}}=N_{\mathrm{GP}}\hspace{0.17em}\left({\eta }_{\mathsf{\Delta }\mathsf{\epsilon }\mathrm{e}*}\hspace{0.17em}{C}_{\mathrm{De}}+{\eta }_{\mathsf{\Delta }\mathsf{\epsilon }\mathrm{ep}*}\hspace{0.17em}\left(C_{\mathrm{De}}+C_{\mathrm{Dp}*}\right)\right)$$

(16)

where η_Δεe* and η_Δεep* are the values analogous to η_Δεe and η_Δεep, respectively, when the integration of the constitutive model is done using the precomputation and interpolation procedure. Their values can, in principle, be different. However, since the discretization of the yield surface is proposed with a very demanding tolerance (10⁻⁴ in the analyses performed above), it is to be expected that η_Δεe* ≈ η_Δεe and η_Δεep* ≈ η_Δεep. In this case, ρ, the ratio between the cost of a traditional explicit integration and the cost of the precomputation–interpolation strategy, is given by the expression

$$\rho =\frac{{\chi }_{\mathrm{PI}}}{{\chi }_{\mathrm{C}}}=\frac{{\eta }_{\mathsf{\Delta }\mathsf{\epsilon }\mathrm{e}}\hspace{0.17em}{C}_{\mathrm{De}}+{\eta }_{\mathsf{\Delta }\mathsf{\epsilon }\mathrm{ep}}\hspace{0.17em}\left(C_{\mathrm{De}}+C_{\mathrm{Dp}*}\right)}{{\eta }_{\mathsf{\Delta }\mathsf{\epsilon }\mathrm{e}}\hspace{0.17em}{C}_{\mathrm{De}}+{\eta }_{\mathsf{\Delta }\mathsf{\epsilon }\mathrm{ep}}\hspace{0.17em}\left(C_{\mathrm{De}}+C_{\mathrm{Dp}}\right)}=\frac{1+\frac{{\eta }_{\mathsf{\Delta }\mathsf{\epsilon }\mathrm{ep}}}{{\eta }_{\mathsf{\Delta }\mathsf{\epsilon }\mathrm{e}}\hspace{0.17em}}\hspace{0.17em}\left(1+\frac{C_{\mathrm{Dp}*}}{C_{\mathrm{De}}}\right)}{1+\frac{{\eta }_{\mathsf{\Delta }\mathsf{\epsilon }\mathrm{ep}}}{{\eta }_{\mathsf{\Delta }\mathsf{\epsilon }\mathrm{e}}\hspace{0.17em}}\hspace{0.17em}\left(1+\frac{C_{\mathrm{Dp}}}{C_{\mathrm{De}}}\right)}$$

(17)

The relative cost is independent of the mesh size. Furthermore, it is necessary to add to the computational costs considered in Equations (15) and (16) the costs associated with the generation of the database {D^p₁, D^p₂, …, D^p_N}. Consequently, it is to be expected that if the mesh is small, the relative importance of the precomputation cost is not negligible. However, the proposed methodology is oriented to large-scale computation, and thus it is a consistent estimator of the efficiency, and for a large computational volume, as shown in Equation (17), it is independent of the mesh size. The ratio η = η_Δεep/η_Δεe expresses the relative importance of the computation of plastic steps in a problem. As expected, given that the precomputation is aimed at optimizing the calculation of the plastic steps, its efficiency is conditioned by the relevance of the plastic strains over the elastic strains, in other words, by the value of η. Using the values of C_De, C_Dp, and C_Dp* from Table 1, the change in ρ with respect to η was plotted and is shown in Figure 8. If η = 0, when the yield is not reached, precomputation is not applied. Only elastic steps are calculated, and both methods are equal (ρ = 1). However, as seen at the end of Section 3, as soon as some plastic steps are required, the precomputation is more efficient (ρ < 1). The larger the volume of the calculation of the plastic steps, the higher the efficiency. When the entire calculation is in the plastic regime (i.e., in normally consolidated conditions), the efficiency tends to a maximum value of 54%.

6. Conclusions

A methodology was presented to parameterise, precompute, and approximate the plastic behaviour of a soil by interpolation. The procedure can be applied to materials whose constitutive models allow the parametrisation and subsequent normalisation of their yield surface as well as the plastic components of the elastoplastic matrix. In this work, the modified Cam Clay model was used as a reference model. After describing the precomputation and interpolation procedure, the efficiency of the method was evaluated by estimating the computational cost of one computational step. Then, a criterion was defined to determine the density of the precomputed database so that the quality of the approximations obtained with the proposed method is like that of solutions based on conventional explicit integration methods. Finally, an inspection was conducted to assess the efficiency of the proposed method in the solution of boundary value problems. It was concluded that the proposed method can reduce the computational burden to as much as 54% of the conventional cost, a significant savings when solving problems in technological applications. This finding is of particular interest in the case of modelling large-scale engineering systems, in sensitivity analysis, or in solving optimisation problems, where a large volume of calculations is required.