# Long-Term Strength of Porous Geomaterials by a Micromechanical Model considering Alternate Wetting and Drying Condition

Authors: Jin Zhang; Shirong Fu; Tao Ni; Bei Han; Chong Shi
Journal: Geofluids
 (2021)

---

## Abstract

This study is devoted to determining the long-term strength of porous geomaterials under alternate wetting and drying condition by statical shakedown analysis. In the framework of micromechanics of porous materials, Gurson’s hollow sphere model with Drucker-Prager solid matrix is adopted as the representative volume element. The effects of alternate wetting and drying are considered as variable water pressure imposed on the inner boundary surface of the unit cell. The cyclic responses are separated as a pure hydrostatic part under compressive/tensive loads and an additional deviatoric part to capture shear effects. The reduction of the long-term strength due to inner water pressure is observed by the illustration of obtained macroscopic criteria with respect to various load parameters. In addition, the accuracy of the analytical solution is also verified by comparing to the results of FEM-based step-by-step computations.


---

## Body

## 1. Introduction

Variable loadings exist widely in engineering structures, such as slops, offshore platform foundations, and pavements [1–3]. The long-term strength of the geomaterials in this condition is visibly reduced comparing to that subjected to a constant load. Actually, experimental observations [4–6] and numerical simulations [7, 8] have already shown that the fracture strains under cyclic loads are significantly lower than those reached in a monotonic way.Additionally, considering the soil or rocks above the water level in the underground engineering, the alternate wetting and drying condition due to the variation of the water level is even more harmful for its long-term stability, which is a common case in natural condition [9]. Consequently, the strength of geomaterials obtained in proportional or constant loading case is accordingly reduced due to the variable loads and alternating wetting condition. To this end, shakedown analysis, a powerful tool to provide the essential information of material under cyclic loads at limit state, is accepted in this research to compute the strength reduction.Two dual theoretical shakedown approaches have been contributed to this subject. The statical one was firstly developed by Melan [10]. He states that the structure shakes down if a so-called time-independent residual stress field can be found, such that the yield condition is always fulfilled in the whole body. The statical approach has been widely used and developed in engineering [11–13] theoretically and numerically, since the obtained limit load is strictly lower than the real one. On the other hand, the kinematical approach is based on a key concept of admissible plastic strain increment, which was firstly introduced by Koiter [14] and developed in [15–17]. Particularly, König and Siemaszko’s extension [18] of Melan’s theorem allows to verify the shakedown condition through critical cyclic loads containing all the vertices in the load domain, instead of any arbitrary load path.Considering the fact that most geomaterials contain microstructures [19, 20] (pores, mineral inclusions, etc.) which induce the inelastic macroscopic behaviors, numerous contributions have been devoted to study the effects of the microstructure on geomaterials within the framework of micromechanics [21–25], especially that of pores [26, 27]. Using Gurson’s hollow sphere model [28], the stress tensors are decomposed as the sum of a pure hydrostatic part under compressive/tensive loads and an additional deviatoric part to capture shear effects [29–31] owing to the symmetric geometry. The first author has already applied the statical shakedown theorem in determining the limit load of porous materials [32–34] with this model subjected to 1 independent cyclic load.In this work, the loading condition is even more complex than in the mentioned works. The hydrostatic and deviatoric parts of the load are no longer related by the imposed macroscopic stress triaxiality. In other words, two independent loads are considered in the current study, for which there are only numerical solutions [35] in the literature. Moreover, the alternate wetting and drying condition is simulated by the variable water pressure applied on the inner boundary surface of the hollow sphere model. The solid matrix is considered obeying Drucker-Prager’s yield law, in order to describe the plastic compressibility and asymmetric behaviors between tension and compression of geomaterials at the macroscale. The obtained macroscopic criterion by homogenization is expected able to predict the long-term strength of porous geomaterials under alternate wetting and drying condition.The paper is organized as follows. In Section2, micromechanic-based formulations and the statical shakedown theorem are briefly recalled. A limit analysis-based macroscopic strength criterion of porous geomaterials proposed in [29] is also provided. The strength reduction due to variable loadings considering alternate wetting and drying condition indicated by the macroscopic shakedown criterion is given in Section 3. The first and second subsections are devoted to the microscopic stress fields and application of the statical theorem, and the variable water pressure is taken into consideration in the third subsection. In Section 4, step-by-step elastoplastic numerical computations are performed for various load cases, and the results are compared with the obtained analytical solution. Conclusions and perspectives are given in the last section.


## 2. Micromechanical Model and Homogenized Macroscopic Criterion for Porous Geomaterials

In this section, a micromechanic-based hollow sphere representative volume element (RVE) of porous geomaterials is firstly introduced. Using this model, we recall also a macroscopic strength criterion proposed by Guo et al. [29], which can describe the plastic compressibility and asymmetric behaviors between tension and compression of the studied materials at the macroscale.


### 2.1. Homogenization of Porous Geomaterials and Hollow Sphere Model

We consider a class of porous geomaterials characterized by a solid phase and pores at the microscale, in which the pores are embedded (see Figure1). The classical Gurson’s hollow sphere model [28] with uniform strain rate boundary condition is adopted to study the macroscopic behaviors of this kind of material by homogenization.Figure 1
Porous geomaterials and hollow sphere model with uniform strain rate boundary condition.The internal and external radii of the chosen hollow sphere RVE are, respectively, noted asa and b. We denote Ω the total volume of the RVE, whereas ω the volume of the voids. Thus, the volume fraction of the void is given by f=ω/Ω=a/b3. The uniform strain rate boundary condition v=D·x is imposed on the outer surface of the hollow sphere, where D is a uniform macroscopic strain rate and x the position vector. The solid matrix is elastoplastic and considered obeying the associated Drucker-Prager plastic law in this work (without considering hardening effects) so as to exhibit the tension-compression asymmetry of geomaterials:
(1)fσ=σeqσ+3ασm−σy≤0,where σeq=3/2σ′:σ′ is the equivalent stress defined from the deviatoric part σ′ of the stress tensor σ and σm=1/3trσ is the mean stress. σy represents the yield stress of the matrix material and α the pressure sensitivity factor related to the friction angle ϕ:
(2)α=tanϕ3.Having the microscopic fields of stressσ and plastic strain rate d, the corresponding macroscopic ones Σ,D of the RVE can be calculated as averages [36]:
(3)Σ=1Ω∫ΩσdV,D=1Ω∫ΩddV.Moreover, the macroscopic stress fieldΣ can be simply obtained by integration on a surface instead of a volume:
(4)Σ=1Ω∫∂Ωσ·n⊗xdS,if the microscopic one σ is statically admissible
(5)σs.t.divσ=0inΩ,σ=0inω,σ·n=0on∂ω,in which n is the unit normal vector outward to the unit cell.According to Hill’s lemma, the following equality relating the dissipation at micro- and macroscales is always fulfilled for admissible fields considering the uniform strain rate boundary condition:(6)1Ω∫∂Ωd:σdV=D:Σ.


### 2.2. Macroscopic Strength Criterion

In the framework of limit analysis, the macroscopic strength criterion of porous geomaterials can be derived from the microscopic plastic yielding criterion of the matrix:(7)Fσ≤0⟶FΣ,f≤0,describing the nonlinear plastic behavior of the homogenized porous geomaterials. Besides, the effects of the void volume fraction can also be reflected. The accuracy of the homogenized criterion depends mainly on the applied stress and strain rate fields of the hollow sphere at microscale. In this subsection, we introduce a homogenized continuum model built around a three-parameter axisymmetric velocity field [29] by Guo et al.In a cylindrical coordinate system in orthonormal frameeρ,eϕ,ez, the following three-parameter velocity field is adopted:
(8)v=A0br3/sρeρ+zez+A1ρeρ+A2zez,where s=1±α for A0≷0. A0, A1, and A2 are parameters.The void growth can be computed as(9)ḟ=3fγ−fA0,in which γ=1−1/s, and the sign of A0 represents void growth (A0>0) or compression (A0<0).By means of limit analysis techniques, an implicit macroscopic criterion of porous geomaterials composed of Drucker-Prager solid matrix is proposed as follows:(10)Σm=131−ffγ−f+α∂Π/∂A0∂Π∂A0,Σeq=fγ−ffγ−f+α∂Π/∂A0∂Π∂De,where ∂Π/∂A0 and ∂Π/∂De can be, respectively, expressed by the Gauss hypergeometric function
(11)∂Π∂A0=ωs/2−1∑0∞1/2n1−s/2n2−s/2n−ω2nn!−f1−s/2∑0∞1/2n1−s/2n2−s/2n−ϖ2nn!,∂Π∂De=∑0∞1/2n−s/2n1−s/2n−ω2nn!−f∑0∞1/2n−s/2n1−s/2n−ϖ2nn!,with an the Pochhammer symbols an=Γa+n/Γa and the pair of ratios ϖ=ω/f1/s, −∞<ω<∞.Noticing that the invariants of the macroscopic stress tensor can be derived similar to their microscopic counterparts,(12)Σm=13TrΣ,Σeq=23Σ′:Σ′,Σ′ being the deviatoric part of Σ.Guo’s homogenized yield criteriaFΣmω,Σeqω=0 mainly defined by Equation (10) for porous geomaterials are plotted in Figure 2 with respect to different values of porosity (red line: f=0.05, green line: f=0.10, and blue line: f=0.20). The classical Drucker-Prager yield criterion (black dash dot line) is also presented in the same figure, in which the effect of porosity of geomaterials is not taken into consideration (f=0). It is interestingly observed that the void volume fraction has a great influence on the closed-form yield surface, which is provided by the macroscopic criterion Equation (10).Figure 2
Macroscopic yield surfaces defined by Guo’s criterion Equation (10) with respect to different values of porosity considering ϕ=15.Different from the phenomenological yield law (e.g., Drucker-Prager model), the dilatant and contractant behaviors of geomaterials in the triaxial compression test can be both reflected by this micromechanic-based model. For instance, within the region of a low confining pressure (−0.5<Σm/σy<0), the plastic dilatancy of geomaterials is easily observed while the contractancy can be obtained under a high confining pressure (Σm/σy<−1.5). Even for the same loading condition, the dilatant and contractant phenomena are different due to the influence of porosity. Under the confining pressure Σm/σy=−1.5, it is located in a dilatant domain for f=0.05 and in a contractant domain for f=0.20, instead. Besides, the compression/tension asymmetric behaviors of geomaterials and the strength reduction at the macroscale due to pores are also observed, according to Guo’s criterion. For more detailed derivation of the above macroscopic criterion, it is referred to [29].


## 2.1. Homogenization of Porous Geomaterials and Hollow Sphere Model

We consider a class of porous geomaterials characterized by a solid phase and pores at the microscale, in which the pores are embedded (see Figure1). The classical Gurson’s hollow sphere model [28] with uniform strain rate boundary condition is adopted to study the macroscopic behaviors of this kind of material by homogenization.Figure 1
Porous geomaterials and hollow sphere model with uniform strain rate boundary condition.The internal and external radii of the chosen hollow sphere RVE are, respectively, noted asa and b. We denote Ω the total volume of the RVE, whereas ω the volume of the voids. Thus, the volume fraction of the void is given by f=ω/Ω=a/b3. The uniform strain rate boundary condition v=D·x is imposed on the outer surface of the hollow sphere, where D is a uniform macroscopic strain rate and x the position vector. The solid matrix is elastoplastic and considered obeying the associated Drucker-Prager plastic law in this work (without considering hardening effects) so as to exhibit the tension-compression asymmetry of geomaterials:
(1)fσ=σeqσ+3ασm−σy≤0,where σeq=3/2σ′:σ′ is the equivalent stress defined from the deviatoric part σ′ of the stress tensor σ and σm=1/3trσ is the mean stress. σy represents the yield stress of the matrix material and α the pressure sensitivity factor related to the friction angle ϕ:
(2)α=tanϕ3.Having the microscopic fields of stressσ and plastic strain rate d, the corresponding macroscopic ones Σ,D of the RVE can be calculated as averages [36]:
(3)Σ=1Ω∫ΩσdV,D=1Ω∫ΩddV.Moreover, the macroscopic stress fieldΣ can be simply obtained by integration on a surface instead of a volume:
(4)Σ=1Ω∫∂Ωσ·n⊗xdS,if the microscopic one σ is statically admissible
(5)σs.t.divσ=0inΩ,σ=0inω,σ·n=0on∂ω,in which n is the unit normal vector outward to the unit cell.According to Hill’s lemma, the following equality relating the dissipation at micro- and macroscales is always fulfilled for admissible fields considering the uniform strain rate boundary condition:(6)1Ω∫∂Ωd:σdV=D:Σ.


## 2.2. Macroscopic Strength Criterion

In the framework of limit analysis, the macroscopic strength criterion of porous geomaterials can be derived from the microscopic plastic yielding criterion of the matrix:(7)Fσ≤0⟶FΣ,f≤0,describing the nonlinear plastic behavior of the homogenized porous geomaterials. Besides, the effects of the void volume fraction can also be reflected. The accuracy of the homogenized criterion depends mainly on the applied stress and strain rate fields of the hollow sphere at microscale. In this subsection, we introduce a homogenized continuum model built around a three-parameter axisymmetric velocity field [29] by Guo et al.In a cylindrical coordinate system in orthonormal frameeρ,eϕ,ez, the following three-parameter velocity field is adopted:
(8)v=A0br3/sρeρ+zez+A1ρeρ+A2zez,where s=1±α for A0≷0. A0, A1, and A2 are parameters.The void growth can be computed as(9)ḟ=3fγ−fA0,in which γ=1−1/s, and the sign of A0 represents void growth (A0>0) or compression (A0<0).By means of limit analysis techniques, an implicit macroscopic criterion of porous geomaterials composed of Drucker-Prager solid matrix is proposed as follows:(10)Σm=131−ffγ−f+α∂Π/∂A0∂Π∂A0,Σeq=fγ−ffγ−f+α∂Π/∂A0∂Π∂De,where ∂Π/∂A0 and ∂Π/∂De can be, respectively, expressed by the Gauss hypergeometric function
(11)∂Π∂A0=ωs/2−1∑0∞1/2n1−s/2n2−s/2n−ω2nn!−f1−s/2∑0∞1/2n1−s/2n2−s/2n−ϖ2nn!,∂Π∂De=∑0∞1/2n−s/2n1−s/2n−ω2nn!−f∑0∞1/2n−s/2n1−s/2n−ϖ2nn!,with an the Pochhammer symbols an=Γa+n/Γa and the pair of ratios ϖ=ω/f1/s, −∞<ω<∞.Noticing that the invariants of the macroscopic stress tensor can be derived similar to their microscopic counterparts,(12)Σm=13TrΣ,Σeq=23Σ′:Σ′,Σ′ being the deviatoric part of Σ.Guo’s homogenized yield criteriaFΣmω,Σeqω=0 mainly defined by Equation (10) for porous geomaterials are plotted in Figure 2 with respect to different values of porosity (red line: f=0.05, green line: f=0.10, and blue line: f=0.20). The classical Drucker-Prager yield criterion (black dash dot line) is also presented in the same figure, in which the effect of porosity of geomaterials is not taken into consideration (f=0). It is interestingly observed that the void volume fraction has a great influence on the closed-form yield surface, which is provided by the macroscopic criterion Equation (10).Figure 2
Macroscopic yield surfaces defined by Guo’s criterion Equation (10) with respect to different values of porosity considering ϕ=15.Different from the phenomenological yield law (e.g., Drucker-Prager model), the dilatant and contractant behaviors of geomaterials in the triaxial compression test can be both reflected by this micromechanic-based model. For instance, within the region of a low confining pressure (−0.5<Σm/σy<0), the plastic dilatancy of geomaterials is easily observed while the contractancy can be obtained under a high confining pressure (Σm/σy<−1.5). Even for the same loading condition, the dilatant and contractant phenomena are different due to the influence of porosity. Under the confining pressure Σm/σy=−1.5, it is located in a dilatant domain for f=0.05 and in a contractant domain for f=0.20, instead. Besides, the compression/tension asymmetric behaviors of geomaterials and the strength reduction at the macroscale due to pores are also observed, according to Guo’s criterion. For more detailed derivation of the above macroscopic criterion, it is referred to [29].


## 3. Strength Reduction due to Variable Loadings considering Alternate Wetting and Drying Condition

This section is devoted to calculate the strength reduction of porous geomaterials subjected to variable loadings, comparing to the homogenized yield criterion (Equation (10)). As mentioned in the first section, the effective load-bearing capacity of materials is also depending on the varying amplitude of applied loads. Similar to pioneering theoretical works [32, 34, 37] on porous media, the stress and strain fields of the hollow sphere RVE are divided into hydrostatic and deviatoric parts, on account of the axisymmetric geometry.It must be remarked that, unlike the previous researches where the variable loadings are imposed proportional during each cycle (T=Σm/Σeq is constant), the hydrostatic parts Σm and Σeq are considered independent in the paper for the purpose of no loss of generality in engineering problems. In other words, the long-term stability of porous geomaterials is concerning a multidimensional shakedown problem instead of one. If the load varies in a certain domain, for which after a transient phase, the material’s response becomes linearly elastic (shakedown). This will ensure the long-term safety of porous geomaterials.


### 3.1. Transformed Variable Loading Problem and Statical Shakedown Theorem

Considering an applied variable loadingΣt in the domain P (see Figure 3(a)), the residual stresses ρ in the considered model can be obtained by subtracting the elastic responses σE from the total stresses σ:
(13)ρx,t=σx,t−σEx,t,where σEx,t is the purely elastic stress responses in the fictitious elastic cell.Figure 3
Transformation of the variable load within a given load domain.
(a)
Original variable load(b)
Equivalent cycles of loadAccording to Melan’s classical shakedown theorem [10], if there exists a time-independent residual stress field ρ¯x, such that
(14)FσEx,t+ρ¯x<0is fulfilled anywhere in the body at any time, the material shakes down. F represents the yield function of the solid matrix (Equation (1)). Then, the dissipated energy is bounded in time, independently on the initial state.Moreover, the key-point time-independent residual stress field belongs to the set of statically admissible fields (Equation (5)) in the present study:
(15)ρ¯∣divρ¯=0inΩ,ρ¯·n=0on∂ω,ρ¯=0inω.Owing to Equation (4), the following special relation relating the residual stresses at macro- and microscales is imposed by the surface integration:
(16)Σr=1Ω∫∂Ωρ¯·n⊗xdS=0,such that the average residual stresses vanish and are identical to the macroscopic elastic stresses in the fictitious elastic body:
(17)Σ=Σr+ΣE=ΣE.Let us introduce the load factorα, which controls the size of the actual load domain P=αP0. Beyond the threshold of the maximum of the admissible load αSD, the body can not keep in the shakedown state under variable loads. As a result, to find the shakedown limit load factor αSD, beyond which the fracture due to fatigue or formation of mechanism (incremental collapse at first cycle), is the main objective of this paper for the long-term stability of porous geomaterials.In practice, the shakedown state under original arbitrary load guaranteed by Melan’s theorem can be transformed owing to König and Siemaszko [18]: if a given structure shakes down over any load path contained within a given load domain P, then it also shakes down over any load path contained within the convex hull of P. In other words, the shakedown condition only needs to be verified over a cyclic load path containing all the vertices of the hyperpolyhedral load domain P (e.g., in Figure 3(b), path 1 ⟶ 2 ⟶ 3 ⟶ 4 ⟶ 1), instead of any arbitrary load path within the same domain. This is also called the convex hull theorem.Considering the notion of the load factorα, the actual load can be decomposed into a combination of several elementary ones
(18)Σt=∑k=1nμktΣk=α∑k=1nμktΣk0,where α is the same factor defined in P=αP0 and Σk0 the reference loads. The load multipliers μkt are independent within the range μk−<μkt<μk+. In this research, the number of independent loads n=2. Consequently, the two elementary loads will construct a 2-dimensional convex load domain having 4 vertices.


### 3.2. Macroscopic Shakedown Criterion of Porous Geomaterials under Variable Loadings

Recalling the previous researches on the shakedown analysis of porous media under one cyclic load, the trial stress fields at the microscale are separated into two parts:(19)σ=σ1+σ2,in order to describe the compressive/tensive and shear effects of the hollow sphere model subjected to pure hydrostatic and deviatoric loads, respectively. It must be pointed out that the construction of the trial stress fields is already provided in the former research [34]. To avoid repetition, we only give the expression of these fields in the following parts.
(1)
As mentioned in Equation (13), the exact stress field subjected to hydrostatic load is considered as the sum of two parts:(20)σ1=σE1+ρ¯1.The elastic stress fieldσE1 in the fictitious body is given in spherical coordinates:
(21)σE1=Σm1−f1+12ar3eθ⊗eθ+eϕ⊗eϕ−2er⊗er.The corresponding residual stress fieldρ¯1 writes as follows:
(22)ρ¯1=C01−ar3γ1+32γar3γeθ⊗eθ+eϕ⊗eϕ−Σm+1−f1+12ar3eθ⊗eθ+eϕ⊗eϕ−2er⊗er.(2)
Similarly, the additional terms of the stress field under deviatoric load is(23)σ2=σE2+ρ¯2,where the fictitious elastic stress field writes as
(24)σE2=−signJ3Σeq31−fa318a2+5r2−5+ν1+3cos2θ2r5−7+5ν−1+3cos2θ2er⊗er+a327a2+5r21−2ν−321a2+5r2−1+2νcos2θ2r5−7+5ν+−1+3cos2θ2eθ⊗eθ+a39a2+25r2−1+2ν−45a2+r2−1+2νcos2θ2r5−7+5ν+1eϕ⊗eϕ+3a312a2−5r21+νsin2θ2r5−7+5ν+3sin2θ2er⊗eθ+eθ⊗er,with ν being Poisson’s coefficient and J3 the third invariant of the macroscopic stress deviator.And the corresponding residualρ¯2 takes the form ρ¯2=ρ¯2a+ρ¯2b, where ρ¯2a is the maximum value of σE2 and ρ¯2b is related to σE2 by
(25)σ′ρ¯2b=C1Krσ′σE2,σmρ¯2b=C1KrσmσE2,where C1 is the constant to be determined and Kr is a function of r, noticing that the detailed derivation of ρ¯2 is firstly provided in [33].Consequently, the full expression of the complete microscopic stress fields (19) has been obtained. Inspiring from former researches [32, 34], the most “dangerous” point of the considered hollow sphere model always locates at the inner boundary (r=a). In this case, the average stress σm and equivalent stress σeq can be readily computed from Equation (19) to Equation (25):
(26)σm=γC0−Σm+−Σm1−f−5ν+13cos2θ−1−7+5νsignJ3+Σeq+−signJ3Σeq31−f+KaC1,(27)σeq2=3γC02−32Σm+−Σm1−f2+ξ2θsignJ3+Σeq+−signJ3Σeq1−f+3KaC12+ξ1θ3γC0−3Σm+−Σm1−fsignJ3+Σeq+−signJ3Σeq1−f+3KaC1,where J3+ (resp., J3−) is the maximum value (resp., minimum value) of the third invariant in the deviatoric space, and
(28)ξ1θ=35ν+57−5ν3cos2θ−12,ξ2θ=225ν2−25ν−1780cos4θ+−375ν2+600ν−375cos2θ+175ν2−325ν−425295ν−72.Owing to König and Siemaszko’s statical theorem [18], only the critical cyclic loading paths containing all vertices of the load domain (Figure 3(b)) need to be verified to guarantee the material’s safety. In this study, the applied loads are considered varying from zero to maximum value:
(29)0<Σm<Σm+,0<Σeq<Σeq+,which is the most common case in engineering. As shown in Figure 4, there are load paths 1, 2, and 3 including 4 vertices: P^10,0, P^20,Σeq+, P^3Σm+,Σeq+, and P^4Σm+,0.Figure 4
Critical cyclic loadings of the considered load domain.Combining the shakedown condition (14) and the microscopic stress field (20), the following equations at the vertices can be obtained at the shakedown limit state (the equality is reached):
(30)94γC02+KaC12ξ2θ+3γC0ξ1θKaC11/2+3αγC0−5ν+13cos2θ−1−7+5νKaC1=σy,(31)94γC0−Σm+1−f2+signJ3+Σeq+31−f+KaC12ξ2θ+3γC0−Σm+1−fξ1θsignJ3+Σeq+31−f+KaC11/2+3αγC0−Σm+1−f−5ν+13cos2θ−1−7+5νsignJ3+Σeq+31−f+KaC1=σy,(32)94γC0−Σm+1−f2+KaC12ξ2θ+3γC0−Σm+1−fξ1θKaC11/2+3αγC0−Σm+1−f−5ν+13cos2θ−1−7+5νKaC1=σy,(33)94γC02+signJ3+Σeq+31−f+KaC12ξ2θ+3γC0ξ1θsignJ3+Σeq+31−f+KaC11/2+3αγC0−5ν+13cos2θ−1−7+5νsignJ3+Σeq+31−f+KaC1=σy.Let us introduce the general stress triaxiality considering the pure elastic responses within the shakedown load domain, which is defined by(34)τ=−Σm+/1−fsignJ3+Σeq+/31−f=γC0KaC1=γC0−Σm+/1−fsignJ3+Σeq+/31−f+KaC1,for the load path connecting P^1 and P^3. Likewise, for vertices P^2 and P^4, one has the following relation:
(35)−τ=γC0−Σm+/1−fKaC1=γC0signJ3+Σeq+/31−f+KaC1.Referring to [34], KaC1 can be eliminated by taking Equation (34) to the first two shakedown conditions (30) and (31). Solving the two equations with respect to the sign of KaC1, one can obtain
(36)signJ3+Σeq+31−fσy=29/4τ2+3τξ1θ+ξ2θ9/4τ2+3τξ1θ+ξ2θ−9α2τ−5ν+13cos2θ−1/−7+5ν2.Taking the general stress triaxiality (34) into account, the macroscopic shakedown criterion of porous geomaterials derived from Equations (30) and (31) is
(37)94Σm+1−f2−3ξ1θΣm+1−fsignJ3+Σeq+31−f+ξ2θsignJ3+Σeq+31−f2−9α2Σm+1−f+5ν+13cos2θ−1−7+5νsignJ3+Σeq+31−f2=2σy94Σm+1−f2−3ξ1θΣm+1−fsignJ3+Σeq+31−f+ξ2θsignJ3+Σeq+31−f2.Similarly, resulting from the last two shakedown conditions (32), Equation (33) and the general stress triaxiality (35), the other part to complete the full macroscopic shakedown criterion is
(38)94Σm+1−f2+3ξ1θΣm+1−fsignJ3+Σeq+31−f+ξ2θsignJ3+Σeq+31−f2−9α2−Σm+1−f+5ν+13cos2θ−1−7+5νsignJ3+Σeq+31−f2=2σy94Σm+1−f2+3ξ1θΣm+1−fsignJ3+Σeq+31−f+ξ2θsignJ3+Σeq+31−f2.According to previous researches, the most “dangerous point” of the adopted hollow sphere model locates at the equator (θ=π/2) or the poles (r=0 or π) depending on the sign of J3+, where
(39)ξ1π2=35ν+525ν−7,ξ2π2=2257ν2−13ν+75ν−72,ξ10=35ν+5−5ν+7,ξ20=225ν2+2ν+15ν−72.The established macroscopic shakedown criterion (Equations (37) and (38)) is illustrated in Figure 5 with porosity f=0.1 and Poisson’s coefficient ν=0.2. The friction angle ϕ of the solid matrix is considered as 30°. The limit surfaces defined by the proposed macroscopic criterion are plotted on a dash dot line (green line for the collapse at equator θ=π/2; blue line for the collapse at poles θ=0). As a result, the effective shakedown domain is bounded by the red solid line.Figure 5
Limit surfaces defined by the proposed macroscopic shakedown criterion of porous geomaterials takingf=0.1, ν=0.2, and ϕ=30.Comparing with the analytic research of porous media under one independent load, the obtained shakedown domain in Figure5 is completely symmetric about the axis Σeq/σy. It can be seen as the “pulsating loading case” in [34], but the hydrostatic part Σm is not related to the deviatoric part Σeq through the stress triaxiality. Moreover, it is obviously seen that the sign of J3+ has no influence on the final shakedown domain (red line). Indeed, this interesting conclusion is in accord with the numerical study [35] by nonlinear optimization based on the interior-point method, which is concerning the shakedown analysis of porous media having von Mises solid matrix.


### 3.3. Long-Term Stability under Alternate Wetting and Drying Condition

As mentioned in the first section, alternate wetting and drying condition occurs in some geoengineering and will reduce its long-term strength. Obviously, an additional stress field needs to be considered for shakedown analysis of porous geomaterials. In this subsection, the effects owing to the alternate wetting and drying condition on the materials are considered as a variable water pressurepw applied at the inner boundary (r=a) of the hollow sphere model.Similar to the construction of stress field (20) under pure hydrostatic load, an additional microscopic stress field is also proposed in the following form:
(40)σw=σEw+ρ¯w,in which the fictitious elastic response and corresponding residual stress field under variable water pressure are
(41)σEw=−pw1−f1+12ar3eθ⊗eθ+eϕ⊗eϕ−2er⊗er,ρ¯1=pw+1−f1+12ar3eθ⊗eθ+eϕ⊗eϕ−2er⊗er−Cw1−ar3γ1+32γar3γeθ⊗eθ+eϕ⊗eϕ,with pw+ being the maximum value of the water pressure and Cw an additional constant. In engineering, the water pressure usually varies in the range 0<pw<pw+.Consequently, resulting from Equations (19) and (40), the total stress at the microscale including the effect of water pressure reads
(42)σ=σ1+σw+σ2,Following the derivation of the shakedown criterion (from Equation (26) to Equation (38)) in the previous subsection, a new macroscopic criterion to predict the long-term stability of porous geomaterials considering alternate wetting and drying condition can be obtained as
(43)94Σm+−pw+1−f2−3ξ1θΣm+−pw+1−fΣeq+31−f+ξ2θΣeq+31−f2−9α2Σm+−pw+1−f+5ν+13cos2θ−1−7+5νΣeq+31−f2=2σy94Σm+−pw+1−f2−3ξ1θΣm+−pw+1−fΣeq+31−f+ξ2θΣeq+31−f2,94Σm+−pw+1−f2+3ξ1θΣm+−pw+1−fΣeq+31−f+ξ2θΣeq+31−f2−9α2−Σm+−pw+1−f+5ν+13cos2θ−1−7+5νΣeq+31−f2=2σy94Σm+−pw+1−f2+3ξ1θΣm+−pw+1−fΣeq+31−f+ξ2θΣeq+31−f2,in which this shakedown condition needs to be fulfilled at θ=0 and π/2 to avoid the fracture due to variable loads in the whole model. The corresponding expressions of ξ1θ and ξ2θ are given by Equation (39).Comparing with the macroscopic criterion in the previous subsection, it can be concluded from the proposed one (43) that the existence of the variable water pressure pw will decrease the effective mean stress and finally reduce the load-bearing capability of the considered material. The amplitude of the variation on the inner water pressure pw depends mostly on the water level in geomaterials. The detailed comparisons will be fully discussed in the next section.


## 3.1. Transformed Variable Loading Problem and Statical Shakedown Theorem

Considering an applied variable loadingΣt in the domain P (see Figure 3(a)), the residual stresses ρ in the considered model can be obtained by subtracting the elastic responses σE from the total stresses σ:
(13)ρx,t=σx,t−σEx,t,where σEx,t is the purely elastic stress responses in the fictitious elastic cell.Figure 3
Transformation of the variable load within a given load domain.
(a)
Original variable load(b)
Equivalent cycles of loadAccording to Melan’s classical shakedown theorem [10], if there exists a time-independent residual stress field ρ¯x, such that
(14)FσEx,t+ρ¯x<0is fulfilled anywhere in the body at any time, the material shakes down. F represents the yield function of the solid matrix (Equation (1)). Then, the dissipated energy is bounded in time, independently on the initial state.Moreover, the key-point time-independent residual stress field belongs to the set of statically admissible fields (Equation (5)) in the present study:
(15)ρ¯∣divρ¯=0inΩ,ρ¯·n=0on∂ω,ρ¯=0inω.Owing to Equation (4), the following special relation relating the residual stresses at macro- and microscales is imposed by the surface integration:
(16)Σr=1Ω∫∂Ωρ¯·n⊗xdS=0,such that the average residual stresses vanish and are identical to the macroscopic elastic stresses in the fictitious elastic body:
(17)Σ=Σr+ΣE=ΣE.Let us introduce the load factorα, which controls the size of the actual load domain P=αP0. Beyond the threshold of the maximum of the admissible load αSD, the body can not keep in the shakedown state under variable loads. As a result, to find the shakedown limit load factor αSD, beyond which the fracture due to fatigue or formation of mechanism (incremental collapse at first cycle), is the main objective of this paper for the long-term stability of porous geomaterials.In practice, the shakedown state under original arbitrary load guaranteed by Melan’s theorem can be transformed owing to König and Siemaszko [18]: if a given structure shakes down over any load path contained within a given load domain P, then it also shakes down over any load path contained within the convex hull of P. In other words, the shakedown condition only needs to be verified over a cyclic load path containing all the vertices of the hyperpolyhedral load domain P (e.g., in Figure 3(b), path 1 ⟶ 2 ⟶ 3 ⟶ 4 ⟶ 1), instead of any arbitrary load path within the same domain. This is also called the convex hull theorem.Considering the notion of the load factorα, the actual load can be decomposed into a combination of several elementary ones
(18)Σt=∑k=1nμktΣk=α∑k=1nμktΣk0,where α is the same factor defined in P=αP0 and Σk0 the reference loads. The load multipliers μkt are independent within the range μk−<μkt<μk+. In this research, the number of independent loads n=2. Consequently, the two elementary loads will construct a 2-dimensional convex load domain having 4 vertices.


## 3.2. Macroscopic Shakedown Criterion of Porous Geomaterials under Variable Loadings

Recalling the previous researches on the shakedown analysis of porous media under one cyclic load, the trial stress fields at the microscale are separated into two parts:(19)σ=σ1+σ2,in order to describe the compressive/tensive and shear effects of the hollow sphere model subjected to pure hydrostatic and deviatoric loads, respectively. It must be pointed out that the construction of the trial stress fields is already provided in the former research [34]. To avoid repetition, we only give the expression of these fields in the following parts.
(1)
As mentioned in Equation (13), the exact stress field subjected to hydrostatic load is considered as the sum of two parts:(20)σ1=σE1+ρ¯1.The elastic stress fieldσE1 in the fictitious body is given in spherical coordinates:
(21)σE1=Σm1−f1+12ar3eθ⊗eθ+eϕ⊗eϕ−2er⊗er.The corresponding residual stress fieldρ¯1 writes as follows:
(22)ρ¯1=C01−ar3γ1+32γar3γeθ⊗eθ+eϕ⊗eϕ−Σm+1−f1+12ar3eθ⊗eθ+eϕ⊗eϕ−2er⊗er.(2)
Similarly, the additional terms of the stress field under deviatoric load is(23)σ2=σE2+ρ¯2,where the fictitious elastic stress field writes as
(24)σE2=−signJ3Σeq31−fa318a2+5r2−5+ν1+3cos2θ2r5−7+5ν−1+3cos2θ2er⊗er+a327a2+5r21−2ν−321a2+5r2−1+2νcos2θ2r5−7+5ν+−1+3cos2θ2eθ⊗eθ+a39a2+25r2−1+2ν−45a2+r2−1+2νcos2θ2r5−7+5ν+1eϕ⊗eϕ+3a312a2−5r21+νsin2θ2r5−7+5ν+3sin2θ2er⊗eθ+eθ⊗er,with ν being Poisson’s coefficient and J3 the third invariant of the macroscopic stress deviator.And the corresponding residualρ¯2 takes the form ρ¯2=ρ¯2a+ρ¯2b, where ρ¯2a is the maximum value of σE2 and ρ¯2b is related to σE2 by
(25)σ′ρ¯2b=C1Krσ′σE2,σmρ¯2b=C1KrσmσE2,where C1 is the constant to be determined and Kr is a function of r, noticing that the detailed derivation of ρ¯2 is firstly provided in [33].Consequently, the full expression of the complete microscopic stress fields (19) has been obtained. Inspiring from former researches [32, 34], the most “dangerous” point of the considered hollow sphere model always locates at the inner boundary (r=a). In this case, the average stress σm and equivalent stress σeq can be readily computed from Equation (19) to Equation (25):
(26)σm=γC0−Σm+−Σm1−f−5ν+13cos2θ−1−7+5νsignJ3+Σeq+−signJ3Σeq31−f+KaC1,(27)σeq2=3γC02−32Σm+−Σm1−f2+ξ2θsignJ3+Σeq+−signJ3Σeq1−f+3KaC12+ξ1θ3γC0−3Σm+−Σm1−fsignJ3+Σeq+−signJ3Σeq1−f+3KaC1,where J3+ (resp., J3−) is the maximum value (resp., minimum value) of the third invariant in the deviatoric space, and
(28)ξ1θ=35ν+57−5ν3cos2θ−12,ξ2θ=225ν2−25ν−1780cos4θ+−375ν2+600ν−375cos2θ+175ν2−325ν−425295ν−72.Owing to König and Siemaszko’s statical theorem [18], only the critical cyclic loading paths containing all vertices of the load domain (Figure 3(b)) need to be verified to guarantee the material’s safety. In this study, the applied loads are considered varying from zero to maximum value:
(29)0<Σm<Σm+,0<Σeq<Σeq+,which is the most common case in engineering. As shown in Figure 4, there are load paths 1, 2, and 3 including 4 vertices: P^10,0, P^20,Σeq+, P^3Σm+,Σeq+, and P^4Σm+,0.Figure 4
Critical cyclic loadings of the considered load domain.Combining the shakedown condition (14) and the microscopic stress field (20), the following equations at the vertices can be obtained at the shakedown limit state (the equality is reached):
(30)94γC02+KaC12ξ2θ+3γC0ξ1θKaC11/2+3αγC0−5ν+13cos2θ−1−7+5νKaC1=σy,(31)94γC0−Σm+1−f2+signJ3+Σeq+31−f+KaC12ξ2θ+3γC0−Σm+1−fξ1θsignJ3+Σeq+31−f+KaC11/2+3αγC0−Σm+1−f−5ν+13cos2θ−1−7+5νsignJ3+Σeq+31−f+KaC1=σy,(32)94γC0−Σm+1−f2+KaC12ξ2θ+3γC0−Σm+1−fξ1θKaC11/2+3αγC0−Σm+1−f−5ν+13cos2θ−1−7+5νKaC1=σy,(33)94γC02+signJ3+Σeq+31−f+KaC12ξ2θ+3γC0ξ1θsignJ3+Σeq+31−f+KaC11/2+3αγC0−5ν+13cos2θ−1−7+5νsignJ3+Σeq+31−f+KaC1=σy.Let us introduce the general stress triaxiality considering the pure elastic responses within the shakedown load domain, which is defined by(34)τ=−Σm+/1−fsignJ3+Σeq+/31−f=γC0KaC1=γC0−Σm+/1−fsignJ3+Σeq+/31−f+KaC1,for the load path connecting P^1 and P^3. Likewise, for vertices P^2 and P^4, one has the following relation:
(35)−τ=γC0−Σm+/1−fKaC1=γC0signJ3+Σeq+/31−f+KaC1.Referring to [34], KaC1 can be eliminated by taking Equation (34) to the first two shakedown conditions (30) and (31). Solving the two equations with respect to the sign of KaC1, one can obtain
(36)signJ3+Σeq+31−fσy=29/4τ2+3τξ1θ+ξ2θ9/4τ2+3τξ1θ+ξ2θ−9α2τ−5ν+13cos2θ−1/−7+5ν2.Taking the general stress triaxiality (34) into account, the macroscopic shakedown criterion of porous geomaterials derived from Equations (30) and (31) is
(37)94Σm+1−f2−3ξ1θΣm+1−fsignJ3+Σeq+31−f+ξ2θsignJ3+Σeq+31−f2−9α2Σm+1−f+5ν+13cos2θ−1−7+5νsignJ3+Σeq+31−f2=2σy94Σm+1−f2−3ξ1θΣm+1−fsignJ3+Σeq+31−f+ξ2θsignJ3+Σeq+31−f2.Similarly, resulting from the last two shakedown conditions (32), Equation (33) and the general stress triaxiality (35), the other part to complete the full macroscopic shakedown criterion is
(38)94Σm+1−f2+3ξ1θΣm+1−fsignJ3+Σeq+31−f+ξ2θsignJ3+Σeq+31−f2−9α2−Σm+1−f+5ν+13cos2θ−1−7+5νsignJ3+Σeq+31−f2=2σy94Σm+1−f2+3ξ1θΣm+1−fsignJ3+Σeq+31−f+ξ2θsignJ3+Σeq+31−f2.According to previous researches, the most “dangerous point” of the adopted hollow sphere model locates at the equator (θ=π/2) or the poles (r=0 or π) depending on the sign of J3+, where
(39)ξ1π2=35ν+525ν−7,ξ2π2=2257ν2−13ν+75ν−72,ξ10=35ν+5−5ν+7,ξ20=225ν2+2ν+15ν−72.The established macroscopic shakedown criterion (Equations (37) and (38)) is illustrated in Figure 5 with porosity f=0.1 and Poisson’s coefficient ν=0.2. The friction angle ϕ of the solid matrix is considered as 30°. The limit surfaces defined by the proposed macroscopic criterion are plotted on a dash dot line (green line for the collapse at equator θ=π/2; blue line for the collapse at poles θ=0). As a result, the effective shakedown domain is bounded by the red solid line.Figure 5
Limit surfaces defined by the proposed macroscopic shakedown criterion of porous geomaterials takingf=0.1, ν=0.2, and ϕ=30.Comparing with the analytic research of porous media under one independent load, the obtained shakedown domain in Figure5 is completely symmetric about the axis Σeq/σy. It can be seen as the “pulsating loading case” in [34], but the hydrostatic part Σm is not related to the deviatoric part Σeq through the stress triaxiality. Moreover, it is obviously seen that the sign of J3+ has no influence on the final shakedown domain (red line). Indeed, this interesting conclusion is in accord with the numerical study [35] by nonlinear optimization based on the interior-point method, which is concerning the shakedown analysis of porous media having von Mises solid matrix.


## 3.3. Long-Term Stability under Alternate Wetting and Drying Condition

As mentioned in the first section, alternate wetting and drying condition occurs in some geoengineering and will reduce its long-term strength. Obviously, an additional stress field needs to be considered for shakedown analysis of porous geomaterials. In this subsection, the effects owing to the alternate wetting and drying condition on the materials are considered as a variable water pressurepw applied at the inner boundary (r=a) of the hollow sphere model.Similar to the construction of stress field (20) under pure hydrostatic load, an additional microscopic stress field is also proposed in the following form:
(40)σw=σEw+ρ¯w,in which the fictitious elastic response and corresponding residual stress field under variable water pressure are
(41)σEw=−pw1−f1+12ar3eθ⊗eθ+eϕ⊗eϕ−2er⊗er,ρ¯1=pw+1−f1+12ar3eθ⊗eθ+eϕ⊗eϕ−2er⊗er−Cw1−ar3γ1+32γar3γeθ⊗eθ+eϕ⊗eϕ,with pw+ being the maximum value of the water pressure and Cw an additional constant. In engineering, the water pressure usually varies in the range 0<pw<pw+.Consequently, resulting from Equations (19) and (40), the total stress at the microscale including the effect of water pressure reads
(42)σ=σ1+σw+σ2,Following the derivation of the shakedown criterion (from Equation (26) to Equation (38)) in the previous subsection, a new macroscopic criterion to predict the long-term stability of porous geomaterials considering alternate wetting and drying condition can be obtained as
(43)94Σm+−pw+1−f2−3ξ1θΣm+−pw+1−fΣeq+31−f+ξ2θΣeq+31−f2−9α2Σm+−pw+1−f+5ν+13cos2θ−1−7+5νΣeq+31−f2=2σy94Σm+−pw+1−f2−3ξ1θΣm+−pw+1−fΣeq+31−f+ξ2θΣeq+31−f2,94Σm+−pw+1−f2+3ξ1θΣm+−pw+1−fΣeq+31−f+ξ2θΣeq+31−f2−9α2−Σm+−pw+1−f+5ν+13cos2θ−1−7+5νΣeq+31−f2=2σy94Σm+−pw+1−f2+3ξ1θΣm+−pw+1−fΣeq+31−f+ξ2θΣeq+31−f2,in which this shakedown condition needs to be fulfilled at θ=0 and π/2 to avoid the fracture due to variable loads in the whole model. The corresponding expressions of ξ1θ and ξ2θ are given by Equation (39).Comparing with the macroscopic criterion in the previous subsection, it can be concluded from the proposed one (43) that the existence of the variable water pressure pw will decrease the effective mean stress and finally reduce the load-bearing capability of the considered material. The amplitude of the variation on the inner water pressure pw depends mostly on the water level in geomaterials. The detailed comparisons will be fully discussed in the next section.


## 4. Numerical Verification by Step-by-Step FEM Computations

This section is devoted to verify the accuracy of the proposed macroscopic criterion of porous geomaterials by comparison to numerical results. The elastoplastic numerical simulations are carried out by Abaqus software in a step-by-step manner to analyze the transient phase before shakedown or collapse.


### 4.1. Implementation of the FEM-Based Computation

Owing to the geometrical symmetry, only a quarter of the hollow sphere model discretized by 1500 axisymmetric elements (see Figure6) is considered, which is refined enough to capture the reliable results. The homogeneous boundary condition v=D·x is imposed on the outer boundary.Figure 6
Initial mesh of the hollow sphere model.As described in Figure4, only 3 critical cyclic load paths containing all the vertices need to be considered to calculate the shakedown limit load, if the variation of inner water pressure is not concerned. On the contrary, under alternate wetting and drying condition, the load domain is modified as shown in Figure 7 located in the first and second quadrants. Consequently, 4 critical cyclic load paths have to be taken into consideration to fulfill the shakedown condition.Figure 7
Critical cyclic loadings of the load domain considering the alternate wetting and drying condition.Nevertheless, the macroscopic stress triaxiality is defined byT3=Σm+/Σeq+ at vertex P^3 describing the shape of the rectangular load domain. Moreover, the stress triaxiality T=Σm/Σeq for all the cyclic load path remains constant in a single cycle, which is implemented in Abaqus by the means of multipoint constraint subroutine (MPC) satisfying the imposed uniform velocity field boundary condition. Notice that such MPC subroutine is firstly provided in the study of polymeric materials [38] and latterly applied in a series of ductile porous media [30, 32, 35] with the same boundary condition.The above computation procedure is performed for an imposed value of stress triaxialityT3 with 100 cycles. The long-term stability of the considered material is evaluated during the transient phase before shakedown or collapse. In practice, the accumulated equivalent plastic strain at the “dangerous” point will be checked in the transient phase, which is supposed to be constant if it is in a shakedown state. By increasing the applied load, the shakedown limit is obtained when the collapse due to fatigue or incremental fracture occurs.The cyclic responses on the inner boundaryr=a at equator θ=π/2 of the hollow sphere at a shakedown state with T3=0.44 are plotted in Figure 8. The porosity f=0.01 and ν=0.2, ϕ=20° of the solid phase are considered. It can be observed that after a transient phase, the plastic response ceases gradually. The strain-stress curve tends to coincide in each cycle, corresponding to 3 critical load paths as shown in Figure 4.Figure 8
Cyclic strain-stress curve (ε11−σ11) under 3 critical loads at r=a, θ=π/2 under shakedown limit load T3=0.44 for porosity f=0.01 with the following parameters of the solid matrix ν=0.2 and ϕ=20°.The described step-by-step computations are performed with the following parameters in Table1. Additionally, different cases of friction angle ϕ=10° and 20° with maximum water pressures pw+=Σm+/4 and 0.28 MPa are also considered.Table 1
Parameters of the material in numerical simulations.
ParametersYoung’s modulusEPoisson’s coefficientνPorosityfYield stressσyValue2 GPa0.20.011 MPa


### 4.2. Comparison and Verification of the Proposed Macroscopic Criterion

The numerical results obtained by the step-by-step elastoplastic computations together with the analytical macroscopic criterion are illustrated in Figures9–11 with respect to the friction angle of the solid matrix ϕ=10° and 20°. The parameters of the material ν=0.2 and f=0.01 are adopted. Owing to the axisymmetry of the shakedown domain (see Figure 5) and the existence of negative stresses in geotechnical engineering, only the part in the fourth quadrant will be compared and discussed.Figure 9
Comparison between the proposed macroscopic criterion and numerical results with respect topw+=0 and 0.25Σm+. The friction angle of the solid matrix ϕ=10°.Figure 10
Comparison between the proposed macroscopic criterion and numerical results with respect topw+=0 and 0.25Σm+. The friction angle of the solid matrix ϕ=20°.Figure 11
Comparison between the proposed macroscopic criterion and numerical results with respect topw+=0 and 0.28 MPa. The friction angle of the solid matrix ϕ=20°.In general (Figures9–11), an accordance between the results obtained by the proposed macroscopic criterion and the step-by-step computations can be observed, indicating that the established criterion has the ability to predict the long-term safety of the considered porous geomaterials under hydromechanical variable loads. More preciously, for the pure hydrostatic load case (Σeq+=0), the analytical solution σm+=31−f/3/2+3α3/2−3α (see [34] for the derivation) fits the value of numerical one, because of the fact that the constructed microscopic stress field is the exact one. On the other hand, with the increase of the shear loads (Σeq+), small differences are remarked due to the approximation of the deviatoric part of the trial stress field. The applied elastic and residual stress fields satisfy the statically admissible condition in an average way, which can be improved by a more accurate solution. Owing to the statical shakedown theorem, a quasi-lower bound is provided by the proposed criterion.Considering 1 single figure (Figures9 or 10), the effective safety load domain is bounded by the solid lines, in which the blue and red lines represent the yield surfaces of the proposed criterion with water pressure pw+=0 and 0.25Σm+, while the green line is the limit surface of the limit analysis-based one [29]. It can be concluded that the long-term strength of geomaterials is obviously reduced due to the variation of the water pressure. The reduction of the “width” of the yield surface can be estimated by the ration of pw+/Σm+−pw+. Besides, around the domain of pure deviatoric load (Σeq+=0), the effective safety domain is defined by the limit analysis-based criterion (10). The collapse of the material is owing to the development of a mechanism (incremental collapse in the first loading cycle) instead of fatigue.Figure11 illustrates the comparison between the analytical and numerical results with the water pressure pw+=0.28 MPa. The maximum value of the water pressure is a constant value, which is different from the two previous cases shown in Figures 9 and 10. The corresponding yield surface has been translated towards the right side, and the safety domain is accordingly reduced.


## 4.1. Implementation of the FEM-Based Computation

Owing to the geometrical symmetry, only a quarter of the hollow sphere model discretized by 1500 axisymmetric elements (see Figure6) is considered, which is refined enough to capture the reliable results. The homogeneous boundary condition v=D·x is imposed on the outer boundary.Figure 6
Initial mesh of the hollow sphere model.As described in Figure4, only 3 critical cyclic load paths containing all the vertices need to be considered to calculate the shakedown limit load, if the variation of inner water pressure is not concerned. On the contrary, under alternate wetting and drying condition, the load domain is modified as shown in Figure 7 located in the first and second quadrants. Consequently, 4 critical cyclic load paths have to be taken into consideration to fulfill the shakedown condition.Figure 7
Critical cyclic loadings of the load domain considering the alternate wetting and drying condition.Nevertheless, the macroscopic stress triaxiality is defined byT3=Σm+/Σeq+ at vertex P^3 describing the shape of the rectangular load domain. Moreover, the stress triaxiality T=Σm/Σeq for all the cyclic load path remains constant in a single cycle, which is implemented in Abaqus by the means of multipoint constraint subroutine (MPC) satisfying the imposed uniform velocity field boundary condition. Notice that such MPC subroutine is firstly provided in the study of polymeric materials [38] and latterly applied in a series of ductile porous media [30, 32, 35] with the same boundary condition.The above computation procedure is performed for an imposed value of stress triaxialityT3 with 100 cycles. The long-term stability of the considered material is evaluated during the transient phase before shakedown or collapse. In practice, the accumulated equivalent plastic strain at the “dangerous” point will be checked in the transient phase, which is supposed to be constant if it is in a shakedown state. By increasing the applied load, the shakedown limit is obtained when the collapse due to fatigue or incremental fracture occurs.The cyclic responses on the inner boundaryr=a at equator θ=π/2 of the hollow sphere at a shakedown state with T3=0.44 are plotted in Figure 8. The porosity f=0.01 and ν=0.2, ϕ=20° of the solid phase are considered. It can be observed that after a transient phase, the plastic response ceases gradually. The strain-stress curve tends to coincide in each cycle, corresponding to 3 critical load paths as shown in Figure 4.Figure 8
Cyclic strain-stress curve (ε11−σ11) under 3 critical loads at r=a, θ=π/2 under shakedown limit load T3=0.44 for porosity f=0.01 with the following parameters of the solid matrix ν=0.2 and ϕ=20°.The described step-by-step computations are performed with the following parameters in Table1. Additionally, different cases of friction angle ϕ=10° and 20° with maximum water pressures pw+=Σm+/4 and 0.28 MPa are also considered.Table 1
Parameters of the material in numerical simulations.
ParametersYoung’s modulusEPoisson’s coefficientνPorosityfYield stressσyValue2 GPa0.20.011 MPa


## 4.2. Comparison and Verification of the Proposed Macroscopic Criterion

The numerical results obtained by the step-by-step elastoplastic computations together with the analytical macroscopic criterion are illustrated in Figures9–11 with respect to the friction angle of the solid matrix ϕ=10° and 20°. The parameters of the material ν=0.2 and f=0.01 are adopted. Owing to the axisymmetry of the shakedown domain (see Figure 5) and the existence of negative stresses in geotechnical engineering, only the part in the fourth quadrant will be compared and discussed.Figure 9
Comparison between the proposed macroscopic criterion and numerical results with respect topw+=0 and 0.25Σm+. The friction angle of the solid matrix ϕ=10°.Figure 10
Comparison between the proposed macroscopic criterion and numerical results with respect topw+=0 and 0.25Σm+. The friction angle of the solid matrix ϕ=20°.Figure 11
Comparison between the proposed macroscopic criterion and numerical results with respect topw+=0 and 0.28 MPa. The friction angle of the solid matrix ϕ=20°.In general (Figures9–11), an accordance between the results obtained by the proposed macroscopic criterion and the step-by-step computations can be observed, indicating that the established criterion has the ability to predict the long-term safety of the considered porous geomaterials under hydromechanical variable loads. More preciously, for the pure hydrostatic load case (Σeq+=0), the analytical solution σm+=31−f/3/2+3α3/2−3α (see [34] for the derivation) fits the value of numerical one, because of the fact that the constructed microscopic stress field is the exact one. On the other hand, with the increase of the shear loads (Σeq+), small differences are remarked due to the approximation of the deviatoric part of the trial stress field. The applied elastic and residual stress fields satisfy the statically admissible condition in an average way, which can be improved by a more accurate solution. Owing to the statical shakedown theorem, a quasi-lower bound is provided by the proposed criterion.Considering 1 single figure (Figures9 or 10), the effective safety load domain is bounded by the solid lines, in which the blue and red lines represent the yield surfaces of the proposed criterion with water pressure pw+=0 and 0.25Σm+, while the green line is the limit surface of the limit analysis-based one [29]. It can be concluded that the long-term strength of geomaterials is obviously reduced due to the variation of the water pressure. The reduction of the “width” of the yield surface can be estimated by the ration of pw+/Σm+−pw+. Besides, around the domain of pure deviatoric load (Σeq+=0), the effective safety domain is defined by the limit analysis-based criterion (10). The collapse of the material is owing to the development of a mechanism (incremental collapse in the first loading cycle) instead of fatigue.Figure11 illustrates the comparison between the analytical and numerical results with the water pressure pw+=0.28 MPa. The maximum value of the water pressure is a constant value, which is different from the two previous cases shown in Figures 9 and 10. The corresponding yield surface has been translated towards the right side, and the safety domain is accordingly reduced.


## 5. Conclusion

In this paper, a macroscopic criterion based on a statical shakedown theorem to predict the long-term strength of porous geomaterials under alternate wetting and drying condition has been established. In the framework of micromechanics of porous materials, Gurson’s hollow sphere model is adopted, and the effects of alternate wetting and drying are considered as a variable inner pressure. The obtained criterion depends on the first and second invariants of the macroscopic stress tensor, Poisson’s coefficient, and the porosity. The reduction of the effective safety domain due to the alternate wetting condition is observed, depending on the amplitude of the variation of water pressure. In the end, the accuracy of the analytical solution has been evaluated by comparing with the results of FEM-based step-by-step computations.In the outlook, special efforts are needed to improve the accuracy of the analytical solution under deviatoric loads. The application of this method on structural computation is also a challenging topic.



---

*Source: 6617453-2021-02-08.xml*