(PDF) An enriched constitutive modelling framework for localised failure of geomaterials

Nguyen, G. (2014). An enriched constitutive modelling framework for localised failure of Nguyen, G. (2014). Obras y Proyectos 15, 33-39 geomaterials. Obras y Proyectos 14, 33-39 An enriched constitutive modelling framework for localised failure of geomaterials Marco para un modelo constitutivo extendido para geomateriales con fallas localizadas Fecha de entrega: 6 de diciembre 2013 Fecha de aceptación: 5 de mayo 2014 Giang D. Nguyen School of Civil, Environmental and Mining Engineering, North Terrace Campus, The University of Adelaide, SA 5005 Australia,

Localised failure in geomaterials is preceded and Fallas localizadas en geomateriales están precedidas y accompanied by intensive deformation and irreversible se dan en conjunto con deformaciones intensas y cambios micro-structural changes of the material in a small but micro estructurales irreversibles del material en una región finite size region. Shear, compaction, and dilation bands de tamaño finito, pero pequeña. Corte, compactación y observed in soils and porous rocks are typical examples bandas de dilatación observables en suelos y rocas porosas of phenomena that lead to localised failure. The width h son ejemplos típicos de fenómenos que conducen a fallas of the localisation band has been experimentally shown localizadas. Ha sido experimentalmente demostrado que el to be a physical quantity related to the microstructure ancho h de la banda de localización es una cantidad física of the material. On the other hand, numerical methods relacionada con la micro estructura del material. Por otro for the solution of boundary value problems usually lado, métodos numéricos para la resolución de problemas de introduce another length scale H, as a result of the valor en la frontera usualmente introducen otra longitud de spatial discretisation of the considered domain into escala H como un resultado de la discretización espacial del smaller ones over which the constitutive response of dominio considerado en partes más pequeñas sobre el cual la the material is defined in terms of incremental stress- respuesta constitutiva del material está definida en términos strain relationships. While h, as a physical quantity, de relaciones incrementales de tensión-deformación. is fixed, H varies with the resolution of the numerical Mientras h como cantidad física está fija, H varía con discretisation. Since h scales with the material la resolución de la discretización numérica. Dado que h microstructure and therefore is usually much smaller escala con la micro estructura del material y por lo tanto than the resolution of the numerical discretisation, es usualmente mucho más pequeño que la resolución de la the case H > h is considered in this study, e.g. failure discretización numérica, el caso H > h es considerado en este behaviour governed by a localisation band of width estudio, por ejemplo el comportamiento en falla gobernado h embedded in an elastic bulk of nominal side H. We por una banda de localización de ancho h inserta en una present a general constitutive modelling framework to masa elástica de lado nominal H. Se presenta un marco de connect these two scales, and corresponding responses modelamiento constitutivo general para conectar estas dos of the materials inside and outside the localisation zone. escalas, y las respuestas correspondientes del material dentro We demonstrate how this approach can help obtain y fuera de la zona de localización. Se demuestra como esta physically meaningful solutions that are independent estrategia puede ayudar a obtener soluciones con significado of the spatial discretisation in numerical analysis. físico que son independientes de la discretización espacial en Numerical analyses of localised failure in quasi-brittle análisis numéricos. Análisis numéricos de falla localizada en materials are used to further highlight the features and materiales cuasi frágiles son además usados para destacar applicability of the proposed approach. las características y aplicabilidad de la estrategia propuesta. Keywords: length scales, constitutive modelling, Palabras clave: escalas de longitud, modelo constitutivo, localised failure, discontinuity, bifurcation, damage, falla localizada, discontinuidad, bifurcación, daño, energía fracture energy de fracturación 33 Nguyen, G. (2014). An enriched constitutive modelling framework for localised failure of geomaterials. Obras y Proyectos 15, 33-39 Sluys, 2001; Samaniego and Belytschko, 2005; Sanborn Introduction and Prévost, 2011). These more sophisticated methods Failure of geomaterials such as the formation of sea ice have been extensively used to address failure modelling leads in the Arctic (Jirasek and Bazant, 1995; Sulsky et al., at large scales, e.g. H > h, and usually idealise the finite 2007), rock fracture in underground mining, etc., usually width localisation zone as a zero thickness surface involves the material behaviour at various scales and stages across which the displacement field is discontinuous (the of deformation. In particular, inelastic deformation and strong discontinuity case). We are also aware of earlier fracture localise in narrow zones, while the surrounding works on weak discontinuity (e.g. only the strain field is bulk, usually of several orders of magnitude larger in extent, discontinuous) using the localisation zone embedded in is unloading elastically. Numerical modelling of such large finite elements (e.g. Belytschko et al., 1988; Sluys and scale failure processes (dimensions of several kilometres) Berends, 1998; and Garikipati and Hughes, 2000). The is computationally challenging, as the behaviour of the key idea of these approaches is to enhance the deformation material inside and outside the localisation zone should be mode of the special finite element so that both inelastic incorporated in a numerical model. In this regard, classical behaviour in the localisation zone and the elastic shrinking continuum models usually lack a length scale to correctly of the bulk continuum can be adequately accounted for. As capture the localised failure of softening materials. While a consequence, all such approaches involve finite element enrichment of such models with an internal length through re-formulations, e.g. modification and/or introduction of the use of nonlocal/gradient regularisation (Pijaudier- shape functions, and hence result in the dependence of the Cabot and Bazant, 1987; Chen and Schreyer, 1987) is a approach on the type of finite element used for the numerical mathematically and probably physically rigorous way, the discretisation. Can we come up with a methodology better application of such enhancements is severely restricted than the smeared crack approach, but less complex than by the available computational resources. This is because the EAS or XFEM? the locations of the failure zones are generally unknown and the considered domain can be of several orders of The key idea described in the present paper is to enhance magnitude bigger than the characteristic width h of the the constitutive behaviour description, rather than the localisation zone, resulting in a very high resolution of the finite element formulation, with a length scale related discretisation, and consequently very large model size. to the width of the localisation zone used for large scale analysis of failure, e.g. h < H. We address the nature of the In the literature, fracture energy regularisation based on localised failure from the constitutive modelling point of the smeared crack concept (Cedolin and Bazant, 1980) is view without having to resort to variational formulations probably the simplest way to cope with softening-related for the discretisation using the finite element method. The issues in the failure analysis of solids/structures. However, interface with any spatial discretisation scheme is taken into this suffers from the drawback that the constitutive account only through the size h of the sub-domain (Figure behaviour must be unphysically scaled with the resolution 1), while connection with any kind of constitutive behaviour of the discretisation to meet the requirement on the energy is specified through its tangent stiffness. The proposed dissipation. In addition, once a coarse spatial discretisation approach is therefore straightforwardly applicable to any is used (e.g. in large scale analysis), the specific fracture material model and any spatial discretisation scheme. energy obtained from the (surface) fracture energy and the resolution of the discretisation (see equation (8) below) becomes smaller than the elastic strain energy at peak (Jirasek and Bazant, 1995), resulting in the inadmissible snapback instability in the constitutive response. Other enhancements to the discretisation scheme include Enhanced Assumed Strain EAS (e.g. Larsson et al., 1996; Oliver, 1996; Borja, 2000; Foster et al., 2007) and XFEM, Figure 1: Localisation zone (shaded) embedded in a volume the eXtended Finite Element Methods (e.g. Wells and (Nguyen et al., 2012) 34 Nguyen, G. (2014). Obras y Proyectos 15, 33-39 Constitutive modelling framework We aim to develop at first a general framework that can fit any constitutive model. The starting point is the localisation the relaxation strain rate of the elastic bulk is: zone of size h embedded in a volume of nominal size H, introduced in relation to the Fracture Process Zone FPZ observed in the failure of quasi-brittle materials, or the shear band in soils. This is the stage beyond the homogeneous deformation of the material (Figure 1), e.g. at the onset or after the bifurcation of the material behaviour (Borja, For elastic unloading of the bulk material with stiffness 2000). It is assumed in this framework that, during further modulus ao, the stress rate of the continuum model is deformation, the material undergoes elastic unloading in (Nguyen et al., 2012): the region outside the localisation zone. While h is related to the inelastic behaviour and failure at the lower scale, the size H of the elastic bulk is merely a numerical feature resulting from the discretisation of the domain under Inside the localisation zone, the constitutive relationship consideration. Grid spacing in the finite difference method rate is assumed of the general form: and the domain size ascribed to an integration point in the finite element method are typical examples of the meaning we attach to the discretisation characteristic size H. As a consequence, H can vary depending on the required where ai is the tangent stiffness. The above two equations resolution of the numerical discretisation, while h can be are linked through the internal equilibrium dictating considered a material property and hence is invariant with the continuity of traction across the boundary of the the numerical discretisation. We view the configuration at localisation zone: failure in Figure 1 as a composite material consisting of two phases: an inelastic localisation zone embedded in an elastic bulk. For this, the linear scaling of the total strain The distinction from a regular constitutive model lies rate applies: in the coupled equations (4-6). From (4), in order to obtain the stress rate from a given strain rate , the constitutive equation (5) dictating the inelastic response of where f = h/H is the volume fraction of the embedded the localisation band, together with the internal equilibrium localisation zone in Ω and subscripts “i” and “o” are used equation (6), are needed. The volume fraction f and the to denote quantities belonging to the behaviour inside inelastic response inside the localisation zone contribute to and outside the localisation zone, respectively (Figure 1). the relaxation of the stress rate in the elastic bulk. Equations Across the boundary of the localisation zone, the internal (4-5) can be worked out to result in the following form of equilibrium in terms of traction continuity must be met, stress-strain relationship (Nguyen et al., 2012): enforcing coupling between the inelastic localisation zone and outside elastic bulk. In other words, the elastic bulk is enhanced by the inelastic behaviour of the localisation zone (Nguyen et al., 2012). Denoting n the normal vector, where is a tensor obtained [ů] the relative velocity between opposite sides of the from the sizes and acoustic tensors of the behaviours localisation band, and adopting the following form for the inside and outside the localisation band. As can be seen, strain rate inside the localisation band (Vardoulakis et al., the overall constitutive behaviour in this case takes into 1978; Kolymbas, 2009; neglecting the small homogeneous account both the responses and sizes of inelastic and elastic term in ): 35 Nguyen, G. (2014). An enriched constitutive modelling framework for localised failure of geomaterials. Obras y Proyectos 15, 33-39 zones, in addition to the anisotropy introduced by the localisation. Size effects are therefore explicitly integrated in the constitutive behaviour. What are the new features? The above approach can be viewed as a two-stage smearing process. The velocity jump between the two sides of the Figure 2: 1D constitutive behaviour and corresponding sizes localisation band is at first smeared over the localisation band with physical size h (equation (2)). Therefore the The relationship between the FPZ size h, specific integration of inelastic constitutive behaviour for stress dissipation g, as the area under the stress-strain curve, and internal variables in the localisation zone is only and fracture energy G, as the energy released due to the related to this physical size and are invariant with respect creation of new surface area can be written as: to the discretisation size H. To obtain the relaxation strain rate and then stress rate of the elastic bulk, a compatibility argument (equation (3)) is utilised, effectively consisting of further smearing of the velocity jump over the discretisation size H. The difference with the traditional The internal equilibrium in this 1D case reads: smeared crack approach is obvious: in the latter, everything while the strain rate inside the FPZ takes the simple form is smeared over the discretisation size H, meaning that . From equations (4-6) we can write: the constitutive response must vary with the resolution of the discretisation to maintain the same amount of energy dissipation. As a consequence, very big H required in large scale modelling due to limited computing resource leads to This results in the strain rate inside the localisation zone unphysical scaling of model parameters which may still be in terms of the total strain rate , as: unable to prevent snapback in the constitutive behaviour (e.g. Jirasek and Bazant, 1995; Sulsky et al., 2007). The current approach, in contrast, gives a direct access to additional degrees of freedom related to the strain within Using (4), we can then write the stress rate in terms of the localisation band of size h, in which the constitutive law strain rate , as: is physical and hence does not snap back. Numerically, this means that the stress return algorithms for rate constitutive equations can always be performed. We illustrate the above points using a simple constitutive In the context of quasibrittle failure, the fracture energy model in 1D setting (Figure 2). This is a softening G (equation (8)) is a material constant representing the behaviour observed in failure of quasibrittle materials like contributions from the release of surface energy due rocks or concrete in which the pre- and post-peak responses to micro-cracking in the FPZ of size h. Therefore the are represented by the slopes ao and ai of the stress-strain inelastic behaviour in the localisation zone characterised curve. The material is homogeneous up to the peak point, by the softening modulus ai (Figure 2) is dependent on G after which the response bifurcates into two branches with and the size h of the FPZ, as intrinsic material parameters. slopes ai and ao corresponding to the inelastic behaviour Expressing everything in terms of G and h, we obtain: inside the Fracture Process Zone FPZ of size h and elastic behaviour in the outside bulk of size H-h (Figure 2). 36 Nguyen, G. (2014). Obras y Proyectos 15, 33-39 bar length H on the overall response, while keeping h fixed. It can be seen that regardless of the response, the area under the stress-displacement curve remains unchanged and is The stress rate in (13) is exactly that obtained from a always equal to the dissipation in the FPZ (see equation smeared crack approach, e.g. smearing G over the size H (8)). The stress strain response inside the FPZ (Figure 3a) of the finite element. Snapback in the overall response, e.g. is invariant with the size H, while that is not the case with , is present if H is sufficiently large to smeared crack approach (Figure 3c) due to the fact that make the denominator positive. However, the constitutive everything is smeared over the bar length H. In short, the response is integrable thanks to the use of an internal single stress component in the smeared crack approach equilibrium equation (6), together with the enrichment requires the variation of the constitutive behaviour with stress σi and strain εi . Regardless of the discretisation size respect to the discretisation to meet the requirement on H, the strain rate of the response inside the FPZ is always the dissipation. The stress-strain behaviour in such cases a positive quantity representing the separation of the just reflects the overall stress-displacement response, e.g. material. Therefore beyond the elastic regime, a negative essentially strain is obtained by smearing the displacement value of , due to snapback in the overall response, over the length H. The coupled stresses in the proposed still results in a positive strain rate for the integration constitutive modelling framework allow us to keep a of inelastic incremental constitutive equations. In the meaningful constitutive response inside the FPZ that is context of numerical failure analysis, εi is monotonically invariant with the discretisation (Figure 3a). Essentially, increasing with failure progression and it can be used as all parameters of the model remain unchanged with respect a control parameter in any indirect displacement control to the resolution of the discretisation, a property that is solution schemes (e.g. the local arc-length control by May missing in traditional smeared crack approach. and Duan (1997) and Yang and Proverbs, (2004)). Physically, snapback in a structural sense (e.g. finite element response) is a consequence of insufficient resolution of the discretisation in large scale failure analysis that, according to the literature, usually requires the enhancement of deformation mode of finite elements, e.g. using EAS or XFEM type enrichments. The proposed framework is another kind of enhancement, relying solely on the constitutive behaviour at integration point level with the use of the enrichment stress σi and strain εi . The advantage over traditional smeared crack is obvious, while the simplicity and practicality advantages over sophisticated enrichment methods like EAS or XFEM can also be seen. We will further illustrate this point in a forthcoming paper. Figure 3: Size effects on constitutive behaviour: a) normalised stress-strain response in the FPZ, b) normalised stress-normali- Numerical examples sed displacement, and c) normalised stress-strain response of a The first numerical example shows the response of a bar of smeared crack model length H, with unit cross sectional area and an embedded In the second example, a 10 mm long bar, clamped at one FPZ of size h (see Figure 2). Only the constitutive response end and free at the other, is discretised using 2D finite is concerned in this example, which is equivalent to the elements under plane strain conditions. Localised failure behaviour of a single linear 1D finite element of length H. is triggered off by weakening the element next to the The fracture energy is taken as 15 times the elastic strain clamp (tensile strength reduced by 10%). We used a linear energy at peak. Figure 3b shows the effects of varying the softening law (see expressions (8-13)) with the following 37 Nguyen, G. (2014). An enriched constitutive modelling framework for localised failure of geomaterials. Obras y Proyectos 15, 33-39 properties: Young’s modulus E = 30000 MPa, Poisson’s ratio v = 0.2, tensile uniaxial strength S = 3 MPa, and fracture energy G = 0.0015 Nmm/mm2. Figure 4a shows that the structural response of the bar is insensitive to the discretisation. The evolution of the displacement profiles in 2 cases, coarse and fine meshes, are also depicted and seen to coincide above the resolution of the coarse mesh. In the third example, a three point bending test is analysed using a simplified version of the constitutive model proposed in Nguyen and Korsunsky (2008), embedded in the above constitutive modelling scheme. The model is implemented in an in-house numerical code based on the Material Point Method MPM (Sulsky et al., 1995). Details on the implementation of the model along with some computational issues related to the application of the proposed scheme will be covered in a forthcoming paper. Figure 5: Three-point bending test: geometry and load-deflection The geometrical data and material properties are taken from response the experimental test of Petersson (1981) and illustrated in Figure 5. As can be seen in Figure 5, the numerical responses are insensitive to the resolution of the discretisation and also Conclusions the integration scheme of the MPM. Integration schemes We developed a new framework that allows the integration using 4 and 16 material points MPs per element result in of a length scale in constitutive models. The development almost same response. In all cases, the numerical predictions is based on the localised nature of failure in geomaterials closely follow the experimental trends. and utilises the internal equilibrium across the localisation band. This results in constitutive models possessing a length scale, and featuring coupled stress behaviour. The fact that both kinematical compatibility and traction continuity are enforced at the constitutive level in the formulation, e.g. at integration points, makes the implementation in any numerical code straightforward. As a consequence, the new approach is applicable to any existing constitutive model and also any discretisation scheme. This is totally different from traditional approaches that always require the unphysical scaling of model parameters with the resolution of the discretisation (e.g. smeared crack approach), or modification of existing finite elements (e.g. EAS or XFEM). Numerical examples in this paper, in the context of quasibrittle failure, demonstrate the essential features and the promising performance of the new approach. Acknowledgements Financial support from the Australian Research Council for project DP1093485 is gratefully acknowledged. Figure 4: A bar under tension: a) load-displacement response, and b) displacement profiles 38 Nguyen, G. (2014). Obras y Proyectos 15, 33-39 Oliver, J. (1996). Modelling strong discontinuities in solid References mechanics via strain softening constitutive equations. Part 1: Belytschko, T., Fish, J. and Engelmann, B.E. (1988). A finite fundamentals. International Journal for Numerical Methods in element with embedded localization zones. Computer Methods Engineering 39, 3575-3600 in Applied Mechanics and Engineering 70, 59-89 Petersson, P.E. (1981). Crack growth and development of Borja, R.I. (2000). A finite element model for strain localization fracture zones in plain concrete and similar materials. Report analysis of strongly discontinuous fields based on standard TVBM-1006, Div. of Build. Mat., Lund Institute of Technology, Galerkin approximation. Computer Methods in Applied Lund, Sweden Mechanics and Engineering 190, 1529-1549 Pijaudier-Cabot, G. and Bazant, Z.P. (1987). Nonlocal damage Cedolin L. and Bazant, Z.P. (1980). Effect of finite element theory. Journal of Engineering Mechanics 113(10), 1512-1533 choice in blunt crack band analysis. Computer Methods in Samaniego, E. and Belytschko, T. (2005). Continuum– Applied Mechanics and Engineering 24(3):305-316 discontinuum modelling of shear bands. International Journal Chen, Z. and Schreyer, H.L. (1987). Simulation of soil-concrete for Numerical Methods in Engineering 62, 1857-1872 interfaces with nonlocal constitutive models. Journal of Sanborn, S.E. and Prévost, J.H. (2011). Frictional slip plane Engineering Mechanics 113, 1665-1677 growth by localization detection and the extended finite element Foster, C.D., Borja, R.I. and Regueiro, R.A. (2007). Embedded method (XFEM). International Journal for Numerical and strong discontinuity finite elements for fractured geomaterials Analytical Methods in Geomechanics 35, 1278-1298 with variable friction. International Journal for Numerical Sluys, L.J. and Berends, A.H. (1998). Discontinuous failure Methods in Engineering 72, 549-581 analysis for mode-I and mode-II localization problems. Garikipati, K. and Hughes, T.J.R. (2000). A variational International Journal of Solids and Structures 35, 4257-4274 multiscale approach to strain localization - formulation for Sulsky D., Zhou S-J. and Schreyer, H.L. (1995). Application of multidimensional problems. Computer Methods in Applied a particle-in-cell method to solid mechanics. Computer Physics Mechanics and Engineering 188(1-3), 39-60 Communications 87, 236-252 Jirasek, M. and Bazant, Z.P. (1995). Particle model for quasibrittle Sulsky, D., Schreyer, H., Peterson, K., Kwok, R. and Coon, fracture and application to sea ice. Journal of Engineering M. (2007). Using the material-point method to model sea ice Mechanics 121, 1016-1025 dynamics. Journal of Geophysical Research 112, C02S90 Kolymbas, D. (2009). Kinematics of shear bands. Acta Vardoulakis, I., Goldscheider, M. and Gudehus, G. (1978). Geotechnica 4, 315-318 Formation of shear bands in sand bodies as a bifurcation problem. Larsson, R., Runesson, K. and Sture, S. (1996). Embedded International Journal for Numerical and Analytical Methods in localization band in undrained soil based on regularized strong Geomechanics 2, 99-128 discontinuity—theory and FE-analysis. International Journal of Wells, G.N. and Sluys, L.J. (2001). A new method for modelling Solids and Structures 33, 3081-3101 cohesive cracks using finite elements. International Journal for May, I.M. and Duan, Y. (1997). A local arc-length procedure for Numerical Methods in Engineering 50, 2667-2682 strain softening. Computers & Structures 64(1-4), 297-303 Yang, Z.J. and Proverbs, D. (2004). A comparative study of Nguyen, G.D., Einav, I. and Korsunsky, A.M. (2012). How numerical solutions to non-linear discrete crack modeling to connect two scales of behaviour in constitutive modelling of concrete beams involving sharp snap-back. Engineering of geomaterials. Géotechnique Letters (special issue on Fracture Mechanics 71, 81-105 Geomechanics Across the Scales) 2(3), 129-134 Nguyen, G.D. and Korsunsky, A.M. (2008). Development of an approach to constitutive modelling of concrete: isotropic damage coupled with plasticity. International Journal of Solids and Structures 45(20), 5483-5501 39