This paper studies the use of the Embedded Finite Element Method (E-FEM) for the numerical modelling of triaxial fracture processes in non-homogeneous quasi-brittle materials. The E-FEM framework used in this study combines two kinematics enhancements: a weak discontinuity allowing the model to account for material heterogeneities and a strong discontinuity allowing the model to represent local fractures. The strong discontinuity features enriched fracture kinematics that allow the modelling of all typical fracture modes in three dimensions. A brief review is done of past work using similar enriched nite element frameworks to approach this problem. The work continues by establishing the theoretical basis of each kind of discontinuity formulation and their superposition through the Hu-Washizu variational principle. Afterwards, two groups of simulations have been done for discussing the performance of this combined E- FEM model: homogeneous simulations and simple heterogeneous simulations. Simple homogeneous material simulations aim to test the capabilities of the strong discontinuity model featuring full 3-D kinematics. Simple heterogeneous simulations show numerical applications of the model to the problem of a single spherical inclusion embedded into a homogeneous matrix. Comparisons will be made with another E-FEM model considering a single local fracture mode approach to discuss the di erences on the representation of fracture physics under all explored conditions. A concluding statement is made on the bene ts and complications identi ed for the E-FEM framework in this kind of applications. Keywords: E-FEM modelling, embedded weak discontinuity, embedded strong discontinuity, enhanced nite elements, incompatible modes, local heterogeneity modelling, local fracture, triaxial fracture, nite element enhancements, mesoscale analysis 1. Introduction The numerical study of fracture phenomena in composite materials requires a detailed consideration of their heterogeneous structure, which contributes to the emergence of complex and unexpected mechanical behaviours. Multiscale analysis approaches have been devised in recent decades to develop mathematical models capable of capturing and predicting their response [1]. Classical approaches are based on homogenisa- tion principles, where it is assumed that it is possible to represent a given composite material as a completely homogeneous domain whose mechanical behaviour is governed by a suciently complex material law taking parameters from studies at smaller scales. As an example, the representative volume element (RVE) is a widely spread approach that has been successfully used both in academic research and industry [2, 3]. Any of these homogenisation approaches require the analysis of composite materials at smaller scales where the e ects of heterogeneities can be assessed in accurate manner, whether through their spatial layout, their material phase individual properties or possible local fractures that might signi cantly in uence the behaviour within a limited domain. It also requires signi cant statistical analysis to ensure that the speci c smaller scale material domain used to pull out the homogenisation procedure is suciently representative with respect to natural variations happening on the heterogeneous structure of the given composite. 1 It is clear that any multiscale approach bene ts from having an ecient yet suciently accurate numerical analysis method at the smaller scales to ful ll all homogenisation requirements in a reasonable time frame. On this matter, the conventional nite element method, being one of the most prominently used techniques for the mechanical analysis of homogeneous materials, presents considerable challenges both in the problem of heterogeneity distributions and in the explicit modeling of local fractures. On one hand, a mesh has to be adapted to represent the regions having di erent material phases. Depending on the natural shapes of these phases and the kind of information that can be measured about them on a 3D setting, intensive adaptive meshing techniques might be required to ensure mesh continuity through these regions while maintaining an accurate depiction of all material interfaces. Even if such mesh con gurations are feasible using currently available meshing software, the mere act of meshing these intricate structures would also bring the question of the resulting quality of the element domains. Element distortions and any other kind of geometrical anomalies will have an impact on the quality of results coming from the numerical analysis as a whole [4]. Excessive mesh density requirements to avoid bad quality issues will translate into long computational solving times. On the other hand, the modelling of a fracture, even if referring to a completely homogeneous domain, already poses numerous critical challenges for a standard nite element approach [5]. The most evident is that the path of the fracture should be explicitly identi ed and delimited beforehand on the mesh itself. Spontaneous failure calculations can also be implemented, but the standard methods do not possess any mathematical features capable of modelling fracture separation and/or propagation other than complete element sti ness damage mechanisms or even direct element removal under certain failure criteria. While these techniques have gained acceptance in both academy and industry for their ease of use, they are well known to produce considerable solution stability problems and modelling limitations [6, 7], deviating predictions from a realistic crack opening and propagation behaviours. As a nal obstacle, the statistical analysis needed for robust homogenisation processes will require considering all these requisites for a large number of numerical simulations at once, only increasing computational demands. In this context, quasi-brittle materials pose speci c modelling challenges that must also be addressed by the numerical approaches proposed to face this analysis problem. The most prominent issue is the com- plexity of their fracture processes featuring multiple cracks happening in a simultaneous fashion, collectively exhibiting a network which, by itself, has clearly identi able spatial characteristics and evolution behaviour depending on both large scale and small scale properties. What in a large scale would be called a de ned fracture , at smaller scales would rather be a complex combined behaviour of multiple individual local cracks giving form to a dominant macro-crack in the material through a coalescence process. This is especially true if the material presents an intricate heterogeneous structure, where material phase interfaces are well-known sources of local stress concentrations that end up producing these individual crack networks. In order to capture this kind of behaviour, the modelling of local fractures should be detailed enough to reveal this fracture kinematics and load triaxiality dependencies. While there is evidence that it suces to model a basic set of local failure modes to accomplish this task [8], local cracks still need considerable exibility to capture fundamental fracture kinematics such as normal separation, tridimensional sliding, closure following sustained compression and even frictional models, among other e ects still under research [9]. Advanced nite element approaches provide with attractive solutions to overcome the aforementioned problems on local scale material modelling, proposing integrating frameworks capable of representing both material heterogeneities and enriched local cracks using non-adaptive meshes. Among them we can nd the Extended Finite Element Methods (X-FEM) [10], the Base Force Element Method (B-FEM) [11] or the Em- bedded Finite Element Methods [12], among others. Many of them are well-endowed with a mathematical exibility capable of modelling di erent kinds of discontinuities to account for a variety of material phases and local fracture modes. They do this through mathematical enhancements on their supporting functions whether in a nodal base (X-FEM), through their elemental mechanical elds (E-FEM) or through new 2 mechanical state variables (B-FEM). Some of them are even mesh-less, in the sense that their supporting functions do not depend on a de nite division of elemental domains in space. The representation capabil- ity of these approaches is greatly enhanced at the expense of increasing operational and implementation complexities. This is due to their deeper mathematical de nitions. Therefore, a full heterogeneity/fracture modelling scheme in 3-D is certainly harder to implement than with any standard approach and therefore harder to track in recent literature. In this sense, the E-FEM framework, being the choice of the authors of this work for the current study, retains a reasonable balance between mathematical complexity and representational capability. This is mainly because it is based on the method of incompatible modes [13, 14], allowing the building of inde- pendent mechanical eld enhancements on each element without ensuring rigorous global continuity. From an operational standpoint, this allows working all heterogeneity and fracture mathematical enhancements in an internal-element fashion, with the possibility to use operator-split methods [15] and thus condensing all element internal e ects before attempting a global displacement solution step, whether linear or nonlin- ear. The E-FEM framework can then be developed in such a way that the global FEM numerical solution engine may be left untouched, enabling the support of a variety of available FEM solution platforms for implementation. Out of the domain of nite element approaches, work has been done on developing many other techniques to face the modelling problem explored in this work, such as discrete elements [16], phase eld modelling [17], strain injection techniques [18] or still standard nite element methods with highly adapted meshes and cohesive models [19]. These techniques have already matured to the point of approaching 3D domains, but there is still lack of some required qualities for the mechanical solution sought by the authors of the current work like, for instance, the representation of objective continuum stress elds for the case of discrete elements. The authors of this work have chosen to focus on the use of the E-FEM framework for approaching the problem of triaxial failure in quasi-brittle materials. In Section 2, a review is made on how the E-FEM framework has recently evolved to approach this problem. In Section 3, the theoretical basis of the E-FEM framework used for this study is described. It makes the integration of a weak and a generalized strong discontinuity model. No detailed development of each of the discontinuity formulations will be done as the bases have already been discussed in other works [12, 20, 21, 22, 23]. Only the essential information characterising the outstanding features of each model will be introduced, as well as the implications of their superposition under the light of the Hu-Washizu variational framework. Section 4 will introduce the numerical simulations done using this generalized E-FEM approach to discuss its performance in di erent problem scenarios allowing to scrutinise each of its enhancements in progressive fashion. Two di erent kinds of simulations will be presented in this work: homogeneous and simple heterogeneous. Homogeneous simulations will show how the strong discontinuity component of this E- FEM framework manages to capture stress triaxiality dependence in a controlled fracture path. Afterwards, simple heterogeneous simulations will feature a single spherical inclusion within an homogeneous matrix material considering di erent load conditions to discuss on the mixed 3D fracture processes unfolding due to the presence of di erent material phases on a simple layout. A comparison is made with respect to another E-FEM model considering a single fracture kinematic mode approach to make an assessment on the resulting behaviour of three-dimensional fracture processes emerging in this problem. Section 5 will close this work with a nal word on the overall e ectiveness of the E-FEM framework, with potential future works. 2. Evolution of the E-FEM framework approaching composite quasi-brittle failure problems From the founding works that established the main base of the E-FEM framework [24, 13], this analysis approach was originally conceived to model localisation phenomena in general materials as an alternative to 3 standard techniques with adapted meshes. It started with the idea of using weak discontinuity enhancements to represent the presence of shear localisation bands, associating speci c damage behaviour laws to all domains falling between two parallel lines having an arbitrary separation (shear band thickness). As the mathematical depth of these developments evolved, the framework began to turn towards the use of strong discontinuities equipped with discrete or regularized crack behaviour laws to gain objectivity in the de nition of strain localisation regions [25, 26]. Applications to quasi-brittle materials began to emerge at this stage, but only for one or two-dimensional problems [27, 28]. The pioneering work by Wells and Sluys [23] started with a full deployment of the approach for 3D problems, pushing the boundaries of the framework and evidencing new theoretical weaknesses in the mathematical structure of strong discontinuity enhancements. They managed to implement the E- FEM framework for a linear tetrahedron including the rigid body displacement crack kinematic modes of normal separation and sliding, making the comparison between variationally symmetric and asymmetric strong discontinuity enhancements. The numerical reproduction of classical concrete fracture test setups in three dimensions like double-notched and single-edge-notched beams were achieved by using these models considering a homogeneous material domain. In all these applications, the local fracture interface is de ned as a plane having from one to three translational degrees of freedom. At the same time, the works on improving mathematical robustness and kinematic consistency in the framework are sustained by the notable work of Jirasek and Oliver [26, 25, 29] in the early 2000s. Non- symmetrical variational schemes were favored to achieve consistent modelling of both statics and kinematics simultaneously. An enrichment of the fracture formulation was explored by implicitly incorporating a ro- tational degree of freedom in the kinematics of the local fracture plane of a constant stress triangle (CST) element [30], establishing the basis for non-uniform and enriched local fracture kinematics. The idea of having fully local enriched fracture modes for improving the consistency of fracture kinematics was later developed in detail on a 2D setting [22]. The authors managed to fully equip fracture kinematics with translational and rotational degrees of freedom in addition to the capability for one of the segregated do- mains of the element to have a simple lateral tension/compression mode. This allowed to reach new levels of variational consistency and better element internal equilibrium conditions, while having internal fracture variables that really described a physically meaningful state of the local fracture. The idea was further used in [31, 32, 33, 34, 35] and nally in the works developed in [20], where a deep theoretical assessment was made for its generalisation on a 3D setting. The present work considers this last development. Concerning the applications to heterogeneous quasi-brittle fracture, the interest of multiscale approaches and the eventual need to incorporate the E-FEM framework into homogenisation procedures began with the notable works by Markovic [36]. In particular, it was Markovic who set the point of departure of heterogeneity modelling by reviving the e ective use of the weak discontinuity formulation taking the foundations of previous works on the X-FEM approach [37]. This line of research gave rise to the rst successful integration of strong and weak discontinuities within the framework for one-dimensional bar elements in [38, 15]. In these works, an idealized 3D heterogeneity layout composed by perfectly spherical inclusions within a homogeneous matrix was used for the modelling of concrete samples. The use of bar elements greatly simpli ed the mathematical complexity inherent to the strong and weak discontinuity enhancements, avoiding all kinematic and variational consistency complications already found in past works. In this approach, a small cubic concrete sample was portrayed as a 3D truss made up of 1D bars, a part of which would be enhanced with weak discontinuity functions to represent material interfaces when these are intersected by borders of di erent material phases. The other pure elements would eventually activate a strong discontinuity enhancement if localisation was attained by means of a given criterion, for then following a discrete traction-separation law. Elements already enhanced with the weak discontinuity are also equipped with the strong one at the same location. While this method was successful in predicting typical tensile resistance values for concrete, the disadvantage was that it did not return satisfactory values for tension/compression resistance ratios. The use of 1D elements also prevented having an objective perspective of local stress and strain eld distributions. 4 This line of research began its rst truly 3D approach in the work done by Roubin et al. [12]. The 3D generalisation of the weak discontinuity enhancement model for heterogeneities was based on the works of Markovic [36], while a strong discontinuity model equipped with a discrete traction separation law was used for local fractures. The strong discontinuity model was capable of representing a single fracture mode of normal separation, activated by a Rankine localisation criterion and followed by an exponential traction- separation law. The representation of heterogeneity distributions was evolved to make use of probabilistic excursion sets, which enhanced the random nature of material phases improving their packaging quality and conserving a degree of smoothness on their interfaces. Despite the limitation of using a single fracture mode kinematics approach, typical tension and compression resistance ratios were found to be more reasonable than the previous works with 1D bars. This was considered as relevant evidence for the hypothesis of having some quasi-brittle materials like concrete to locally exhibit predominant mode I fracture modes, collectively yielding to more complex fracture processes following the interfaces set by inclusion layouts. Having at this point approached 3D problems where material continuum can be modelled in a more direct fashion through 3D nite elements, one of the most discussed characteristics of these E-FEM models was the lack of cross-elemental continuity for the emerging cracks. As the mathematical foundations of these enhancements are intrinsically based on the method of incompatible modes [13, 14], there is no guarantee to keep continuity of a local crack to neighboring elements once one of them starts to develop. This leads to the simultaneous creation of multiple local cracks that, by the mere mechanism of sti ness damaging and internal force redistribution, will tend to favour a spontaneous coalescence phenomenon giving rise to a larger scale fracture without explicitly aiming for it. While nonlocal crack continuity within the E- FEM framework has already been well studied and implemented in multiple applications for homogeneous material simulations [39, 31, 40], the line of research currently discussed ([38, 12, 41, 8, 42]) intentionally keeps the locally discontinuous nature of the E-FEM models to favour this behaviour in the context of complex heterogeneous materials. This decision is supported by multiple studies that favour the hypothesis of multi-cracks at smaller scales contributing to the emergence of larger fracture processes in such materials [43, 44]. Hauseux [41] later achieved the simulation of fracture sliding modes (mode II) with the same framework in the context of simulations of excavations of geomaterials. The strong discontinuity model was still based on a single fracture mode, but it considered mode II local fracture kinematics activated by means of a Mohr- Coulomb localisation criterion, followed by a discrete exponential traction-separation law. Even if these excavation simulations were done considering homogeneous materials, Hauseux also experimented with the simulation of heterogeneous cubic samples, considering compressive loads under precon nement conditions. Heterogeneities in this case were still represented by means of mathematical idealisations, such as packed spheres containing di erent material phases having di erent size ratios. Typical Mohr-Coulomb resistance behaviour of some rock samples was captured with this model. The works of Stamati [45, 8] added a new dimension to the E-FEM framework experimental validation by means of test setups having in-situ measurements in concrete samples that allowed direct comparisons with numerical simulations at the mesoscale. The data acquisition included accurate 3D descriptions of the heterogeneity distributions (largest pores and aggregates) within the samples in question by means of X-ray tomography scans. A linear quasi-static displacement loading pro le is imposed in the sample, stopping the loading process at determined points to make a complete 360scan of it, assessing its entire state at the given load step. The process continues until breaking the sample completely, where a last scan is retrieved for capturing its terminal fracture distribution. No signi cant evolution on the mathematical structure of the E-FEM framework was registered at these stages. Simulations of these samples using the same single fracture mode E-FEM framework allowed performing very direct validations for the predictions of both tridimensional crack networks and global resistances. Finally, this review will mention the works done by Sun on including crack reclosure and healing within the discussed E-FEM models [46]. In such works, mode I and mode II strong discontinuity models were 5 further enhanced with full crack reclosure mechanisms. Comparisons were made between models having exclusively a mode I fracture modeling approach or a mode II. This work evidences that, while single local fracture mode formulations are able to successfully predict a variety of large-scale phenomena under speci c load conditions, these are restrained by a physical representation limit and will not be able to approach other local phenomenology that remains relevant for the modelling of tridimensional fracture processes in quasi-brittle heterogeneous materials. In order to have a fully functional E-FEM framework, the strong discontinuity formulation should be generalized to gain enough exibly to handle all local requirements at once, thus acquiring a more universal character on its application to all types of mechanical problems. The current work continues with this line of research, attempting such generalisation on the strong discontinuity model and its integration with a weak discontinuity model. The next section will detail the reasons that led the authors to revise in considerable detail the theoretical bases of each part of the models coming from the E-FEM framework. 3. Double enhancement E-FEM formulation with generalized fracture kinematic modes This section will start by providing a concise description of the weak and strong discontinuity displace- ment eld enhancements. This work will not present a detailed derivation of all mathematical operators and relations associated to these basic eld enhancement de nitions and related variational analysis within the framework. The reader can refer to [20] and [21] to have a more complete theoretical background on these topics. The fundamental modelling of an embedded discontinuity within a nite element starts by having a body representing an elemental domain ehaving boundaries @ is crossed by a surface  dthat establishes the presence of a strain (weak) or displacement (strong) discontinuity, partitioning the complete domain in two subregions +and . The subregion interface  dis characterized at each of its points by normal and parallel unit vectors ^n,^tand ^m, respectively. These are commonly referred to as the local frame on the E-FEM context. A typical scheme in two dimensions can be appreciated in Figure 1 Figure 1: Basic schematic of an embedded discontinuity (whether weak or strong) in 2D From now on, while all theoretical premises presented in this study are applicable to a wide variety of nite element geometries and con gurations, this study will mostly particularise all works to a speci c base 3D nite element: a linear tetrahedron. The interface  drepresenting whether the local element rupture or the division between di erent material phases is assumed to be planar, characterized by a constant and unique local frame ( ^n,^t,^m) and segmenting a given tetrahedron in two subvolumes VandV+. 3.1. Weak discontinuity formulation The weak discontinuity formulation within the E-FEM framework is used to model the presence of di erent material phases within a single element domain. It introduces a jump on the strain eld to account for the presence of di erent materials while still retaining displacement continuity (thus the reason of calling it weak discontinuity). Its basic construction assumes that a general displacement eld ucan be decomposed 6 as the sum of a homogeneous standard displacement uand an enhanced displacement eld euin the following general form: u=u+eu (1) where the enhancement euhas the role of adding a displacement slope discontinuity representing the change in material characteristics, i.e., the jump on the strain eld. A typical particularisation to a tetrahedron element can be appreciated in Figure 2, assigning di erent linear elastic material properties E;to each of the subdomains . Figure 2: Basic schematic of a weak discontinuity in 3-D for the modelling of material heterogeneities within a tetrahedral element The model taken for all developments and simulations presented in this study is the one typically used on the line of research followed in [12] and analysed in detailed in [21]: Weak discontinuity enhancement eu= ^n(xxd) ["]n^n+ ["]t^t+ ["]m^m (2)  =( +=V Vx2 + =V+ Vx2 (3) where ^n(xxd) is the normal distance from a given point xto the nearest location xdof the surface d. The internal variables [ "]n;["]t;["]mrepresent the magnitude of the strain jump when crossing  d. As studied in [21], this weak discontinuity formulation proposal entails a permanent variational incon- sistency depending on element geometry and its alignment with respect to the interface  d. However, if the model features a mesh with good aspect ratios ( <3:0) and if the ratio between linear elastic parameters falls within the order of 5 or less, it has been demonstrated that this formulation grants a reasonable represen- tation of the heterogeneous sti ness response and the expected strain eld distributions [21]. Considering this, the fact that the enhancement bene ts from unconditional numerical stability and that it also allows for the use of a completely symmetric solver, its use for the current study is justi ed. 3.2. Strong discontinuity formulation The role of the strong discontinuity in this work is to model the displacement jump associated to an element local fracture. It is de ned by having a general displacement eld uexpressed as a continuous base elduadded to a Heaviside function Hd, whose trigger location is that of the discontinuity surface  d. It is described as: u=u+Hd[juj]; (4) 7 where [juj] is a vector that describes the displacement jump in each direction, This study takes advantage of a non-uniform de nition of [ juj] eld components at the interface to represent fully enriched fracture kinematics describing the motion of the +domain with respect to . In general, rigid body translations [juj]n;[juj]t;[juj]m, rigid body rotations n;t;mand simple tensile-compressive strains in each direction n;t;mcan be all identi ed on a general linear de nition for each component of [ juj]. Figure 3 illustrates an example of this enriched set of fracture kinematic modes for a domain +on a tetrahedron. The centroid of dis taken as the origin for de ning all rigid body modes. Figure 3: Illustration of an enriched fracture kinematic mode set including rigid body translations and rotations for the + domain with respect to . The e ect of simple strains n;t;mis not included for the sake of clarity in the gure. Eq. 4 does not allow a coherent imposition of boundary conditions, since those that are prescribed on pertain the vector u, while the remaining that belong to +have their information contained in u+Hd[juj]. To overcome this, a typical technique in the framework [25] has been to introduce a new auxiliary eld 'function allowing for the following change of variables: ^u=u+'[juj]; (5) where ^uis a uni ed vector that has its node displacements completely consistent with uwhile retaining the strong discontinuity information continuously spread inside thanks to the 'function. The only mathematical constraint'is required to satisfy is the following: '=( 1xi2 + 0xi2 i= 1;2;3;4 (6) wherexidenotes a nodal position. Introducing this variable change along with the generalised de nition for fracture kinematics through a kinematic modes vector allows reaching an enriched de nition for the displacement eld accounting for a strong discontinuity: Strong discontinuity enhancement u=^u+ (H')J (7) where Jis an interpolation matrix based on local ;; coordinates operating on . The latter contains all 8 the unknown fracture kinematic modes of translation, rotation and simple strain: J=2 41 0 0 0   0 0 0 1 0000 0 0 10 0 0 3 5 (8a) T=[juj]n0[juj]t0[juj]m0ntmnnttmmT(8b) A robust de nition of 'is made to avoid variational and kinematic inconsistencies as much as possible ([23, 20]). The 'function will have a cubic, piece-wise ( +; domains), individual de nition per each local component ( n;t;m ) of [juj]. Using the local coordinates ( ;; ), a compound 'can be expressed using a complete cubic polynomial base P3: '=2 4' n 0 0 0' t 0 0 0' m3 5 ' j=PT 3  jj=n;t;m P3=1      222 2 22 22 2 333T = + n  n + t  t + m  mT(9) This extended de nition for 'allows the usage of  jcoecients as free parameters to aid the satisfaction of numerous consistency constraints that are helpful on building a solid strong discontinuity formulation structure [20]. 3.3. Variational integration of both discontinuity enhancements 3.3.1. Enhancement superposition To model both heterogeneities and local element fractures in the same element domain, the aforemen- tioned de nitions have to be integrated within a unique mathematical framework. The main assumption by doing this superposition is that both the material interface and the fracture interface coincide perfectly at the same surface during all the analysis. Based on the previous de nitions, the following can be reached by making a linear superposition: u=^u+eu+ (H')J (10) The general strain eld "associated to this displacement eld can be reached by applying a symmetric gradient operator rs() =1 2h r()T+r()i : "=rsu="+e"+^" (11a) "=rsu (11b) e"=rseu=e"= h ["]n(n ^n)s+ ["]m(n ^m)s+ ["]t n ^tsi (11c) ^"=rs[(H')J] = HrsJ'rsJrs'J |{z } bounded+[^n (J)]s |{z} unbounded=^"b+^"u (11d) Here, the resulting strain eld has been expressed as a sum of the di erent respective strain elds coming from both weak (e") and strong ( ^") discontinuity enhancements, as well as the standard strain eld ". The strain eld coming from the strong discontinuity enhancement is further classi ed as bounded ( ^"b) and unbounded (^"u), depending on the presence of the Dirac . The strain eld coming from the weak discontinuity eld is expressed using tensor products ( ) between the di erent unit vectors conforming the local frame. 9 3.3.2. Variational principle and its associated discretisation strategy The Hu-Washizu framework has been chosen as the variational principle in this work. It treats dis- placement, strain and stress elds ( real elds ) (u;";) as well as their respective variations ( virtual elds ) (u;";) as completely independent of each other. This allows a considerable exibility on the discretisa- tion strategy for each of the elds. Expressed already in a Voigt vector format, the equation system entailed by this variational framework is the following one: Z @utdVZ utfbdVZ @ uttdA= 0 (12a) Z et(@u")dV= 0 (12b) Z e"t(("))dV= 0 (12c) wheretandfbare the boundary traction and body force vectors. The @operator is the equivalent of the symmetric gradient operator rsin a Voigt format. An important distinction is made between the stress eld(") that is calculated from a constitutive law taking the real strain eld "and the real stress eld , which is independent. In the same way, the real strain eld "is not necessarily equal to the symmetric gradient of the displacement eld @u. The discretisation strategy for each of the real ( and virtual) elds must be able to capture the physics sought by the study allowing at the same time for mathematical simplicity and numerical eciency. This work will keep with the strategy already presented for both strong and weak discontinuity enhancements presented in [20, 21]. The displacement eld uand its variation uare discretized by only considering the standard base displacement eld ^uwithout the addition of any enhancements: u=^u=Nd (13a) u=Nd (13b) withNas a standard interpolation matrix and dthe standard nodal displacement vector. dis the cor- responding variation. This strategy means that only the eld ^uis used for describing node positions and imposing boundary conditions. The strain eld "and its variation "conserve all kinematics description terms as stated in Eqs. 11a-11d. Thus, it is this eld that has the main role of representing the physics of both heterogeneities and potential local fractures and to spread it within the variational framework. The discretisation of the weak discontinuity model considers the grouping of its internal strain jump variables ["]n, ["]t, ["]mwithin a weak discontinuity vector [j"j] =["]n["]t["]mTalong with its variation [j"j]. The realenhancement eld in the strong discontinuity model, as seen in Eq. 11d, has already a frac- ture kinematics modes vector de ned having all pertinent interpolation information. On the other hand, it has been preferred to discretize the virtual version of the strong discontinuity eld by only considering a three-mode constant virtual displacement jump eld [juj] =[juj]n[juj]t[juj]mT. The resulting discretisation can be stated as: "=Bd+G w[j"j] +G0 s (14a) "=Bd+G w[j"j] +G s[juj] (14b) 10 where weak enhancement operators G whave been introduced, as well their strong counterparts G0 s;G s. The apostrophe in these latter indicate an interpolation based on an enriched 9 fracture kinematic modes set to di erentiate with respect to a three-mode interpolation without apostrophe, like in the virtual operator G s. The operator for real and virtual weak discontinuity enhancements can be retrieved by putting Eq. 11c on a proper Voigt format, giving rise to an operator Hw: Gw= Hw;Hw=2 6666664n2 x nxmx nxtx n2 y nymy nyty n2 z nzmz nztz nxny+nynxnxmy+nymxnxty+nytx nzny+nynznzmy+nymznzty+nytz nznx+nxnznzmx+nxmznztx+nxtz3 7777775(15) where the components from each local frame vector ^n;^t;^mare explicitly used in Hw. For the strong discontinuity, Eq. 11d can used directly to obtain the real matrix operator: G0 s=H@J'@J@'J|{z} bounded+HsJ|{z} unbounded=G0 sb+G0 su (16) where a distinction has been made again between bounded and unbounded expressions. Hsis identi ed as a projection operator in Voigt format on the ^ndirection, with the following structure: Hs=2 6666664nx0 0 0ny0 0 0nz nynx0 0nzny nz0nx3 7777775(17) G sis actually de ned as the classical three-mode component matrix operator used in previous works using the EAS (enhanced assumed strain) approach [23, 12, 20]: G s=A VeHs+Hs=G sb+G su (18) The stress discretisation is approached by adopting an independent interpolation style for real and virtual elds (;): =Ss (19a) =Ss (19b) where S;Sare the respective interpolation matrices and s;sare the real and virtual stress vectors. The interpolation matrix Sremains a free choice, which determines the possible shape (and order) of the de nitive stress eld on the element. The matrix Sremains rather constrained by the structure of both weak and strong discontinuity enhancements. In this work, a uniform behaviour for the stress eld will be chosen, having Sas the identity matrix and sas a six component vector. The stress eld (") calculated from a constitutive law will be de ned by following the Discrete Strong Discontinuity Approach (DSDA), where only the bounded part of the real strain eld (Eq. 14a) is contem- plated. From Eqs. 11a-11d, it is seen that the only singular term on the strain de nition is identi ed as ^"u. A linear elastic constitutive law having a second order tensor Cwill be considered in this study, whose linear parameters depend on the material phase present on each domain +; . This is an exigence coming from the weak discontinuity model, which respects each material phase with its own material law. Having all this in mind, (") can be set as: (") =C"b=C Bd+G w[j"j] +G0 sb (20) 11 where: C=2 6666664c 1c 2c 20 0 0 c 2c 1c 20 0 0 c 2c 2c 10 0 0 0 0 0 c s0 0 0 0 0 0 c s0 0 0 0 0 0 c s3 7777775(21a) c 1=E(1) (1 +) (12); c 2=E (1 +) (12); c s=E 2 (1 +)(21b) The unbounded expression, which in uences the interface directly, is compensated by a traction-separation law, discussed later in this section. 3.3.3. Summary of variational analysis Substitution of all discretised elds into Eqs. 12a-12c and their subsequent treatment yield the main equations governing the behaviour of all discontinuity enhancement internal variables. Working the second equation from the variational system (Eq. 12b) and keeping the independence between the internal variables [ j"j] and, two di erent orthogonality constraints are found for the virtual stress interpolation matrix S: Z STG wdV=0 (22a) Z STG0 sdV=0 (22b) These relations only demand to prove that, for a chosen de nition of the real eld enhancements, a matrix Scapable of being simultaneously orthogonal to both enhancement operators exists. It can be shown that a matrix Scan be always found having polynomial components with high enough order to yield the necessary free parameters to satisfy both Eqs. 22a, 22b [20]. Other than this veri cation, Eqs. 22a, 22b represent no impact to any structure within this E-FEM approach. It is the last variational equation (Eq. 12c) that grants the most relevant foundations for the model structure. Taking the arbitrary nature of all independent vectors ( d;[j"j];[juj]) involved in ", three equations are derived: Z BTdV=Z BT(")dV (23a) Z GT wdV=Z GT w(")dV (23b) Z GT sdV=Z GT s(")dV (23c) The rst one (23a), as it is, serves to simplify the internal-external force balance. It is also the only direct relation that helps to calculate the real stress eld as a function of the constitutive stress eld (") since Bis xed for a given nite element. If the stress elds had the same polynomial order, a strong equality could established. It is clearly not the case for this study. However, a volume-weighted average can be devised: =1 VeZ (")dV (24) Eq. 23b is used to obtain the governing expression for the weak discontinuity law. It is a custom to assume that both sides of Eq. 23b are equal to zero. Having zero in the left side will allow the imposition 12 of the patch test constraint that nishes the de nition of G wthrough  (which has been provided to the reader beforehand for the sake of intelligibility). The patch test demands the expression to be satis ed at least for the case of a constant real stress eld. Therefore, choosing already constant in the rst place simpli es the process. The right hand in Eq. 23b, having particularised G w, yields the governing relation for the weak discontinuity model: Z GT wC Bd+G w[j"j] +G0 sb dV= 0 (25a) !Kwbd+Kww[j"j] +Kws= 0 (25b) Kwb=V+V VeHT w C+C B (25c) Kww=V+V V2eHT w VC++V+C Hw (25d) Kws=V+V VeHT w C+G0+ sbCG0 sb (25e) where subdomain integrations have been done through +; , introducing the de nition of speci c sti ness matrices Kwb,KwwandKws. For the case of strong discontinuity operators, averaged versions G0+ sb,G0 sb were de ned per subdomain as: G0+ sb=1 V+Z +G0 sbdV; G0 sb=1 VZ G0 sbdV (26) The governing relation coming from the strong discontinuity model provides with the means to calculate the average traction vector Ton the crack surface  dwhich, in this case, having chosen a constant real stress eld, is equal to the actual traction vector Titself. Again, left and right hands of Eq. 23c are simultaneously set to zero as a deliberate assumption in the framework. First, in an analogous process to the weak discontinuity model, the zero at the left side of Eq. 23c provides information to set a de nition for GT s, which has been already shown to the reader (Eq. 18). GT sis also assumed as a constant matrix (the simplest form) having the same bounded-unbounded constitution as its real counterpart. Its unbounded term is de ned exactly as with the real enhancement, only without the matrix J(three modes instead of nine). All integrations related to it can be done using basic properties of the Dirac function . The bounded term is then determined from this left-hand relation, again, considering the satisfaction of a patch test constraint for a constant  eld. The equality between both left and right hands in Eq. 23c allows nding an expression for the disconti- nuity traction vector Tas a function of the traction vector Tebased on the constitutive stress ("). Both traction vectors arise when working the unbounded terms, which are multiplied by the operator Hswhich is nothing but a projection operator onto the discontinuity surface  d. Using then the equivalence between stress elds in Eq. 24, one can reach the convenient conclusion that Te=T=T. Ultimately, the zero in the right hand of Eq. 23c establishes the governing relation for the traction through  das a function of Te=T: T=Z GT sC Bd+G w[j"j] +G0 sb dV (27a) !T=Ksbd+Ksw[j"j] +Kss (27b) Ksb=1 VeHT s V+C+VC B (27c) Ksw=V+V V2eHT s C+C Hw (27d) Kss=1 VeHT s V+C+G0+ sbVCG0 sb (27e) 13 where additional sti ness matrices Ksb,Ksw,Ksshave been de ned. From both governing equations 25b, 27b, the internal state variable [ j"j] can be eliminated to express the traction vector Tas a function of standard nodal displacements dand the fracture kinematic modes vector : [j"j] =K1 ww(Kws+Kwbd) (28a) !T=Te+M (28b) Te=Ked= KsbKswK1 wwKwb d (28c) M=KssKswK1 wwKws (28d) Eq. 28b is the standard form for presenting the traction vector associated with the local crack on the current E-FEM framework. Eq. 28b determines the state of crack traction depending on loads and the kinematic state of the crack. On one hand, the Tevector represents the traction imposed by the current load in the element, depending on the standard displacement vector d. It is a variable mechanical demand that will guide the evolution of the fracture process within the element. The Mmatrix, referred to in this work as the fracture sti ness , accounts for the mechanical impact of having a kinematic state of fracture. It does not depend on load but only on a mixture between element geometry characteristics and basic properties of the surface  d. This composition of load e ects and the crack state yield a current state for T. In turn,Tis to be controlled (damaged) by a traction separation law. From now on, the description of all remaining mathematical structures in the model will be done on the local frame ^n;^t;^m(with the respective ;; coordinates). For instance, the traction vector Tis now expressed asT=TnTtTmT. The same will be done with other vectors and matrices in later sections. At last, the internal-external elemental force balance can be established from the rst relation of the Hu-Washizu system (Eq. 12a), and using Eq. 23a as it is to express the real stress eld as a function of the constitutive stress eld (") directly: fe int=Z Bt(")dV=Z NtfbdV+Z NttdA=fe ext (29) where the classical distinction of internal and external elemental forces fe int,fe exthas been brought upon. Special interest is placed on the internal forces fe int, where Eq. 20 can be used to further develop: fe int=Kbbd+Kbw[j"j] +Kbs (30a) Kbb=Bt V+C++VC B (30b) Kbw=V+V VeBt C+C Hw (30c) Kbs=Bt V+C+G0+ sb+VCG0 sb (30d) where, again, intermediary sti ness matrices Kbb;Kbw;Kbshave been de ned. The use of intermediary Kmatrices such as in these latest de nitions allows a compact presentation of the main relations governing the behaviour of internal variables [ j"j];and nodal displacements d(Eq. 25b, 27b and 30a), and will help to present the global solution process in Section 3.7. 3.4. Localisation criteria A localisation criterion is required to designate a local failure within an element and to start introducing the impact of fracture mechanics through the expressions developed earlier plus a set of traction separation laws (Section 3.5). It should be extensive in the sense that it should cover all possible stress limit states within a given element. 14 Localisation also has the secondary role in the case of homogeneous elements of determining the ori- entation and location of the fracture surface  d. The base element explored in this work remains a linear tetrahedron whose real stress eld has been set as uniform. Thus, a localisation criterion can objectively determine the orientation but not the exact location of the local fracture. On this matter, the authors of this study have decided to take the approach followed in [12, 41] and make  dto pass through the centroid of the element reaching localisation. For all heterogeneous elements, it has been already stated that the crack plane will be made perfectly coincident with the material interface plane. This work has taken the traction vector components Tn;Tt;Tmto articulate the localisation criterion on the normal-shear stress space ( n; ) on d. A completely closed criterion has been built by composing a piece-wise function with three di erent sub-criteria: •ARankine criterion meant to set the limit for predominantly tensile state of stresses, represented by a straight vertical line on the ( n; ) plane with a positive abscissayR. •Acompression limit criterion intended to limit the amount of compression undertaken by an element. This is also represented by a straight vertical line on the ( n; ) plane with a negative abscissayC. •AMohr-Coulomb criterion for covering all remaining possible limit combinations of shear-tensile and shear-compressive state of stresses on the same plane. This is represented by a couple of symmetric lines (with respect to the naxis) described with the equation =C(tan )n This piece-wise curve remains very simple (4 material parameters) and ensures that all possible stress failure states are covered. Figure 4 illustrates this compound criterion. Figure 4: The envelope localisation criterion on the ( n; ) plane used in this work showing notable intersecting points. A proportional load is shown beginning on a zero state up to reaching an intersection with the curve at ( y; y) . To detect if a given load has breached this localisation criterion, a di erent assessment is done whether a fully homogeneous or heterogeneous element is examined. Figure 5 shows the approach for each case. A heterogeneous element already has a potential  dorientation de ned, so the only task is to take the e ective traction on this plane, given directly by the components of Te(Eq. 28c) to calculate the current location (ni; i) on the (n; ) plane: ni=Tn; i=q T2 t+T2m (31) Note that the failure is in a certain way forced to emerge on the plain given by the material phase interface, regardless if the criterion has already been violated by considering another di erent orientation for the traction plane. Having a current traction location point (A), a line can be constructed beginning from the origin assuming a proportional load behaviour, and an intersection (B) can be found between this line and a localisation criterion. The actual distance ABfrom the envelope curve tothe current load point through 15 this line (outwards is positive) de nes the localisation criterion  wfor this case, depending on the type of intersection: w=8 >>>< >>>:yR ni1 jTj; Intersection with Rankine lineyC ni1 jTj; Intersection with compression line C i+nitan 1 jTj;Intersection with Mohr-Couloumb line(32) where it has been recognized that jTj=q 2 i+2ni. The element is considered in a localised state if  w is equal or larger than zero, establishing the intersection point ( y; y) as the initial yield parameters for operating the respective traction-separation laws (Section 3.5): (y; y) =8 >>>< >>>: yR;yR ni i Rankine Intersection yC;yC ni i Compression Intersection C i+nitan ni;C i+nitan i Mohr Coulomb Intersection(33) Note that while the normal direction ^nof the plane is fully determined, the parallel direction remains arbitrary. In this work it is assumed that the direction of the resultant shear traction component iis assigned with the ^tunit vector exclusively. Thus, at localisation, Tmis set to zero. Figure 5: Illustration for the assessment of the di erent localisation criteria whether if treating a heterogeneous element (left) or a homogeneous element (right), which uses the classical 3D Mohr circle theory. The compression line and any related calculations with it are not shown for the sake of clarity in this illustration. For a homogeneous element, the state of stresses on the ( n; ) plane can be portrayed by drawing the biggest of its three Mohr stress circles with the help of the principal stresses 1;2;3. The distances from the circle to the nearest points on each sub-criterion line are monitored for detecting localisation, producing three criteria  R;MC;Cfor Rankine, Mohr-Coulomb and compression lines, respectively: R=1YR C=YC3 MC=r(Ctan ) cos (34) Sign conventions are the same as with heterogeneous calculations. Figure 5 does not show the distance Cfor the sake of clarity. When the Mohr circle grows suciently to touch one of the lines, the predicted intersection values ( y; y) are set as initial yield parameters: (y; y) =8 >< >:(yR;0) Rankine Intersection (yC;0) Compression Intersection cos2 (+Ctan );cos2 (Ctan ) Mohr Coulomb Intersection(35) 16 The direction ^nof the normal to the plane corresponding to the maximized criterion can be calculated by considering the angle (Figure 5) between the radial line containing the intersection and the horizontal axis. From classical 3D Mohr circle theory, this represents the double of the rotation angle that the actual stress plane has, expressed in principal stress coordinates, around the 2axis on the principal stress reference frame. For the case of localisations occurring with Rankine or compression criteria, it is clear that ^nwill be whether the direction associated to 1or3, respectively ( = 0;). For the case of a collision with the Mohr-Coulomb line, can be expressed as a function of using basic geometry. The cosine directors (v1;v2;v3) for ^nin the principal stress frame of reference can then be calculated as: v1= cos 2 v2= 0 v3= sin 2(36) The directors v1;v2;v3should then be rotated into the current frame of reference to continue the analysis. Again, as with the heterogeneous element case, parallel directions cannot be fully speci ed. For the case of Rankine or compression localisations, ^tis just set as the direction of the following the next eigenvector from the one selected for ^n. For the Mohr-Coulomb criterion, the direction of the resultant shear stress (knowing ^n) on dis set as ^t. In all cases, ^mis calculated by a product vector between ^nand^t. The collection of traction-separation laws explained later in this work also involve the crack surface stress components not explicitly included in the traction vector T, as it is just a projection of an entire stress tensor (a vector if presented in a Voigt format) on the ^ndirection. In the local frame, three components of the stress tensor jdassociated to the crack plane  dcan be already identi ed as Tn;Tt;Tmdirectly: jd=2 4nn=Tnnt=Ttnm=Tm tn=Tttttm mn=Tmmtmm3 5 (37) Thus, it is pertinent to establish a yield value for ytt;ytm;ymmfor further damage treatment. An easy option would be to take these stress components from the current stress state whose traction vector just breached the localisation criteria. But, this traction point could not be necessarily placed close enough to the criteria lines, so that the ytt;ytm;ymmvalues would not be correct. This is especially true when managing constant load steps that could breach a localisation criterion leaving big gaps in between. As the aforementioned calculations only consider the ( n; ) plane, there is no way to assess these stress components accurately. As estimation of their state at the calculated yield point ( y; y) must be made. Whether in the case of a homogeneous or a heterogeneous element, a unique estimation approach is proposed: a linear scale factor is calculated based on the current load point ( ni; i) and the established yield intersection (y; y). Then, a complete tensor jdis scaled using this factor, for then retrieving a projected value for ytt;ytm;ymmassociated to the state ( y; y). The goal is to be as most consistent as possible with the state of stresses at localisation, as it will determine the initial behaviour for all fracture kinematic modes () in the traction-separation laws described in the next section. 3.5. Traction-separation laws and fracture mechanics A traction separation law has the unique role of describing the behaviour of fracture mechanics directly once the element has reached a localisation state. It talks entirely from the perspective of the fracture, without any regard to native element characteristics. Its de nition starts by setting an expression for relating an equivalent stress eqassociated to the discontinuity surface  dand a post-localisation behaviour functionq, whereqshall involve fracture state variables ( ) along with intrinsic fracture physical properties: m=eqq=0 (38) where a mcriterion is always assumed as zero for de ning physically admissible states for eq. Note that Eq. 38 has been cast in a vector mode since multiple traction-separation relations can be de ned. 17 3.5.1. Main traction-separation relations As already mentioned before, the strong discontinuity part of the framework is already expressed as a function of the traction vector T. Therefore, it is intuitive to build traction separation laws around it, and seteqsimply asTderiving one traction-separation equation for each traction component Tn;Tt;Tm: eqq=Tq=Te+Mq=0 (39a) Tn=Ten+9X kMnkk=qn (39b) Tt=Tet+9X kMtkk=qt (39c) Tm=Tem+9X kMmkk=qm (39d) where the vector Thas been expanded on its load-driven (elastic) vector Teand fracture mode-driven component Mfeaturing the 39 fracture sti ness matrix. Consequently, three di erent qjfunctions will control the evolution for each of the traction components. While the qjfunctions can have a versatile de nition to represent a variety of hardening phenomena, in this work these functions will exclusively have the role of damage drivers, decreasing the traction values associated to the crack as it grows through the evolution of its kinematic modes ( ),e.g., as the crack separates ([ juj]n0) or slides more ([ juj]t0, [juj]m0). For simplifying the approach, it has been decided not to involve all the nine possible kinematic modes in theqj functions. Emphasis will be placed on the rigid body kinematic modes [ juj]n0;[juj]t0;[juj]m0to control the evolution of the damaging process for all traction vector components. qis then de ned as: q=2 4qn qt qm3 5=2 4yn yt ym3 5eyn GfI[juj]n0+yt GfII[juj]t0+ym GfII[juj]m0 (40) whereyjare the initial yield values set from the localisation calculations (Section 3.4). The GfIand GfIIconstants are de ned as internal fracture physical parameters that can be perceived as fracture surface energies for mode I and mode II fracture types coming from the classical theory of fracture mechanics, respectively. The choice of an exponential damage law allows to consider the compound contribution of each fracture rigid body mode to develop the crack completely until reaching a point in which the local crack no longer o ers any resistance on any direction in a smooth fashion. Such state is referred to as a terminal separation condition . As theoretically these conditions will strictly hold at in nite separation, the authors has decided that a crack progression returning a damaging factor of at least e4(0.0183) is enough for all following discussions. 3.5.2. Alternate traction-separation relations To complete the damaging process on the stress components ytt;ytm;ymmnot explicitly present in T (refer to the discussion in the 3.4), three additional relations of the same style are added to the system. These relations arise from the projections of on the other remaining directions ^tand ^m, completely unrelated to that happening on the ^ndirection. Indeed, it should be noted that satisfaction of the Eqs. 39b-39d during all the crack development process does not guarantee a decrease of ytt;ytm;ymm. If no conditions are imposed on these components, the model will inevitably present stress locking, which has been already shown as problematic [23, 20]. This is actually the main motivation to enrich fracture kinematics beyond rigid body translations: the liberalisation of additional kinematic modes is needed in order have enough parameters to establish consistent internal equilibrium within an element. To build relations analogous to Eqs. 39b-39d for the ^tand ^mdirections, the corresponding projections are extracted by proposing alternate 18 Hst,Hsmprojection operators (in the local frame) instead of Hsin Eq. 27a-27e: Hst=2 66666640 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 03 7777775Hsm=2 66666640 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 03 7777775(41) Following the structure of Eq. 27a, this would produce two di erent traction vectors T0andT00corresponding the^tand ^mdirections, respectively: T0=1 VeZ HT stC Bd+G w[j"j] +G0 sb =T0 e+M0 (42a) T00=1 VeZ HT smC Bd+G w[j"j] +G0 sb =T00 e+M00 (42b) where corresponding load driven vectors T0 e,T00 eand sti ness matrices M0,M00come along. For the ^t projection, only two traction components T0 t;T0 mare needed to capture the stress components tt;tm, while for the ^mprojection, only one is ultimately required ( T00 m) for addressing mm. The following system of three alternate traction-separation relations can then be reached: T0 m=T0 et+9X kM0 tkk=qtt (43a) T0 t=T0 em+9X kM0 mkk=qtm (43b) T00 m=T00 em+9X kM00 mkk=qmm (43c) 2 4qtt qtm qmm3 5=2 4ytt ytm ymm3 5eyn GfI[juj]n0+yt GfII[juj]t0+ym GfII[juj]m0 (43d) where new damaging functions qtt;qtm;qmmhave been de ned, related to yield parameters ytt;ytm;ymm coming from the localisation calculations (estimated at the end of Section 3.4). Note that the exponential decay function in Eq. 43d remains the same as with the main traction-separation equation system (Eq. 40). Satisfaction of Eqs. 39b-39d, 43a-43c will guarantee a proper and complete damage process for any element regardless of its geometric properties, fracture surface characteristics and load evolution patterns. 3.5.3. Crack reversibility and a phenomenological fracture mechanics approach As this framework has liberalised multiple degrees of freedom for local fractures, a variety of physical situations can virtually happen to a given element having a certain load, geometry and crack con gurations. One could think about a crack relaxation and its potential closure, or a change in sliding direction, among others. If these local phenomena are not correctly assessed from the modelling standpoint, they will induce non-physical and instability-prone behaviours in a number of elements. For instance, allowing negative values for any crack rigid body modes [ juj]j0on any of the exponential expressions de ning the qjfunctions will return damage factors over 1.0. If it is chosen in this case to just freeze such modes for a predicted reversible behaviour, Eqs. 39b-39d, 43a-43c will respond to this arbitrary overconstraint by returning physically incoherent values for other kinematic modes and/or the spurious growth of traction vector components. All 19 these issues promote the emergence of numerical disturbances in the model that will prevent the successful solution of a particular problem. To prevent this, such o ending elements could be killed during numerical simulation runtime by enforcing an arbitrary damage mechanism on them, thus disabling the propagation of these abnormal e ects. This, however, would mean that the model would deliberately neglect local physical phenomena in favour of numerical hygiene and, depending on the element count, it could go up to the point of in uencing the response of a given system considerably, preventing a meaningful study. For this reason, the authors of this work have made their best to incorporate the modelling of certain local crack physics in a mostly basic but ecient style. The reader might nd these additions to lie on the edge of overburdening the E-FEM framework and might probably nd a hard time justifying their existence or even their accurate representation. In any way, a reasonable compromise has been found between the introduction of new physics along their parameters, the associated computational costs and the quality of modelling results. Going forward on this matter, one of the most practical features of the model is that, by the way the main traction-separation laws (Eqs. 39b-39d) were de ned, they can be used to model direct fracture mechanics with a variety of e ects. Indeed, the system conformed by Eqs. 39b-39d can be perceived as a cohesion force balance on each of the ^n;^t;^mdirections normalized by the crack surface area . As such, they represent the direct mechanics of +as a physical body with respect to : X iTni=qn (44a) X iTmi=qt (44b) X iTti=qm (44c) where the damaging functions qn;qt;qmrepresent a cohesion force (per surface area) between the bodies. Atterminal separation conditions , these cohesion forces are driven to zero, which completely renders the +body independent from . However, Eqs. 44a-44c can be taken as a base for representing further e ects adding the presence of additional forces in the corresponding directions and adapting the expressions as necessary, keeping physical coherence in the model. In this sense, the authors of this work have added the possibility to model crack reclosure, crack compression and crack sliding with frictional resistance. Crack reclosure on this kind of E-FEM framework has been addressed in [46], and the ideas managed in this speci c research are taken by the authors of this work to incorporate a basic closing mechanism. This modelling approach implies the de nition of irreversible crack separation/sliding states that determine the current state of cohesion damage (qn;qt;qm), while having an actual crack separation/sliding states that describe the current physical state of the crack (kinematic state of +). In the works done in [38, 12, 41], allirreversible crack kinematic modes are made permanently equal to actual crack kinematic modes, thus reclosure modelling is not possible in such models. Again, it is only in [46] that this distinction has been made, albeit only for a single kinematic mode formulation. In this work, as only crack rigid body modes [juj]n0;[juj]t0;[juj]m0are involved in the qn;qt;qmdamaging functions, only these can serve as a base for irreversible crack mode de nitions. These are named as jujn,jujtandjujm, respectively. They shall only have positive change rates: _jujj0 (45) Whenever the state of load imposes growing conditions for all crack rigid body modes further from their accumulated/irreversible values jujn,jujtandjujm, Eqs. 44a-44c are solved exactly as shown considering nonlinear exponential decay functions qn;qt;qm. If the load state is such that it promotes a crack rigid body 20 response less than the one already registered by jujn,jujt,jujmin any direction, then the current cohesion forceqjbecomes a simple linear spring model, allowing the decreasing kinematic mode [ juj]j0to return through a straight line down to a certain value, possibly reaching closure ([ juj]j0= 0). Figure 6 shows this basic reversibility behaviour for both normal separation and sliding cases from points 1 to 3. Figure 6: Illustration for reversible and irreversible behaviours for crack normal separation (left) and sliding kinematic modes (right). The path from 1 to 2 is irreversible, while 2-3 is depicted as a reversible path forced by an insucient load that fails to continue crack evolution to point 2a. If the load is reversed enough, the crack can will be eventually closed (3). Further load reversibility induces compression for the left plot and sliding on the opposite side of the crack for the right plot (3-4). Eventually, one can even reach irreversible behaviour zones once again (4-5). There is an inherent asymmetry between the normal and parallel models since the normal model is the only one capable of exhibiting compression states. The sti ness of the linear spring driving this reversible behaviour varies depending on the current values ofjujn,jujtandjujm. For a given mode that is currently on a decreasing state, its balance equation changes to the following form: X iTji=" yj jujjeyn GfIjujn+yt GfIIjujt+ym GfIIjujm# [juj]j0(46) where the right hand is the equation of a straight line with its slope value enclosed in brackets, and having an intersection with the origin [0 ;0]. Note that the slope depends on irreversible values jujn,jujtandjujm, which are assumed to be completely determined from a previous load step, so that Eq. 46 continues to be solved for an actual kinematic state [juj]j0, resulting in a value less than the currently stored jujj. For the speci c case of having [ juj]n0reversing enough to reach a value of zero, the crack is considered to be in a state of compression , and normal contact force between the +; will start to develop ( i.e., a negative value of Tn). The value of this normal force will depend on the amount of compression (negative [juj]n0) and a proposed contact "sti ness". This is actually what in contact analysis literature would be called an interface sti ness , in this case regularised by a characteristic length due to the nature of the normalized force balance. In this work, it has been estimated using the young modulus Ecof the least sti of the materials found in +; divided by the characteristic length lcof the whole element in question. This brings the following behaviour for the ^ndirection under compression: Tn=Ec lc [juj]n0(47) The summation symbol has been removed as from now on no further forces/e ects will be accounted for the^ndirection. Note also that the slope of the straight line continues to be strictly positive. Indeed, for negative values of Tn, Eq. 47 will return a negative value for [ juj]n0. The sti ness of the straight line 21 in Eq. 47 is naturally much larger than that used for crack closing in Eq. 46. Again, Figure 6 (left) shows this behaviour for Tnup to point 4. This contact interface cannot reach in nite compression values, as a compaction limit has been set. This value remains precisely the limit stress under compression yC already de ned for localisation. Reaching this compression limit will reach a state of material compaction, temperingTnto a nal value TCin a permanent fashion. The authors acknowledges that more physically consistent compaction models that consider microporosity properties exist in current literature. The present simple approach was just deemed sucient to complete the physical behaviour for the model at this stage of development. This is nally depicted by point 5 in Figure 6 for Tn. For the crack sliding directions ^j,^m, there is no special consideration for any sign change, and the reversibility relation (Eq. 46) can be still used to predict the actual values for [ juj]t0and [juj]m0. However, their behaviour will fall back to irreversible decay if their respective absolute values go beyond jujtand jujm,i.e., if the current sliding values go negative enough to attain the current accumulated sliding on the opposite sense. This means that sliding damage is considered as symmetric, and it can be accumulated whether triggering from the positive side of the ^j,^mdirections or their negative sides, as well. This is portrayed also on the right side of Figure 6. The last important e ect to account for on the ^j,^mdirections is the friction induced by normal compres- sion. Whenever a negative value of Tnis calculated, an associated resisting traction Tt,Tmwill emerge, and will totally prevent or moderately diminish further sliding motion between +; . This dissipation is introduced as a simple Coulomb friction model as follows: X iTti=TtTt=Tet+9X kMtkkTt=qt (48a) X iTmi=TmTm=Tem+9X kMmkkTm=qm (48b) Ttjmax=tjTnj; Tmjmax=mjTnj (48c) Ttqt; Tmqm (48d) wheretcharacterises the friction on the sliding direction ^t(recalling that it was chosen during localisation analysis as the direction of the shear traction vector component resultant, Section 3.4) and where m corresponds to the friction model on the perpendicular sliding direction ^m. Also, ahas been introduced to select the proper sign corresponding to oppose the driving motion by TtorTm. Independent Coulomb friction coecients t;mhave been introduced. Note that, while Eqs. 48a, 48b have been written for the case of nonlinear decay behaviour, these also apply for crack closure states . The conditions in Eq. 48d imply that the friction model is to be applied once the respective damaged cohesion normalized forces qt;qmare weakened enough to fall below the generated friction. This avoids having a duplicated resistance to sliding motion, as the cohesion forces associated to the original traction separation law already imply some physics due to surface roughness while developing the crack. The application of friction in Eqs. 48a, 48b does not exactly work as an addition of a residual constant to the entire traction curve, but rather as a resistance threshold preventing a current equilibrium position to move further whether in the main crack direction or backwards. Figure 7 has been conceived to better illustrate the intended use of Coulomb friction. Departing from a current equilibrium point at [ juj]t0i, the current load step will attempt to move from the current state, whether further through the nonlinear curve by means of a positive traction load change  T+ ton the crack sti ness line with xed slope Mttor to retract on the opposite direction by means of a negative traction load change  T t. If friction is present, this available traction load has to overcome this amount of friction, otherwise the crack will remain on its current equilibrium state. This threshold given by Ttis the same for both sliding directions. If external load e ects manage to overcome friction, the e ective sliding will be provided by a diminished traction T tTt, to end whether in a state [ juj]+ t0or [juj] t0depending on load direction. 22 Figure 7: Illustration of the application of the friction model to resist further sliding on the ^tdirection. Even after all measures taken concerning the incorporation of these local mechanics e ects, there are still some rare cases of extremely unexpected element behaviour. At this point, where the o ending element count represented a substantially tiny proportion of the active global model, the authors of this work nally rendered the act of deliberately killing these elements as totally justi ed. 3.5.4. Closing relations, a complete view of the system At this point, even if enough relations have been de ned to determine (and properly damage) the full state of the traction separation vector Tand all remaining stress components ytt;ytm;ymmassociated to the crack surface, the complete algebraic system described so far remains open. For the authors of this work, the existence of three extra fracture kinematic modes is regarded as an opportunity: the framework allows to de ne three additional auxiliary relations for any kind of mathematical or phenomenological representation improvements. Three nal closing relations are de ned to fully uncouple rigid body translations by driving the compound e ect of all other (non-translational) kinematic modes to zero on each of the main traction balance equations. This way, the three main traction equations can be written having their linear left hands based on a single variable. This way, three nal relations are added to the system, closing it completely. A summary can be described as follows: Main traction-separation laws (Tn;Tt;Tm) - Ten+Mnn[juj]n0=8 >>>>< >>>>:qn [juj]n0>jujn (qn)i1 jujn[juj]n00[juj]n0jujn Ec lc[juj]n0yCTn0 qC TnyC(49a) Tet+Mtt[juj]t0Tt=( q t [juj]t0 >jujt qtj jujt[juj]t0 [juj]t0 <jujt(49b) Tem+Mmm[juj]m0Tm=( q m [juj]m0 >jujm qmj jujm[juj]m0 [juj]m0 <jujm(49c) 23 Alternate traction-separation laws (ytt;ytm;ymm) - T0 et+9X kM0 tkk=qtt (50a) T0 em+9X kM0 mkk=qtm (50b) T00 em+9X kM00 mkk=qmm (50c) Closing relations - 9X k=4Mnkk= 0 (51a) 9X k=4Mtkk= 0 (51b) 9X k=4Mmkk= 0 (51c) In this summary, the fact that all the non-diagonal terms from the rst three columns of the M matrix have been nulli ed by making the correct choice for the '-related functions [20], leaves only the Mnn;Mtt;Mmmterms (M11;M22;M33) as slope coecients for the crack rigid body translations [ juj]n0;[juj]t0and [juj]m0, respectively. Note that the closing relations (Eqs. 51a-51c) have been also implicitly considered in the rst three traction-separation equations (Eqs. 49a-49c) to isolate rigid body translations completely. So, basically, these three rst equations remain completely independent, for then solving the behaviour of the remaining rotation and simple axial strain modes in the rest of the system. Once again, the choices taken in this work have led the crack rigid body translation modes to get the highest modelling priority. 3.6. Solution of the nonlinear traction-separation law equation system The last discussion in this section will address the solution of the traction-separation equation system depicted by Eqs. 49a-51c. As a rst step it should be recognised that the rst three equations 49a-49c conform a system completely independent by itself. This means that a solution for [ juj]n0;[juj]t0;[juj]m0can be obtained from this system, which can be later used to solve the remaining six equations. There is, however, a nonlinear multivariate exponential linking the three equations at once. It should be remembered that there could be the case in which reversibility is detected on one or more directions. Should that be the case, the corresponding rigid body translation [ juj]j0can be solved and used immediately on the other equations. In the worst case scenario of having three nonlinear relations, an e ective approach is to condense to a single univariate equation. This can be done by nding linear relations between each pair of variables by suppressing the common exponential factor. An example, a reduction to [ juj]n0can be done as follows: [juj]t0=yt ynMtt Ten+Mnn[juj]n0 Tet Mtt(52a) [juj]m0=ym ynMmm Ten+Mnn[juj]n0 Tem Mmm(52b) Once with that, the equation treating the ^ndirection is taken to the general form a+b[juj]n0=ce(d[juj]n0+f)(53) 24 which, following the typical approach in [12], has a closed solution expressed with the rst branch of the LambertW0function: [juj]n0=1 AW0(C)B A(54a) A=d B =d ba C =d bce(f+d ba)(54b) Whether for nonlinear or reversible behaviour cases, the solution of the system is viewed as the intersection of each of the lines represented by the left hand of each of Eqs. 49a-49c with the behaviour law functions (right hand of the equations). The ordinate of the intersection represents the current value for the entire traction vector Tjfor each corresponding direction. The straight lines have xed slopes coming from their Mnn;Mtt;Mmmcrack sti ness coecients. Because of the nature of the strong discontinuity models and the way the'was designed, these crack sti ness values are strictly negative . The lines can only vary their position depending on the load demand implied by the components of the vector Te. By the way the traction-separation phenomenological laws were conceived, there will always exist a solution for the subsystem on the [ juj]n0;[juj]t0;[juj]m0crack rigid body translations. Figure 8 illustrates the intersection for multiple positions of the representative straight lines just discussed. Having [ juj]n0;[juj]t0;[juj]m0, the remaining equations of the system are entirely linear, and thus the solution for the remaining kinematic modesn;t;m;n;t;mcan be readily obtained by any conventional linear system solution method. Figure 8: Solution of individual equations for crack rigid body translations as intersections between a straight line and the behaviour laws proposed in Section 3.5.3. Three exemplary lines are introduced for analysis in the normal ^nand parallel ^t directions. The position of all these lines is regulated by the load driven traction vector Te, which determines the intersection with the vertical axis. The slope of the lines never change, and it is always negative. 3.7. Linearisation of the model and global solution The linearisation process for the application of nonlinear global iterative solutions can be accomplished by departing from the main relations governing each of the internal discontinuity variables [ j"j],and the nodal displacements d(Eqs. 25b, 30a and the system 49a-51c). The process starts by stating expressions for linear increments of key dependent variables in these governing equations as a function of linear increments 25 of a set of independent variables in the following fashion: Y=nvX i=1@Y @XiXi (55) where general vectors X,Yare designated as independent and dependent variables, respectively. In this work, di erent dependent variable increments  Yare de ned to aid the calculation of the complete state of the global mechanical system through Eqs. 25b, 30a and 49a-51c, while the increments of the main inde- pendent variables  Xin the framework remain always as  d,  [j"j] and following all the developments so far. The rst and most standard case is the global internal-external force balance, in which the increment Yis taken as the global residual of internal and external nodal forces from successive iterations kandk+1 for a set of assembled element force balance equations. For this, the linearisation of Eq. 30a is assembled repeatedly ( Aoperator) for all neelements currently iterating during the application of a load step p+ 1 coming from an already solved step p. This yields the following: nel A( Kbbd (k+1) p+1+Kbw [j"j] (k+1) p+1+Kbs (k+1) p+1) =nelA( fe int (k+1) p+1fe ext p) (56) In the E-FEM framework, which is based on the principle of incompatible modes [13], no global continuity is guaranteed for any of the internal variables characterising both weak and strong discontinuities. Hence, none of the linearisations associated to Eqs. 25b, 49a-51c go through any conventional assembly process. Indeed, there will be one relation to be satis ed for each element in an independent fashion. For instance, the weak discontinuity governing relation (Eq. 25b) can also be linearised by de ning its balance through a dependent variable wat element level, which is always sought as zero during the nonlinear solution process. Taking increments in this expression leads to the following: Kwbd (k+1) p+1+Kww [j"j] (k+1) p+1+Kws (k+1) p+1= w p+1=0w k p+1(57) where Eq. 57 is never assembled in a global format. Linearisation of the strong discontinuity model governing relations is not evident since di erent sets of traction separation laws were incorporated with di erent considerations (Eqs. 49a- 49c, 50a-50c) and there are also the closing relations coming from purely algebraic considerations (Eqs. 51a-51a). In [12], a cascade linearisation approach was followed, in which a criterion (as the one de ned in Eq. 38) is established as the base dependent variable. Then, a rst linearisation step is done with respect to the increment of an intermediate variable  Talong with  . Finally, the increment  Tis expressed as increments of all main independent variables, merging any common expressions at the end. In this work, such approach will be taken to address Eqs. 49a- 49c, 50a-50c. The rst three main traction separation equations can be grouped having mas their main dependent variable, considering exactly the vectorTas an intermediate independent variable to express the linearisation as: m=@m @TT+@m @ (58a) m=Ksbd+Ksw[j"j] + Kss@q @  (58b) The next three traction separation relations (Eqs. 50a-50c) are treated the same way by setting a zero balance aas their dependent variable. A three component traction vector Ta= T0 tT0 mT00 mTis 26 designated this time as the intermediate independent variable and the cohesion force expressions are grouped asqa=qttqtmqmmT, yielding: a=Ksbd+Ksw [j"j] + Kss@qa @  (59a) where double star ( ) notations have been introduced to de ne sti ness matrices Ksb;Ksw;Kssfor each of these alternate traction components coming from other projections di erent than ^n. Finally, linearisation of the closing Eqs. 51a-51c is done by de ning a dependent variable cand directly taking an increment with respect to : c=@ @2 43X j=19X k=4Mjkk3 5 (60) Once the nine relations have been linearised, they can be assembled to build a single matrix block-based equation to pack down the nine relations into a single linearised expression de ning  0as a collective dependent increment: Ksbd (k+1) p+1+Ksw [j"j] (k+1) p+1+Kss (k+1) p+1= 0=00 k p+1(61) where the matrices Ksb,Ksw,Ksscontain the aforementioned small-scale assembly of linearisations of all traction separation laws and closing equations previously developed. Having all these linearisations, the same global solution approach in [12] is followed. A complete linearised system is rst built in a block-matrix format as: 2 4KbbKbw Kbs KwbKww Kws KsbKswKss3 5(k) p+12 4d  [j"j] 3 5(k) p+1=2 64nelAffe intfe extg w 03 75(k) p+1(62) Considering, once again, that the weak and strong discontinuity models only have immediate in uence at a local level, the typical global strategy solution for a nonlinear iteration k+ 1 is to take the solution of the standard displacement vector dsolved on the previous iteration kto internally update the actual state of the weak and strong discontinuity variables [ j"j];using the nonlinear equation system (Eqs. 49a-51c) and Eq. 25b. With this, the linear system 62 is condensed in the  [ j"j];increments to reach a nal expression involving exclusively the increment of the standard nodal vector  dand a condensed elemental sti ness matrix Ksc, which is to one to be assembled and eventually used to execute a numerical method of choice by a standard nite element global solution engine. The nal condensed global FEM equilibrium equation system has the following form: nelA( Ksc (k+1) p+1d (k+1) p+1) =nelAffe intfe extg (k) p+1(63a) Ksc (k+1) p+1=KbbKbwKbsKww Kws KswKss1Kwb Ksb (k) p+1(63b) As a nal step, for the calculation of the residuals at the right hand of Eq. 63a, element internal forces fe intcan be easily retrieved each iteration kby using whether Eq. 29 or 30a depending if real stresses  have been already calculated during program runtime: fe int (k) p+1=VeBt (k) p+1=Kbbd+Kbw[j"j] +Kbs (k) p+1(64a) 27 4. Numerical Simulations In this section, an implementation done within the FEAP program will be tested under a variety of conditions to study the response of the generalised E-FEM formulation and its capability for predicting di erent kinds of fracture processes. This implementation was done on an instance of the FEAP program [47] adapted as a component for multi-scale and coupled multi-physics simulations (coFEAP) through the software component framework CTL [48]. Within the FEAP program, the implementation only consisted in introducing an element code routine and auxiliary local routines associated to the present model, without altering absolutely anything of the global solution process routines, whether linear or iterative non-linear. Using this nalised program, two types of simulations will be presented to the reader: •Homogeneous sample simulations , which are primarily a test for the strong discontinuity model consistency. These setups will feature completely homogeneous elements without the presence of any weak discontinuities. Torsional load simulations will take place to discuss the emergence of controlled tridimensional fracture processes. By controlled , the authors of this work imply that a speci c pattern for a tridimensional fracture process is already expected for a given simulation setup knowing its geometry and load conditions. In this sense, a pure torsion case should qualitatively return a well- known tridimensional fracture process. •Basic heterogeneous sample simulations , in which the development of fracture processes around a de nite region having material interfaces will be studied. A cubical homogeneous matrix material domain having a spherical inclusion in the center is prescribed with di erent load con gurations. This allows to test the weak discontinuity and its in uence on the initiation and growth of fractures in this model. Tension and compression load cases are presented to the reader, tracking the development of the local crack networks at various points of the evolution of fracture processes. All these simulations will be done side-by-side with a numerical model having the single mode I E-FEM framework as proposed in [12, 49]. The comparison will help to study the added value implied by the introduction of the fracture mode generalisation and also to better determine the reliability limits for the single mode I framework. It is important to remark that is not meant to execute a formal and extensive validation procedure for the generalisation built in this E-FEM framework. The authors of this work would rather like to present this section to the reader as an illustrative showcase of the new model capabilities, triggering ideas that will make a solid case for future works. 4.1. Homogeneous simulations The following homogeneous material torsion simulations are meant to introduce the reader to the process of tracking the evolution of a tridimensional fracture process. The proposed numerical model is a cubic domain having a torsion load on its upper (+ z) face. The basic description of the model is shown in Figure 9. The length Lcof the cube is 10 mm, and it features a non-structured tetrahedral mesh having a characteristic length of 0.1 mm. To imitate pure torsional load as much as possible while simultaneously maintaining a displacement-controlled scheme, a quick dedicated sub-routine has been implemented. An angular torsion value is set as the main load reference, for then calculating a set of polar coordinates ri;ifor each node iat the upper face of the cube taking the center of the cube as the origin. Based on each node's radial and angular position, a magnitude and direction of displacement on the x;ydirection is assigned (uxi;uyi), which is consistent with the reference angular torsion . The lower face of the cube is kept completely xed. For this simulation, a nal torsional angle of 0.03was prescribed (equivalent to a sample shear strain of 3.7 104mm/mm at the outermost bers). A special feature of this setup is that it is has a weakened material band with a height  hat the center of the model. Here, the rst occurrence of crack initiation is promoted by letting the remaining material bulk (above and below the band) to have slightly increased strength values ( yR,C). This bulk has exactly 28 Figure 9: Model descriptions for homogeneous torsion simulations having a 80% weakened material band and a torsion angle of 0.03. The weak material band in the model is depicted in red. the same linear material parameters. This is done to avoid any spurious concentrations near the vicinity of the regions having boundary conditions to avoid a direct bias on fracture initiation at the lower and upper faces. For these simulations, the middle band has been set to  h= 0:8Lc. Material properties for both formulations can be found in Tables 1 and 2. For the remaining bulk material, the only di erence is the increased yR;Cparameters (by 10%) in both formulation cases. Generalised modes formulation input parameters: weakened band E+= 14000 MPa += 0:2E= 14000 MPa = 0:2 yR= 9 MPa C= 12 MPa tan = 0:6yC=500 MPa GfI= 0:0005 MPammGfII= 0:001 MPamm t= 0:6 m= 0:6 Table 1: Input material parameters used for the generalised modes formulation for the weakened material band in the homo- geneous torsion simulations. Single mode formulation input parameters: weakened band E+= 14000 MPa += 0:2E= 14000 MPa += 0:2 yR= 9 MPa GfI= 0:00001 MPamm Table 2: Input material parameters used for the single modes formulation for the weakened material band in the homogeneous torsion simulations. These simulations aim to track the formation and evolution of the crack networks for the two formulations being compared. Speci cally, it is of interest to study how the fracture process switches its nature during the global post-localisation stage. Di erences are expected within the single mode and generalised mode formulations when shear and fracture sliding demands become more signi cant during the load application. To make a clear point on this, Figures 10, 11 and 12 are proposed to illustrate the initiation and consolidation of the local crack networks for both models. Figure 10 focuses on the initiation of the fracture process, xed at the rst load step where the global behaviour of the model becomes nonlinear. Here, it can be observed how both models start to promote the emergence of classical 45fracture processes for a brittle or quasi-brittle material under torsion, forming an angular periodic pattern of local crack networks. Each of the cube model faces presents this initiation as these are the outermost surface layers of the body under torsion, which are well-known to be the rst regions to reach high levels of stress following the classical theory of mechanics of materials. The entirety of these local crack networks have been localised under the mode I criterion, i.e., all in tension failure (shown 29 Figure 10: Elements conforming the local crack networks within the cube domain for a torsion load just right after beginning global nonlinear behaviour. All blue elements correspond to a tensile (mode I) failure model. Each network develops a tilted classical brittle fracture process with no di erence between single and generalised models. The cube shows also a translucent contour of resultant nodal displacements having red as maximums and red as minimums. Figure 11: Elements conforming the local crack networks within the cube domain for a torsion load almost halfway through load application. All blue elements correspond to a tensile (mode I) failure model, and red elements are associated to a mode II failure model. The networks have evolved to converge to the center of the cube. The generalised model starts to exhibit the emergence of a shift in the fracture process towards a mode II behaviour. as blue elements within the transparent cube). Almost no di erence can be perceived between the networks predicted by the single mode and generalised modes formulations. Figure 11 portrays a state where the fracture process related to mode I failure and kinematics has completely developed in both models. The four local crack networks have converged in a slightly helical 30 Figure 12: Elements conforming the local crack networks within the cube domain for a torsion load at the end of load application. All blue elements correspond to a tensile (mode I) failure model, and red elements are associated to a mode II failure model. The fracture process has not signi cantly evolved for the single mode formulation, whereas the generalised model predicts a mode II behaviour for the conclusion of the fracture process. fashion, all reaching the core of the cubical sample. Again, the single mode and the generalised modes formulations make the same overall prediction for the shape of this consolidation of local crack networks. The generalised model, however, starts exhibiting small regions containing mode II related failures and kinematics (shown as red elements). These are non-consolidated networks and they still do not have any in uence in global mechanical behaviour. Finally, Figure 12 captures the nal state of both numerical simulations at the end of torsion application. This is where a critical divergence in behaviour can be acknowledged between both formulations. The single mode formulation is not able to predict any signi cant evolution of the fracture process but only a slight expansion of the already predicted shape. On the other hand, the generalised formulation allows a complete change of nature for the fracture process: it starts to focus exclusively on mode II failure and kinematics, expanding through the entirety of upper and lower borders of the weakened material band. This allows to nalise the entire fracture process for the sample and to release all the remaining stored energy therein. This simple setup helps to illustrate the reader that a single mode E-FEM formulation will be properly capable of portraying a full 3D fracture process while the general state of stresses promotes the emergence of the corresponding single mode within the domain of the simulation model. Once the current state of the structure along with a given load demands the fracture process to change its nature due to its nonlinear behaviour, a single mode formulation will fail to capture subsequent fracture evolution. 4.2. Basic Heterogeneous Simulations The next set of simulations aims to study the e ect of weak discontinuities in the formation and propa- gation of a tridimensional fracture process considering a simple heterogeneity distribution and having simple load cases. The main simulation setup will focus on a homogeneous matrix domain having the presence of a spherical inclusion made of a stronger material. Fracture processes are expected to start at the material interfaces and to evolve in particular ways depending on each load case. This gives hint on basic failure mechanisms happening at the mesoscale for many quasi-brittle materials. 31 The geometry of the heterogeneous sample is a cube with a side dimension of 10 mm having a perfectly spherical inclusion of 6 mm of diameter. This material phase distribution is projected onto a mesh conformed, once again, by regular (good aspect ratios) but unstructured tetrahedral elements. The average characteristic length of the elements is 0.15 mm. As a consequence, a number of these elements have been cutby the inclusion boundaries, generating a number of material interfaces modelled through the weak discontinuity approach. Figure 13 shows the general mesh built and already projected with the material phases for all these simulations, as well as the generated weak discontinuity interfaces for the sphere boundary. Figure 13: General model mesh description for the basic heterogeneous setups featuring a single spherical inclusion. On the left, the mesh of all elements touched by the inclusion sphere within a translucent mesh of the remaining bulk of the material matrix domain. In the middle, a mid-cut section of the entire mesh, distinguishing exclusive mortar material elements (in blue), exclusive aggregate material elements (in gray) and bi-phase elements (in red) generally found at the sphere border. On the right, a sphere reconstructed by taking the resulting planar interface surfaces  dcoming from each bi-phase element, along with their respective normal vector. The latter usually stands as a cross check for correct inclusion boundary construction. The parameters given for both single and generalised modes formulations can be found in Tables 3 and 4. Note that the interface failure parameters take the resistance limits coming from the weakest material phase (mortar). This is a general assumption done in this work. Pure tensile and pure compressive load cases Single Mode Formulation Input Parameters Mortar-only elements E+= 14000 MPa += 0:2E= 14000 MPa = 0:2 yR= 9 MPa GfI= 0:00001 MPamm Aggregate-only elements E+= 70000 MPa += 0:2E= 70000 MPa = 0:2 yR= 45 MPa GfI= 0:00005 MPamm Interface elements (- for mortar, + for aggregate) E+= 70000 MPa += 0:2E= 14000 MPa = 0:2 yR= 9 MPa GfI= 0:00001 MPamm Table 3: Input material parameters used for the single mode formulation in the simple heterogeneous simulations in this work. have been contemplated for the following simulations. As in previous sections, each load is exerted until a full tridimensional fracture process is developed in the entire domain, reaching a terminal point in which the cohesion between the mortar and the inclusion can be considered as completely broken , rendering the model no longer capable of o ering any kind of mechanical resistance for further progress in displacement application. 32 Generalised Modes Formulation Input Parameters Mortar-only elements E+= 14000 MPa += 0:2E= 14000 MPa = 0:2 yR= 9 MPa C= 36 MPa tan = 0:6yC=500 MPa GfI= 0:0001 MPammGfII= 0:6 MPamm t= 0:6 m= 0:6 Aggregate-only elements E+= 70000 MPa += 0:2E= 70000 MPa = 0:2 yR= 45 MPa C= 180 MPa tan = 0:6yC=500 MPa GfI= 0:0001 MPammGfII= 0:6 MPamm t= 0:6 m= 0:6 Interface elements (- for mortar, + for aggregate) E+= 70000 MPa += 0:2E= 14000 MPa = 0:2 yR= 9 MPa C= 36 MPa tan = 0:6yC=500 MPa GfI= 0:0001 MPammGfII= 0:6 MPamm t= 0:6 m= 0:6 Table 4: Input material parameters used for the generalised modes formulation in the simple heterogeneous simulations in this work. 4.2.1. Tensile simulations In a rst model, a pure tensile load will be imposed exactly as in the previous homogeneous material simulations, xing the face corresponding to the lower zvertical position and having a displacement on the opposite face on the positive zdirection. Again, two adjacent lateral faces are again restricted in the x andydirections to avoid rigid body modes. As a reference, the maximum vertical displacement imposed in both setups is of 0.006 mm. Based on previous discussions, no signi cant di erences are expected between both E-FEM formulations, as the single-mode version is known to perform in a very similar manner than the generalised model given that the global load is purely tensile. Figure 14 shows the global response for each model concerning total vertical reactions. Indeed, both models attain approximately the same global maximum resistance, having some slight di erences in post-localisation behaviour. Figure 14: Global stress-strain curves under a tensile load for the single spherical inclusion setups for both single mode and generalised modes formulations. The sample stress and sample strain are calculated exactly the same way as in previous sections. 33 The direct development of the tridimensional fracture process is described through a series of pictures depicting the micro-crack networks at multiple steps of the load application process, done exactly as in Section 4.1. The state of local fracture networks is rst captured at the moment of maximum structural resistance for both formulations 15, and then a state portraying completely post-localisation behaviour is retrieved 16. All these images are taken by making a plane cut at the middle of the cube model, revealing what happens around and within the spherical inclusion with detail. Figure 15: Elements conforming the local crack networks for a global strain corresponding to the maximum sample resistance. All blue elements correspond to a tensile (mode I) failure model. The fracture process initiates at the top and bottom of the inclusion for both formulations. For the single-mode formulation, there is only one possible progression state for a local crack, which is opening. Thus, the coloured contours for this case only show a single blue colour, indicating that the elements concerned did break on a Rankine criterion, exhibiting a certain amount of normal separation exclusively, without any possibility for crack reversal, local compression or sliding. For the generalised model, di erent colours are assigned to elements that broke under a determined criterion (a blue palette for Rankine, a red palette for Mohr-Coulomb or compaction limit), and also by their current local crack state: opening, relaxing or closed, the last necessarily implying a certain amount of compression and a certain frictional resistance to slide. As anticipated, the local crack network initiation and eventual coalescence for a fully developed fracture process was practically the same for both E-FEM formulations. The only perceptible di erence is that the generalised model is able to provide slightly more insight about elements exhibiting relaxation (lighter blue coloured elements). Absolutely no element was observed having any kind of local mode II failure. The interface between mortar and spherical inclusion is clearly segregated through pure tensile planes. Again, it is striking to see this uniform fault plane predictions having no global crack tracking algorithms. However, the lack of priority to develop a single fracture process is also a trait that drives the presence of two simultaneous failure bands, which is not intuitive from a physically practical standpoint. It should be recalled that this simulation setup remains an unusual representation featuring a perfectly centered sphere within an equally perfect cubic homogeneous matrix domain. 34 Figure 16: Elements conforming the local crack networks for a global strain corresponding to a completely developed fracture process. All elements are still related to mode I local failures, reporting some crack relaxations (lighter blue colours) for the generalised formulation. The fracture process has a very de nite double band shape for both formulations, observing no signi cant di erences. 4.2.2. Compression simulations In this case, a pure compression load is imposed by just reverting the displacement direction in the setups presented in the previous section, having exactly the same surrounding boundary conditions. There is a signi cant di erence in both global resistance and ductility (de ned hereby as the sample strain correspond- ing to maximum resistance "max) between single mode and generalised modes formulations. Thus, di erent maximum compression distances have to be used to fully developed their respective fracture processes: -0.05 mm and -0.1 mm, respectively. The outcome of these simulations is presented in the same style as in the tensile load case. Figure 17 shows the behaviour of globally de ned stress-strain curves for both models. This time, the evolution of the fracture process will be presented in a separate fashion for the single and generalised formulations, as they have signi cant di erences in multiple stages through the development of local crack networks. Each step as identi ed with a strain value relative to the sample ductility ( "max), covering states before and during the maximum resistance load steps as well as complete global post- localisation behaviours. Figure 18 shows a four step description (( "1"4)) for the case of the single mode formulation, containing naturally only blue element networks corresponding to mode I local failures. The single mode formulation presents a double band initiation at the lower and upper hemispheres of the inclusion ( "1), the fracture process then propagates completely around it ( "2). At maximum resistance, the fracture process begins to extend out of the sphere vicinity in a rather uncontrolled and irregular fashion ("3). At post-localisation, the single mode formulation exhibits a sudden fall in resistance, and the fracture process seems to evolve in an arbitrary manner, extending to most of the mortar domain ( "4). Cohesion between the inclusion and matrix material is completely broken, nishing ending the evolution of the fracture process. Figures 19 and 20 describe six evolution steps ( "1"6) for the case of the generalised model. A two band network initiation and an eventual takeover of the inclusion surface is also observed for this formulation. The elements in this network are predominantly associated to a mode I failure ( "1), but it can be observed a small 35 Figure 17: Global stress-strain curves under a compressive load for the single spherical inclusion setups for both single mode and generalised modes formulations. amount of regularly di used mode II elements in the upper and lower regions ( "2). The behaviour starts to diverge from the single mode formulation from the third step ( "3), where shear-compression zones initiate at the top and bottom faces of the mortar domain. These regions contain exclusively elements associated to a mode II local failure, and also in a crack closure state, presenting frictional behaviour (the reddest colour in the palette denotes this state). These zones extend vertically until fusing with the inclusion surface fracture process ("4). At maximum resistance, a middle shear band emerges and surrounds the inclusion surface completely ( "5), and nally the fracture extends from this band to all the mortar domain in a vertical fashion, nishing the tridimensional fracture process. ( "6) Note that most of the mode II-related elements within these local networks remain in a crack closure/compression state. These results suggest that these predominant state serves to provide more ductility and resistance to the sample, and the authors continues to state that the modelling of this local fracture phenomenon is of utter importance for the prediction of compressive failure at larger scales in quasi-brittle materials. Overall, it is observed once again that both formulations share certain stages of the behaviour of the tridimensional fracture process, but the generalised model is able to grant a larger insight of further steps where mode II local failures and compression having frictional sliding become signi cantly important. Again, this is re ected on the global stress-strain response. Indeed, it appears to be evident that the modelling of mode II-related phenomena cannot be neglected under certain load cases. 5. Conclusions This work has introduced the reader to an integration of a generalised strong discontinuity model with a weak discontinuity model to conform a numerical analysis approach capable of delivering meaningful representations of tridimensional fracture processes for quasi-brittle materials. A theoretical background is provided, describing the fundamental mathematical aspects of the generalised fracture mode kinematics for the strong discontinuity and the weak discontinuity. A clear description of localisation criteria has been provided, along with the traction separation equation system that bene ts from the generalised structure of the fracture kinematic modes to incorporate more robust local crack physics. It is worth remembering that the authors of this work gave priority to the in uence of crack rigid body displacements modes (separation and sliding) for the damage-regulating de nitions. 36 Figure 18: Fracture process evolution for the compressive load case considering the single mode formulation. Four di erent stages are reported at di erent strains relative to the maximum resistance strain "max A number of numerical simulations have been done to attest to general capabilities of the generalised models for di erent situations. The reader should always keep in mind that this work was never intended to be a formal validation and application of the model to the context of a speci c material study. Despite this fact, the authors recognises that many of the results found through Section 4 are quite revealing on the behaviour of tridimensional fracture processes. In particular, there is the participation of mode II mechanics on the emergence of di erent global and local crack network behaviours. Results suggest that a model capable of considering local sliding, compression and frictional mechanics will be able to gain more insight at the fracture processes under highly compressive or shear demands at larger scales, speci cally those featuring larger ductility and resistance through more energy dissipation mechanisms. This could serve as an evidence for stating that mode I local failures are not suciently enough to tell the completely story of tridimensional fracture processes at small scales for quasi-brittle materials. 37 Figure 19: Fracture process evolution for the compressive load case considering the generalised modes formulation (four stages). Colours having a blue palette denote those related to mode I local failures, while the red palette corresponds to those related to mode II. The reddest colour speci cally indicates elements under crack closure compression. The present work has shown that it is still possible to remain with an E-FEM framework having a strictly local mathematical structure, bene ting of enough simplicity without having to make a signi cant incursion on a global FE numerical solution platform. It should be reminded that the essence of this line of research was to achieve a suciently predictive advanced FE numerical model keeping its computational burdens and implementation complexities to a minimum. In this sense, the authors of this work believe that the E-FEM framework is still an attractive option worth exploring, o ering unique advantages over other advanced FE techniques.