(PDF) A constitutive modelling framework

Engineering Fracture Mechanics xxx (2013) xxx–xxx Contents lists available at ScienceDirect Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure Giang D. Nguyen a,⇑, Alexander M. Korsunsky b, Itai Einav c a School of Civil, Environmental and Mining Engineering, The University of Adelaide, Adelaide, SA 5005, Australia b Department of Engineering Science, The University of Oxford, Parks Road, Oxford OX1 3PJ, UK c School of Civil Engineering, The University of Sydney, Sydney, NSW 2006, Australia a r t i c l e i n f o a b s t r a c t Article history: We propose a constitutive modelling framework with enhanced kinematics to capture Received 13 May 2013 localised mode of deformation. The total strain is decomposed into two components to Received in revised form 4 September 2013 reﬂect an inelastic localisation band embedded in an elastic bulk. This is the usual case Accepted 1 November 2013 in numerical analysis of localised failure in geomaterials, when the size of the localisation Available online xxxx band is very small compared to an element of the discretised domain under consideration. The proposed framework takes into account the sizes and corresponding behaviours of the Keywords: two inelastic and elastic zones and hence gives derived constitutive models a length scale. Length scales Constitutive modelling This is an essential feature in dealing with size effect issues as a consequence of localised Localised failure failure in geomaterials. The proposed framework is applied to a constitutive model for the Discontinuity failure analysis of quasi-brittle materials. The implementation algorithms are developed Quasi-brittle and novel features are illustrated through numerical examples. Fracture process zone Ó 2013 Elsevier Ltd. All rights reserved. Softening Damage 1. Introduction Localised failure in geomaterials involves the inelastic loading inside a very thin, but ﬁnite thickness, band, while the material outside the band usually unloads elastically. The dissipation properties of the material are therefore fully governed by what actually happens inside this very thin localisation zone. The thickness of this band is related to the microstructure of the material and is a physical quantity that brings size effect features to the material response [1,2]. Shear/compaction/dila- tion bands observed in soils, porous rocks, shear bands in metallic materials and Fracture Process Zone (FPZ) in quasi-brittle materials are typical examples of this localised failure. In such cases the homogeneity of the continuum is lost, with discon- tinuity of kinematic ﬁelds across the boundary of the localisation band. This means the deﬁnitions of stress and strain over a volume element containing both the bulk elastic and an embedded inelastic localisation zone are not physically meaningful anymore. Given the above aspects of localised failure, it is essential that constitutive models take all of them into account by representing the size and behaviour of the material inside the localisation band as well as those of the surrounding bulk elas- ticity. Failure to do so, which is the case in most classical constitutive models, leads to the loss of physical meaning of the macro stress and strain, as these quantities in such cases are deﬁned over a non-homogenous volume element. ⇑ Corresponding author. Tel.: +61 8 83132259. E-mail addresses:

[email protected]

(G.D. Nguyen). 0013-7944/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 2 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx Nomenclature GF mode I fracture energy gF speciﬁc fracture energy; gF = GF/h h effective width of the fracture process zone H effective size of the spatial discretisation f volume fraction Lc critical length ft0 uniaxial tensile strength E Young modulus m Poisson’s ratio ao isotropic elastic stiffness tensor of the bulk material ai isotropic elastic stiffness tensor of the material inside the localisation zone aTi tangent stiffness tensor of the material inside the localisation zone Di isotropic elastic compliance tensor of the material inside the localisation zone Ao acoustic tensor associated with ao Ai acoustic tensor associated with ai ATi acoustic tensor associated with aiT Ki elastic stiffness of the localisation zone idealised as a surface KTi tangent stiffness of the localisation zone idealised as a surface n normal vector of the localisation band Xo volume of the material outside the localisation zone Xi volume of the material inside the localisation zone X total volume of the sub-domain; X = Xi + Xo D scalar damage variable [u] displacement jump e total average strain ei total strain of the material inside the localisation zone eo total strain of the material outside the localisation zone epi plastic strain of the material inside the localisation zone eei elastic strain of the material inside the localisation zone r average stress tensor r+ ‘‘positive’’ part of the stress tensor r from the eigenvalue decomposition ri stress tensor of the material inside the localisation zone ro stress tensor of the material outside the localisation zone rm the mth principal stress ti traction at the boundary of the localisation zone pm the unit vector of the mth principal direction Q tensor deﬁned by the ratio between plastic strain rate and damage rate W Helmholtz free energy potential of a unit volume gi Gibbs free energy potential for a unit volume inside the localisation zone Wi Helmholtz free energy potential for a unit volume inside the localisation zone Wo Helmholtz free energy potential for a unit volume outside the localisation zone W C Helmholtz free energy potential for the localisation zone idealised as a surface U dissipation potential (total rate of dissipation) for a unit volume U i dissipation potential for a unit volume inside the localisation zone Uo dissipation potential for a unit volume outside the localisation zone /D damage component deﬁning the explicit form of the dissipation potential /p plasticity component deﬁning the explicit form of the dissipation potential UD damage dissipation (damage rate of dissipation) Up plastic dissipation (plastic rate of dissipation) vD damage energy (thermodynamically associated with the damage variable D) F⁄(D, ri) function controlling the overall rate of dissipation F(D) function controlling the evolution of damage process x parameter controlling the partition of fracture energy Ep parameter controlling the initial slope of the softening stress–strain curve n parameter controlling the slope change of the softening stress–strain curve dk non-negative plasticity-like multiplier y loading function, or combined damage-yield function in stress/strain space y⁄ loading function, or combined damage-yield function in mixed energy–stress space d50 median diameter for which half of the sample is ﬁner dmax maximum aggregate diameter Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx 3 Mathematically it leads to the ill-posedness of the Boundary Value Problems (BVPs) and consequently the dependence of the numerical solutions on the spatial discretisation [2]. Several approaches have been proposed to tackle the modelling of localised failure in the analysis of solids and structures. In the context of constitutive modelling, they range from simple and practical one such as the smeared crack approach [3], to mathematically sophisticated regularisations such as Cosserat [4] and nonlocal/gradient theories [5–12]. The key point of these approaches is to introduce a length scale related to the localised failure into the constitutive equations. In the smeared crack approach [3], it is the size of the ﬁnite element that governs the constitutive equations. Essentially the constitutive behaviour must be scaled with the resolution of the discretisation. This scaling is used on the basis of meeting requirement on the energy dissipation, e.g. forcing cracked elements to reproduce a dissipation invariant with the size of the element. However, this is unphysical as there is no unique set of parameters for the same material, due to the variation of the spatial resolution over the computational domain. In addition, once a coarse spatial discretisation is used (e.g. in large scale anal- ysis), the speciﬁc fracture energy as the ratio between the (surface) fracture energy and the characteristic size of the discret- isation becomes smaller than the elastic strain energy at peak, resulting in the inadmissible snapback instability in the constitutive response [13]. Enrichment of constitutive models with an internal length through the use of nonlocal/gradient regularisation [5–12] is a mathematically and probably physically rigorous approach to account for the size of the localisation zone. However, these models require that the discretisation be ﬁner than the physical width of the localisation zone. Since locations of the failure zones are generally unknown and the considered domain can be of several orders of magnitude bigger than the characteristic width of the localisation zone, very high resolution of the discretisation is usually required resulting in very large model size and hence the applications of such enhancements are severely restricted by the available computational resources. Examples include failure of geomaterials, such as the formation of sea ice leads in the Arctic [13,14], and rock fracture in underground mining. Enrichments to the discretisation scheme such as the Enhanced Assumed Strain (EAS) method (e.g. [15–19]) or the eXtended Finite Element Method (XFEM) [20–23] are popular to address the localised failure in which the width h of the failure zone (Fracture Process Zone in failure of quasi-brittle materials, or shear/compaction bands in soils) is small com- pared to the solid/structure under consideration. If H is the resolution of the spatial discretisation (e.g. size of ﬁnite element, grid size in ﬁnite difference), then it is the case h < H. In these kinds of enhancement, the ﬁnite element shape functions are usually enriched to improve the kinematics of deformation the element can handle. This is due to the embedding of the local- isation zone, usually idealised as a strong discontinuity across the element. The richness of the deformation that includes both the discontinuity/crack opening/shearing and the shrinking of the elastic bulk can then be appropriately captured. This also helps to remove unphysical snapbacks in the constitutive response usually encountered in large scale analysis using smeared crack approach. However a separate constitutive law for the discontinuity, disconnected from the bulk continuum behaviour, is required to model the crack opening/shearing or behaviour inside the localisation band. Other similar approaches to improve the kinematics of ﬁnite elements include that by Belytschko and co-authors [24], Sluys and Berends [25], Garikipati and Hughes [26,27], and more recently Huespe et al. [28] in the context of ductile failure. In Belytschko et al. [24] and Sluys and Berends [25] the localisation zone appears as a ﬁnite width band embedded in a ﬁnite element. While Belytschko et al. [24] built their approach on the Hu-Washizu three-ﬁeld variational principle, which is sim- ilar to the case of the EAS mentioned above, Garikipati and Hughes [26,27] developed their own ﬁnite element with embed- ded zone starting from a variational multiscale approach. In contrast with the use of a strong discontinuity (via cohesive models) in the EAS and XFEM, the ﬁnite elements with embedded localisation zone [24,25,28] introduce weak discontinuity (e.g. only the strain ﬁeld is discontinuous) to the kinematic ﬁelds. In accordance with this weak discontinuity, a continuum model is used to describe the inelastic material behaviour inside the localisation zone. Can we come up with a constitutive modelling framework better than the smeared crack approach and independent of the type of Finite Element (FE), but less complicated than the EAS or XFEM? This is because all the above addressed ap- proaches dealing with the case the h < H (width h of the localisation band smaller than the resolution H of the spatial dis- cretisation) involve ﬁnite element re-formulations, e.g. modiﬁcation and/or introduction of shape functions, and hence result in the dependence of the approach on the type of ﬁnite element used for the numerical discretisation. The key idea described in the present paper is to enhance the constitutive behaviour descriptions with an additional kinematics mode and a corresponding length scale related to the width of the localisation zone to correctly describe localised failure of mate- rials. As a consequence there is no need to enrich the kinematics of existing ﬁnite elements and this facilitates the imple- mentation of the approach in existing ﬁnite element codes. We keep the usual interface with existing ﬁnite element formulation through the deﬁnition of averaged stress, which could be considered as stress on the boundary of a volume ele- ment containing an embedded localisation band. This stress and corresponding averaged strain in fact represent the com- posite behaviour of the continuum idealised as an elastic bulk with an embedded localisation band. Localised failure once triggered after a diffuse failure stage will lead to the bifurcation of constitutive response into inelastic behaviour inside the localisation band and elastic unloading outside it (see Fig. 1 later). It is noted that this bifurcation requires two descrip- tions following difference branches, inelastic loading and elastic unloading, of the same constitutive model. In this case inter- nal variables characterising the inelastic behaviour are deﬁned inside the localisation band, but not smeared out over the whole volume element like in smeared crack approach. An internal equilibrium condition in terms of traction continuity across the boundary of the localisation band is utilised to connect the two sets of stress–strain relationships. The resulting constitutive model consists of two coupled parts: elastic behaviour of the bulk continuum, and inelastic response inside the Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 4 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx Fig. 1. (a) Numerical discretisation and a localisation zone (darkened) (after [29]). (b) Corresponding material responses inside and outside the localisation zone. localisation band driven by the physics of the material failure, e.g. FPZ zone due to micro-cracking in quasibrittle failure, compaction band due to grain crushing in cemented rocks, shear band in metals. This paper follows our earlier publication [29], where some fundamental relationships have been initiated and presented. In this paper we provide more physical ground, implementation algorithms and practical applications to complete the pic- ture of the proposed modelling strategy and provide a basis for further developments. In this sense, Section 2.1 based on Nguyen et al. [29] sets a background for the modelling. The energy aspects of the approach are then addressed in Section 2.2, followed by an illustration of the features of constitutive models derived from this approach in the context of quasi-brit- tle failure. In Section 3 we provide details of a damage model to be used with the proposed approach. The numerical imple- mentation algorithms are then described in Section 4 and the combination of this damage model with the proposed modelling for analysing quasi-brittle failure is presented in Section 5. 2. A framework for constitutive models with embedded localisation zone 2.1. The basis We know that traditional constitutive models are deﬁned over a unit volume in which the deformation is usually as- sumed to be homogeneous. This assumption is however not always valid, especially in modelling materials that exhibit localised failure. Cracking in quasi-brittle materials, shear/compaction bands in soils and soft rocks are typical examples of this localised failure, when there is a strong variation of strain across a small band width. Continuum assumption breaks down and the deﬁnitions of stress, strain or any physical quantities as averages over the volume crossed by the localization band are not appropriate. We aim to develop at ﬁrst a general constitutive modelling framework that can be used for any constitutive models for materials exhibiting localised failure. The focus is the post-localisation behaviour (Fig. 1), as it is where all the above men- tioned issues come. In this sense, the proposed approach is switched on when localised failure is detected, e.g. via the loss of positive deﬁniteness of the determinant of the acoustic tensor [22,30], while the preceding diffuse failure can be treated in a classical way using the same constitutive model (see Fig. 1b). With the above issues of localised failure in mind, it can be anticipated that a single stress/strain measure over the volume element crossed by the localisation band is not adequate. Besides, internal variables of continuum models must be correctly deﬁned over the localisation zone where inelastic defor- mation takes place. These are prerequisites that classical constitutive models do not meet, as all quantities and measures in these models are usually deﬁned over the whole unit volume. This results in the scaling of the structural behaviour with respect to the volume occupied, which in turn leads to severe physical and numerical issues in the modelling. A well known issue in using classical constitutive models without regularisation for analysis of quasi-brittle structures is the convergence to elastic structural snap back (corresponding to zero dissipation) when the mesh size tends to zero [2]. Treatment of the issues, e.g. the smeared crack/deformation approach, is also incomplete if the physical phenomenon behind localised failure is not properly captured. In such a case, the unphysical constitutive snapping back of the inelastic model is encountered if the size of the element is large enough due to the smearing and scaling of the energy across the whole sub-domain (e.g. ﬁnite element), not the physical region where dissipation takes place. Taking the case of a localisation band crossing a volume element, as a simpliﬁed and generalised conﬁguration of localised failure for large scale modelling, we can consider this as a composite material featuring two phases: inelastic behaviour in- side the localisation band of width h, and elastic material outside this band. The deﬁnitions of sizes h and H and correspond- ing volume fraction f = h/H are best illustrated in Fig. 1, and the term ‘‘large scale’’ here has a relative meaning, addressing the case h < H in modelling localised failure. We can see that this is general enough to encompass quasi-brittle failure at the structural/laboratory scale, and even the true large scale modelling of sea ice failure (e.g. in [13,14]) when the localisation band in fact consists of several other lower scale bands. Of course further development of the framework should allow the nested scales of localised failure, as we stressed out in [29], or the multiple intersecting localisation bands/cracks for non-monotonic loadings. While h is related to the inelastic behaviour and failure at the lower scale, the size H of the elastic bulk can be considered as a numerical size measure resulting from the discretisation of the domain under consideration. Grid spacing in the ﬁnite difference method and the domain size ascribed to an integration point in the ﬁnite element method are Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx 5 typical examples of the meaning we attach to the discretisation characteristic size H. As a consequence, H can vary depend- ing on the required resolution of the numerical discretisation, while h is considered as a material property and hence is invariant with the numerical discretisation. We are also aware of the fact that physically h can vary during the deformation process to represent the diffuse to localised failure, and this has been tackled by several researchers (e.g. [31,32]). This effect is not included in the current study, but preserved for future developments. We view the volume element crossed by a localisation band (Fig. 1) as a composite material consisting of two phases. Following mixture theory, the volume averaged total strain rate can be expressed as [29]: e_ ¼ f e_ i þ ð1 f Þe_ o ð1Þ where subscripts ‘‘i’’ and ‘‘o’’ are used for the strain inside and outside the localisation region, respectively. In line with the above strains, in principle it can be expected to have three different measures of stress, ri, ro and r, correspondingly, for this volume element (Fig. 1). However we need to keep in mind that for this conﬁguration of a composite material, the Hill-Man- del condition [33] should be met, stating that the work produced by the macroscopic stress r and strain rate e_ should equate the volume-averaged value due to the stresses and strain rates inside and outside the localised region. Utilising the deﬁnition of the strain rate e_ i (Eq. (2), below), the Hill–Mandel condition, and the continuity of traction across the boundary of the localisation zone, it can be proved that the stresses ro and r coincide. Details on this can be found in Nguyen et al. [29]. This result, subjected to the above three conditions, is in fact a sub set of a more general one. However we will show that this sub set is rich enough to enable us to adequately describe the constitutive behaviour of materials exhibiting localised failure without having to resort to post treatments of the constitutive equations (e.g. using nonlocal/gradient or viscous regularisa- tions) or element enrichments (e.g. using the EAS or XFEM). Following Vardoulakis et al. [34] and Kolymbas [35] the strain rate e_ i measured from the relative velocity ½u _ between opposite sides of the localisation band can take the following form: 1 1 _ s¼ e_ i ¼ ðn ½uÞ _ þ ½u ðn ½u _ nÞ ð2Þ h 2h where n denotes the normal vector of the band (Fig. 1) and is the tensor product. In the above expression, the contribution from the homogeneous strain rate e_ o is neglected, given the fact that its magnitude is negligible compared to the term 1 h _ s inside the band. Further details and reasoning can be found in [29]. ðn ½uÞ Regarding the stress–strain responses, general assumptions about the behaviour are used to make the development appli- cable to any constitutive models. It is reasonable to state that the material inside the localisation band undergoes inelastic deformation while that outside the band is unloading. This unloading process is assumed to be elastic in this current devel- opment. While this is an acceptable assumption in the context of quasi-brittle failure that we are dealing with (e.g. concrete and rock fracture), we need to keep in mind that it may not always be the case, e.g. in shear/crushing failure of granular materials the unloading process could be inelastic. In general, only the tangent stiffnesses of the corresponding behaviours are needed here in the following formulation. The material inside the localisation is loading with tangent stiffness aTi , while that in the outside elastic bulk is assumed to have (elastic) stiffness ao (Fig. 1). Due to these kinds of responses and the vol- ume averaging of the total strain in (1), we can write the following rate constitutive relationships for the inelastic behaviour inside the localisation zone: r_ i ¼ aTi : e_ i ð3Þ and the elastic unloading for the outside bulk: 1 r_ ¼ r_ o ¼ ao : e_ o ¼ ao : ðe_ f e_ i Þ ð4Þ 1f As stated above these stress rates cannot be independent, but are connected through the continuity of traction across the boundary of the localisation band: ðr_ o r_ i Þ n ¼ 0 ð5Þ Using Eqs. (1)–(4), we can write the above internal equilibrium condition (5) of the constitutive model in the following form: 1 f 1 ao : _ s n ¼ aTi : ðn ½uÞ e_ ðn ½uÞ _ s n ð6Þ 1f h h The symmetry of the strain tensor is taken to obtain: f ð1 f Þ T 1 ð1 f Þ T ðao : e_ Þ n ¼ Ao þ _ ¼ Ai ½u Ao þ _ Ai ½u ð7Þ h h H h where AiT = n aTi n and Ao = n ao n are the acoustic (or localisation) tensors associated with the tangent stiffnesses inside and outside the localisation band, respectively. The above gives the velocity jump ½u _ in terms of the total strain rate e_ as: Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 6 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx 1 1 ð1 f Þ T _ ¼ ½u Ao þ Ai ðao : e_ Þ n ¼ C1 ðao : e_ Þ n ð8Þ H h The above expression of ½u_ can then be substituted into Eq. (4), together with the form (2) of the inelastic strain e_ i , to obtain the stress rate r_ as: 1 1 s r_ ¼ ao : e_ ðn ðC1 ðao : e_ Þ nÞÞ ð9Þ 1f H This stress rate r_ and strain rate e_ , as averages deﬁned over the whole volume X, are needed to interface with existing constitutive modelling structure in the numerical implementation of the model. As seen in Eq. (9), their relationship is gov- erned by both the behaviour and corresponding sizes of the material inside and outside the localisation band. In other words, while at the end a relationship r_ e_ is required, the underlying structure of the constitutive equations in this case is different from that of a classical one. We will show the characteristics of the new approach through an example in the coming section. At the moment, there is no other condition on h and H besides the h < H requirement, which is valid if the spatial discret- isation is sufﬁciently larger than the physical size of the localisation zone. We also know that the other case h > H is more suitable for other kinds of enhancement using nonlocal or gradient theories, which is more computationally demanding. Although physically there is no h = 0, (e.g. the localisation zone, although very thin in some cases, must physically have a ﬁnite width) the case h ? 0 and how it can be derived from or connected with the above formulation deserves a discussion, as it is widely encountered in the literature. From the traction at the boundary of the localisation zone, we have 1 1 1 _ s n ¼ n aTi n ½u t_ i ¼ r_ i n ¼ aTi : ðn ½uÞ _ ¼ ATi ½u _ ð10Þ h h h As can be seen in the above equation, the term 1h ATi ¼ KTi has the physical meaning of the stiffness across the localisation band. Since h is physically non-zero, this stiffness KTi must admit a ﬁnite value. A mathematical argument, following Simo et al. [36], can also be used here to arrive at the same condition for KTi . This is based on the non-trivial solution of Eq. (8) when h ? 0. In such a case, it is required that 1h ATi ¼ KTi be bounded so that (8) can give a non-null velocity jump ½u. _ In the context of modelling localised failure, the case h ? 0 is relevant to quasi-brittle failure, when the physical width of the fracture process zone (FPZ) is usually much smaller than the resolution of the numerical discretisation. As an example, h/dmax = 1–4 [37,38] (dmax is the maximum aggregate diameter) has been experimentally observed in concrete, and the con- tinuum assumption in such cases is only valid at a scale much bigger than the aggregate size. Alternatively, as mentioned in the introduction, due to limitation on computational power in large scale modelling of failure, the resolution of the numer- ical discretisation is usually of higher order of magnitude than the physical size of the failure zone. This renders the assump- tion h ? 0 (or more precisely h H) reasonable and favours correspondingly a discrete description of the localisation band using cohesive zone models. The link between discrete and continuum descriptions of material behaviour inside the local- ization zone will be elaborated in the next Sections. 2.2. Energy aspects We show in this section the general thermo-mechanical implications of the above development. Considering the conﬁg- uration of a composite material consisting of two phases, the Helmholtz free energy of the system can be written as: W ¼ f Wi eei ; a þ ð1 f ÞWo ðeo Þ ð11Þ where for isothermal processes, Wi eei ; a and Wo(eo) represent the speciﬁc elastic strain energies of the material inside and outside the localisation zone, respectively. A generalised internal variable a is used here to represent micro-structural pro- cesses occurring in the localisation zone (for simplicity, subscript ‘‘i’’ on a is dropped due to the assumption that inelastic behaviour takes place only inside the localisation zone). Example of a includes the damage variable representing micro- cracking processes, and the breakage internal variable and porosity characterising grain crushing and associated pore col- lapse, respectively, in crushable granular materials [39,40]. These are the processes that directly link the stiffness of the material to the micro-structural changes. For the case h ? 0, the ﬁrst term in the above expression can be written as a surface energy potential WC. In particular, using [u]e as the recoverable part of the total displacement jump [u] and assuming the elastic strain eie takes the same form as the total strain, we have: 1 eei ¼ ðn ½ue Þs ð12Þ h the ﬁrst term in the right hand side of (11) can then be expressed for h ? 0 as: h 1 1 s 1 s f Wi ¼ ðn ½ue Þ : ai ðaÞ : ðn ½ue Þ ð13Þ H2 h h Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx 7 In the above expression, the internal variable a is assumed to directly govern the elastic stiffness ai of the material inside the localisation zone. In the context of damage mechanics [41] it is the effect of the damage variable on the elastic stiffness. This is physically reasonable and does not affect the generality of the results in this development. The above expression can be written as: 1 1 e 1 f Wi ¼ ½u Ki ðaÞ ½ue ¼ WC ð½ue ; aÞ ð14aÞ H2 H or, using f = h/H (see Section 2.1): 1 Wi ¼ WC ð½ue ; aÞ ð14bÞ h where Ki ðaÞ ¼ 1h ½n ai ðaÞ n is the elastic stiffness of the very thin layer idealised as a zero thickness interface. It is consistent with the deﬁnition of the tangent stiffness KiT in the previous section. The Helmholtz energy WC appearing in Eq. (10) rep- resents the total energy per unit area stored in the localisation zone idealised as a zero thickness interface. For this idealised case of h ? 0 (or h H), the speciﬁc Helmholtz free energy potential takes the following form: 1 W¼ WC ð½ue ; aÞ þ Wo ðeo Þ ð15Þ H We can see that the surface quantity WC is smeared over the size H of the domain X to make a meaningful expression for the deﬁnition of W over a unit volume. The above mathematical structure is similar to that assumed by Armero in [42], ex- cept the fact that the unbounded term in [42] appears as 1/H. The current expression (11), however, is inspired by the mix- ture theory without having to resort to the more classical reasoning using bounded and unbounded terms introduced by Simo et al. [36] and later adopted by several researchers (e.g. [15,16,18]). We make a branching at this point to address two related but different cases: ﬁnite and (mathematically) inﬁnitesimal width h of the localisation band. 2.2.1. The case of ﬁnite h For rate-independent isothermal process, the speciﬁc dissipation rate U, due to irreversible processes can be expressed as: f W Ue ¼ W _ P0 ð16Þ where W f is the rate of work for a unit volume in X. Taking into account the composite nature of the material occupying X (see Fig. 1), we have: f ¼ f ri : e_ i þ ð1 f Þro : e_ o W ð17Þ _ ei Substituting (11) and (17) and the incremental relationship e_ i ¼ e þ e in (16) results in: _ pi @ Wi e @ Wi @ Wo Ue ¼ f ri : e_ ei þ e_ pi þ ð1 f Þro : e_ o f : e _ i þ : a _ ð1 f Þ e : e_ o ð18Þ @ eei @a @ eo For reversible processes, there is no dissipation (U ¼ 0) and the rates of internal variables representing inelastic/dissipa- tive processes therefore vanish. The above expression reduces to: @ Wi @ Wo Ue ¼ f ri : e_ ei þ ð1 f Þ ro e : e_ o ¼ 0 ð19Þ @ eei @ eo Since this must hold for any strain rates, and any arbitrarily small volume X, we obtain the elasticity laws: @ Wi ri ¼ ð20Þ @ eei and @ Wo ro ¼ ð21Þ @ eeo The speciﬁc dissipation rate then reduces to: @ Wi Ue ¼ f ri : e_ pi : a_ P 0 ð22Þ @a Therefore for the domain occupied by volume X, the total dissipation rate is: X e ¼ f X ri : e_ p @ Wi : a_ ¼ h AH ri : e_ p @ Wi : a_ ¼ Ah ri : e_ p @ Wi : a_ ¼ XU ð23Þ i i i @a H @a @a Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 8 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx where A is the projected area of the idealised localisation zone (see Fig. 1). We can also write the above in terms of surface quantities, using the relationship (14a): X 1 1 @ WC @ WC _ p ¼ Ah t ½u _ p : a_ ¼ A t ½u : a_ ð24Þ h h @a @a 2.2.2. The case of inﬁnitesimal h While the energy balance equation keeps the same form (16), the speciﬁc Helmholtz free energy is of the form: 1 W¼ WC ð½ue ; aÞ þ Wo ðeo Þ ð25Þ H The corresponding rate of work per unit volume is: f ¼ 1 ti ½u W _ þ ro : e_ o ð26Þ H This leads to the following form of the speciﬁc dissipation rate: e¼ 1 1 @ WC e @ WC @ Wo U _ þ ro : e_ o ti ½u ½ _ u þ : a _ e : e_ o ð27Þ H H @½ue @a @ eo e _ which, after making use of the decomposition of the velocity jump into recoverable (or elastic) and irrecoverable parts, ½u _ p , respectively, can then be rearranged to become: and ½u 1 @ WC 1 1 @ WC @ Wo e¼ U ti e _ e þ ti ½u ½u _ p : a_ þ ro e : e_ o ð28Þ H @½u H H @a @ eo Again, following the same thermo-mechanical argument above for the case of totally reversible deformation, U ¼ 0, we ob- tain the elasticity law for both the traction ti and bulk stress ro: @ WC ti ¼ ð29Þ @½ue @W ro ¼ eo ð30Þ @ eo The dissipation rate per unit volume then reduces to: 1 1 @ WC e¼ U _ p ti ½u : a_ ð31Þ H H @a The total dissipation rate in the domain X is then given as: X 1 1 @ WC @ WC e ¼ AH ¼ XU _ p ti ½u _ p : a_ ¼ A ti ½u : a_ ð32Þ H H @a @a As can be seen in (23), (24) and (32), the total dissipation rate for a domain occupied by X in this case is independent of X, and is solely governed by the inelastic behaviour and the size of the localisation zone (see Eq. (23)) or, in extreme case of h ? 0, the inﬁnitesimally thin surface A and its associated behaviour. In other words, there is no size H in these expressions and the dissipation does not scale with the size associated with the numerical discretisation of the domain under consideration, but with the physical size h of the embedded inelastic zone. This is an important characteristic that classical constitutive models lack. 2.3. Key characteristics of the approach A one-dimensional (1D) illustration is used to show the characteristics of the proposed framework in the context of dam- age mechanics modelling of fracture processes. The main purpose is to demonstrate the importance of a physical length scale in the constitutive model. The smeared crack is used as an example for a physical implication behind its regularisation using a numerical length scale. Two equivalent descriptions, discrete (h ? 0) and continuum (ﬁnite h) presented in the preceding section are depicted in Fig. 2. For quasi-brittle failure considered here, the fracture energy GF (Fig. 2) governs the traction- displacement jump response of the material model. The corresponding speciﬁc dissipation that governs the behaviour of the continuum model is gF = GF/x, as the area under the stress–strain curve (Fig. 2), where x is the length scale involved. For the smeared crack approach, everything is smeared over the size of the element, e.g. x = H (Fig. 3), and the tangent stiffness of the inelastic material can be worked out from Fig. 2: , ! 2GF ao aTi ¼ ao 1 02 ð33Þ Hf t Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx 9 Fig. 2. Correspondence between discrete and continuum models of the localisation zone. Fig. 3. 1D illustrations for the classical smeared crack and the proposed approach. The stress strain relationship is therefore: , ! 2GF ao r_ ¼ ao e_ 1 02 ð34Þ Hf t It can be seen that the constitutive response in this case scales with the size H of the element (coinciding with bar length in this example), as a simple yet effective way to keep the dissipation invariant with the numerical resolution. This avoids the drawback of classical continuum models for fracture/damage processes in having the dissipation linearly dependent on the size H of the discretisation. However this regularisation does not come without a penalty. Firstly, the physical meaning of constitutive model is lost in this case, e.g. the constitutive behaviour depends on the numerical discretisation but not the real physical length scale. Secondly and more signiﬁcantly, when H is large enough to make the denominator in (34) positive: 2GF ao H > Lc ¼ ð35Þ ft02 the constitutive behaviour exhibits snapping back, which in classical plasticity theory leads to a negative plasticity multiplier preventing the integration of the rate constitutive equations. The critical length Lc above imposes a limit on the resolution of the discretisation, beyond which the rate constitutive equations are not integrable. This restriction is a consequence of ignor- ing the mechanisms of localised failure in constitutive modelling. For the proposed approach, x = h, resulting in: , ! 2GF ao aTi ¼ ao 1 02 ð36Þ hf t Since h is a physical quantity measuring the width of the FPZ, the denominator of (36) is always negative, e.g. the speciﬁc fracture energy is always greater than the speciﬁc strain energy at peak stress (Fig. 2), and consequently aTi can never get above zero. The material is homogeneous up to the peak point, after which the response bifurcates into two branches with slopes aTi and ao corresponding to the inelastic behaviour inside the fracture process zone (FPZ) of size h and elastic behav- iour in the outside bulk of size H h (Fig. 3). The steps described in Section 2.1 are followed in this 1D case to obtain an incremental stress–strain relationship. The strain rate inside the FPZ can be expressed as e_ i ¼ ½u=h, _ and the internal equilib- rium (5) becomes r_ o ¼ r_ i , resulting in: 1 ao ðe_ f e_ i Þ ¼ aTi e_ i ð37Þ 1f The strain rate e_ i and corresponding stress rate can be obtained from the above equation as: ao e_ 1 aTi ao e_ e_ i ¼ T ; and r_ ¼ ao ðe_ f e_ i Þ ¼ ð38a; bÞ ð1 f Þai þ fao 1f ð1 f ÞaTi þ fao Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 10 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx Using (36), we can rewrite the above equations as: !, ! , ! 2GF ao 2GF ao 2GF ao e_ i ¼ e_ 1 02 1 02 ; and r_ ¼ ao e_ 1 02 ð39a; bÞ hf t Hf t Hf t While the overall incremental stress–strain relationship should be the same as that in smeared crack approach, the inter- esting part lies behind that stress–strain relationship. The strain ei in this case represents the enhanced mode of the consti- tutive behaviour, similar to the enhanced kinematics of ﬁnite elements in the EAS- or XFEM- based approaches. It is a monotonically increasing quantity irrespective to the resolution of the discretisation represented by H. Snapping back, although still existing (see Eq. (39a)), is no longer associated with the integration of the rate constitutive equation, but trans- ferred to a more structural behaviour. In such cases (H > Lc or the denominator in (39) is positive), the strain rate in the unloading elastic bulk is large enough to dominate the response of the constitutive model, making e_ < 0, but the strain rate e_ i representing the tension of the material inside the localisation zone is always a positive quantity. 3. A constitutive model for the localisation zone We have just described a general approach for a constitutive model featuring two scales of behaviour. What is left is a speciﬁc inelastic model to complete the descriptions of the approach and to further illustrate its capabilities. In the context of quasi-brittle failure, a damage mechanics model is used for the inelastic behaviour of the material inside the localisation zone. The elastic response just follows the elastic branch of this model. 3.1. Formulation The thermodynamic formulation of a coupled damage-plasticity model is brieﬂy described in this section, following our earlier developments [12,29]. The model is formulated for a unit volume of the material inside the localisation zone. Denot- ing D, ei, eip, and eie = ei eip the isotropic (scalar) damage variable, total, plastic and elastic strains, respectively, the Gibbs speciﬁc energy potential gi can be written as: ri : Di : ri gi ¼ ri : epi ð40Þ 2ð1 DÞ where Di is the (isotropic) elastic compliance tensor, expressed in terms of the Young modulus E and Poisson’s ratio m (dij is the Kronecker delta): ð1 þ mÞ 2m ðDi Þijkl ¼ dij dkl þ dik djl þ dil djk ð41Þ 2E 1þm From the energy potential gi, the strain and speciﬁc damage energy vD are: @g i Di : ri ei ¼ ¼ þ epi ð42Þ @ ri ð1 DÞ @g ri : Di : ri vD ¼ i¼ ð43Þ @D 2ð1 DÞ2 For the strong coupling between damage and plasticity, the dissipation potential is assumed of the form: qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Ui ¼ /2D þ /2p ð44Þ in which the contributions from damage (/D) and plasticity (/p) are, respectively: pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ F ðD; ri ÞvD /D ¼ dD ð45Þ cos x ﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ F ðD; ri Þri p /p ¼ pﬃﬃﬃﬃﬃﬃ : dei ð46Þ sin x vD In the above expressions, x is the parameter governing the coupling between damage and plasticity, while F⁄(D, ri) is a function controlling the overall rate of dissipation (see [12,29] for more details). From the dissipation potential, the stress ri and damage energy vD are: @ Ui @ Ui @/p /p @/p ri ¼ ¼ ¼ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð47Þ @dep @/p @depi 2 / þ/ 2 @de p i D p @ Ui @ Ui @/D /D @/D vD ¼ ¼ ¼ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð48Þ @dD @/D @dD 2 / þ/ 2 @dD D p Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx 11 The loading function y⁄ in mixed stress-energy space (vD and ri) is obtained as a result of a degenerate Legendre trans- formation of the dissipation potential (44): !2 vD r : ri i y ¼ @/D þ @/p @/ 160 ð49Þ @dD p @dei : @depp i The above loading function gives the following ﬂow rules (with dk being a non-negative multiplier having stress/energy unit in this model): @y cos2 x dD ¼ dk ¼ 2dk ð50Þ @ vD F ðD; ri Þ 2 @y v sin x dep ¼ dk ¼ 2dk D ri ð51Þ @ ri F ðD; ri Þri : ri In pure stress space, the loading function y is obtained by substituting (45) and (46) into (49): ri : C : ri y ¼ vD F ðD; ri Þ ¼ F ðD; ri Þ 6 0 ð52Þ 2ð1 DÞ2 As shown in Nguyen and Korsunsky [12], from the above loading function, the symmetry of response in tension and com- pression in stress space can be overcome to capture correctly the weaker behaviour in tension in quasi-brittle materials. In such a case, function F⁄(D, ri) can be chosen as follows: ri : Di : ri FðDÞ F ðD; ri Þ ¼ rþi :Di :rþi ð53Þ 2ð1 DÞ2 2ð1DÞ2 where r+ denotes the positive part of a stress tensor r, which is decomposed into positive and negative parts using the eigen- value decomposition [43,44]: X 3 rþij ¼ hrm ipm m i pj ð54Þ m¼1 In the above expression, hi is the Macaulay bracket, pm is the unit vector of the mth principal direction and rm is the mth principal stress. Further details and some properties of the decomposition can be found in [43,44]. Therefore the loading function, after simpliﬁcation, becomes: rþi : Di : rþi y¼ FðDÞ 6 0 ð55Þ 2ð1 DÞ2 Function F(D) in the above loading condition governs the inelastic behaviour of the material model. Following our earlier works on modelling quasi-brittle failure [12,45], it can take the form: 2 ft02 E þ Ep ð1 DÞn FðDÞ ¼ n ð56Þ 2E Eð1 DÞ þ Ep ð1 DÞ where ft0 is the uniaxial tensile strength, and Ep and n are two parameters governing the damage evolution, the identiﬁcation and roles of which can be found in the above papers. 3.2. Tangent stiffness tensor The formulation of the continuum tangent stiffness aTi is presented in this section for the implementation in the next sec- tion. From (42), the stress can be written in incremental form as: ri ri dri ¼ ð1 DÞ D1 i : dei depi dD ¼ ð1 DÞ ai : dei depi dD ð57Þ ð1 DÞ ð1 DÞ Taking the consistency condition of the loading function, we can write: @y @y dy ¼ : dri þ dD ¼ 0 ð58Þ @ ri @D Substituting the stress increment (57), and ﬂow rules (50) and (51) in the above equation and solving it for dD, we obtain: ð1 DÞ @@yr : ai : dei dD ¼ i ¼ Pi : dei ð59Þ ð1 DÞ @@yri 1 : ai : Q i þ ð1DÞ @y @ ri @y : ri @D Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 12 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx where depi ri : Di : ri ri Qi ¼ ¼ tan2 x ð60Þ dD 2ð1 DÞ2 ri : ri and ð1 DÞ @@yr : ai P¼ i ð61Þ ð1 DÞ @@yri 1 : ai : Q i þ ð1DÞ @y @ ri @y : ri @D The incremental stress–strain relationship and tangent stiffness aiT can then be obtained by substituting (59) and (60) in (57): ( " # ) ri dri ¼ ð1 DÞ ai ðai : Q i Þ þ Pi : dei ¼ aTi : dei ð62Þ ð1 DÞ2 4. Numerical implementation 4.1. Stress return algorithm An algorithm is presented in this section for the stress update of the constitutive model, given the input total strain incre- ment de, sizes h of the localised band and H of the spatial discretisation, and the complete knowledge of the constitutive behaviour inside and outside the band (Fig. 1). The coupled stresses in the model require an algorithm different from that for classical constitutive model, in the sense that both stresses should be stored and updated at any time increment. While the steps for inﬁnitesimal increments described in Section 2, equivalent to the forward Euler algorithm, can also be used in the numerical implementation, its accuracy is an issue due to the ﬁnite step size of the increment in numerical analysis. An improvement using Runge–Kutta 4th order method is therefore developed here. The total stress increment dr is of the form: 1 f s dr ¼ ao : de ðn d½uÞ ð63Þ 1f h where the increment of displacement jump is obtained from (see Eq. (8)): d½u ¼ C1 ðao : deÞ n ð64Þ with C dependent on both sizes (h and H) and the corresponding behaviours of the material inside and outside the localisa- tion zone: 1 ð1 f Þ T C¼ Ao þ Ai ð65Þ H h While Ao = n ao n depends on the pristine elastic stiffness ao and therefore is constant in this case, ATi ¼ n aTi n is gov- erned by the tangent stiffness aTi which is a function of the stress ri and other internal variables of the inelastic model (see Eq. (62)). The solution to (64) is therefore implicit, due to the fact that this stress and corresponding internal variables are also dependent on ei, and hence on [u] (see Eq. (2)). Denoting C1 = F(ri ([u]),D([u])) = F([u]), and [u]0 the displacement jump at the previous converged step, we can rewrite (64) in the following form: d½u ¼ Fð½uÞ ðao : deÞ n ð66Þ A Runge–Kutta like solution for (66) is as follows: 1 d½u ¼ ðd½u1 þ 2d½u2 þ 2d½u3 þ d½u4 Þ ð67Þ 6 where d½u1 ¼ Fð½u0 Þ ðao : deÞ n ð68aÞ 2 0 1 d½u ¼ Fð½u þ d½u =2Þ ðao : deÞ n ð68bÞ 3 0 2 d½u ¼ Fð½u þ d½u =2Þ ðao : deÞ n ð68cÞ d½u4 ¼ Fð½u0 þ d½u3 Þ ðao : deÞ n ð68dÞ For a given d[u]k (k = 1..4), the strain increment deki of the material inside the localisation zone can be calculated from (2), as: Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx 13 Fig. 4. Performance of stress return algorithms. (the strain increment is mDe, with De = 5 106). 1 s deki ¼ _ kÞ ðn d½u ð69Þ h This is then followed by the update of stress rki , damage Dk and plastic strain epk i using an implicit stress return algorithm for classical constitutive models. Details on this algorithm can be found in our previous works (e.g. [46,47]). Once the ﬁnal d[u] is obtained from (67), the same stress return procedure is carried out to obtain the stress ri and damage D for the mate- rial inside the localisation zone. The overall stress increment dr is then calculated from (63). Of course there are several different and possibly better algorithms that we can implement and test. This is however out of the scope of this paper and is therefore left for a future investigation. The performance of the above Runge–Kutta (R–K) algorithm is satisfactory for our analysis of quasi-brittle failure in this paper. A comparison of the performance of the forward Euler and R–K like algorithm above is shown in Fig. 4. This is a tension test on concrete, using our proposed approach de- scribed in Section 2, and damage model presented in Section 3. The model parameters are the same as those for the mixed mode test in Example 5.2. As can be seen the R–K like algorithm performs much better than the Euler one and is almost insensitive to the size of the increment. 4.2. Interface with ﬁnite elements The last two points we want to address in the implementation are the determination of the normal n of the localisation band, and the discretisation size H involved in the formulation. Generally for soils, the loss of positive deﬁniteness of the acoustic tensor ATi can be used for the determination of the band orientation. In such cases the normal n is based on the min- imum determinant of ATi (e.g. [29,30]). For quasi-brittle fracture of concrete and rocks, the orientation determined by the ﬁrst principal stress is appropriate (e.g. [23,48]) and hence used in this study. In addition, for small strain assumption in this study, the rigid body rotation of the localisation band is not taken into account. The localisation band with speciﬁed width h is ‘‘inserted’’ to the description of constitutive behaviour once the onset of inelastic deformation governed by the criterion (55) is met. In the interface with ﬁnite elements, this orientation of the localisation band is required to determine the effec- tive size H to be used in the constitutive equations. We take a look at the following situation of a 2D ﬁnite element of area X with an embedded FPZ of width h (Fig. 5). While a single integration point in this situation brings no issue, as the sizes H and h can be straightforwardly ﬁgured out from Fig. 1, the use of multiple integration points to represent a crack requires a more systematic way to interface with ﬁnite ele- ments. At a ﬁrst sight, assigning a fraction of the width h to individual integration points seems to be a solution. This how- ever suffers from the fact that there is no general way to calculate this fraction of h, given different integration rules and generally different ﬁnite element types. More seriously in such a case, h, a physical quantity, is no longer independent from the discretisation. A solution to this is to use an effective size H determined from the ﬁnite element, not from the area/volume occupied by the integration points. The width h is ﬁxed in this case. The orientation of the FPZ is at ﬁrst used to determine the effective length L (Fig. 5) of the crack at an integration point, assuming the traction and displacement jump [u] are uniform along this Fig. 5. Determination of the effective discretisation size H. Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 14 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx Fig. 6. Insensitivity of the response with respect to the number of integration points. length. The effective size H can then be determined from L and the area X of the ﬁnite element, as depicted in Fig. 5. Brieﬂy speaking, in the case of multiple integration points per element, the constitutive behaviour scales with the size H associated with the element (see Eq. (39b)), not with its integration point. This helps keep the same h for all integration points inside the element and guarantees the correct amount of energy dissipation. A similar derivation of the effective size H, but applied to a totally different approach, can be found in Schreyer et al. [49]. A test using a single 2D linear quadrilateral ﬁnite element (Fig. 6) shows the effectiveness of the proposed algorithm. The response of the element in this case is insensitive to the num- ber of integration points used. A simpler alternative way to interface with ﬁnite element is to use an effective H determined from the area (AFE) or vol- pﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃ ume (VFE) of the ﬁnite element, e.g. H ¼ AFE in 2D or H ¼ 3 V FE in 3D. This makes the implementation in commercial FE code easier, as the user deﬁned material routine does not need to know anything about the FE, except a length representing the element size. Of course the aspect ratios of FEs in this case must not be so irregular to give an appropriate measure of the effective size H. Numerical examples in Section 5 will demonstrate the effectiveness of proposed implementation algorithms. It is also noted that these proposed effective size H are just ﬁrst degree approximations for interfacing the new constitutive modelling scheme with ﬁnite elements. Large deformation issues will need more complicated algorithms to take into ac- count the rigid body rotation of the band together with distortion of ﬁnite elements. These have not been taken into account yet in the current developments. Although the results are satisfactory for benchmark tests (see Section 5 below), further investigations are still required for a better way to couple this proposed constitutive modelling scheme with the spatial discretisation. 5. Numerical examples The above approach, model and algorithms are implemented in our in house numerical code based on the Material Point Method (MPM; [50,51]). The use of the MPM, in its quasi-static version, at the moment is just for convenience in the devel- opment and the research activities of our group. The proposed model can be implemented in any numerical methods de- voted to solving BVPs in solid mechanics. For the quasi-brittle fracture tests under small strain assumptions in this paper, there is no signiﬁcant difference between the results from the Finite Element Method (FEM) and the MPM, as the material points in the latter almost stay in the same elements during the failure process. In all examples, linear quadrilateral elements and structured meshes are employed for the spatial discretisation. The tolerances for the stress return algorithm and the incremental-iterative scheme for global equilibrium based on Newton–Raphson method are set to 105 and 104, respec- tively. We note here that for the use of continuum tangent stiffness (62), quadratic convergence is not achieved and the iter- ation scheme is in fact quasi-Newton method. The numerical examples are to highlight the insensitiveness of the numerical solutions with respect to the resolution of the discretisation, the performance of the numerical implementation algorithms and the anisotropy introduced by the localisa- tion band. The anisotropy in this case is induced by an oriented localisation band embedded in the constitutive behaviour. While the inelastic model inside the localisation band (Section 3) is isotropic, the overall constitutive response is anisotropic thanks to the orientation dependent of this band. In all examples, pure damage behaviour corresponding to x = 0 is employed, as the coupling between damage and plasticity in the underlying inelastic model (Section 3) is not a focus in this case. 5.1. Three-point bending test This is the experiment carried out by Petersson [52]. The geometry of the beam and loading condition are depicted in Fig. 7. The following material properties are obtained from the experiment [52]: Young modulus E = 30,000 MPa, Poisson’s ratio m = 0.2, uniaxial tensile strength ft0 ¼ 3:33MPa, and fracture energy GF = 0.124 Nmm/mm2. The procedures proposed in [12,45] are used to determine model parameters corresponding to the above material properties. We note that in the case of fracture modelling, the width h of the localisation band is related to the fracture energy GF of the material. In other words, both h and the constitutive behaviour in the localisation zone affect the amount of dissipation and should be properly determined from the fracture property of the material (see [12,45] for details). Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx 15 Fig. 7. Three point bending test: geometry and loading. For the given GF above, two different values of h are used here to obtain two corresponding sets of model parameters, and to show the invariance of the numerical solution (Table 1). In the numerical analysis, thanks to symmetry, only half of the beam is modelled and the fracture energy corresponding to this conﬁguration is half the original one, to avoid doubling the total energy dissipation (due to 2 cracks modelled). Three different ﬁnite element meshes (coarse, medium and ﬁne) consisting of uniform 4-node quadrilateral elements are used, with element sizes being 20 10 mm2, 10 ⁄ 5 mm2, and 6.7 ⁄ 3.3 mm2, respectively. The results can be seen in Fig. 8a, showing the insensitiveness of the numerical solutions with respect to mesh density. A comparison with an isotropic damage model (Section 3) using smeared crack regularisation is also depicted (Fig. 8b). The observed stress locking of the proposed approach comes directly from its anisotropy. In this case, the ﬁrst principal stress rotates from the horizontal direction when moving away from the axis of symmetry. As a consequence, the crack orienta- tions, determined from the direction of the ﬁrst principal stress, deviate from the vertical one due to the effects of shear stress at points off the axis of symmetry. This facilitates the transfer of horizontal stress rxx across the cracks even at high damage level. In contrast, for isotropic damage models, the whole stress tensor is affected by damage, resulting in no locking. Similar numerical observations have been documented in [53]. In contrast with efforts to remove the locking at the end of the failure process by switching to isotropic model (e.g. [53]), we believe that allowing secondary cracking in the model is a physically better strategy. In reality this locking should be relieved due to the activation of new cracks, but the numerical prediction currently does not allow this activation because it takes into account only a single oriented FPZ at a material point (see Fig. 1). The remediation of this stress locking issue can be done by more advanced algorithms to take into account more than one cracks at a material point. This is a subject of on-going research. Fig. 9a highlights the inter-dependency of h and the constitutive response to reproduce the correct amount of dissipation. The objectiveness of the solutions with respect to the numerical integration scheme, e.g. the number of material points (MPs) per element, is illustrated in Fig. 9b. This conﬁrms the effectiveness of the algorithm to determine the discretisation size H in Section 4b. 5.2. Mixed mode test To examine further the responses of the proposed approach and model in the numerical analysis of quasi-brittle failure, the mixed mode cracking of a double edge notched (DEN) specimen is numerically simulated in this example. The experi- mental tests of the DEN specimens examined were carried out by Nooru-Mohamed et al. [54] and the corresponding data are extracted here. The geometrical data of the specimen along with the boundary conditions are shown in Fig. 10 above. Only one loading path was considered in this example: biaxial loading (path 2a, [54]), in which the axial tensile P and lateral compressive shear load Ps were applied to the specimen so as to keep the ratio dn/ds unchanged throughout the test (dn/ds equals to 1.0 in load path 2a). The average deformations used to control the loading process are calculated as follows 0 0 0 [54]: dn ¼ ðdM dM þ dN dN Þ=2, and ds ¼ dPs dPs . However in the numerical test, incorporation of the above deformation controls is difﬁcult [55] and an alternative way is adopted in this study. Following Jefferson [55], the rigid movements of Table 1 Model parameters for three-point bending test. Parameters (see Eq. (56)) h = 0.1 mm h = 0.5 mm Ep (MPa) 35.813 179.926 n 0.149 0.179 Table 2 Model parameters for mixed-mode test. Parameters (see Eq. (56)) h = 0.05 mm h = 0.1 mm Ep (MPa) 8.184 16.372 n 0.128 0.137 Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 16 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx Fig. 8. (a) Mesh independent response. (b) Comparison (using medium mesh) with isotropic model of Section 3 enhanced by smeared crack regularisation. Fig. 9. Insensitiveness of the numerical solution (using medium mesh) with respect to: (a) different values of width h and (b) number of material points per element. the upper left edge and top edge of the specimen are prescribed and used as ds and dn, respectively, while the lower right edge and bottom edge are kept ﬁxed in both directions. The material properties provided by experiments (E ¼ 32; 000 MPa;m ¼ 0:2; ft0 ¼ 3:0 MPa; and GF = 0.11 Nmm/mm2) result in the model parameters listed in Table 2. We use three different uniform meshes with ﬁnite element sizes of 2.5 ⁄ 2.5 mm2 (coarse), 1.25 ⁄ 1.25 mm2 (medium) and 0.625 ⁄ 0.625 mm2 (ﬁne). The numerical results in Fig. 11 follow the trend observed in experiments. Under prediction of the elastic stiffness is seen, as a result of using the simpliﬁed displacement control described above. The over prediction of the maximum load and post-peak branches are mostly due to the use of a simple constitutive model in this framework. This is however also a trend in other previous numerical studies (e.g. [55,56]) using totally different models and approaches. It is noted that the simple damage model presented in Section 3 is used to complete the description and illustration of fea- Fig. 10. Double-edge notched test: geometry and loading. Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx 17 Fig. 11. Convergence upon mesh reﬁnement. (using h = 0.05 mm; 9 material points per element). Fig. 12. Performance of the proposed approach (anisotropic) and the isotropic model of Section 3 (using coarse mesh, h = 0.05 mm; 9 material points per element). Fig. 13. Failure pattern. tures of the proposed scheme. More advanced constitutive models for the material behaviour inside the localisation band can be developed to yield better matches with experiments. This issue of a particular constitutive model is however not the focus of this paper. The convergence of the numerical results upon reﬁning the discretisation is shown in Fig. 11. A comparison with isotropic damage model of Section 3, enhanced by smeared crack regularisation, is also carried out and shown in Figs. 12 and 13. It is clear that isotropic models fail to capture the correct response in this mixed mode test and this observation is consistent with our previous results using nonlocal isotropic damage model [47]. In the proposed approach, although the underlying inelastic model is isotropic, the anisotropy introduced by the oriented FPZ and the use of coupled stresses help capture correctly the stiffening shear response. Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 18 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx Fig. 14. Insensitiveness of the results with respect to h and number of MPs per element (while GF is kept ﬁxed). Fig. 15. Insensitiveness of the results with respect to 2 different ways to determine effective size H of ﬁnite elements. The complex H is based on Fig. 5, pﬃﬃﬃﬃﬃﬃﬃ while the simple one is H ¼ AFE (using medium mesh with 9 MPs/element andh = 0.05 mm). The performance of the implementation algorithms is demonstrated in Figs. 14 and 15. We observe the independence of the numerical solutions with respect to integration rule and two different ways to determine effective size H of ﬁnite ele- ments. At the structural scale, the relationship between the width h of the FPZ and the inelastic behaviour inside the FPZ (Section 2.3) is also respected, resulting in the insensitiveness of the numerical solution with respect to h since the fracture energy GF is kept ﬁxed. 6. Conclusions We addressed the issues of correctly capturing localised failure in quasi-brittle materials with a new constitutive mod- elling approach. A length scale related to the width of the localisation band is introduced in the development of constitutive models to yield a material model that captures the size effect. The derived constitutive models in this case possess two stres- ses/strains corresponding to the behaviour in the localisation zone and the bulk elastic. These stresses are coupled via an internal equilibrium equation enforcing the traction continuity across the boundary of the localisation band. We believe that this is a simpler, more direct and physical way to tackle localised failure at the constitutive level, compared to existing ap- proaches using smeared crack and nonlocal/gradient regularisation. Implementation algorithms were developed to provide an interface with the existing structure of ﬁnite element codes. The numerical results in this study show the promising fea- tures of the new approach. This is only an initial step towards a better strategy for the modelling of material failure. The localisation zone has been idealised as a thin band across the volume, although reality is much more complicated, e.g. FPZ at the continuum scale in fact consists of many micro-cracks with different orientations. Further developments are under way to take into account this nested scale nature of localised failure for constitutive model with better predictive capabilities. Acknowledgements This research was supported under the Australian Research Council’s Discovery Projects funding scheme (project number DP1093485). Support from the University of Sydney through the Research Support Grant 2013-00053 for Giang Nguyen is also gratefully acknowledged. Comments from Drs. Irene Guiamatsia and Rachel Gelet in revising the manuscript are appreciated. Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx 19 References [1] Jirásek M. Objective modeling of strain localization. Revue française de génie civil 2002;6:1119–32. [2] Jirásek M, Bazant ZP. Inelastic analysis of structures. John Wiley & Sons; 2002. [3] Cedolin L, Bazant ZP. Effect of ﬁnite element choice in blunt crack band analysis. Comput Methods Appl Mech Engng 1980;24:305–16. [4] Vardoulakis I. Shear-banding and liquefaction in granular materials on the basis of a Cosserat continuum theory. Ing Arch 1989;59:106–13. [5] Bazant ZP. Why continuum damage is nonlocal: micromechanics arguments. J Engng Mech 1991;117:1070–87. [6] Chen Z, Schreyer HL. Simulation of soil–concrete interfaces with nonlocal constitutive models. J Engng Mech 1987;113:1665–77. [7] Pijaudier-Cabot G, Bazˇant ZP. Nonlocal damage theory. J Engng Mech 1987;113:1512–33. [8] De Borst R, Mühlhaus H-B. Gradient-dependent plasticity: formulation and algorithmic aspects. Int J Numer Meth Engng 1992;35:521–39. [9] Peerlings RHJ, De Borst R, Brekelmans WAM, De Vree JHP. Gradient-enhanced damage for quasi-brittle materials. Int J Numer Meth Engng 1996;39:3391–403. [10] Jirásek M. Nonlocal models for damage and fracture: comparison of approaches. Int J Solids Struct 1998;35:4133–45. [11] Hien Poh L, Swaddiwudhipong S. Over-nonlocal gradient enhanced plastic-damage model for concrete. Int J Solids Struct 2009;46:4369–78. [12] Nguyen GD, Korsunsky AM. Development of an approach to constitutive modelling of concrete: isotropic damage coupled with plasticity. Int J Solids Struct 2008;45:5483–501. [13] Jirásek M, Bazant ZP. Particle model for quasibrittle fracture and application to sea ice. J Engng Mech 1995;121:1016–25. [14] Sulsky D, Schreyer H, Peterson K, Kwok R, Coon M. Using the material-point method to model sea ice dynamics. J Geophys Res 2007;112. [15] Borja RI. A ﬁnite element model for strain localization analysis of strongly discontinuous ﬁelds based on standard Galerkin approximation. Comput Meth Appl Mech Engng 2000;190:1529–49. [16] Foster CD, Borja RI, Regueiro RA. Embedded strong discontinuity ﬁnite elements for fractured geomaterials with variable friction. Int J Numer Meth Engng 2007;72:549–81. [17] Larsson R, Runesson K, Sture S. Embedded localization band in undrained soil based on regularized strong discontinuity—theory and FE-analysis. Int J Solids Struct 1996;33:3081–101. [18] Oliver J. Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 1: fundamentals. Int J Numer Meth Engng 1996;39:3575–600. [19] Wells GN, Sluys LJ, de Borst R. Simulating the propagation of displacement discontinuities in a regularized strain-softening medium. Int J Numer Meth Engng 2002;53:1235–56. [20] Borja RI. Assumed enhanced strain and the extended ﬁnite element methods: a uniﬁcation of concepts. Comput Methods Appl Mech Engng 2008;197:2789–803. [21] Samaniego E, Belytschko T. Continuum–discontinuum modelling of shear bands. Int J Numer Meth Engng 2005;62:1857–72. [22] Sanborn SE, Prévost JH. Frictional slip plane growth by localization detection and the extended ﬁnite element method (XFEM). Int J Numer Anal Meth Geomech 2011;35:1278–98. [23] Wells GN, Sluys LJ. A new method for modelling cohesive cracks using ﬁnite elements. Int J Numer Meth Engng 2001;50:2667–82. [24] Belytschko T, Fish J, Engelmann BE. A ﬁnite element with embedded localization zones. Comput Methods Appl Mech Engng 1988;70:59–89. [25] Sluys LJ, Berends AH. Discontinuous failure analysis for mode-I and mode-II localization problems. Int J Solids Struct 1998;35:4257–74. [26] Garikipati K, Hughes TJR. A study of strain localization in a multiple scale framework—the one-dimensional problem. Comput Methods Appl Mech Engng 1998;159:193–222. [27] Garikipati K, Hughes TJR. A variational multiscale approach to strain localization – formulation for multidimensional problems. Comput Methods Appl Mech Engng 2000;188:39–60. [28] Huespe AE, Needleman A, Oliver J, Sánchez PJ. A ﬁnite thickness band method for ductile fracture analysis. Int J Plast 2009;25:2349–65. [29] Nguyen GD, Einav I, Korsunsky AM. How to connect two scales of behaviour in constitutive modelling of geomaterials. Géo Lett 2012;2:129–34. [30] Das A, Nguyen GD, Einav I. Compaction bands due to grain crushing in porous rocks: a theoretical approach based on breakage mechanics. J Geophys Res 2011;116:B08203. [31] Comi C, Mariani S, Perego U. An extended FE strategy for transition from continuum damage to mode I cohesive crack propagation. Int J Numer Anal Meth Geomech 2007;31:213–38. [32] Moës N, Stolz C, Bernard PE, Chevaugeon N. A level set based model for damage growth: the thick level set approach. Int J Numer Meth Engng 2011;86:358–80. [33] Hill R. Elastic properties of reinforced solids: some theoretical principles. J Mech Phys Solids 1963;11:357–72. [34] Vardoulakis I, Goldscheider M, Gudehus G. Formation of shear bands in sand bodies as a bifurcation problem. Int J Numer Anal Meth Geomech 1978;2:99–128. [35] Kolymbas D. Kinematics of shear bands. Acta Geotechnica 2009;4:315–8. [36] Simo JC, Oliver J, Armero F. An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. Comput Mech 1993;12:277–96. [37] Bazˇant ZP, Pijaudier-Cabot G. Measurement of characteristic length of nonlocal continuum. J Engng Mech 1989;115:755–67. [38] Otsuka K, Date H. Fracture process zone in concrete tension specimen. Engng Fract Mech 2000;65:111–31. [39] Einav I. Breakage mechanics–part I: theory. J Mech Phys Solids 2007;55:1274–97. [40] Rubin MB, Einav I. A large deformation breakage model of granular materials including porosity and inelastic distortional deformation rate. Int J Engng Sci 2011;49:1151–69. [41] Lemaitre J. A course on damage mechanics. Springer Verlag; 1992. [42] Armero F. Large-scale modeling of localized dissipative mechanisms in a local continuum: applications to the numerical simulation of strain localization in rate-dependent inelastic solids. Mech Cohes-Frict Mat 1999;4:101–31. [43] Ladeveze P. Sur une theorie de l’endommagement anisotrope. Rapport interne No 34: LMT Cachan, France; 1983. [44] Ortiz M. A constitutive theory for the inelastic behavior of concrete. Mech Mater 1985;4:67–93. [45] Nguyen GD, Houlsby GT. Non-local damage modelling of concrete: a procedure for the determination of model parameters. Int J Numer Anal Meth Geomech 2007;31:867–91. [46] Nguyen GD, Einav I. A stress-return algorithm for nonlocal constitutive models of softening materials. Int J Numer Meth Engng 2010;82: 637–70. [47] Nguyen GD. A thermodynamic approach to constitutive modelling of concrete using damage mechanics and plasticity theory [DPhil thesis]. University of Oxford; 2005. [48] Alfaiate J, Simone A, Sluys LJ. Non-homogeneous displacement jumps in strong embedded discontinuities. Int J Solids Struct 2003;40: 5799–817. [49] Schreyer HL, Sulsky DL, Munday LB, Coon MD, Kwok R. Elastic-decohesive constitutive model for sea ice. J Geophys Res: Oceans 2006;111:C11S26. [50] Sulsky D, Chen Z, Schreyer HL. A particle method for history-dependent materials. Comput Meth Appl Mech Engng 1994;118:179–96. [51] Sulsky D, Zhou S-J, Schreyer HL. Application of a particle-in-cell method to solid mechanics. Comp Phys Commun 1995;87:236–52. [52] Petersson PE. Crack growth and development of fracture zones in plain concrete and similar materials. Report TVBM-1006: Div. of Build. Mat. Lund Institute of Technology, Lund, Sweden; 1981. [53] Jirásek M, Zimmermann T. Rotating crack model with transition to scalar damage. J Engng Mech 1998;124:277–84. Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006 20 G.D. Nguyen et al. / Engineering Fracture Mechanics xxx (2013) xxx–xxx [54] Nooru-Mohamed MB, Schlangen E, van Mier JGM. Experimental and numerical study on the behavior of concrete subjected to biaxial tension and shear. Advan Cem Based Mater 1993;1:22–37. [55] Jefferson AD. Craft—-a plastic-damage-contact model for concrete. II. Model implementation with implicit return-mapping algorithm and consistent tangent matrix. Int J Solids Struct 2003;40:6001–22. [56] Zi G, Rabczuk T, Wall W. Extended meshfree methods without branch enrichment for cohesive cracks. Comput Mech 2007;40:367–82. Please cite this article in press as: Nguyen GD et al. A constitutive modelling framework featuring two scales of behaviour: Fundamentals and applications to quasi-brittle failure. Engng Fract Mech (2013), http://dx.doi.org/10.1016/j.engfracmech.2013.11.006

(PDF) A constitutive modelling framework featuring two scales of behaviour: fundamentals and applications to quasi-brittle failu