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

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.


Introduction
Variable loadings exist widely in engineering structures, such as slops, offshore platform foundations, and pavements [1][2][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][5][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][12][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][16][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][22][23][24][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][30][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][33][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 Section 2, 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, stepby-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.

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.

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 Figure 1). 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. The internal and external radii of the chosen hollow sphere RVE are, respectively, noted as a 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/bÞ 3 . 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: where σ eq = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð3/2Þσ′ : σ′ q is the equivalent stress defined from the deviatoric part σ′ of the stress tensor σ and σ m = ð1/3ÞtrðσÞ is the mean stress. σ y represents the yield stress of the matrix material and α the pressure sensitivity factor related to the friction angle ϕ: 2 Geofluids 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]: Moreover, the macroscopic stress field Σ can be simply obtained by integration on a surface instead of a volume: if the microscopic one σ is statically admissible 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: 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: 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 frame fe ρ , e ϕ , e z g, the following three-parameter velocity field is adopted: where s = 1 ± α for A 0 ≷ 0. A 0 , A 1 , and A 2 are parameters.
The void growth can be computed as in which γ = 1 − ð1/sÞ, and the sign of A 0 represents void growth (A 0 > 0) or compression (A 0 < 0). By means of limit analysis techniques, an implicit macroscopic criterion of porous geomaterials composed of Drucker-Prager solid matrix is proposed as follows: where ∂Π/∂A 0 and ∂Π/∂D e can be, respectively, expressed by the Gauss hypergeometric function with ðaÞ n the Pochhammer symbols ðaÞ n = ðΓða + nÞÞ/ðΓðaÞÞ and the pair of ratios ϖ = ω/f 1/s , −∞ < ω < ∞.
Noticing that the invariants of the macroscopic stress tensor can be derived similar to their microscopic counterparts, Σ ′ being the deviatoric part of Σ. Guo's homogenized yield criteria FðΣ 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).
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 3 Geofluids 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].

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.

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 σ: where σ E ðx, tÞ is the purely elastic stress responses in the fictitious elastic cell. According to Melan's classical shakedown theorem [10], if there exists a time-independent residual stress field ρðxÞ, such that is 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: Owing to Equation (4), the following special relation relating the residual stresses at macro-and microscales is imposed by the surface integration: such that the average residual stresses vanish and are identical to the macroscopic elastic stresses in the fictitious elastic body: Let us introduce the load factor α, which controls the size of the actual load domain P = αP 0 . 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 4 Geofluids 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 where α is the same factor defined in P = αP 0 and Σ 0 k the reference loads. The load multipliers μ k ðtÞ are independent within the range μ − k < μ k ðtÞ < μ + 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.

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: 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: The elastic stress field σ E ð1Þ in the fictitious body is given in spherical coordinates: The corresponding residual stress field ρ ð1Þ writes as follows: (2) Similarly, the additional terms of the stress field under deviatoric load is where the fictitious elastic stress field writes as with ν being Poisson's coefficient and J 3 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 σ Eð2Þ and ρ ð2bÞ is related to σ Eð2Þ by where C 1 is the constant to be determined and KðrÞ 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): where J 3+ (resp., J 3− ) is the maximum value (resp., minimum value) of the third invariant in the deviatoric space, and 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: 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 1 ð0, 0Þ,P 2 ð0, Σ eq+ Þ,P 3 ðΣ m+ , Σ eq+ Þ, andP 4 ðΣ m+ , 0Þ.
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):  Figure 4: Critical cyclic loadings of the considered load domain.

Geofluids
Let us introduce the general stress triaxiality considering the pure elastic responses within the shakedown load domain, which is defined by for the load path connectingP 1 andP 3 . Likewise, for verticeŝ P 2 andP 4 , one has the following relation: Referring to [34], KðaÞ C 1 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 KðaÞ C 1 , one can obtain Taking the general stress triaxiality (34) into account, the macroscopic shakedown criterion of porous geomaterials derived from Equations (30) and (31) is 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 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 J 3+ , where 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.
Comparing with the analytic research of porous media under one independent load, the obtained shakedown domain in Figure 5 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 J 3+ 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.

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 pressure p w 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: in which the fictitious elastic response and corresponding residual stress field under variable water pressure are with p w+ being the maximum value of the water pressure and C w an additional constant. In engineering, the water pressure usually varies in the range 0 < p w < p w+ . Consequently, resulting from Equations (19) and (40), the total stress at the microscale including the effect of water pressure reads 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 Geofluids 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 p w 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 p w depends mostly on the water level in geomaterials. The detailed comparisons will be fully discussed in the next section.

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.

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 Figure 6) is considered, which is refined enough to capture the reliable results. The homogeneous boundary condition v = D · x is imposed on the outer boundary. As described in Figure 4, 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.
Nevertheless, the macroscopic stress triaxiality is defined by T 3 = Σ m+ /Σ eq+ at vertexP 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 triaxiality T 3 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 boundary r = a at equator θ = π/2 of the hollow sphere at a shakedown state with T 3 = 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 6: Initial mesh of the hollow sphere model.

Geofluids
The described step-by-step computations are performed with the following parameters in Table 1. Additionally, different cases of friction angle ϕ = 10°and 20°with maximum water pressures p w+ = Σ m+ /4 and 0.28 MPa are also considered.

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 Figures 9-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.
In general (Figures 9-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+ = ð3ð1 − 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 (Figures 9 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 p w+ = 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 p w+ /ðΣ m+ − p w+ Þ. 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. Figure 11 illustrates the comparison between the analytical and numerical results with the water pressure p w+ = 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.

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.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors state that they do not have any financial or nonfinancial conflict of interests.