Numerical analysis of the indentation of a poro-elastoplastic half space by a rigid cylinder

Poro-elasto-plastic cylindrical indentation as a means for material characterization is explored numerically in this work. When a rigid cylinder is being pressed into a semi-infinite, fully saturated poroelastic medium via step displacement loading, undrained and drained elastic constants are reflected in the responses at the early and late times, while hydraulic diffusivity is manifested in the transient behavior. Such an indentation test can therefore be used in principle to identify material parameters. For geomaterials, plastic deformation may occur due to the indentation action. How the hydromechanically coupled indentation response is affected by plastic deformation is investigated here using finite element analysis. We consider the problem to be plane strain with frictionless contact. Drucker-Prager failure criteria with and without a cap are considered. We show that if there is only limited plastic deformation, the concept of using the transient response for determining hydraulic diffusivity may be extended from poroelasticity to poro-elasto-plasticity.


Introduction
Poroelastic material characterization is critical for hydromechanically coupled problems such as reservoir production for hydrocarbon recovery and fluid injection-induced seismicity, among many others.In general, an ensemble of undrained, drained, and unjacketed tests based on multiaxial setups is needed to identify the full set of material parameters required for linear poroelasticity analysis.These tests are particularly time-consuming for low-permeability materials.
Alternatives of using contact mechanics-based approaches such as spherical and cylindrical indentation have been suggested in the literature [1,2] and realized in laboratory experiments for soft materials such as polymeric gels [3 − 5].Feasibility of these methods lies in the fact that undrained and drained elastic constants of the porous medium are reflected in the responses at the early and late times, while hydraulic diffusivity or the coefficient of consolidation is manifested in the transient behavior only.In the indentation process, the characteristic length scale is associated with the contact length or radius.Therefore, both size of the indenter tip and the depth of penetration could be used to adjust the characteristic time for pore pressure diffusion.
For poroelastic cylindrical indentation, Hui et al. [1] derived a theoretical solution using the method of Muskhelishvili [5] assuming that the rigid indenter is subjected to a Heaviside step displacement and the constituents of the poroelastic half space are incompressible.The indentation action is modeled by prescribing the normal displacement over the contact area, the size of which remains constant over time.Meanwhile, comprehensive theoretical studies on poroelastic spherical indentation have recently been conducted by the co-authors [7 -10].Master curves that can be readily used for experimental interpretation have been rigorously constructed for various loading and surface drainage conditions.
For geomaterials, plastic deformation may occur in the indentation process.It is therefore necessary to understand how the indentation response will be affected by plastic deformation.It is shown in [8] using a Drucker-Prager failure criterion for the case of spherical indentation with step displacement loading that if the material is of sufficient strength such that plastic yielding occurs only at the undrained limit, the effect of plastic deformation on the normalized indentation force response in the transient phase may be neglected.This conclusion makes it feasible to use the spherical indentation test to determine hydraulic diffusivity for geomaterials.
In the present work, we examine the indentation of a poro-elasto-plastic half space by a rigid cylinder using the finite element method in ABAQUS/CAE 2018.The particular case when the indenter is subjected to a step displacement loading is investigated.The problem is assumed to be plane strain and the surface is assumed to be impermeable under the indenter, but permeable everywhere else.Frictionless contact between the rigid indenter and the porous medium is explicitly modeled.Drucker-Prager criterion and the modified Drucker-Prager/Cap criterion are used to examine the effects of both shear and volumetric plastic yielding.

Theoretical consideration
For cylindrical indentation in a linear elastic material, the relation for the nominal half length of contact, ܽ, depth of penetration, ݀, and indenter radius, ܴ, is purely geometrical [10], The elastic solution for the indentation force [10] is applicable for poroelastic indentation at the undrained ‫ݐ(‬ = 0) and drained ‫ݐ(‬ → ∞) limits if appropriate elastic constants are used, where ‫ܩ‬ and ‫ݒ‬ are the shear modulus and the drained Poisson's ratio and ߶ is a poroelastic parameter that is a function of Biot's coefficient ߙ, storage coefficient ܵ, and the drained bulk and shear moduli ‫ܭ‬ and ‫.ܩ‬For a general linear poroelastic material with compressible constituents, the undrained Poisson's ratio can be expressed as ‫ݒ‬ ௨ = 1 − ߶/2.
A workflow to determine the hydraulic diffusivity can be constructed as follows.After obtaining measurement of the indentation force as a function of time, we can use the early and late time force asymptotes to determine ߱ and to normalize the indentation force.Hydraulic diffusivity ܿ ௩ can be determined by matching the experimental data against an appropriate solution of ‫ܨ‬ ‫ݐ(‬ * ) for a given ߱.

Numerical model setup
A numerical model of width ‫ݓ‬ = 40 mm and height ℎ = 50 mm is used to represent half of the domain.The left side of the model ‫ݔ(‬ = 0) represents the contact axis.Roller support is applied at ‫ݔ‬ = 0 and ‫ݖ‬ = ℎ, i.e., the bottom of the domain, while the lateral boundary at ‫ݔ‬ = ‫ݓ‬ is traction free.The drainage boundary condition is set as impermeable at ‫ݔ‬ = 0, zero pore pressure at ‫ݔ‬ = ‫ݓ‬ and ‫ݖ‬ = ℎ .The surface ‫ݖ(‬ = 0) is impermeable within the contact area (0 ≤ ‫ݔ‬ ≤ ܽ) , and fully permeable everywhere else ‫ݔ(‬ > ܽ).
The mesh is constructed using eight-node plane strain quadrilateral, biquadratic displacement, and bilinear pore pressure elements (CPE8RP), which satisfies the so-called LBB (LadyzhenskayaBabuška-Brezzi) condition and does not have numerical issues of pore pressure oscillation at the undrained limit [8].Frictionless contact between the indenter and the domain is explicitly modeled.A cylindrical indenter of radius ܴ = 10 mm at a depth of penetration ݀ = 0.1 mm corresponds to a nominal contact half length ܽ = 0.76 mm .The following material properties based on the Gulf of Mexico shale [11] are used: Young's modulus ‫ܧ‬ = 1.853GPa,Poisson's ratio ‫ݒ‬ = 0.22, porosity ݊ = 0.3 and permeability ߢ = 1 × 10 ିଵଽ m ଶ .We assume that fluid viscosity ߤ = 10 ିଷ Pa ⋅ s and the solid and fluid phases are incompressible.This set of material parameters corresponds to ߱ = 0.56.

Drucker-Prager failure criterion
A linear Drucker-Prager (D-P) failure criterion can be expressed as, ݂ = ඥ3‫ܬ‬ ଶ − ‫ܫ‬ ଵ tan ߚ/3 − ݇ = 0, where ‫ܫ‬ ଵ and ‫ܬ‬ ଶ are the stress invariants and ݇ and ߚ are the yield stress and the friction angle, respectively.Note that ‫ݍ‬ = ඥ3‫ܬ‬ ଶ is the so-called Mises equivalent stress.For geomaterials, since material properties are generally obtained in the laboratory based on a Mohr-Coulomb (M-C) criterion, cohesion and friction angle, ܿ and ߮, of a M − C criterion are used as inputs here.Assuming the M-C yield surface inscribes the D-P failure envelope, we have [12], The corresponding plastic potential can be expressed as ݃ = ඥ3‫ܬ‬ ଶ − ‫ܫ‬ ଵ sin ߰/ඥ3 + sin ଶ ߰, where ߰ is the dilatancy angle (0 ≤ ߰ ≤ ߮).The flow rule is associative when ߮ = ߰.
It follows from the poroelastic simulation that the critical cohesion below which incipient plastic failure occurs is ܿ = 30.5MPaif we assume ߮ = ߰ = 20 ∘ .Therefore, we choose three representative cases, ܿ = 17,20 and 30MPa, to analyze the effect of plastic deformation.
The numerical results show that at ܿ = 30MPa, since plastic deformation is very limited, the overall response is nearly identical to the poroelastic case, as evidenced by a comparison in Figure 1 for the effective horizontal stress ߪ ௫ distributions along the contact axis between the poro-elastoplastic cases with ܿ = 17,20,30MPa and the poroelastic case at both the undrained and drained limits.Tensile stress is positive here.The points where the stress distribution is not smooth mark the elasto-plastic interface.Effect of plastic deformation becomes notable at ܿ = 17 and 20MPa.At the undrained limit, plastic deformation reduces the tensile horizontal stress within the plastic zone.As a result of stress redistribution, the maximum tensile stress now occurs at the lower boundary of the plastic zone at both the undrained and drained limits.These three cases with ܿ = 17,20,30MPa share a few common characteristics in the evolution of plastic deformation.Plastic yielding occurs instantaneously at the undrained limit and there is no further accumulation of plastic deformation during the transient phase.It can be observed from the contours of the equivalent plastic strain in Figure 2 that the plastic zone is contained beneath the contact area.Size and shape of the plastic zone as well as the location and magnitude of the maximum equivalent plastic strain remain unchanged over time.
As cohesion decreases from 20 to 17MPa, plastic deformation intensifies and size of the plastic zone increases as expected.At ܿ = 20MPa, the maximum equivalent plastic strain is ߳ ୫ୟ୶ = 0.52% and is located on the contact axis.At ܿ = 17MPa, ߳ ୫ୟ୶ = 1.188%, more than twice that of ܿ = 20MPa.However, its location is now away from the contact axis.
Additional simulations with various cohesion values allow us to categorize the plastic deformation mechanism into three distinct types.When 30.5 > ܿ ≥ 19.7MPa, the plastic zone is contained within the matrix undergoing elastic deformation and the location of ߳ ୫ୟ୶ is at the contact axis.In contrast, when 19.7 > ܿ ≥ 17MPa, while the plastic zone is still contained, the location of ߳ ୫ୟ୶ has moved to the side, away from the contact axis.When ܿ < 17MPa, plastic yielding now occurs in the transient phase and the plastic zone could extend to the free surface and become uncontained.As a result, numerical convergence becomes a rather challenging issue.As far as the pore pressure distribution is concerned, the normalized pore pressure contours at ‫ݐ‬ * = 0,0.1,1and 10 for ܿ = 20MPa, see Figure 3, is very similar to those from ܿ = 30MPa and the poroelastic case, even though the maximum and minimum pore pressure values are different.
At early time, pore pressure is the largest right underneath the indenter.Overtime, location of the maximum pore pressure moves downwards along the contact axis as a result of surface drainage.At ‫ݐ‬ * = 10, the pore pressure contours show clear evidence of the drainage boundary effect on the right side.At all times, there is a rather small negative pore pressure zone near the edge of the contact area due to strong shear induced dilatancy.As cohesion decreases to ܿ = 17MPa, the contours at ‫ݐ‬ * = 0 are clearly affected by the expanded size of the plastic zone, see Figure 4.However, shape of contour plots at ‫ݐ‬ * = 0.1,1 and 10 does not exhibit marked differences from those of ܿ = 20MPa.Increases in the plastic zone size and in the magnitude of the plastic strain result in a stronger effect of shear induced dilatancy.Consequently, magnitude of the maximum pore pressure at the same dimensionless time is smaller at ܿ = 17MPa.
The most important finding from this set of simulations is that the normalized force relaxation response is only slightly affected by the plastic deformation, even though the magnitude of the indentation force at the early and late time has changed, see Figure 5(a).It can be seen that ‫ܨ‬ ‫ݐ(‬ * ) relaxes only slightly slower when cohesion decreases from 30 to 17MPa with the discrepancy mostly within ‫ݐ‬ * ∼ 1 − 100.Note that the simulation is terminated at ‫ݐ‬ * = 10 ସ .The indentation forces at ‫ݐ‬ * = 0 and 10 ସ are chosen as the two force asymptotes.

Modified Drucker-Prager/Cap model
Since the indentation action could result in large compressive stress underneath the contact area, we use the modified Drucker-Prager/cap model here to investigate the potential effect of compaction failure.The yield surface now has two additional components, a cap and a transition surface that connects the cap smoothly to the Drucker-Prager surface.Parameters that describe these two surfaces are: ‫‬ -the hydrostatic compression yield stress, ܴ -the cap eccentricity, and ߙ -a parameter for the transition surface.Interested readers are referred to the ABAQUS manual for detailed descriptions.In this study, we choose ܴ = 0.2 and ߙ = 0.01.In this model, the flow rule is associative only in the cap and the transition part of the yield surface.For the shear failure surface, the flow rule is defined by a plastic potential similar to the transition yield function.As a result, even if we have the same shear strength parameters, results from the cap model will be different from those based on the linear Drucker-Prager model.
In the cap hardening model, the evolution of the hydrostatic yield stress ‫‬ is expressed as an adhoc exponential growth function of the plastic volumetric strain Δߝ ௩ . Two cap models are employed in this study, see Table 1.The volumetric plastic strain Δߝ ௩ from the cap model is an increment to the state when onset of compaction failure occurs.Between the two models, ‫‬ is smaller in Cap 2 and its magnitude is comparable to the critical cohesion value ܿ = 30.5MPa.
We should point out that the choice of Cap 2 model, though somewhat arbitrary, is however not unreasonable.
Plastic deformation from the two cap models with shear strength parameters kept as ܿ = 20 MPa and ߮ = 20 ∘ are shown in Figure 6.In theory, plastic strain could have contributions from all three components of the yield surface.However, with Cap 1 model, there is no plastic strain resulted from the transition and the cap surface and plastic deformation due to shear failure occurs only at the undrained limit.Denote the equivalent plastic strain from the Drucker-Prager yield surface as ߳ .In this case, characteristics of plastic deformation are the same as the one without the cap, see Figure 6(a).However, due to the non-associative flow rule, the size of the plastic zone and the magnitude of the equivalent plastic strain become much larger, ߳ ୫ୟ୶ = 0.966%.When Cap 2 model is used, shear plastic deformation still occurs only at the undrained limit, but there is plastic strain from the transition surface at ‫ݐ‬ * > 0.1.Compaction failure occurs right underneath the contact area, see Figures 6(c) -(d).The maximum values are ߳ ௧୫ୟ୶ = 0.112% at ‫ݐ‬ * = 1 and ߳ tmax = 0.167% at ‫ݐ‬ * = 10, an order of magnitude smaller than that due to shear failure, ߳ ୫ୟ୶ = 1.083%, located now away from the contact axis, see 6(b).Compared with the case without a cap, effect of plastic deformation on the pore pressure distribution is relatively minor.The maximum and minimum values from the two cap models at the four dimensionless times ‫ݐ‬ * = 0,0.1,1and 10 are: Cap 1, ‫‬ ୫ୟ୶ = 77.1146,39.9892, 10.7083, 2.6196MPa, and ‫‬ ୫୧୬ = −9.5185,−0.1569, −0.0758, −0.0113MPa; Cap 2, ‫‬ ୫ୟ୶ = 74.4157,41.2200,10.9975,2.5922MPa, and ‫‬ ୫୧୬ = −9.5709,−0.1538, −0.0745, −0.0111MPa .Decrease in the hydrostatic yield stress in Cap 2 causes the maximum pore pressures to decrease slightly.The most notable difference between the results with and without a cap is the rather large negative pore pressure, ‫‬ ୫୧୬ , at ‫ݐ‬ * = 0 in the cap models.The negative pore pressure occurs near the contact edge.However, the size of this negative pore pressure zone remains small and the magnitude decreases quickly with time.
A comparison of ‫ܨ‬ ‫ݐ(‬ * ) from the three simulations with ܿ = 20MPa is shown in Figure 5 The results suggest that for this particular cohesion value, the cap models do not have a strong influence on the force relaxation response.Differences from the three simulations are minor.

Conclusion
A finite element analysis is conducted to examine cylindrical indentation in a poro-elasto-plastic medium with step displacement loading.Effects of cohesion and plastic yielding due to compaction failure are investigated using the Drucker-Prager failure criterion and a modified Drucker-Prager/Cap model.The numerical results suggest that if plastic deformation is limited, namely, if shear plastic failure occurs only at the undrained limit and is contained, the normalized transient force response is only slightly affected by the plastic deformation and can therefore be approximated by a poroelastic solution.As such, the concept of using the transient response to determine hydraulic diffusivity for geomaterials may be extended from poroelasticity to poroelastoplasticity.The findings here are consistent with our previous work on spherical indentation [9].However, since failure due to tensile crack growth is not considered here and the analysis on the effect of cap model is not exhaustive, further work will be needed in order to evaluate the applicability of this method for a particular material.

Figure 1 .Figure 2 .
Figure 1.Comparison of the effective horizontal stress ߪ ௫ distributions along the contact axis at the undrained and drained limits between the poro-elasto-plastic cases with ܿ = 17,20,30MPaand the poroelastic case.

Table 1 .
Two cap hardening models