A micro-mechanical constitutive model for heterogeneous rocks with non-associated plastic matrix as implicit standard materials

Abstract

In this work, we shall propose a new micro-mechanical constitutive model for the estimation of effective elastic-plastic behaviors of heterogeneous rocks. A bi-potential based incremental variational (BIV) approach is developed in order to take into account non-uniform local strain fields of constituents. The studied materials are composed of a non-associated and pressure sensitive plastic matrix, elastic inclusions and/or voids. For clarity, the local behavior of matrix is first described by an elastic perfectly-plastic model. Based on the bi-potential theory to dealing with non-associated plastic flow, the solid matrix is considered as pertaining to implicit standard materials (ISMs). The effective incremental bi-potential and macroscopic stress tensor are then estimated through an extension of the incremental variational method initially established for generalized standard materials(GSMs). The accuracy of the BIV model is verified by comparing the model’s predictions with the reference results obtained from direct finite element simulations. Furthermore, by assuming that the general formulation obtained for the perfectly plastic matrix remains valid for each loading increment, the BIV model is extended to considering that the solid matrix exhibits an isotropic hardening by using an explicit algorithm. The accuracy of the extended BIV model is also validated by a series of comparisons with the reference solutions obtained by direct finite element simulations for both inclusion-reinforced composites and porous materials. Both local and macroscopic responses are compared. As an example of application, the extended BIV model is finally applied to estimating the mechanical responses of typical claystone and sandstone under different loading paths.

Introduction

Rocks are usually regarded as typical composites, which are used in a very wide range of engineering constructions. These materials contain different kinds of heterogeneities at different scales. Pores and inclusions are two main families of heterogeneities. Furthermore, these materials are composed of several mineral phases of different properties. The mineral compositions may significantly vary in space, for instance with geological depth. Laboratory studies have shown that the macroscopic physical and mechanical properties of these materials are affected by heterogeneities and mineral compositions. So far, different kinds of macroscopic models, mainly elastic-plastic and damage models have been developed. Directly fitted from laboratory tests, these models are able to correctly reproduce the main features of mechanical behaviors of those materials. However, they are not able to properly consider the effect of heterogeneities and mineralogical compositions on the macroscopic mechanical responses.

Based on linear homogenization techniques, micro-mechanical models have first been established during the last decades for modeling induced damage in brittle rocks (Zhu et al., 2008, Zhu et al., 2016, Zhao et al., 2018, Zhang et al., 2019). Important advances have also been obtained on micro-mechanical modeling of plastic deformation in ductile and porous rocks by using nonlinear homogenization methods. For instance, clayey rocks have been characterized as composites constituted of a plastic clay matrix in which calcite and quartz grains are embedded (Guéry et al., 2008, Jiang et al., 2009). In some multi-scale models, the microstructure of clayey rocks has further been enriched by considering the clay matrix as a porous material at the microscopic scale (Shen et al., 2012). The effective inelastic behavior of the porous clay matrix has been estimated by using the Hill incremental method (Hill, 1965). As for metallic composite materials, it was found that the use of the original Hill’s incremental method produced too stiff mechanical behaviors (Suquet, 1996, Chaboche et al., 2005). The main reason is the fact that uniform local strain fields are assumed in constituents of composites in the Hill’s method. In order to improve the numerical performance of this method, artificial techniques, such as isotropization of tangent elastic-plastic stiffness tensor, have been proposed. This correction technique has also been applied to clayey rocks (Guéry et al., 2008, Jiang and Shao, 2009, Shen et al., 2012). However, all those correction techniques are generally not based on any physical background.

Meanwhile, advanced nonlinear homogenization techniques have been developed for composite materials considering non-uniform local fields in constituent phases (Castañeda, 1991, Castañeda, 1992, Castañeda, 2002, Lahellec and Suquet, 2007a, Lahellec and Suquet, 2007b, Lahellec and Suquet, 2013, Boudet et al., 2016, Brassart et al., 2011, Brassart et al., 2012, Danas and Castañeda, 2012), just to mention some representative ones. In particular, variational principles based on the use of a ”linear comparison composite (LCC)” were proposed for the mean field homogenization method of nonlinear elastic composites (Castañeda, 1991, Castañeda, 1992, Castañeda, 2002), and used to generate improved bounds and more generate estimates for the nonlinear elastic-plastic composites (Castañeda and Suquet, 1997, Danas and Castañeda, 2012). By extending these previous works, a new incremental variational method has been established (Lahellec and Suquet, 2007a, Lahellec and Suquet, 2007b) for modeling effective nonlinear properties of viscoelastic composites without local threshold or hardening. In this new method, equivalent interval variables (EIV) are introduced to capture the non-uniform local plastic strain fields. Further, the same authors have proposed a rate variational model (RVP) by considering a non-uniform field of plastic strain rate (Lahellec and Suquet, 2013). More recently, the EIV method has been extended to modeling elastic-(visco) plastic composites with local threshold and isotropic and/or linear kinematic hardening (Boudet et al., 2016). On the other hand, based on the variational principle established in Ortiz and Stainier (1999), alternative incremental variational models have been proposed in Brassart et al., 2011, Brassart et al., 2012 for studying elastic-(visco) plastic composites with local isotropic hardening. The EIV method has further been extended to the description of geological materials with a pressure-dependent Drucker-Prager plastic matrix (Zhao et al., 2019). However, all these previous models have been developed in the scope of Generalized Standard Materials (GSMs) (Halphen and Nguyen, 1975) with an associated plastic flow rule.

Extensive experimental results have clearly shown that for most rocks, a non-associated plastic flow rule is required for correctly modeling the coupling between shear and volumetric strains. These materials cannot be considered as Generalized Standard Materials. As a first approximation, the microstructure of these materials at a selected relevant length scale, for instance micrometer, can be characterized by the representative unit cell shown in Fig. 1. Several sets of elastic inclusions (mineral grains in rocks) are embedded in a plastic matrix (clay matrix in clayey rocks). The local behavior of the matrix is generally described by a non-associated and pressure sensitive plastic model. The incremental variational methods developed for the GSMs cannot be directly used to estimate the effective mechanical behaviors of rocks.

In order to generalize the incremental variational principles to heterogeneous rocks, the idea here is first to transform these non-GSMs into a class of implicit standard materials (ISMs). This is done with the help of the bi-potential theory initially developed for macroscopic elastic and plastic behaviors of non-GSMs (De Saxcé and Feng, 1991, De Saxcé and Bousshine, 1998, De Saxcé, 1995). This theory has been successfully used for modeling soils and rock-like materials with non-associated plastic models (Bodovillé and De Saxcé, 2001, Bodovillé, 2001, Hjiaj et al., 2003, Berga, 2012). Moreover, the bi-potential theory is naturally suitable for developing a variational approach of constitutive modeling.

With the help of the bi-potential theory, the aim of this work is to develop a new incremental variational method for estimating the effective elastic-plastic behavior of heterogeneous rocks composed of a non-associated and pressure sensitive plastic matrix. This is based on the construction of an incremental elastic-plastic bi-potential for ISMs by using an implicit time-discretization scheme.

On the other hand, ductile and porous heterogeneous rocks generally exhibit plastic hardening. For instance, in the case of an isotropic hardening, the internal friction or cohesion can evolve during plastic deformation. In the case of heterogeneous rocks idealized in Fig. 1, plastic hardening occurs in the plastic matrix. This mechanism should be taken into account. However, the formulation of an incremental variational model for materials with a pressure-sensitive plastic matrix with plastic hardening may becomes mathematically very complex. By taking the incremental nature of the approach, a simplified explicit method is proposed in this paper. The new bi-potential base incremental variational model (BIV) is first developed by considering a perfectly plastic matrix. Then at the end of each loading increment, the plastic properties are updated but frozen for next loading increment. The plastic matrix is then considered as a material without hardening during the current increment.

The proposed new BIV model is validated by comparing model’s predictions and numerical results issued from direct finite element simulations for both perfectly plastic and plastic with hardening cases. Finally, the new BIV model is applied to estimating the effective mechanical responses of typical claystone and porous sandstone in various loading paths.

Throughout this paper, the following notions of tensorial products of any second order tensors $\mathit{A}$ and $\mathit{B}$ will be used: ${\left(\mathit{A}\otimes \mathit{B}\right)}_{\mathit{ijkl}}=A_{\mathit{ij}}{B}_{\mathit{kl}}$ and $\mathit{A}:\mathit{B}=A_{\mathit{ij}}{B}_{\mathit{ij}}$. Fourth order tensors are denoted by blackboard bold characters, and one can define ${\left(\mathbb{C}:\mathit{B}\right)}_{\mathit{kl}}=C_{\mathit{ijkl}}{B}_{\mathit{kl}}$. The symbol $‖\mathit{A}‖=\sqrt{\mathit{A}:\mathit{A}}$ is used to denote the norm of any second order tensor $\mathit{A}$. With the second order identity tensor $\mathit{\delta }$, usually used fourth order isotropic identity tensor $\mathbb{I}$ and fourth order hydrostatic projects $\mathbb{J}$ are expressed in the components form as $I_{\mathit{ijkl}}=\frac{1}{2}\left({\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}}+{\delta }_{\mathit{il}}{\delta }_{\mathit{jk}}\right)$ and $J_{\mathit{ijkl}}=\frac{1}{3}{\delta }_{\mathit{ij}}{\delta }_{\mathit{kl}}$, respectively. The fourth order deviatoric projects $\mathbb{K}=\mathbb{I}-\mathbb{J}$ is then obtained. Moreover, the fourth-order tensors $\mathbb{J}$ and $\mathbb{K}$ have the properties: $\mathbb{J}:\mathbb{J}=\mathbb{J},\mathbb{K}:\mathbb{K}=\mathbb{K},\mathbb{J}:\mathbb{K}=\mathbb{K}:$ $\mathbb{J}=\mathbf{0}\text{.}$.