# Time homogenization for clays subjected to large numbers of cycles

Aurélie Papon, Zhenyu Yin, Yvon Riou, Pierre-Yves Hicher

► To cite this version:

Aurélie Papon, Zhenyu Yin, Yvon Riou, Pierre-Yves Hicher. Time homogenization for clays subjected to large numbers of cycles. International Journal for Numerical and Analytical Methods in Geomechanics, 2013, 37 (11), pp.1470-1491. 10.1002/nag.2092 . hal-01007385

HAL Id: hal-01007385

https://hal.science/hal-01007385

Submitted on 19 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Time homogenization for clays subjected to large numbers of cycles

A. Papon¹,², Z.-Y. Yin¹,³, Y. Riou⁴ and P.-Y. Hicher¹

¹Research Institute in Civil and Mechanical Engineering, UMR CNRS 6183 - Ecole Centrale Nantes, University of Nantes, 1 rue de la Noë, BP 92101, 44321 Nantes cedex 3, France

²The University of Queensland, Earth Systems Science Computational Centre (ESSCC), School of Earth Sciences, St Lucia, QLD 4072, Australia

³Department of Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

⁴MMGC Department, Ecole Centrale Nantes, 1 rue de la Noë, BP 92101, 44321 Nantes cedex 3, France

This paper discusses the reliability and the efficiency of a time homogenization method employed to reduce the computational time during cyclic loading for two common geotechnical tests and two elastoplastic models for clays. The method of homogenization is based upon splitting time into two separate scales. The first time scale relates to the period of cyclic loading and the second to the characteristic time of the fatigue phenomenon. The time homogenization method is applied to simulate an undrained triaxial test (homogeneous stress state) and a pressuremeter test (nonhomogeneous stress state) under one-way cyclic loading on normally consolidated clay. This method is coupled with two elastoplastic models dedicated to cyclic behavior of clay (a bounding surface plasticity model and a bubble model). Both linear and nonlinear elasticities are considered. The difficulty encountered when applying this method to models introducing nonlinear elasticity and kinematic hardening is pointed out. The performance of time homogenization related to the main parameters is numerically investigated by comparison with conventional finite element simulations.

KEY WORDS: time homogenization; cyclic loading; clay; constitutive models; pressuremeter test; finite element method

## 1. INTRODUCTION

Fatigue in soils may be a significant phenomenon for many structures such as wind power plants, offshore installations, and railways [1]. As permanent strains accumulated during cyclic loading are likely to generate damages on such structures, their design requires the reliable prediction of strains and stresses in soils. To do this, a classical finite element (FE) simulation can be selected. However, it often leads to prohibitive computational time. Moreover, Wichtmann et al. [2] indicated that such simulations are only suitable for a relatively low number of cycles (inferior to 50). Indeed, for higher numbers of cycles, the classical FE simulation may generate systematic numerical errors of the order of the physical strain accumulation.

In order to reduce computational time and sometimes to limit numerical errors, different numerical strategies have been developed. They are classified according to the constitutive model used. Two types can be distinguished: “heuristic” or “phenomenological” models describing only the average behavior without the oscillating part, and models called in this paper “conventional”, based on the theory of plasticity and describing the whole behavior of soil during each cycle.
Based on experimental observations, the pseudo-creep method [3] involves a heuristic model. Indeed, an elastoviscoplastic model is used to represent the average behavior of soil subjected to cyclic loading by simulating a creep test. In this case, the viscous constitutive parameters have to be modified to account for the mapping between the viscous and the cyclic behaviors, or in other words, for the mapping between time and cycles.

Wichtmann and coworkers [2, 4] proposed to use alternatively one conventional and one phenomenological models during one simulation. The first two cycles of the simulation are simulated with the conventional model. Then, based on the results provided by this model, a high-cycle accumulation model is used to predict the strain accumulated during a certain number of cycles without tracing the oscillating strain path. Regularly during the simulation, some control cycles are simulated with the conventional model to update, if necessary, the strain amplitude (for more details see [2, 4]).

Conventional models can be involved in different methods. We refer, in this article, to four main methods related in literature:

* One of them is the cycle jump method, which proceeds in two steps. The first cycle (or more) is simulated and, according to the results, the value of one (or more) variable is extrapolated over a number of cycles initially fixed or updated at each increment. Different versions of this method have been proposed in the literature by Van Paepegem et al. [5], Kiewel et al. [6], and Cojocaru and Karlsson [7].

* Another method was proposed by Lions et al. [8] and later by Farhat and Chandesris [9] using a time-parallel algorithmic scheme to take advantage of parallel processors.

* The rapid cycling creep method introduced by Ponter [10, 11] assumes the decomposition of the total stress into two parts: a cyclic elastic stress $\sigma^{el}$, which would be developed in the case of a purely elastic constitutive model, and a self-equilibrating stress $\rho$ (also called residual stress) caused by the inelasticity introduced by the current constitutive model. If the loading cycle period is small enough compared with the characteristic material time, the residual stress $\rho$ can be considered as constant. Using this assumption, Ponter and Brown [12] proposed an FE scheme completed by Spiliopoulos [13] to solve rapid cycling creep problems.

* Originally developed for multiscale spatial problems by Bensoussan et al. [14] followed by Sanchez-Palencia [15], the asymptotic homogenization method has been adopted by many authors [16–21] to solve multiscale time problems.

The aim of this paper is to test a numerical technique, which is likely to overcome simulation problems related to cyclic loading and could be involved in design computation of geotechnical structures. The asymptotic time homogenization (hereafter referred to as time homogenization) is adopted because it provides a mathematical justification in a general framework and can be applied to many different problems including coupling problems. Compared with the applications proposed in the literature, the present study, dedicated to the cyclic behavior of clay, introduces two geotechnical configurations and two models introducing different constitutive mechanisms (isotropic and kinematic hardenings, and linear and nonlinear elasticities), which are analyzed separately. This systematic study will allow the authors to discuss the use of the time homogenization method for cyclic loading in the area of geotechnical engineering.

In this paper, the principle of this method is first described, so that main theoretical conclusions can be drawn. Then this technique is coupled with two different constitutive models for simulating an undrained triaxial test and a pressuremeter test under one-way cyclic loading. A sensitivity study with respect to the main parameters related to the time homogenization technique is carried out, and the efficiency of the method is tested. This study appears as a precondition for discussing the interest of time homogenization for the design of geotechnical structures.

## 2. PRINCIPLE OF TIME HOMOGENIZATION

Apart from the time parallelism method, all methods mentioned in the introduction of this paper assume a slow evolution of the properties of the material compared with the rapid cyclic loading.
This hypothesis is justified by experimental observations. The time homogenization technique uses this assumption for defining two time scales: the cyclic load needs to be described by an accurate time scale, whereas a coarser time scale is enough for the material evolution under large numbers of cycles.

## 2.1. Definition of two time scales

The aim of the time homogenization is to separate the effects that take place at two time scales. To do so, we define the time characteristic of the cycles $\tau_0$ (the period) and the time characteristic of the fatigue phenomenon $t_r$. The time $t_r$ corresponds qualitatively to the duration necessary for distinguishing the fatigue effects. Therefore, it depends on the material properties. To illustrate the concept of the time characteristic of a phenomenon, we refer to a definition often used in control engineering, when considering a variable $y$ as a function of the time $t$, which evolves according to the relation $y(t) = y_\infty(1 - e^{-t/t_r})$ (see Figure 1). On a graph representing the evolution of the variable $y$ related to fatigue as a function of the time $t$, $t_r$ is defined as the time corresponding to the intersection between the line of the asymptotic behavior and the tangent at the origin.

The ratio of $\tau_0$ over $t_r$ is a small positive value noted as $\zeta$. A second time scale $\tau$ describing the loading can be derived from this ratio $\zeta$ and from the natural time scale $t$, describing the long-term fatigue:

$$ \tau = \frac{t}{\zeta} \qquad (1) $$

The time homogenization method, as explained further, consists in approximating displacements and stresses by using asymptotic expansions with respect to the ratio $\zeta$. Because we further consider only the first order of the asymptotic expansions (see section 2.2), the smaller the ratio $\zeta$, the better the results.

## 2.2. Setting of equations

We consider a spatial domain $\Omega$, whose boundary is noted as $\Gamma$, subjected to volumetric forces $b$. Boundary conditions are applied on, respectively, $\Gamma_f$ and $\Gamma_u$, corresponding to the portions of $\Gamma$ where the stress and displacement conditions, are set, respectively. The displacement conditions are set to zero. We also assume the small strain hypothesis and the decomposition of the total strain tensor $\varepsilon$ into the elastic part $\varepsilon^{el}$ and the inelastic one $\varepsilon^{in}$. The constitutive model is assumed to be linear elastic viscoplastic and the flow rule is given by the operator $B$ according to Guennouni's terminology [17]. We consider an isotropic hardening, whose variable is noted as $\alpha$ and operator is noted as $D$. Note that, in the applications, we will introduce nonlinear elasticity and kinematic hardening.

The loading $f$ can be described by a $\tau_0$-periodic function with respect to $\tau$, and depending on the two time scales $t$ and $\tau$. All the variables of the problem, displacements, stresses, strains, hardening, named

<figure><img src="image_1.png" /><figcaption>Figure 1. Definition of the time characteristic $t_r$ of the fatigue phenomenon.</figcaption></figure>
$\phi$, are assumed to be $\tau_0$-periodic with respect to $\tau$, and depending on $t$ and $\tau$. We also consider that the operators $B$ and $D$ do not depend explicitly on time and that the loading $f$ is given by:

$$f_i(x, t, \tau) = \lambda(t, \tau) f_i^*(x) \quad (2)$$

In order to simplify the following relations, we mention the time and spatial dependencies only when it can be confusing otherwise. Moreover, we use the conventional notation:

$$\phi^{\zeta}(x, t) = \phi(x, t, \tau) \quad (3)$$

Considering this, the problem can be summarized as follows:

$$
\left\{
\begin{array}{l}
\text{Equilibrium equation: } \sigma_{ij,j}^{\zeta} + b_i^{\zeta} = 0 \text{ on } \Omega \\
\text{Constitutive equation: } \sigma_{ij}^{\zeta} = C_{ijkl} : (\dot{e}_{kl}^{\zeta} - \epsilon_{kl}^{\zeta}) \text{ on } \Omega \\
\text{Flow rule: } \dot{\epsilon}_{ijkl}^{\zeta} = B^{\zeta}(\sigma, \alpha) \\
\text{Hardening rule: } \dot{\alpha}^{\zeta} = D^{\zeta}(\sigma, \alpha) \\
\text{Kinematic equation: } \dot{e}_{ij}^{\zeta} = \frac{u_{i,j}^{\zeta} + u_{j,i}^{\zeta}}{2} \text{ on } \Omega \\
\text{Initial conditions: } u_{i}^{\zeta}(t=0) = \tilde{u}_{i} \text{ on } \Omega_{u} \text{ (so that the boundary conditions are satisfied)} \\
\qquad \sigma_{ij}^{\zeta}(t=0) = \tilde{\sigma}_{ij} \text{ on } \Omega_{f} \text{ (so that the boundary conditions are satisfied)} \\
\text{Boundary conditions: } u_{i}^{\zeta} = 0 \text{ on } \Gamma_{u} \\
\qquad \sigma_{ij}^{\zeta} \cdot n_{j} = f_{i}^{\zeta} \text{ on } \Gamma_{f}
\end{array}
\right.
\quad (4)
$$
where $C_{ijkl}$ is the elastic tensor, $u_i$ the components of the displacement vector, and $\sigma_{ij}$ the components of the stress tensor. Note that summation convention is followed.

Note that, in this section dealing with the principle of time homogenization, we do not consider, for sake of simplicity and clarity, the excess pore pressure. Therefore, we do not distinguish the effective stress and the total stress. The following sections dedicated to geotechnical applications will make this distinction.

The variables $\phi$ of the problem are sought as an asymptotic expansion in power of $\zeta$. In fact, it is enough to assume the expansions of the displacement $u$, the inelastic strain $\dot{e}^{\text{in}}$, and the hardening variable $\alpha$ as made by Guennouni [17]. Thanks to kinematic and constitutive equations, it may be possible to conclude for the remaining variables.

$$\phi^{\zeta}(x, t) = \sum_{M=0,1,...}^{\infty} \zeta^M \phi_M(x, t, \tau) \quad (5)$$

where $\phi_M (M=0,1,...)$ are $\tau_0$ periodic with respect to $\tau$.

Once the expansions are inserted into the equations of the problem summarized in Equation (4), we determine the equations corresponding to different orders of $\zeta$ by identifying the coefficients of the polynomials in $\zeta$. This identification can be justified by the fact that all the equations are fulfilled whatever the loading period is, as long as the loading is applied to the same material.

We consider especially a constant volumetric force $b$ and a sinusoidal loading $f$. Therefore, the expansions of the volumetric force and the boundary conditions include only one term (of order $\zeta^0$). Note that the time derivation is given by the chain rule:
$$ \dot{\phi}^{\zeta} = \frac{\partial \phi}{\partial t} + \frac{1}{\zeta} \frac{\partial \phi}{\partial \tau} \qquad (6) $$

The identification of the coefficients of the terms $\zeta^{-1}$, $\zeta^0$, respectively, leads to the systems of Equations (7) and (8), respectively.

$$ \left\{ \begin{array}{l} \frac{\partial \epsilon_{0kl}^{in}}{\partial \tau} = 0 \\ \frac{\partial \alpha_0}{\partial \tau} = 0 \end{array} \right. \qquad (7) $$

$$ \left\{ \begin{array}{l} \sigma_{0ij,j} + b_i = 0 \\ \sigma_{0ij} = C_{ijkl} : (\epsilon_{0kl} - \epsilon_{0kl}^{in}) \\ \frac{\partial \epsilon_{1kl}^{in}}{\partial \tau} + \frac{\partial \epsilon_{0kl}^{in}}{\partial t} = B(\sigma_0, \alpha_0) \\ \frac{\partial \alpha_1}{\partial \tau} + \frac{\partial \alpha_0}{\partial t} = D(\sigma_0, \alpha_0) \\ \varepsilon_{0ij} = \frac{(u_{0i,j} + u_{0j,i})}{2} \\ u_{0i}(t = \tau = 0) = \tilde{u}_i \\ \sigma_{0ij}(t = \tau = 0) = \tilde{\sigma}_{ij} \\ u_{0i} = 0 \quad \text{on} \quad \Gamma_u \\ \sigma_{0ij}.n_j = f_i \quad \text{on} \quad \Gamma_f \end{array} \right. \qquad (8) $$

The system of Equations (7) shows that the $\zeta^0$ order terms of the inelastic part of the strain and of the hardening variable are $\tau$ independent. It means that the inelastic effects are noticeable only at the long-term scale.

We define in Equation (9) the operator $\langle \phi \rangle$, which averages a function $\phi$ over a period $\tau_0$. It allows us to introduce the decompositions given by Equation (10).

$$ \langle \phi \rangle = \frac{1}{\tau_0} \int_{0}^{\tau_0} \phi \cdot d\tau \qquad (9) $$

$$ \begin{aligned} \rho_{0ij}(x, t, \tau) &= \sigma_{0ij}(x, t, \tau) - \langle \sigma_{0ij}(x, t, \tau) \rangle \\ \delta_{0ij}(x, t, \tau) &= \varepsilon_{0ij}(x, t, \tau) - \langle \varepsilon_{0ij}(x, t, \tau) \rangle \\ \chi_{0ij}(x, t, \tau) &= u_{0ij}(x, t, \tau) - \langle u_{0ij}(x, t, \tau) \rangle \end{aligned} \qquad (10) $$

Let us apply this operator to the system of Equations (8) considering the property $\langle \frac{\partial \phi}{\partial \tau} \rangle = 0$, if the function $\phi$ is $\tau_0$ periodic with respect to $\tau$. It leads to a new problem (Pb*) (Equation (11)), called global problem, whose initial conditions will be further specified. Note that the averaging operator has no effect on the terms $\epsilon_{0kl}^{in}$ and $\alpha_0$, because they are $\tau$ independent. It appears that the flow and hardening rules of the time equivalent material correspond to the time averaged rules of the initial material.
$$
\left\{
\begin{aligned}
\langle \sigma_{0ij} \rangle_j + b_i &= 0 \\
\langle \sigma_{0ij} \rangle &= C_{ijkl} : (\langle \epsilon_{0kl} \rangle - \epsilon_{0kl}^{in}) \\
\frac{\partial \epsilon_{0kl}^{in}}{\partial t} &= \langle B(\sigma_0, \alpha_0) \rangle \\
\frac{\partial \alpha_0}{\partial t} &= \langle D(\sigma_0, \alpha_0) \rangle \\
\langle \epsilon_{0ij} \rangle &= \frac{\langle u_{0i,j} \rangle + \langle u_{0j,i} \rangle}{2} \\
\langle u_{0i} \rangle &= 0 \quad \text{on } \Gamma_u \\
\langle \sigma_{0ij} \rangle \cdot n_j &= \langle f_i \rangle \quad \text{on } \Gamma_f
\end{aligned}
\right.
\qquad (11)
$$

The subtraction of the system of Equations (8) from the equations of the problem (Pb*) (Equation (11)) brings to light a second problem (Pb**) (Equation (12)). It is a linear elastic problem independent of (Pb*).

$$
\left\{
\begin{array}{l}
\rho_{0ij,j} = 0 \\
\rho_{0ij} = C_{ijkl} : \delta_{0kl} \\
\langle \delta_{0ij} \rangle = \frac{\langle \chi_{0i,j} \rangle + \langle \chi_{0j,i} \rangle}{2} \\
\chi_{0i} = 0 \quad \text{on } \Gamma_u \\
\rho_{0ij} n_j = f_i - \langle f_i \rangle \quad \text{on } \Gamma_f
\end{array}
\right.
\qquad (12)
$$

Considering the special expression of the loading (Equation (2)), the problem (Pb**) can be solved:

$$
\begin{align}
\chi_{0i}(x, t, \tau) &= (\lambda(t, \tau) - \langle \lambda(t, \tau) \rangle) \cdot \chi_{0i}^*(x) \nonumber \\
\rho_{0ij}(x, t, \tau) &= (\lambda(t, \tau) - \langle \lambda(t, \tau) \rangle) \cdot \rho_{0ij}^*(x) \tag{13}
\end{align}
$$

The following initial conditions, deduced from the solution of the problem (Pb**), can be specified
and complete the formulation of the problem (Pb*).

$$
\begin{align}
\langle u_{0i} \rangle (t=0) &= \tilde{u}_i - \chi_{0i} (t=\tau=0) \nonumber \\
\langle \sigma_{0ij} \rangle (t=0) &= \tilde{\sigma}_{ij} - \rho_{0ij} (t=\tau=0) \tag{14}
\end{align}
$$

Note that the subtraction of the system of Equations (8) from the equations of the problem (Pb*)
(Equation (11)) leads to two additional equations corresponding to the flow and hardening rules.
These equations introduce the terms $\epsilon_{1kl}^{in}$ and $\alpha_1$. In this study, only the terms of order $\zeta^0$ are
sought. Note that the consideration of the terms of order superior to $\zeta^0$ would allow us to get more
accurate results. But considering usual geotechnical uncertainties, it does not seem relevant to take
into account terms of superior orders. However, further practical applications, which are beyond
the scope of this paper, should be carried out to confirm this hypothesis.

2.3. *Theoretical and numerical consequences*

The former development leads to the decomposition of the initial problem into two sub-problems:

* A global problem (Pb*), describing the long-term behavior. The loading corresponds to the time averaged loading, and the constitutive model is said to be time homogenized (i.e., the plastic strain rate and the hardening variable rate correspond to the time average of the initial ones).

* A local problem (Pb**), describing the short-term behavior. This linear elastic problem can be solved independently of the rest of the problem, and the loading corresponds to the oscillating part of the initial loading.
To solve numerically the global problem (Pb*), we can use "big increments" equal to several cycles. For the sake of simplicity, we consider only big increments corresponding to multiple full cycles. However, before solving it, we need to determine the homogenized constitutive model, that is, to determine the rate of plastic strain and hardening variable averaged over a cycle from the initial constitutive model. It leads to the use of hundreds of increments ("small increments") over the first cycle of a big increment. The resolution of the linear elastic problem (Pb**) is independent of the increments used, so that we may select the same coarse partition as the one used for the global problem.

The resolution of the global problem consists, firstly, in determining the plastic strain by linear extrapolation and, secondly, in solving the remaining equations in Equation (11), and especially the equilibrium equation. In section 3.2, we consider an undrained triaxial test under one-way cyclic loading. Using the homogeneity of the test and the undrained condition, the determination of the total strain and of the pore pressure is straightforward. In section 3.3, the FE simulation of a pressuremeter test under one-way cyclic loading with the time homogenization implies a necessary modification of a standard FE code. The structure of the FE code and the FE resolution of the conventional elastoplastic problems (according to the "initial stress method" [22]) have facilitated this numerical work, because the force associated with the plastic strain is determined. Knowing the plastic strain by linear extrapolation and, therefore, the associated force, the final displacement and the final pore pressure are calculated from the elastic problem considering the associated force and the linear elastic relationship. The treatment of nonlinear elasticity is detailed in section 3.2.

The analysis presented in section 2.2 considers the loading period $\tau_0$ and the time characteristic of the fatigue phenomenon $t_r$ as the relevant units of time for the whole domain $\Omega$. However, in the presence of significant plasticity, these values may not be identical at all locations in a boundary value problem. According to section 2.2, the accuracy of the time homogenization technique is based on the value of the ratio $\zeta = \tau_0/t_r$, the largest value of $\zeta$ becoming the relevant one for the whole domain. However, as long as the order of magnitude of $\zeta$ remains the same, the homogenization technique can be applied. This is assumed in the application to the pressuremeter test. Note that phase difference among the different locations does not affect the application of the homogenization technique as long as the ratio $\zeta$ is of the same order of magnitude.

In this section, the principle of the time homogenization technique has been shown. Note that this technique has been developed in a mathematical framework in order to reduce the computational time. It does not consider any rheological aspects. Therefore, if the initial constitutive model does not take into account a given physical phenomenon, then the use of the time homogenization will not include this phenomenon.

### 3. APPLICATIONS

The aim of this section is to test the efficiency of the time homogenization method compared with the conventional FE simulation. Note that analytical solutions would have been more suitable because they would not introduce systematic numerical errors. Unfortunately, the constitutive relationships studied here are too complicated to derive analytical solutions.

Firstly, we describe the two cyclic constitutive models used in this study. Then, we apply the time homogenization method coupled with these models to an undrained triaxial test under one-way cyclic loading. Finally, we simulate a pressuremeter test under one-way cyclic loading using this method implemented in the FE software CESAR-LCPC.

#### 3.1. Constitutive models

In order to validate the time homogenization method, we selected simple constitutive models based on elastoplasticity. Even if these two constitutive models do not account for such phenomena as viscosity, they propose with a relative reduced number of parameters a first approximation of the cyclic behavior of clay. Moreover, both models are based on usual concepts introduced in cyclic constitutive models, like bounding surface and kinematic hardening [23]. It makes their structures close to many other more sophisticated models developed in the literature, for example, [24, 25].
However, even if the selected models are relatively simple, the reasoning proposed in section 2.2 is dedicated to elastic-viscoplastic models and needs therefore to be adapted for elastoplastic models. Indeed for elastoplasticity, the plastic multiplier is obtained through the consistency equation, whereas for elastoviscoplasticity following Perzyna's overstress theory [26], it is obtained through the constitutive definition of the viscous nucleus [27]. However, Guennouni [17] proposed to consider an elastoplastic problem as the limit of a series of elastoviscoplastic problems and established mathematically the same conclusions as the ones obtained for an elastoviscoplastic model. His conclusions are directly applied to the elastoplastic models used in this study.

The elastoplastic models do not introduce any time dependency. Therefore, the application of the time homogenization method requires the consideration of a mapping between time and cycle. Note that we keep using the same terminology (time, long-term, slow, etc.).

### 3.1.1. Bounding surface plasticity model.
The first model refers to the theory of bounding surface plasticity introduced in geotechnics by Dafalias and Hermann [28]. It is a bounding surface plasticity derived from the Modified Cam-Clay model; therefore, it assumes two surfaces. The introduction of a second surface gives a better representation of the cyclic behavior of clays, compared with a classical elastoplasticity with isotropic hardening. Indeed, the bounding surface plasticity represents better the reloading part by allowing for a smooth transition between the elastic behavior and the elastoplastic behavior defined by the classical elastoplasticity.

An intermediate domain is defined between the elastic domain (inside the sub-yield surface, also referred to as the loading surface) and the classical elastoplastic domain (on the yield surface, also referred to as the bounding surface). The two surfaces are homothetic with respect to a projection center $\sigma'_0$. Note that the following $\sigma$ denotes the total stress tensor and $\sigma'$ the effective stress tensor. Figure 2 represents the two surfaces, and the projection center corresponds to the origin of the coordinates. The current stress $\sigma'$ defines the loading surface. An “image” stress $\bar{\sigma}'$ is defined as the projection of the current stress $\sigma'$ on the bounding surface with respect to the projection center (Figure 2).

The plastic strain increment is given by the flow rule:

$$d\epsilon^p = \left\langle \frac{1}{K_p} \frac{\partial F}{\partial \sigma'} d\sigma' \right\rangle \frac{\partial F}{\partial \sigma'} \quad (15)$$

where $\langle x \rangle = x$ if $x \ge 0$, otherwise 0. $\epsilon^p$ is the plastic strain tensor, $F$ the function representing the bounding surface, and $K_p$ the plastic modulus, which will be further specified.

Equation (15) shows that a loading increment pointing to the inside of the loading surface leads to a purely elastic behavior. When the current stress $\sigma'$ is on the bounding surface, the continuity of the plastic strain is guaranteed by Equation (16). In this case, the constitutive model is reduced to the same form as that of the Modified Cam-Clay model.

<figure><img src="image_2.png" /><figcaption>Figure 2. Representation of the bounding and loading surfaces in (p', q) plan.</figcaption></figure>
$$ \frac{1}{K_p} \frac{\partial F}{\partial \sigma'} d\sigma' = \frac{1}{\bar{K}_p} \frac{\partial F}{\partial \bar{\sigma}'} d\bar{\sigma}' \quad (16) $$

where $\bar{K}_p$ is the plastic modulus, when the current stress state $\sigma'$ is located on the bounding surface.

$\bar{K}_p$ can be determined on the basis of the consistency relation applied to the image stress state $\bar{\sigma}'$. The plastic modulus $K_p$ at the current stress state $\sigma'$ depends on the plastic modulus $\bar{K}_p$ at the image stress state $\bar{\sigma}'$ and on the distance between the current and image stress states. Based on the relation proposed by Manzari and Nour [29], the plastic modulus is given by:

$$ K_p = \bar{K}_p + H_0 \cdot \frac{1 + e_0}{\lambda - \kappa} \cdot p_c^{\prime 3} \cdot (\beta - 1) \quad (17) $$

where $H_0$ is a model constant controlling the rate of decay of stiffness with strain for the current stress state within the bounding surface (reloading path): a big value of $H_0$ results in small values of plastic strains (see Equation (15) and [29]). In [29], the order of magnitude considered for $H_0$ ranges from 10 to 10,000. $p_c$ is the isotropic hardening variable defining the size of the bounding surface. $\beta$ is equal to the ratio $p_c/p_{c\sigma}$ (see Figure 2). When the distance $\beta \ge 1$ between the two surfaces decreases, the plastic modulus decreases progressively. Therefore, there is a progressive transition controlled by the cyclic parameter $H_0$ between the elastic behavior and the elastoplastic behavior defined by the classical elastoplasticity, as mentioned previously. The parameters $e_0$, $\lambda$, and $\kappa$, respectively, initial void ratio, slope of the normal consolidation line, and slope of the swelling line, are similar to the Modified Cam-Clay model.

The elastic volumetric behavior of this constitutive model follows the common equation of the Modified Cam-Clay model (Equation (18)), and we assume a constant Poisson's ratio $\nu$.

$$ d\epsilon_v^e = \frac{\kappa}{1 + e_0} \frac{dp'}{p'} \quad (18) $$

where $\epsilon_v^e$ is the volumetric elastic strain and $p'$ the mean effective stress. Equation (18) implies elastic modulus (bulk modulus, shear modulus) depending on the mean effective stress, which leads to energy generation/dissipation in a closed “elastic” cycle. Obviously, this can affect the results obtained for a large number of cycles. In spite of this, following previous studies (e.g., [30–32]), we consider this formulation because when we investigate a Hooke-type linear elasticity, it allows us to have a constant Poisson’s ratio and to derive easily a constant Young’s modulus (more details in section 3.2).

### 3.1.2. Bubble model.
In the previous model, only elastic strains occur during unloading, whereas experiments underline also a plastic behavior. Compared with the previous model, the second model also assumes two surfaces to guarantee a smooth transition between the elastic behavior and the elastoplastic behavior defined by the classical elastoplasticity, but it introduces a kinematic hardening to allow for a better representation of the real behavior (i.e., elastoplastic behavior during unloading). The loading surface in the previous model becomes the kinematic yield surface in this model (Figure 3). The present model is based on the two-surface model, also called bubble model,

<figure><img src="image_3.png" /><figcaption>Figure 3. Representation of the bounding and kinematic yield surfaces in (p', q) plan.</figcaption></figure>
developed by Al-Tabbaa [33]. The two surfaces have the same shape and the ratio between the sizes of the kinematic surface and of the bounding surface noted R is considered as a constitutive parameter. The current stress is located either inside (elastic domain) or on (elastoplastic domain) the kinematic surface.

The bounding surface corresponds to the Modified Cam-Clay yield surface:

$$F_b = \left( p' - \frac{p'_c}{2} \right)^2 + \frac{3(s : s)}{2M^2} - \frac{p'_c^2}{4} \quad (19)$$

where $p'_c$ is the isotropic hardening variable representing the size of the bounding surface, s is the deviatoric stress tensor, and M is the slope of the critical state line.

The kinematic surface has the same shape as the bounding surface, but with different size and location:

$$F_y = (p' - p'_a)^2 + \frac{3((s - s_a) : (s - s_a))}{2M^2} - \frac{R^2 p'_c^2}{4} \quad (20)$$

where $p'_a$ is the mean effective stress and $s_a$ the deviatoric stress tensor at the center of the kinematic yield surface.

The bounding surface evolves with the same isotropic hardening rule as in the Modified Cam-Clay model:

$$\frac{dp'_c}{p'_c} = \frac{1 + e_0}{\lambda - \kappa} d\epsilon_v^p \quad (21)$$

where $\epsilon_v^p$ is the volumetric plastic strain.

The kinematic yield surface evolves according to a combination of isotropic and kinematic hardening rules, as follows:

$$\begin{aligned} dp'_a &= \frac{dp'_c}{p'_c} p'_a + S \left[ \frac{p' - p'_a}{R} - \left( p' - \frac{p'_c}{2} \right) \right] \\ ds_a &= \frac{dp'_c}{p'_c} s_a + S \left[ \frac{s - s_a}{R} - s \right] \end{aligned} \quad (22)$$

where the scalar S is obtained by using the consistency condition for the kinematic yield surface.

The plastic flow rule is assumed to be associated and thus the plastic strain is given by:

$$d\epsilon^p = \left\langle \frac{1}{A_0 + A_1} \frac{\partial F_y}{\partial \sigma'} d\sigma' \right\rangle \frac{\partial F_y}{\partial \sigma'} \quad (23)$$

where $A_0$ is the plastic modulus given by the Modified Cam-Clay model and $A_1$ is the plastic modulus given by an interpolation rule.

Al-Tabbaa [33] specified that the interpolation rule is not unique. Based on the works of Hashiguchi [34], Al-Tabbaa proposed a plastic modulus $A_1$ depending on the measure of the distance b between the kinematic yield surface and the bounding surface. Therefore, when the two surfaces are in contact, the value $A_1$ is equal to 0 and only the plastic modulus $A_0$ is mobilized. In this study, the plastic modulus $A_1$ is given by:

$$A_1 = \frac{4(1 - e_0)}{\lambda - \kappa} \psi \left( \frac{b}{b_{\text{max}} - b} \right)^{\xi} \left( \frac{p'_c}{2} \right)^3 R^2 \quad (24)$$

where $\psi$ is a model constant controlling the rate of decay of stiffness with strain for the stress state on the kinematic yield surface and within the bounding surface (unloading and reloading paths), similar to $H_0$ in Equation (17). $\xi$ is a model constant, superior to 1, controlling the variation of the plastic
hardening modulus for the stress state on the yield surface and within the bounding surface (unloading and reloading paths): the higher the value of $\xi$, the smoother the transition between the behaviors on the kinematic yield surface within and on the bounding surface. $b_{max}$ is the maximal value of $b$ equal to:

$$b_{max} = p'_{c}(1 - R) \qquad (25)$$

According to Al-Tabbaa [33], the measure of the distance $b$ between the kinematic yield surface and the bounding surface is given by the dot product of the outward normal $n$ to the kinematic surface at the current stress state, and the vector $\beta$, which links the current stress state (point A in Figure 3) to the stress state on the bounding surface with the same outward normal (point B in Figure 3).

As noted for the first model (Equation (17)), the definition of the plastic modulus (Equation (24)) guarantees the smooth transition between the elastic behavior and the elastoplastic behavior defined by the classical elastoplasticity. The elastic behavior of this model is the same as in the first model.

The model parameters can be separated into two categories: the parameters related to the monotonic behavior ($v, \kappa, \lambda, e_0, M, p'_{c0}$) and those related to the cyclic behavior ($R, \psi, \xi$). Note that $p'_{c0}$ is the initial value of $p'_{c}$.

For the implementation of this constitutive model, we follow the concept of changing the active yield surface introduced by Grammatikopoulou et al. [24].

In the following section, we present successively various applications of increasing complexity, so that the difficulties related to the use of the time homogenization method can be detected in a systematic way.

## 3.2. Simulations of undrained triaxial test under one-way cyclic loading

The simulation of cyclic triaxial tests does not cause any major computational time problem. However, it seems important to consider such a homogeneous test for a first estimation of the efficiency of the time homogenization.

Even though the aim of this paper is purely numerical, the choice of loading and the determination of the constitutive parameters are based on experimental data [35]. Undrained triaxial tests carried out on normally consolidated clays under one-way cyclic loading (maximal deviatoric stress $q_{max} = 65$ kPa and initial mean effective stress $p' = 200$ kPa) are simulated. The parameters related to the monotonic behavior are kept constant and their values are given in Table I. In the special case of one-way cyclic test, the evolutions of the axial permanent strain and of the permanent excess pore pressure are characteristic of the fatigue phenomenon. The permanent strain, respectively the permanent excess pore pressure, is defined as the strain, respectively the excess pore pressure, remaining after a cycle.

The analysis presented in section 2.2 specifies a linear elasticity, whereas both constitutive models introduce a nonlinear elasticity. In order to analyze the effects of this nonlinearity on the time homogenization method, we propose to study successively both cases. Firstly, the elastic model is linearized and transformed into a Hooke's model with constant Poisson's ratio and Young's modulus. The elasticity tensor therefore remains constant and is set to its initial value, depending on the initial value of $p'$ and on the values of $\kappa, v$ and $e_0$ (see Equation (18)). Secondly, the initial elastic model is used and the constitutive equation of the global problem (Pb*) is modified as:

<table><caption>Table I. Modified Cam-Clay parameters for the simulation of triaxial tests.</caption><thead><tr><th>v</th><th>κ</th><th>λ</th><th>e<sub>0</sub></th><th>M</th><th>p'<sub>c0</sub> (kPa)</th></tr></thead><tbody><tr><td>0.3</td><td>0.08</td><td>0.23</td><td>1.15</td><td>1</td><td>200</td></tr></tbody></table>
$$ \langle \sigma_{0ij} \rangle = \langle C_{ijkl} \rangle : (\langle \epsilon_{0kl} \rangle - \epsilon_{0kl}^{in}) \qquad (26) $$

According to Equation (26), the elasticity tensor is determined after a preliminary simulation of the first cycle of each big increment: the elasticity tensor for simulations with homogenization is equal to its average over the first cycle of the big increment (see section 2.2).

The use of time homogenization requires the crucial determination of the number of cycles included in a big increment used for the resolution of the global problem (Pb*) (see section 2.2), in other words the “size of the cycle jump”. This integer is noted in the following $N_c$ and is considered as a measure of the computational time gain, because the calculations with homogenization are about $N_c$ times faster than the ones without. The determination of $N_c$ therefore results from a compromise between the computational time and the accuracy of the results.

Another crucial parameter of the time homogenization is the ratio $\zeta$ of the loading period $\tau_0$ over the time $t_r$ characteristic of the fatigue phenomenon. As mentioned in section 2.1, the smaller the $\zeta$, the more accurate the results. Given that the constitutive models are time independent, we expressed both parameters $\tau_0$ and $t_r$ as number of cycles. The period $\tau_0$ remains certainly equal to one cycle, so that the study of the ratio $\zeta$ comes to the study of $t_r$.

Consequently, the sensitivity of the two parameters $t_r$ and $N_c$ towards the accuracy of the resolution by homogenization, compared with the conventional FE resolution, is investigated in the following subsections 3.2.1 and 3.2.2.

Note that the first 10 cycles were conventionally simulated because the first cycles can be considered as irregular according to the important nonlinear rate of the permanent strain.

### 3.2.1. Coupling of homogenization with the bounding surface plasticity model

Cyclic constitutive parameters and characteristic time $t_r$. The bounding surface plasticity model introduces one new parameter in comparison with the Modified Cam-Clay model, namely the parameter $H_0$, which characterizes the behavior during reloading. This parameter influences the value of the time characteristic $t_r$ as detailed in the following paragraph. Therefore, even if the sensitivity study concerns explicitly the parameter $H_0$, it also deals implicitly with the characteristic time $t_r$. To determine the reference value $H_0^{ref} = 3000$, we used experimental results [35] and consider two additional values: 1000 and 10,000. Figure 4 (respectively 5) represents the evolution of the axial permanent strain (respectively the evolution of the permanent excess pore pressure) with the number of cycles. Note that a big value of $H_0$ results in small values of plastic strains, as mentioned in subsection 3.1.1.

The time characteristic is determined by estimating the time at which a change appears in the behavior of the material: from a rapid evolution of the fatigue phenomenon to a slow one. Figure 5 shows that the evolution of the permanent excess pore pressure stabilizes after a certain number of cycles and it seems that the asymptotic value is the same for the three different values of $H_0$.

<figure><img src="image_4.png" /><figcaption>Figure 4. Evolution of the permanent axial strain as a function of the number of cycles assuming the nonlinear elastic bounding surface plasticity model during the simulation of a one-way cyclic undrained triaxial test (conventional simulation).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 5. Evolution of the permanent excess pore pressure as a function of the number of cycles assuming the nonlinear elastic bounding surface plasticity model during the simulation of a one-way cyclic undrained triaxial test (conventional simulation).</figcaption></figure>

Therefore, the evolution of the permanent excess pore pressure allows for a reliable determination of the time characteristic $t_r$, which corresponds to the intersection between the tangent at the origin and the asymptote. Graphically, we observe similar characteristic times for both versions of the model (linear and nonlinear elasticity). The value of $t_r$ is equal to 400, 800, and 1800 cycles for $H_0$ = 1000, 3000, and 10,000, respectively. According to section 3.2, the highest value of $t_r$ leads to the smallest value of $\zeta$, and, therefore, to the most accurate solution.

**Linear elasticity.** Figures 6 and 7 allow us to compare the evolutions of the axial permanent strain and of the permanent excess pore pressure in the case of simulations with and without homogenization for the three values of $H_0$. The parameter $N_c$ is set to 200. The relative error between conventional simulations and simulations with homogenization tends to decrease during the simulation. In the unfavorable case ($H_0$ = 1000, $\zeta$ = $2.5 \times 10^{-3}$), the relative error on the axial permanent strain is equal to 94% at 200 + 10 cycles, to 14% at 2000 + 10 cycles, and to 1.5% at 30,000 + 10 cycles. The axial permanent strain is overestimated, which is in agreement with what should be expected. Indeed, in this model, plastic strain occurs only during loading and the permanent strain rate decreases during the simulation. Therefore, because the time homogenization linearizes the evolution of the plastic strain, it overestimates automatically the permanent strain. The permanent excess pore pressure is also overestimated; it is due to similar causes.

Simulations with and without homogenization for different values of $N_c$ are represented in Figures 8 and 9. The parameter $H_0$ is set to 1000 (unfavorable case, $\zeta$ = $2.5 \times 10^{-3}$). As expected, the smaller the $N_c$, the smaller the relative error between conventional simulations and simulations with homogenization. Figures 8 and 9 highlight also the linear extrapolation of the plastic strain specified by the development in section 2.2 and the progressive correction of the plastic strain increment.

<figure><img src="image_6.png" /><figcaption>Figure 6. Comparison between the simulations of a one-way cyclic undrained triaxial test without and with homogenization ($N_c$ = 200): evolution of the permanent axial strain assuming the linear elastic bounding surface plasticity model.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 7. Comparison between the simulations of a one-way cyclic undrained triaxial test without and with homogenization ($N_c = 200$): evolution of the permanent excess pore pressure assuming the linear elastic bounding surface plasticity model.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 8. Comparison between the simulations of a one-way cyclic undrained triaxial test without and with homogenization ($H_0 = 1000$): evolution of the permanent axial strain assuming the linear elastic bounding surface plasticity model.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 9. Comparison between the simulations of a one-way cyclic undrained triaxial test without and with homogenization ($H_0 = 1000$): evolution of the permanent excess pore pressure assuming the linear elastic bounding surface plasticity model.</figcaption></figure>

**Nonlinear elasticity.** The same study was performed using the previous model with nonlinear elasticity. It could be observed as previously that the results converge towards a limit solution. However, in the case of nonlinear elasticity, the limit solution does not correspond to the solution without homogenization. Figure 10 representing the evolution of the axial permanent strain with the number of cycles shows a bias between conventional simulations and simulations with homogenization. This gap is related to the incorrect capture of the nonlinear elastic behavior by the time homogenization. This error is inevitable whatever the solving method used. Even for $N_c = 1$, a
<figure><img src="image_10.png" /></figure>

Figure 10. Comparison between the simulations of a one-way cyclic undrained triaxial test without and with homogenization: evolution of the permanent axial strain assuming the nonlinear elastic bounding surface plasticity model.

bias appears and seems to stabilize after 1000 cycles, because of the stabilization of the elasticity matrix with the mean effective stress.

*Rule of thumb.* The aim of this subsection is to propose a rule of thumb easy to apply by the user. Knowing the value of ζ, the user can refer to this rule to set the value of N<sub>c</sub> given the relative error estimated at t<sub>r</sub> cycles. Note that this kind of rule of thumb is a priori applicable only for this model. Figures 11 and 12 represent the relative error on the permanent axial strain at t<sub>r</sub> cycles as a function of the product ζ N<sub>c</sub><sup>2</sup> with linear and nonlinear elasticity, respectively. Note that this product

<figure><img src="image_11.png" /></figure>

Figure 11. Relative error on the permanent axial strain at t<sub>r</sub> cycles as a function of the product ζ N<sub>c</sub><sup>2</sup> assuming a linear elasticity.

<figure><img src="image_12.png" /></figure>

Figure 12. Relative error on the permanent axial strain at t<sub>r</sub> cycles as a function of the product ζ N<sub>c</sub><sup>2</sup> assuming a nonlinear elasticity.
corresponds to the ratio of $N_c$ over the number of increments required to reach $t_r$. The evolution of the relative error can be considered as linear. The difference between the Y-intercepts is caused by the bias related to the nonlinearity of the elasticity. Note that the Y-intercept for the linear elastic case is only slightly different from zero (1.7%), which is satisfactory given the basic rule proposed.

### 3.2.2. Coupling of homogenization with the bubble model.
The aim of this section is to investigate if a different constitutive model can influence the efficiency of the time homogenization method. We consider the bubble model first with a linear elasticity and then with a nonlinear elasticity. In section 2.2, a kinematic hardening is not considered in the mathematical development of the homogenization technique. A linear extrapolation of the center of the kinematic surface is not possible, because the theory of plasticity requires that the stress state remains within or on the yield surface, which is not guaranteed if the stress state and the hardening variables are independently updated. In this work, the stress state at the end of each cycle is assumed to be on the surface (plasticity) and is used to update the center of the kinematic surface (for more details see [36]). Note that for triaxial tests, this assumption is true, so that the simulations without homogenization and with homogenization for $N_c = 1$ are identical.

*Cyclic constitutive parameters and time characteristic t<sub>r</sub>.* The bubble model introduces four cyclic parameters. However, in this study, we examined only the sensitivity of the model’s response with respect to parameter $\psi$, whose influence is similar to the parameter $H_0$ in the bounding surface plasticity model as shown Figure 13. The values of the three other cyclic parameters are given in Table II and are in agreement with experimental results [35]. For the sensitivity study with respect to $\psi$, the reference value is set to $10^{10}$, and we consider two additional values $10^9$ and $10^{11}$. Note that the value of $R$ is in agreement with values found in the literature [24].

The rule given in the subsection 3.2.1 for the estimation of the characteristic time $t_r$ is applied: the values of $t_r$ corresponding to the different values of $\psi$ are equal to 175, 250, and 375 cycles for both linear and nonlinear elasticity and, therefore, are smaller than the values obtained with the first model. To set the value of $N_c$, we use the rule of thumb determined from the results given by the first model as a basic approximation. $N_c$ is set to 100, so that the relative error on the permanent axial strain at $t_r$ cycles would remain inferior to 35% for linear elasticity and to 80% for nonlinear elasticity.

*Linear elasticity.* Figure 14 represents the evolution of the axial permanent strain for three values of $\psi$. The axial permanent strain given by the homogenization technique is alternatively overestimated and underestimated so that the relative error on the axial permanent strain between conventional simulations and simulations with homogenization decreases, then increases. This error is equal to

<figure><img src="image_13.png" /><figcaption>Figure 13. Evolution of the permanent axial strain as a function of the number of cycles assuming the nonlinear bubble model during the simulation of a one-way cyclic undrained triaxial test (conventional simulation).</figcaption></figure>

<table><caption>Table II. Cyclic parameters of the bubble model for the simulation of triaxial tests.</caption><thead><tr><td>ψ</td><td>ξ</td><td>R</td></tr></thead><tbody><tr><td>10<sup>9</sup>, 10<sup>10</sup>, 10<sup>11</sup></td><td>8</td><td>0.1</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>Figure 14. Comparison between the simulations of a one-way cyclic undrained triaxial test without and with homogenization: evolution of the permanent axial strain assuming the linear elastic bubble model.</figcaption></figure>

66% at 100+10 cycles, to 13% at 2000+10 cycles, to -5% at 10,000+10 cycles and to -9% at 30,000+10 cycles. This numerical observation may be explained by the fact that the bubble model allows plastic strain to occur during both loading and unloading (contrary to the bounding surface plasticity model). Indeed, because of the linear extrapolation of the plastic strain, the permanent strain is overestimated during loading and underestimated during unloading. If $N_c$ is too high, the balancing process caused by the time homogenization is relatively instable: if the first overestimation is too large, then it leads to an underestimation. According to the importance of each stress path, the permanent strain can be overestimated or underestimated. For the design of structures by means of time homogenization and for security reasons, the user has to be aware of this possible underestimation of permanent strains when using a more sophisticated model and, therefore, he or she has to select carefully the value of $N_c$.

**Nonlinear elasticity.** The original bubble model combines a nonlinear elasticity and a kinematic hardening so that we can expect a bias caused by the nonlinear elastic behavior and a possible overestimation or underestimation of the axial permanent strain. Figure 15 shows a better approximation given by the time homogenization method for $N_c$ = 100 with a nonlinear elasticity than with a linear elasticity (Figure 14). Therefore, it seems that the two kinds of errors balance each other: overestimation caused by the nonlinearity of the elasticity and underestimation of the permanent strain caused by the kinematic hardening. Therefore, without any additional numerical experiences, it seems difficult to predict the quality of the approximation.

**Rule of thumb.** In this subsection, we set the value of $N_c$ considering the rule of thumb derived from the first model. Given the previous results related to the second model, we can check a *posteriori* in what extend this *a priori* rule of thumb can be used. The evolution of the relative error on the permanent axial strain at $t_r$ cycles as a function of the product $\zeta N_c^2$ assuming a linear and a nonlinear elasticity is represented in Figures 11 and 12, respectively. The evolution can be

<figure><img src="image_15.png" /><figcaption>Figure 15. Comparison between the simulations of a one-way cyclic undrained triaxial test without and with homogenization: evolution of the permanent axial strain assuming the nonlinear elastic bubble model.</figcaption></figure>
considered linear only for $\zeta N_c^2$ superior to 25. Considering this range of values, we observe a similar slope for both models with linear and nonlinear elasticity (about 0.006). For $\zeta N_c^2$ inferior to 25, the relative error tends to zero for both versions of the bubble model, whereas for the bounding surface plasticity model, only the case with linear elasticity tends to zero. This fact is in agreement with the previous comments. Indeed, for linear elasticity, the time homogenization for $N_c = 1$ leads to the same results as the ones given by the conventional method. A nonlinear elasticity introduces a necessary bias between the simulations without homogenization and with homogenization for $N_c = 1$. For the first model, this bias remains, even for $\zeta N_c^2 << 1$, whereas for the second model, it is balanced by the introduction of the kinematic hardening.

### 3.2.3. Comments.
Four applications to undrained triaxial tests with different complexities have been presented in this section. Because they do not lead to prohibitive computational time, these applications allow us to test the efficiency of the time homogenization method and to point out the major difficulties encountered.

In simulations with isotropic hardening and linear elasticity, the time homogenization technique is very effective and can be used without any specific recommendations. The introduction of nonlinear elasticity leads to a systematic bias, and the application to the bubble model shows possible overestimation and underestimation of permanent strains.

To guarantee a certain reliability of the results, the user has to carry out in every instance a sensitivity study with respect to the parameter $N_c$. Note that this sensitivity study can be reduced by establishing rules of thumb, which can be improved with the increasing number of experiences of the user. Such study is shown in Figures 14 and 15. Eventually, the user has to find the best compromise between the accuracy of the approximation given by the homogenization method and the computational time. Note that the latter depends not only on the parameter $N_c$ but also certainly on the total amount of cycles simulated, which can be determinant.

## 3.3. Simulations of pressuremeter test under one-way cyclic loading

### 3.3.1. Numerical modeling.
Following the previous systematic study, we consider in this section the application to a cyclic *in situ* test, for which computational time may be prohibitive. Because the pressuremeter test is a boundary value problem, which requires the use of an FE software, the time homogenization resolution has been implemented in the FE software CESAR-LCPC as explained in section 2.3. The implementation has been validated by simulating undrained triaxial tests and by comparing the results with the ones presented in the previous section.

To investigate the effects of the time homogenization in a non-homogeneous test, we selected the bubble model. Indeed, the introduction of kinematic hardening makes this model more interesting for design purposes. However, in the framework of this study, we decided to linearize the elasticity (see section 3.2), which allows us to avoid the systematic bias caused by the nonlinear elasticity and to focus on the resolution of a non-homogeneous problem. Moreover, by doing so, we prevented a possible balancing effect between the kinematic hardening and the nonlinearity of the elasticity highlighted in subsection 3.2.2.

Considering the necessary computational time using a conventional FE resolution, only one set of constitutive parameters was considered. Table III summarizes the values of the parameters. Moreover, for the same reason, the reference simulation without homogenization is limited to 100 cycles. However, even if the comparison between simulations with or without homogenization is reduced to 100 cycles, it allows us to draw some preliminary conclusions. A one-way cyclic pressuremeter test without any water flow on a normally consolidated clay was simulated. Owing to the axisymmetry of the problem, a two-dimensional FE model was sufficient. Moreover, plain-strain

<table><caption>Table III. Parameters of the linear elastic bubble model for the simulation of the cyclic pressuremeter test.</caption><thead><tr><th>v</th><th>E (kPa)</th><th>λ</th><th>e<sub>0</sub></th><th>M</th><th>p'<sub>c0</sub> (kPa)</th><th>R</th><th>ψ</th><th>ξ</th></tr></thead><tbody><tr><td>0.3</td><td>1612.5</td><td>0.23</td><td>1.15</td><td>1</td><td>50</td><td>0.1</td><td>10<sup>10</sup></td><td>8</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>Figure 16. Finite element model of the pressuremeter test.</figcaption></figure>

condition was assumed in the vertical direction. Eight-node rectangular elements were used, so that the model contained 891 nodes and 252 elements (see Figure 16). The ratio of the outer diameter to the inner diameter of the model was taken sufficiently high to model the condition of infinite medium [37]. The loading was stress controlled and the maximal increase of the probe pressure $Δp$ was equal to 25 kPa.

According to the results given by the triaxial tests, the characteristic time $t_r$ is equal to 250 cycles ($\zeta = 4 \times 10^{-3}$). Because the number of simulated cycles is limited to 100 ($< t_r = 250$ cycles), we cannot use directly the rule of thumb defined previously to determine $N_c$. Considering the results given in Figure 14, we first set $N_c$ to 25. This first attempt led to instabilities at the beginning of the simulation, showing that the correction by the homogenization is more sensitive in a non-homogeneous simulation. Therefore, we set $N_c$ to 10, which means that the calculations with homogenization are about 10 times faster than the ones without.

### 3.3.2. Results.
Figure 17 represents the evolution of the ratio $u_{\text{perm}}(a)/a$, where $a$ is the radius of the probe and $u_{\text{perm}}(a)$ the displacement of the wall at the end of a cycle, for conventional simulations and simulations with homogenization ($N_c = 1$ and $N_c = 10$).

In this case, the simulations without and with homogenization for $N_c = 1$ are not identical. The relative error on the ratio $u_{\text{perm}}(a)/a$ is maximum (7%) at the 100th cycle. The gap appears at the second increment solved by homogenization and is related to the update of the center of the kinematic surface (see subsection 3.2.2). Indeed, the update of this center assumes a plastic behavior at the end of each cycle for each integration point, whereas this condition is not satisfied because of the non-homogeneity of the problem (more details are given in [36]).

To compare this result with the previous ones, we consider the simulation of the triaxial test assuming the bubble model with the smallest value of $\zeta N_c^2$ (i.e., 3.6), knowing that the simulation of the pressuremeter test leads to $\zeta N_c^2 = 0.4$. For the simulation of the triaxial test, the error on the permanent axial strain at 100 cycles is equal to 5%. Given the important ratio between the two values of $\zeta N_c^2$, we would expect a relative error on the ratio $u_{\text{perm}}(a)/a$ inferior to 5% for the pressuremeter test. The higher error observed (9.7%) is caused by the non-homogeneity of the problem as explained previously. Therefore, the rules of thumb established for triaxial tests cannot at the present state be directly extrapolated to pressuremeter tests and a fortiori to geotechnical structures. However, further numerical and experimental studies dedicated to different non-

<figure><img src="image_17.png" /><figcaption>Figure 17. Comparison between the simulations of a one-way cyclic undrained pressuremeter test without and with homogenization: evolution of the ratio $u_{\text{perm}}(a)/a$ assuming the linear elastic bubble model.</figcaption></figure>
homogeneous tests could allow the user to improve these rules and, thus, to increase the reliability of
the simulations with homogenization.

4. CONCLUSIONS

Considering the significant consequences of cyclic loading on many geotechnical structures, the
design of such structures requires reliable and efficient numerical methods. Because of the
prohibitive computational time related to FE simulations of millions of cycles, numerical
strategies have been developed. In this study, the authors considered the time homogenization
method, because of its mathematical justification in a general framework and its possible
application to many different problems. Two constitutive models dedicated to cyclic behavior
of clay (a bounding surface plasticity model and a bubble model) were selected and applied for
modeling two common geotechnical tests. The results point out the sensitive parameters related
to the time homogenization technique and the difficulties encountered when applying this
method to particular constitutive models. The time homogenization process has been clearly
defined by Guennouni [17] if the constitutive model introduces pure isotropic hardening and
linear elasticity. The introduction of nonlinear elasticity in both models and the consideration
of kinematic hardening in the bubble model require new strategies for the update of some
variables. Especially, the elasticity matrix is updated by its average over the first cycle of the
big increment and the center of the kinematic surface is updated so that the stress state remains
on the kinematic surface. The rapidity and the accuracy of the simulations using the
homogenization method were compared with the conventional FE simulations.

Firstly, an undrained triaxial test under one-way cyclic loading was simulated with and without
homogenization. The parameter Nc controls the computational time gain, as the calculations with
homogenization are about Nc times faster than the ones without. The value of Nc depends on the
time characteristic of the fatigue phenomenon tr and on the accuracy of the approximation expected
after a certain number of cycles. As long as the constitutive model introduces pure isotropic
hardening and linear elasticity, the time homogenization method leads in a short computational time
to very satisfactory results. However, the introduction of nonlinear elasticity in both models and the
consideration of kinematic hardening in the bubble model generate disturbances, even if the results
remain acceptable considering the reduction of the computational time.

Secondly, the time homogenization method was implemented in FE software. This numerical
work was facilitated by the similarities between the resolution by homogenization and the initial
stress method. To stay close to design purposes, a one-way cyclic pressuremeter test was simulated,
using the linear elastic bubble model. The results are affected by the strategy used for the updating
of the kinematic variables, which is less efficient for the pressuremeter test than for the triaxial test.
Therefore, further studies on the determination of the kinematic variables have to be performed in
order to reduce the error caused by the time homogenization.

To help the user applying the time homogenization technique and especially to set the value of Nc,
rules of thumb could be established. In this study, we derived a rough rule for each model from the
results given by the triaxial tests. These rules have to be extended by considering different stress
levels as well as two-way cyclic loading, and to be developed for different materials in order to
finally be able to deal with complex boundary value problems. Only the comparison with
experiments currently allows us to estimate the numerical error introduced by the homogenization
technique. Considering the results given by this first study, it seems difficult up to now to predict if
a sensitivity study with respect to Nc and the use of the rules of thumb will be sufficient to reduce
efficiently the computational cost. The progressive increase of the parameter Nc during the
simulations would improve the efficiency of the method. However, it has been disregarded in this
first approach because of the necessary introduction of an additional parameter. A simple and
systematic study has been preferred.

In light of the conclusions of this study, we believe that the time homogenization method warrants
further investigation. The joint development of constitutive models and of numerical strategies may
become a necessary condition for an efficient implementation. This work represents only the first
stage of the study of structures subjected to large numbers of cycles. The reliability of the time homogenization method in geotechnics can only be improved by more feedback derived from numerical and experimental comparisons.

## ACKNOWLEDGEMENTS

The authors would like to thank the French institute of sciences and technology for transport, development and networks for making available the development source of the FE software CESAR-LCPC, as well as the Pays de la Loire Region for financial support and all the partners of the French research Project SOLCYP (behavior of piles submitted to cyclic loading), granted by the French National Agency for Research (ANR), in which the present study has been developed.

## REFERENCES

1. Andersen KH. Bearing capacity under cyclic loading offshore, along the coast, and on land. *Canadian Geotechnical Journal* 2009; **46**:513-535.

2. Wichtmann T, Niemunis A, Triantafyllidis T. Strain accumulation in sand due to cyclic loading: drained triaxial tests. *Soil Dynamics and Earthquake Engineering* 2005; **25**:967-979.

3. Hicher PY. Comportement mécanique des argiles saturées sur divers chemins de sollicitations monotones et cycliques. Application à une modélisation élastoplastique et viscoplastique. *Thèse de Doctorat d'Etat*, Université Paris 6, 1985.

4. Niemunis A, Wichtmann T, Triantafyllidis T. A high-cycle accumulation model for sand. *Computers and Geotechnics* 2005; **32**:245-263.

5. Van Paepegem W, Degrieck J, De Baets P. Finite element approach for modelling fatigue damage in fibre-reinforced composite materials. *Composites* 2001; **Part B 32**:575-588.

6. Kiewel H, Aktaa J, Munz D. Application of an extrapolation method in thermocyclic failure analysis. *Computer Methods in Applied Mechanics and Engineering* 2000; **182**:55-71.

7. Cojocaru D, Karlsson AM. A simple numerical method of cycle jumps for cyclically loaded structures. *International Journal of Fatigue* 2006; **28**:1677-1689.

8. Lions JL, Maday Y, Turinici G. Résolution d'edp par un schéma en temps « pararéel ». *Comptes Rendus de l'Académie des Sciences*, Paris 2001; **332**(1):661-668.

9. Farhat C, Chandesris M. Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid-structure applications. *International Journal for Numerical Methods in Engineering* 2003; **58**:1397-1434.

10. Ponter ARS. Deformation, displacement, and work bounds for structures in a state creep and subject to variable loading. *Journal of Applied Mechanics* 1972; **39**:953-958.

11. Ponter ARS. The analysis of cyclically loaded creeping structures for short cycle times. *International Journal of Solids and Structures* 1976; **12**:809-825.

12. Ponter ARS, Brown PR. The finite element solution of rapid cycling creep problems. *International Journal for Numerical Methods in Engineering* 1978; **12**:1001-1024.

13. Spiliopoulos KV. Simplified methods for the steady state inelastic analysis of cyclically loaded structures. In Inelastic Analysis of Structures under Variable Loads. Theory and Engineering Applications, Weichert D, Maier G (eds). Kluwer Academic Publishers: Dordrecht, Boston, London, 2000; 213-232.

14. Bensoussan A, Lions JL, Papanicolaou G. Asymptotic Analysis for Periodic Structures. Studies in Mathematics and its Application. North-Holland: Amsterdam New York Oxford, 1978.

15. Sanchez-Palencia E. Non-Homogeneous Media and Vibration Theory. Lecture Notes in Physics, 127. Springer-Verlag: Berlin Heidelberg New York, 1980.

16. Kodaissi E, Aubry D, Meimon Y. Homogenization of cyclic loading on viscoplastic offshore foundations. In International Conference on Numerical Methods in Engineering: Theory and Applications, Swansea, Middleton J, Pande GN (eds). Balkema: Rotterdam Boston, 1985; 465-473.

17. Guennouni T. Sur une méthode de calcul de structures soumises à des chargements cycliques : l'homogénéisation en temps. *Modélisation mathématique et Analyse numérique* 1988; **22**:417-455.

18. Shahroul I, Meimon Y. Calculation of marine foundations subjected to repeated loads by means of homogenization method. *Computers and Geotechnics* 1995; **17**:93-106.

19. Boutin C, Wong H. Study of thermosensitive heterogeneous media via space-time homogenisation. *European Journal of Mechanics A/Solids* 1998; **17**:939-968.

20. Pecker A, Prevost JH, Dormieux L. Analysis of pore pressure generation and dissipation in cohesionless materials during seismic loading. *Journal of Earthquake Engineering* 2001; **5**(4):441-464.

21. Yu Q, Fish J. Temporal homogenization of viscoelastic and viscoplastic solids subjected to locally periodic loading. *Computational Mechanics* 2002; **29**:199-211.

22. Zienkiewicz OC, Taylor RL. The Finite Element Method, Vol. 2 (5th edn). Butterworth-Heinemann: Oxford, 2000.

23. Muir Wood D. Approaches to modelling the cyclic stress-strain response of soils. In Cyclic Loading of Soils: From Theory to Design, O'Reilly MP, Brown SF (eds). Blackie and Son: Glasgow and London, 1991; 18-69.
24. Grammatikopoulou A, Zdravkovic L, Potts DM. General formulation of two kinematic hardening constitutive models with a smooth elastoplastic transition. *International Journal of Geomechanics* 2006; **6**:291–302.

25. Sivasithamparam N, Kamrat-Pietraszewska D, Karstunen M. An anisotropic bubble model for soft clays. In Seventh European Conference on Numerical Methods in Geotechnical Engineering, Trondheim, Norwegian, Benz T, Nordal S (eds). CRC Press/Balkema: Rotterdam, 2010; 21–26.

26. Perzyna P. Fundamental problems in viscoplasticity. *Advances in Applied Mechanics* 1966; **9**:243–377.

27. Yin ZY, Hicher PY. Identifying parameters controlling soil delayed behaviour from laboratory and *in situ* pressuremeter testing. *International Journal for Numerical and Analytical Methods in Geomechanics* 2008; **32**(12):1515–1535.

28. Dafalias YF, Herrmann LR. Bounding surface formulation of soil plasticity. In Soil Mechanics – Transient and Cyclic Loads, Pande GN, Zienkiewicz OC (eds). Wiley and Sons: New York, 1982; 253–282.

29. Manzari MT, Nour MA. On implicit integration of bounding surface plasticity models. *Computers and Structures* 1997; **63**(3):385–395.

30. Muir Wood D. Soil Behavior and Critical State Soil Mechanics. Cambridge University Press: Cambridge New York Melbourne, 1990.

31. Dafalias YF, Manzari MT, Papadimitriou AG. SANICLAY: simple anisotropic clay plasticity model. *International Journal for Numerical and Analytical Methods in Geomechanics* 2006; **30**(12):1231–1257.

32. Yin ZY, Karstunen M, Chang CS, Koskinen M, Lojander M. Modeling time-dependent behavior of soft sensitive clay. *ASCE Journal of Geotechnical and Geoenvironmental Engineering* 2011; **137**(11):1103–1113.

33. Al-Tabbaa A. Permeability and stress–strain response of speswhite kaolin. *Ph.D. Thesis*, Cambridge University, 1987.

34. Hashiguchi K. Two- and three-surface models of plasticity. In Fifth International Conference on Numerical Methods in Geomechanics, Nagoya, Kawamoto T, Ichikawa Y (eds). Balkema: Rotterdam Boston, 1985; 285–292.

35. Hicher PY. Contribution à l'étude de la fatigue des argiles. *Ph.D. Thesis*, Ecole Centrale des Arts et Manufactures, Paris, 1979.

36. Papon A. Modélisation numérique du comportement des sols. Homogénéisation temporelle et identification des paramètres sous très grands nombres de cycles. *Ph.D. Thesis*, Ecole Centrale de Nantes, 2010.

37. Rangeard D, Hicher PY, Zentar R. Determining soil permeability from pressuremeter tests. *International Journal for Numerical and Analytical Methods in Geomechanics* 2003; **27**:1–24.