An explicit formulation of the macroscopic strength criterion for porous media with pressure and Lode angle dependent matrix under axisymmetric loading

Jincheng Fan, Laurence Brassart, Wanqing Shen, Xiurun Ge

a Institute of Deep Earth Sciences and Green Energy, College of Civil and Transportation Engineering, Shenzhen University, Shenzhen, 518060, China

b School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai, 200240, China

c Department of Engineering Science, University of Oxford, Oxford, OX1 3PJ, UK

d Department of Materials Science and Engineering, Monash University, Clayton, VIC, 3800, Australia

e Univ. Lille, CNRS, Centrale Lille, UMR 9013 - LaMcube - Laboratoire de Mécanique, Multiphysique, Multi-échelle, 59000, Lille, France

Keywords:Porous materials Unified strength criterion Homogenization Gurson criterion

ABSTRACT This paper aims to propose an explicit formulation of the macroscopic strength criterion for porous media with spherical voids.The matrix is assumed rigid and perfectly plastic with yield surface described by the three-parameter strength criterion, which is Lode angle and pressure dependent and capable of accounting for distinct values of the uniaxial tensile strength, uniaxial compressive strength (UCS) and equal biaxial compressive strength(eBCS). An exact upper bound of the macroscopic strength is derived for porous media subjected to purely hydrostatic loading. Besides, an estimate of the macroscopic strength profile of porous media under axisymmetric loading is obtained in parametric form.Moreover,a heuristic strength criterion in explicit form is further developed by examining limit cases of the parametric strength criterion.The developed strength criteria are assessed by finite-element based numerical solutions. Compared with the parametric strength criterion which involves cumbersome functions, the heuristic one is convenient for practical applications.For specific values of the matrix’s strength surface,the proposed heuristic strength criterion can recover the well-known Gurson criterion.The present work also addresses the effect of the ratio of matrix’s eBCS to UCS on the macroscopic strength of porous media.For matrix with distinct values of eBCS and UCS,neglecting the difference between eBCS and UCS would result in an underestimation of the macroscopic strength, especially when the pressure is large.? 2021 Institute of Rock and Soil Mechanics, Chinese Academy of Sciences. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction

Strength prediction of porous media has been one of the research focuses recently. A number of theoretical investigations have been conducted to determine the macroscopic strength of porous materials for various local yield criteria. Gurson (1977)considered hollow sphere and cylinder geometries to represent porous materials with von Mises matrix, and derived the effective yield function using kinematic limit analysis. Following this pioneering work, Cazacu and Revil-Baudard (2017) proposed an analytic yield criterion for porous solids with pressure insensitive matrix described by a plastic potential depending on the second and third invariants,and compared it with the macroscopic criteria for porous Mises and Tresca materials. The study of macroscopic strength of porous material with pressure-sensitive matrix obeying Drucker-Prager (DP) criterion was conducted by Jeong and Pan(1995) and later extended by Guo et al. (2008). Recently, Shen et al. (2020) reviewed the macroscopic strength criteria for porous materials with DP matrix and evaluated them through numerical solutions. Based on an improvement of the strength predictions for purely shear loading, Shen et al. (2020) proposed one heuristic strength criterion, which could give better predictions of the macroscopic strength. For materials with local yield function depending on the third invariant of stress, one can refer to recent investigations in the context of Tresca matrix(Cazacu et al.,2014),Mohr-Coulomb (MC) matrix (Anoukou et al., 2016) and matrix described by the strength criterion proposed by Bigoni and Piccolroaz (Brach et al., 2018). As an alternative method, statical limit analysis has been adopted to formulate the macroscopic strength criterion for porous materials with von Mises matrix(Cheng et al.,2014;Shen et al.,2015)and Hill orthotropic matrix(El Ghezal et al., 2017). Effect of the void shape on the macroscopic strength of porous media was discussed in Shen et al.(2017)for DP matrix and Monchiet et al. (2008) for Hill-type matrix.

An accurate description of the macroscopic strength of a particular porous material is contingent on the accuracy of the local plasticity model(Revil-Baudard and Cazacu,2014;Cazacu and Revil-Baudard,2017).Strength of geomaterials,such as rocks and concrete are pressure and Lode angle dependent(Chemenda and Mas,2016).Neither the von Mises criterion nor the DP criterion can account for the effectof the thirdinvariantof stress.Incontrast,the Tresca model depends on the third invariant,however it is pressure-independent and hence does not apply to pressure-dependent materials.Although the MC criterion has been widely adopted in the mechanics community to describe the shear dominated failure mechanism (Jaeger et al., 2009; Labuz and Zang, 2012), it neglects the effect of the middle principal stress and fails for materials with distinct values of uniaxial compressive strength (UCS) and equal biaxial compressive strength (eBCS). Besides, the established strength criterion for porous materials with MC matrix under axisymmetric loading in Anoukou et al.(2016)is in parametric form,and the parametric functions are cumbersome.In fact,in the present work we also developed an explicit formulation of the macroscopic strength criterion for porous MC materials under axisymmetric loading,which is more applicable for practical applications.

Recent studies have revealed that the eBCS of geomaterials such as concrete and rocks is larger than their UCS (Kupfer et al.,1969;Brown, 1974; Amadei and Robison, 1986; Hussein and Marzouk,2000). Lee et al. (2004) experimentally investigated the failure of plain concrete and showed the eBCS was about 17%higher than the UCS. Sirijaroonchai et al. (2010) performed equal biaxial compression tests on high-performance fiber-reinforced cement material and concluded that the eBCS was about 1.5 and 1.6 times the UCS respectively for hooked and spectra fiber-reinforced cement material. Guo (2014) reviewed the main experimental results for the eBCS of concrete and pointed out the ratio of eBCS to UCS of concrete lay in 1.15-1.35. Experimental study of the eBCS and UCS of concrete at early ages has been conducted by Dong et al.(2016),and experimental results revealed the ratio of eBCS to UCS approximately decreased from 3.5 to 1.2 within the first 7 d and remained at a value of 1.15 up to the age of 28 d. Yun et al. (2010) observed that the difference between the BCS and UCS of granite increased dramatically with the confining pressure up to some threshold value, beyond which the BCS dropped, but still remained higher than the UCS.In that study,the eBCS of granite was about 1.13 times the UCS on average.Although the MC criterion can account for the strength difference under uniaxial tensile and compressive loading,it cannot account for different eBCS and UCS values.In this respect,the established macroscopic yield functions in the literature are not expected to be accurate for such porous materials.

Lately, the three-parameter strength criterion has been considered as an appropriate candidate to describe the yield surface of geomaterials (Yu, 2017, 2018). Through appropriate choice of parameter values in the expression of the yield function,the threeparameter strength criterion can retrieve classical strength criteria,such asthe MC and Tresca criteria.It can also linearlyapproximate the well-known Matsuoka-Nakai criterion (Matsuoka and Nakai,1974)and Willam-Warnke criterion(Willam,1975).Thanks to its ability to account for distinct eBCS and UCS, the three-parameter strength criterion is promising in view of describing the behavior of a wide classofmaterials.Motivatedbytheseconsiderations,hereweattempt to derive macroscopic strength functions for porous media whose local behavior obey the three-parameter strength criterion.

The paper is organized as follows. In Section 2, the kinematic limit analysis-based homogenization is briefly reviewed.In Section 3, the plasticity model of the matrix is presented and the plastic dissipation function for the three-parameter strength criterion is derived.Section 4 delivers the exact solutions of limit stresses and velocity fields in a hollow sphere subjected to isotropic loading with plastic behavior obeying the three-parameter strength criterion,which is used to construct the trial velocity field in the kinematic limit analysis. An exact solution of the upper bound of the macroscopic strength of the hollow sphere under purely hydrostatic loading is also given. In Section 5, an estimate of the macroscopic strength for the axisymmetric loading is proposed in parametric form, followed by an improvement of the established criterion to allow the exact upper bound solution to be fully retrieved in the case of purely hydrostatic loading. Numerical assessment of the parametric strength criteria is performed in Section 6.For practical applications,we develop a heuristic strength criterion in explicit form for porous media subjected to axisymmetric loading in Section 7.The concluding remarks are given in Section 8.

2. Theoretical background

We consider a porous material occupying a domain Ω consisting of traction-free pores, ω0, embedded in a rigid, perfectly-plastic matrix ω1. The strength of matrix is specified by a yield function F（σ） defining a convex set C of plastically admissible stresses:

We assume that the local flow behavior complies with Hill’s maximum dissipation principle (Hill, 1950), which states that,among all plastically admissible stresses, the actual stress corresponding to a given strain rate d maximizes the plastic dissipation:

This equation identifies π(d) as the support function of the convex set C.The maximum dissipation principle implies normality(associativity) of the plastic flow rule:

where ˙γ is the plastic multiplier, such that ˙γ ≥0 and ˙γF = 0(Kuhn-Tucker conditions).When the yield surface is not regular,the condition generalises as d ∈NC(σ),where NCis the cone of outward normals to C at σ (Maugin,1992).

The domain Ω is subjected to prescribed velocity boundary conditions on its external boundary ?Ω corresponding to a given macroscopic strain rate D:

Here we assume that the pores do not intersect the domain boundary. The affine velocity boundary conditions ensure that the macroscopic strain rate equals the average microscopic strain rate:

where n is the outward unit normal to the external boundary ?Ω.The above expression based on surface integral was adopted because the strain rate is undefined within the pore volume. For non-porous material,it is equivalent to the classical volume average strain rate.On the other hand, the macroscopic stress Σ is defined as the volume average of the microscopic stress field:

where the stress vanishes within the pores.

We introduce the set K（D） of kinematically and plastically admissible velocity fields:

Kinematically admissible velocity fields are continuous within the matrix and satisfy the affine boundary conditions (Eq. (4)).Plastically admissible strain rates are compatible with the plastic flow rule (Eq. (3)). Plastically admissible conditions corresponding to the three-parameter strength criterion will be given in the next section. The set of statically admissible stress fields is given by

where ?ω0represents the pore surface. The consistency between macroscopic and microscopic work rates is expressed by the Hill-Mandel lemma (Hill,1967):

which holds for any admissible couple (d,σ*), where d and σ* are not necessarily related by the flow rule.

Let us introduce the macroscopic dissipation function Π(D) = 〈π(d)〉. From Eq. (2) together with Hill’s lemma, it follows that (Suquet,1985):

By definition, the macroscopic dissipation Π(D) is the support function ofIt follows that the macroscopic flow rule also complies with a normality principle,and that the macroscopic strength can be obtained from the macroscopic dissipation as

For a given macroscopic strain rate D,the macroscopic dissipation is the solution to the following variational problem(Lubliner,2008):

The optimal velocity field v is the solution of the mechanical boundary value problem, and the corresponding stress field satisfies the local flow rule (Eq. (3)) and the mechanical equilibrium condition ?·σ =0.

In practice,an exact analytical solution to the problem(Eq.(13))is often beyond reach. Upper bounds on the dissipation function can then be obtained by considering a trial field

We follow the approach proposed by Guo et al. (2008) and Anoukou et al. (2016), and consider a class of kinematically admissible trial fields defined by a limited number of adjustable parameter. Such trial fields thus belong to a subset of the set of kinematically admissible fields,Optimal values of the parameters (leading to the lowest upper bound) are then obtained by minimizing ^Π:

3. The three-parameter strength criterion

3.1. Yield function

Shear stress-based strength criteria have been extensively developed to describe failure dominatedbyshear-slip mechanisms.For any stress state,there exist three principal shear stresses,among which only two are independent.Different from the classical Tresca and Von Mises criteria, which respectively take into account the maximum shear stress and all the three principal shear stresses,Yu(1983)proposed a twin-shear yield criterion that assumes that yielding starts when the sum of the two larger principal shear stresses reaches a critical value.Later,considering the different effects of the two larger principal shear stresses on the yielding of materials,Yu and He(1992)proposed a unified yield criterion through a simple and unified mathematical formula,which has great flexibility for reproducing a series of strength-domain shapes.The unified yield criterion can fully retrieve the Tresca criterion and linearly approximate the Von Mises criterion.To account for the effect of normal stresses on the yielding,Yu et al.(1999)further developed the unified strength theory,which is capable of retrieving the unified yield criterion and MC criterion.A complete review of the development of Yu’s strength theory can be found in Yu(2018).As a generalization of the unified strength criterion,the three-parameter strength criterion is superior in characterizing distinct values of uniaxial tensile strength(UTS),UCS and eBCS(Yu,2017).

The three-parameter strength criterion can be expressed in terms of principal stresses as follows (Yu, 2017):

where the three parameters M, N and Q can be identified from classical experimental tests, as explained below. The additional variable b ∈[0,1]represents the relative effect of the intermediate principal stress σ2.By varying b from 0 to 1,a series of yield surfaces(i.e. distinct material models) can be obtained (see Fig.1).

The above yield functions can also be written in terms of stress invariants (,) and Lode angle (θσ):

By evaluating the yield condition corresponding to the yield functions (Eqs. (17) and (18)) in the specific cases of uniaxial tension,uniaxial compression and equal biaxial compression,the three parameters M, N and Q can be expressed in terms of the UTS (σt),UCS(σc)and eBCS(σcc).Through introducing α=σt/σcand β=σcc/σc, the following relations are obtained:

Experimental values for media usually give α ≤1 and β ≥1. It can be verified that these assumptions on the strength ratios yield 0 ≤M ≤1, N ≥0 and Q > 0. If the material exhibits the same strength under uniaxial compression and equal biaxial compression(i.e.β=1),we have N=0 by Eq.(21)and the three-parameter strength criterion reduces to the unified strength theory(Yu et al.,1999).

Fig.1. Various yield surfaces of the three-parameter strength criterion corresponding to different values of the parameter b.The yield surface is represented in the deviatoric plane (here the principal deviatoric stresses, and are not ordered).

3.2. Maximum plastic dissipation

The maximum plastic dissipation function is essential for the derivation of the homogenized strength criterion for porous materials. The maximum dissipation functions for classical yield criteria, such as Tresca, Mises, DP and MC have been developed.However, to the best of our knowledge, the maximum plastic dissipation for the three-parameter yield criterion has not yet been proposed.

The maximum dissipation function is given by π(d) = σ : d,where σ satisfies F（σ） = 0 and is related to d by the flow rule.For all three regimes and arbitrary value of b ∈[0,1], it can be calculated that:

4. Macroscopic strength criterion under purely hydrostatic loading

In this section, we will establish the macroscopic strength criterion for porous media subjected to purely hydrostatic loading.Following Gurson (1977),we consider the simple hollow sphere model as an idealization of the porous structure (see Fig. 2 in which ridenotes the inner radius and reis the external radius). The porosity is given by f =/. The inner surface is traction free.

Fig. 2. The hollow sphere model.

4.1. Construction of the trial velocity field for a hollow sphere under isotropic loading

To construct the trial velocities used in the kinematic limit analysis,the limit velocity field in the hollow sphere subjected to isotropic loading is firstly carefully analyzed, followed by the development of the macroscopic strength criterion for porous media under purely hydrostatic loading. Since the plastic flow rule depends on the stress states (see Eq. (3)), the static analysis of the stress field in the hollow sphere subjected to isotropic loading is presented before the kinematic analysis of the velocity field.

4.1.1. Static analysis

Due to the spherical symmetry of the geometry and loading conditions,the only non-zero components of stress in the spherical coordinate system represented in Fig.2 are σr(r),σθ(r)and σφ(r).In the absence of body force,static equilibrium requires that:

In general, we may have either σr> σθ= σφor σθ= σφ> σr,depending on the sign of the applied pressure on the external boundary. In any case, it is possible to express the relationship between the limit stress components in the following form:

where the constants A and B can be identified from the corresponding yield function. Inserting Eq. (24) into the equilibrium equation in Eq. (23) yields an ordinary differential equation for σr(r).Integrating along the radius and using the inner free boundary condition,σr(ri) = 0, we successively obtain

(1) Case σr>σθ=σφ

In this case, the yield function (Eq. (17)) applies, and the two constants A and B in Eq. (24) are identified as

where ι=M+4 N/[3(1+b)]=1-2α/[β(1+α)]<1.It follows that A < 0 and B < 0. According to Eqs. (25) and (26), it follows that σθ= σφ≤σr< 0. Thus, this case corresponds to a hollow sphere under isotropic compressive loading.

(2) Case σθ=σφ>σr

This case corresponds to the yield function (Eq. (18)), and the two constants A and B are calculated as

Obviously,A>0 and B>0.With reference to Eqs.(25)and(26),we can then show that σθ=σφ≥σr>0.This case thus corresponds to a hollow sphere under isotropic tensile loading.

Inserting Eq. (27) or (28) into Eqs. (25) and (26), the limit stresses can be expressed as

where Aε= 2 [(1 + b)M + N]/[(1 + b) (ι + ε)], in which ε = 1 for isotropic tensile loading and ε = -1 for isotropic compressive loading. In particular, the limit pressure p applied on the external surface of the sphere is easily obtained as p=-σr(re)(p is assumed positive for compressive loading).

4.1.2. Kinematic analysis

Due to the spherical symmetry of the geometry and loading conditions, the velocity field is purely radial, i.e. v = vrer, and the non-zero strain rate components are given by It can be verified from the definition of the Lode angle that the case σr> σθ= σφcorresponds to θσ= 0 and the case σθ= σφ> σrcorresponds to θσ=π/3.These stresses are located on edges of the yield surface of three-parameter strength criterion except for b=1.The plastic flow at such stresses should therefore comply with the generalized flow rule:

where F+and F-are the yield functions of the two intersecting faces, and γ+and γ-are the corresponding non-negative plastic multipliers.

(1) Case σr>σθ=σφ

Setting σ1=σrand σ2=σθ=σφ=σ3,straightforward evaluation of the flow rule on the edge θσ= 0 gives

According to Eq. (31), the strain components in Eqs. (34) and(35) are equal, which results inThus, Eqs. (33)-(35) can be rewritten as

Substituting Eqs. (36) and (37) into Eq. (31), we obtain the following ODE for the plastic multiplier:

Integrating along the radius then leads to

Inserting Eq. (39) into Eq. (37) and using the fact that vr= rdθ,the radial velocity is obtained:

For prescribed velocity boundary condition on the external surface,vr(re)=ve,and the unknown plastic multiplier on the pore surfacecan be identified.One also readily verifies from the latter equation that the surface velocity veis negative (remember that ι<1).In other words,the case σr>σθ=σφcorresponds to isotropic compression, consistently with the static analysis.

(2) Case σθ=σφ>σr

Setting σ1=σθ=σφ=σ2and σ3=σr,evaluation of the flow rule on the θσ=π/3 gives

From the equality dθ = dφ, it also follows thatThe integration procedure to obtain the plastic multiplier and the radial velocity is identical to the previous case, and therefore we only provide the final result:

Evaluation of the above expression for r = reallows then to expressin terms of the surface velocity.In this case,the velocity ve=v(re)is positive,and this case corresponds to isotropic tension.Eqs. (40) and (44) can be rewritten in a unified way:

where ε=1 if ve>0(isotropic tensile loading)and ε=-1 if ve<0(isotropic compressive loading).

4.2. Formulation of the macroscopic strength criterion under purely hydrostatic loading

The macroscopic strength criterion for porous media under purely hydrostatic loading is established in the following. In the spherical coordinate framework(r,θ,φ),the macroscopic strain rate corresponding to purely hydrostatic loading can be expressed as

For purely hydrostatic loading, the plastically admissible velocity in the matrix is exactly the limit velocity under isotropic loading(i.e.Eq.(45)).With reference to Eq.(4),the trial velocity can be expressed in terms of the mean macroscopic strain rate as

where ζ = (3 -εM)/[3 (1 + ει)].

According to Eq. (31), the non-zero strain rate components corresponding to the velocity in Eq. (47) can be computed as

The trace of such strain rate is expressed as

Inserting Eq. (49) into Eq. (22) yields

The macroscopic dissipation function under purely hydrostatic loading can be calculated as

According to Eq. (12), the macroscopic strength under purely hydrostatic loading is

5. Macroscopic strength criterion under axisymmetric loading

Since the present work concentrates on the macroscopic strength criterion for porous materials which are subjected to axisymmetric loading, we focus on the case where b = 0, which corresponds to the innermost yield locus among all possible values for b (Fig.1). In principle, the case 0 < b ≤1 could be addressed following a similar procedure, at the expense of heavier calculations. As mentioned in Section 4, the hollow sphere model is adopted.

The macroscopic strain rate is chosen in the following form:

where eρ, eφand ezare the unit vectors of the local basis. The corresponding mean and equivalent macroscopic strain rates are

In the following,we propose a set of trial velocity fields defined by one adjustable parameter and derive the corresponding conditions of plastic admissibility.Next,we calculate the optimal value of the adjustable parameter by addressing the minimum problem(Eq.(16)) and obtain the corresponding limit stress by differentiation.To this end, note the following equalities:

5.1. Trial velocity field under axisymmetric loading

Since the hollow sphere model is subjected to axisymmetric loading, the velocity field in the matrix should be in the following form:

Following Guo et al. (2008), Anoukou et al. (2016) and Brach et al. (2018), we consider trial velocity fields as

The first term vhomis defined as

where C1and C2are two constants.The superscript“hom”refers to the fact that the strain rate resulting from this velocity field is homogeneous. The second term on the right-hand side of Eq. (57) is chosen of the same form as the exact velocity in the hollow sphere under isotropic loading derived in the previous section, here expressed in the cylindrical coordinate system:

where w = C0（re/r）3ζ, ζ = (3 -εM)/[3 (1 + ει)], in which ε = 1 if C0>0 and ε=-1 if C0<0.Also note that ζ<1 if ε=1 and ζ>1 if ε=-1.The superscript“het”is introduced to denote that the strain rate resulting from this velocity field is heterogeneous.We note the radial coordinate is related to the polar coordinates byand sin θ = ρ/r, cos θ = z/r. Note that the adjustable parameter C0does not coincide with the surface velocity of the sphere under axisymmetric loading.

Kinematically admissible conditions require the trial velocity field(Eq.(57))to be compatible with the affine boundary condition(Eq. (4)), which leads to

Thus,among the three adjustable parameters C0,C1and C2,only one of them is independent. We consider C0as the only independent parameter.

5.2. Plastically admissible conditions

In the kinematic limit analysis,the local strain rate derived from the trial velocity field should be plastically admissible,i.e.it should be compatible with the flow rule. Since the yield surface is not smooth,plastically admissible conditions must be derived for three possible regimes, depending on whether the stress lies on a face,edge or apex of the yield surface. In the case where b = 0, the plastically admissible condition can be summarized as |d1|+ |d2|+|d3| ≤（3+ N）dm/（M+ N）. In principle, the relation between the independent variable C0and the macroscopic strain rate parameters(Dmand Deq)should exactly comply with such local plastically admissible condition to make the support function finite in the matrix. Due to the complex relation between the principal strain rates and C0,Dmand Deq,it is not possible to obtain a deterministic relation accounting for all the local conditions. Following the approach adopted in Guo et al. (2008), Anoukou et al. (2016) and Brach et al. (2018), instead of satisfying the plastically admissible conditions locally, the parameters are required to obey the plastically admissible conditions in an average sense:

where〈dm〉1represents the volume average of the mean trial strain rate in the matrix domain ω1and can be calculated as 〈dm〉1=Dm+ （f - f1-ζ）C0/（1 - f）, ω = 2ζ|C0|/Deq, J（ω） is provided in Appendix A, and H is a function of Deqand C0.

5.3. Formulation of the parametric strength criterion under axisymmetric loading

Inserting Eq. (22) into Eq. (16), we can reformulate the minimum problem as follows:

which is subjected to the average plastically admissible constraint(Eq.(61)).

The constrained minimum problem is to be solved using the Lagrangian method combined with the Karush-Kuhn-Tucker(KKT)condition, with the Lagrangian given by

where λ is the Lagrange multiplier.

The solution to the problem is given by λ > 0 and C0which satisfy the following conditions:

Eq. (64) shows that plastic admissibility requires that the strict equality holds in Eq.(61).Eq.(65)can be solved for λ,giving λ as a function of C0and the macroscopic strain rate:

where I0is introduced to ensure the continuity of the derived macroscopic strength function and given in Appendix A.

However, in general, it is not possible to solve Eq. (64) analytically and obtain C0as an explicit function of D.Therefore,analytical expression of the macroscopic dissipation function in terms of the optimal value of the adjustable parameter C0cannot be found,making calculation of the macroscopic stress by analytical differentiation impossible.Therefore,in the following,an alternative(but equivalent) approach is pursued.

First, the variation of Eq. (62) is calculated as

Besides, the variation of the plastic admissibility condition (Eq.(64)) yields

where

in which Ieq= I0/（1-ζ）is introduced to ensure the continuity of the macroscopic strength.

Inserting Eq. (68) into Eq. (67), we obtain

From Eq. (69), the macroscopic mean and equivalent limit stresses can be identified:

In these expressions,the value of the parameter C0needs to be calculated numerically by solving Eq. (64).

5.4. Improvement of the parametric strength criterion (Eqs. (70)and (71))

Due to the fact that the purely hydrostatic loading can be considered as a limiting case of axisymmetric loading, the macroscopic strength of porous media subjected to purely hydrostatic loading can be estimated by the parametric strength criterion(Eqs.(70) and (71)). Note that we have given the exact solution of macroscropic hydrostatic strength (Eq. (52)) in Section 4. The estimated strength predicted by the parametric strength criterion will be compared with the exact solution in the following.

As the general state of axisymmetric loading goes to the limit state of macroscopic purely hydrostatic loading, the trial velocity components vhetin Eq. (57) dominates the velocity field in the matrix. Remember that ω = 2ζ|C0|/Deqand Deq= 2|C1-C2|/3.Thus,ω →+∞for purely hydrostatic loading,which is followed by

where η1= (M + N)/(3 - M) and η2= (M + N)/(3 + 2 N + M).

Inserting Eq. (72) into the parametric strength criterion (Eqs.(70) and (71)) results in

which is the estimated macroscopic strength under purely hydrostatic loading obtained by the parametric strength criterion. We note that the estimated macroscopic hydrostatic tensile strength(MHTS) corresponds to ε = 1 and the estimated macroscopic hydrostatic compressive strength (MHCS) corresponds to ε = -1.

The exact macroscopic hydrostatic strength(Eq.(52))in the case where b = 0 can be calculated as

Further analysis shows that the hydrostatic strength predicted by Eq. (73) is not always equal to the exact value calculated by Eq.(74).To illustrate it,we compare the estimated hydrostatic strength with the exact solution in Fig.3.As β=1,the estimated hydrostatic strength coincides with the exact value. However, as β = 1.15, the MHCS is underestimated by Eq.(73),especially when the ratio α of matrix’s UTS to UCS is lower. Thus, the parametric strength criterion cannot always recover the exact solution for purely hydrostatic compression, which is due to the relaxed plastic admissibility condition (Eq. (61)) used for the derivation of the parametric strength criterion.The averaged constraint is prone to broaden the space of admissible trial velocity, and thus the minimum problem(Eq. (62)) might be solved with a broadened constraint space,which results in an underestimation of the macroscopic strength.

To solve this issue, it is assumed that η satisfies the following equation:

Solving Eq. (75), we can obtain

which means that only if η in Eq.(72)satisfies the above equation,the parametric strength criterion could exactly recover the macroscopic strength under purely hydrostatic loading.

In fact, as3= 1, η in Eq. (72) automatically complies with the condition (Eq. (76)), and thus the estimated MHTS coincides well with the exact solution(see Fig.3).However,as3=-1,η in Eq.(72)does not necessarily satisfies such condition unless in the case where N = 0 (β = 1), and the estimated MHCS is lower than the exact solution as β>1.As a result,an improved parametric strength criterion is proposed by revising η2as

6. Assessment of the macroscopic strength criteria with finite element solutions

In this section, the macroscopic strength criterion for porous media under purely hydrostatic loading and the parametric strength criteria for porous media subjected to axisymmetric loading are assessed with finite element (FE) based numerical solutions. The thick-walled spherical shell (see Fig. 2) is numerically investigated with the FE software ABAQUS.Due to the symmetry,a quarter of the shell is investigated by an axisymmetric model with a discretization of 1050 biquadratic axisymmetric quadrilateral elements CAX8, as illustrated in Fig. 4 for the case at f = 0.1. For comparison, the solid matrix is assumed elastic-perfectly plastic,obeying the three-parameter strength functions (Eqs. (17) and(18)), and the plastic deformation obeys the associated flow rule.Since the routine for solving elasto-plastic behavior governed by the three-parameter strength criterion is not included in the software ABAQUS, a User Material Subroutine was developed. In this routine, the normality rule was used wherever the normal to the yield surface of the three-parameter strength criterion is unambiguously defined. At the singular point of the yield surface, the approach developed by Koiter(1953)was applied to define the flow direction. For state update, the implicit elastic predictor/return mapping scheme introduced in de Souza Neto et al. (2009) was applied.Following the numerical procedure developed by Faleskog et al.(1998)and Cheng and Guo(2007),which enables the loading to be applied with a constant macroscopic stress triaxiality (Σm/Σeq),a user subroutine Multi-Points Constraints was used to apply the uniform strain rate boundary condition on the outer surface.

Fig. 3. Comparison of the macroscopic strength estimated by the original parametric strength criterion (Eq. (73)) with the exact upper bound (Eq. (74)) under purely hydrostatic loading at different values of α: (a) β = 1, f = 0.2; and (b) β = 1.15, f = 0.2. Σcm: MHCS; Σtm: MHTS.

Fig. 4. Meshes used in the FE solution for f = 0.1.

First, the macroscopic strength of porous media under purely hydrostatic loading is studied. The comparison of the hydrostatic strength computed by Eq.(74)with the numerical results is given in Table 1. The FE solutions agree well with the exact values. The maximum deviation reaches 0.93% for MHCS and 0.02% for MHTS,which can serve as a cross-validation of these approaches.

Table 1 Macroscopic strength of porous media under purely hydrostatic loading at α = 0.33.

Next, the yield profile of the porous material with α = 0.59,β=1.15 and f=0.1,loaded by axisymmetric loading is investigated respectively by the original parametric strength criterion, the improved parametric strength criterion and the FE analysis. The strength profiles are depicted in the (normalized) stress space(Σm/σt, (Σ11- Σ33)/σt) (see Fig. 5). Indeed, for axisymmetric loading, the Cauchy stress is expressed as Σ = Σ11(e1?e1+ e2?e2) + Σ33(e3?e3), and the third invariant of such deviatoric stress is= - 2（Σ11-Σ33）3/27. For axisymmetric stress such that Σ11< Σ33,>0, and otherwise,<0. For the parametric criteria,the solutions to case A (3=1, sgn（Dρ-Dz） =1)and case B(3=1,sgn（Dρ-Dz） = -1)respectively correspond to the strength loci in the first quadrant (Σm> 0,<0) and the fourth quadrant (Σm> 0,>0). The solutions to case C (3= -1,sgn（Dρ-Dz） = - 1) and case D (3= -1, sgn（Dρ-Dz） = 1)respectively correspond to the strength loci in the third quadrant(Σm< 0,>0) and the second quadrant (Σm< 0,<0). Irrespective of the parametric criteria and the FE approach, the presence of voids in the matrix leads to the closing of the yield profile, i.e. the macroscopic strength of porous materials is no longer infinite when they are subjected to loading with large pressure. The yield locus obtained by the improved parametric strength criterion is in a good agreement with the FE solution,whereas the strength predicted by the original parametric criterion is obviously lower than the FE solution as Σm/σt<-2,which shows the superiority of the improved parametric criterion over the original parametric criterion.

In the following,we apply the improved parametric criterion to study the yield loci of porous materials.As β=1,i.e.the UCS equals the eBCS, the solid matrix obeys the MC criterion. For illustration,we take for example porous MC materials with frictional angles of φ = 30°, 25°and 20°, which have been investigated in Anoukou et al. (2016) and Pastor et al. (2016). Correspondingly, the ratios of UTS to UCS are α=0.33,0.41 and 0.49.As seen in Fig.6,φ has a strengthening effect on MHCS and weakening effect on MHTS.As φ goes from 20°to 30°, the MHCS increases from 10.8c to 37.2c,

Fig. 5. Comparison of yield profiles predicted by the original parametric strength criterion (dash line), the improved parametric strength criterion (solid line) with FE solution (points) in the (Σm/σt, (Σ11 - Σ33)/σt) space (α = 0.59, β = 1.15, f = 0.1).

Fig. 6. Comparison of yield profiles predicted by the improved parametric strength criterion (solid line) with FE solution (points) at different values of φ in the (Σm/c,(Σ11 - Σ33)/c) space (β = 1, f = 0.1).

whereas the MHTS decreases from 1.5c to 1.1c. Irrespective of the improved parametric criterion and the FE analysis, φ has a strengthening effect on the macroscopic strength as Σm/c<-5.For stress states corresponding to Σm/c>-5,the strengthening effect decreases with increasing φ. As Σmapproaches the MHTS, φ will have a weakening influence on the macroscopic strength. The predicted macroscopic strength by the improved parametric criteiron is in a good agreement with the numerical solution except for that in the third quadrant,where the predicted strength locates inside the FE solution and the difference between the strength predicted by the improved parametric criterion and FE analysis tends to become larger as φ increases.

7. Heuristic strength criterion under axisymmetric loading

The improved parametric strength criterion for porous media is not convenient for practical use due to the involving cumbersome functions (see Appendix A). Besides, the estimated macroscopic strength by the improved parametric strength criterion may exhibit a large difference from the FE solution in the quadrant (Σm< 0,>0)(see Fig.6).Following the strategy adopted by Jeong and Pan(1995), we will derive an explicit equation of the macroscopic criterion for porous media under axisymmetric loading based on the specific macroscopic strength in the following limiting cases: (1)purely hydrostatic loading,Σeq=0;(2)deviatoric loading,Σm=0;and (3) porosity, f →0.

In the first case under purely hydrostatic loading, the macroscopic strength (Eq. (74)) can be written in an equivalent way as

In the second case, as the hollow sphere is subjected to deviatoric loading, the trial velocity components vhomdominates the velocity field and C0→0. Thus, the macroscopic strength under deviatoric loading can be obtained with the improved parametric criterion by letting ω →0:

and it is noted that under axisymmetric loading, Σeq=- sgn

In the third case, as the porosity goes to zero, the material behaves like the matrix and the macroscopic strength can be approximated by the matrix’s yield function (Eq. (17) and (18)).Therefore, under axisymmetric loading, the macroscopic strength of materials with vanishing pores can be computed as

which can be rewritten in the quadratic form as

Based on such specific macroscopic strength(Eqs.(78)-(81)),a general equation of the heuristic strength criterion is assumed to take the following form:

where G is a dimensionless parameter that satisfies

Obviously, the developed heuristic strength criterion (Eq. (7))recovers the specific strength in the limiting cases. It should be noted that Σmand sgnrespectively represent the dependence of macroscopic strength on pressure and Lode angle. In contrast with the parametric form of the macroscopic strength criterion,the heuristic strength criterion gives an explicit and simple relation between Σmand Σeq. Hereafter, we also call the heuristic strength criterion as the explicit form of the macroscopic strength criterion or the explicit strength criterion.

Using Eq.(21),the explicit strength criterion can be expressed as

where

As the eBCS of the matrix equals UCS, i.e.β = 1, the yield criterion of the matrix reduces to the Mohr-Coulomb criterion. The tensile strength σtand the ratio α of the matrix’s UTS to UCS can be expressed in terms of the frictional angle φ and cohesion c as

Inserting Eq. (85) and β = 1 into Eq. (84), we can obtain an explicit strength criterion for porous materials with Mohr-Coulomb matrix under axisymmetric loading:

where

As φ →0,the MC criterion retrieve the Tresca criterion. Due to Eq. (86), there exist

According to the L′Hospital’ s rule, we can obtain

Inserting Eqs.(87)and(88)into Eq.(86)leads to the formulation of explicit strength criterion for porous materials having Tresca matrix under axisymmetric loading:

which is consistent with the Gurson model.

Assessment of the heuristic strength criterion in explicit form(Eq.(84))is conducted with FE analysis in the remaining part.The function g in Eq. (84) is assumed to take the following form:

The macroscopic strength profiles of porous MC materials with a fixed porosity (f = 0.1) and different values of friction angle are investigated again in Fig.7.Satisfyingly,the poor predictions of the parametric strength criterion in the quadrant(Σm<0,>0)(see Fig. (6)) can be well corrected by the heuristic strength criterion.The effect of β on the macroscopic strength of porous media is illustrated in Figs. 8-10. It appears that the strength profile of porous material is strongly shrank in the compressive domain when β decreases.For matrix with distinct values of eBCS and UCS,neglecting the difference between eBCS and UCS would underestimate the macroscopic strength, especially when the applied pressure is large.

Fig.7. Comparison of yield profiles predicted by the heuristic strength criterion(dash line)with FE solution(points)at different values of φ in the(Σm/c,(Σ11-Σ33)/c)space(β = 1, f = 0.1).

8. Conclusions

In this work,we proposed macroscopic yield criteria for porous media with matrix described by the three-parameter strength criterion.The hollow sphere model subjected to axisymmetric loading is considered. The maximum plastic dissipation function and plastic admissibility conditions for the three-parameter strength criterion are derived.To construct the trial velocity field,the exact limit stresses and velocity fields in the hollow sphere are analyzed under isotropic loading.Based on these work,an exact upper bound of the macroscopic strength of porous medua under purely hydrostatic loading is obtained. Besides, an estimate of macroscopic strength of porous media subjected to axisymmetric loading is proposed in parametric form through kinematic limit analysis with a relaxed plastic admissibility condition. The relaxed constraint is prone to broaden the admissible trial velocity space and thus results in an underestimation of the macroscopic strength.To be able to fully recover the exact upper bound for purely hydrostatic loading, the established parametric criterion is further modified.Moreover,for practical applications,a heuristic strength criterion is also developed in explicit form by examining the limit cases of the improved parametric strength criterion. The dependence of the macroscopic strength of porous media on the Lode angle is reflected through the sign of the third invariant. Numerical assessment of the proposed strength criteria is performed. The macroscopic strength given by the improved parametric strength criterion is shown to be in a good agreement with the numerical solution as the porous material is subjected to triaxial compression or the mean part of the applied axisymmetric stress is positive,whereas the heuristic strength criterion could provide satisfactory predictions of the macroscopic strength except for the macroscopic stress states of triaxial compression with large pressure. The present work also highlights the ratio of eBCS to UCS significantly impacts the macroscopic yield profile. Neglecting the difference between eBCS and UCS would underestimate the macroscopic strength especially when the pressure is large, such as the hydrostatic compressive strength.

Fig.8. Comparison of yield profiles predicted by the heuristic strength criterion(dash line) with FE solution (points) at different values of β in the (Σm/σt, (Σ11 - Σ33)/σt)space (α = 0.33, f = 0.1).

Fig.9. Comparison of yield profiles predicted by the heuristic strength criterion(dash line) with FE solution (points) at different values of β in the (Σm/σt, (Σ11 - Σ33)/σt)space (α = 0.33, f = 0.15).

Fig.10. Comparison of yield profiles predicted by the heuristic strength criterion(dash line) with FE solution (points) at different values of β in the (Σm/σt, (Σ11 - Σ33)/σt)space (α = 0.33, f = 0.2).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This research was partially supported by the National Natural Science Foundation of China(Grant No.51804203)and this support is gratefully appreciated.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2021.03.013.