This paper discusses the consistency of the theoretical basis behind the kinematic models of strong disconti- nuity methods for local fracture simulations using the Embedded Finite Element Method (E-FEM). A brief review is made on the elemental enhancement functions from the current E-FEM literature and how previous works managed to model mode I (normal) and mode II (parallel) fracture kinematics in multiple dimensions. Further analysis is made on how these approaches also introduce unintended mesh dependencies and basic kinematic inconsistencies of the fracture model with respect to its hosting element. Notable work from a few authors discussing some of these issues and their contributions to resolve them is reviewed as well. Based on this analysis, a new proposal of strong discontinuity enhancement functions is introduced to ensure a broader kinematic coherence within the element and to avoid the observed theoretical faults. This is done by making a more extensive use of the exibility granted by the Hu-Washizu variational principle and by in- troducing new algebraic constraints that will ensure more correct fracture kinematics without compromising the acknowledged simplicity of the whole E-FEM framework. Element-level simulations are done to compare the outputs within a group of selected formulation approaches, including the novel proposal. Simulations show that the new element formulation grants a wider level of basic kinematic coherence between the local fracture outputs and element kinematics themselves, demonstrating an increase in robustness that might drive the usefulness and competitiveness of E-FEM techniques for fracture simulations to a higher level. Keywords: E-FEM modelling, embedded strong discontinuity, enhanced nite elements, incompatible modes, local fracture modelling, nite element enhancement functions 1. Introduction One of the core features of the E-FEM modelling approach is the ability to simulate local material fractures by the introduction of strong discontinuity enhancement functions to elemental displacement elds. These functions are driven by internal variables representing the fracture kinematics in the normal and parallel directions to a given fracture plane. This kind of internal variable fracture representation is based on the theory of incompatible modes [1], which allows contiguous elements in the same mesh to have internal kinematic states that do not respect rigorous continuity between them. This gives rise to ecient and transparent solution processes where element fracture kinematic variables can be calculated in a completely internal solution process. The update of all elemental contributions for a nonlinear global solution step can then be done using a traditional FEM solver framework, requiring no modi cations to it unless the framework is deliberately enriched with a non-local analysis features such as, for example, the use of a crack path eld method to insure continuity between elements [2]. Email address: alejandro.ortegalaborin@3sr-grenoble.fr (A. Ortega) This is a preprint August 24, 2021 The E-FEM approach has been praised for its simple, yet powerful capacity for the accurate modelling of local and global fracture processes of quasi-brittle materials over other methods such as X-FEM or the Partition of Unity [3]. These latter have a more robust and deeper de nition of their kinematics attacking directly to the de nition of nodal interpolation functions and their support. This allows a natural continuity across elements and a more enriched description of a fracture in general [4]. However, this additional complexity works as a double-edged sword, bringing more challenges to their numerical implementation. Despite this fact, the E-FEM approach has not become as popular as other alternative FEM options as X-FEM to gain much ground into commercial FEM packages such as ABAQUS or LS-DYNA [5, 6, 7]. While the simplicity of the E-FEM approach has been acknowledged many times, there are some numerical issues that prevent this formulation to gain a consensual solid stance as a reliable and accurate structural numerical simulation approach for material fractures. These have been recognized by many authors working on the eld: mesh dependencies, severe stress locking, unintended coupling between fracture displacement modes, among others [8, 9, 10]. The intent of this work is to make an analysis on how the basic de nitions of strong discontinuity kinematics have been introduced so far on the E-FEM framework (Section 2 and 3) to identify the roots of its potential theoretical faults and related numerical issues (Section 4). Some authors gave their valuable insights on these challenges and propose theoretical enhancements as a workaround [11, 12, 13, 8]. This work attempts to consolidate these e orts to establish a basic understanding that allows for a new set of wiser E-FEM proposals in a general 3D format, avoiding current issues and bringing more robustness and modelling capabilities for the framework. Following the previous considerations, this paper proposes a new speci c E-FEM strong discontinuity enhancement scheme (Section 5). Simulations at element level are done to compare the results to older enhancement proposals (Section 6), with encouraging conclusions (Section 7). It is important to remark that the present work will focus on the branch of the E-FEM framework related to SKON (Statically and Kinematically Optimal Non-symmetric) formulations, which is the one that produces the most physically consistent models, both kinematically and statically [14, 15, 16, 17]. The authors believe it is worth continuing the theoretical works in this line since its physical correctness sets a solid foundation for modelling further complex crack phenomena like reclosure, local compression and explicit local friction, as well as the possibility to consider other kind of discontinuities in a simultaneous fashion. 2. Analysis of basic internal strong discontinuity kinematics The starting point of every fracture representation done within the E-FEM framework is to enrich the basic displacement eld of a continuous body with the direct addition of a mathematical discontinuity function representing the displacement jump from one split region of the body to the other: u=u+Hd[juj]; (1) where the theoretical displacement of the body uis set as the composition of a regular displacement uand a jump normally introduced by means of a Heaviside function having the location of the fault surface  d as its main trigger. An internal fracture kinematics vector [ juj] is used to represent the general motion of one of the fractured regions with respect to the other. Figure 1 and Figure 2 are very common illustrations used by authors in the eld to visualize the discontinuous enrichment for the fracture of a body split in two regions +and . The origin, shape and orientation of the surface  drepresenting the fracture remains completely inde- pendent from this analysis, and will depend on the speci c strain localisation criterion chosen for studying a given material. As the general philosophy of the E-FEM framework tends to aim for formulation simplicity, 2 Figure 1: Basic schematic of an embedded strong discontinuity in 2D the most common fracture surface of choice is typically a plane or a line. This work assumes this simple geometry as the fracture surface from now on. We will also assume that the analysis is done on a general 3D domain unless speci c remarks for 2D or 1D cases are made. While Eq. 1 already allows for the calculation of a valid strain eld directly through a symmetric gradient operatorrs() =1 2[r() +rt()], there is a couple of important basic remarks worth making at this point. 2.1. Kinematic consistency of boundary condition imposition A discretized element whose kinematics is governed by Eq. 1, does notallow a consistent imposition of boundary conditions on any traditional global FEM solution engine. This is because the displacements on the borders of the domain of the element ( i.e., the nodes on it) are directly reachable through the vector u, while the displacements on the borders of +are rather a composition of uand [juj]. We can use the vectorualone to impose boundary conditions on , but not on +since the actual displacement in this region remains a composition of two di erent vectors u+ [juj]. Figure 2: Example of the behavior of the eld ufor a 1D element (conformed by nodes N1andN2) with a constant displacement jump [ juj], as well as a corresponding ^ u This issue has been well identi ed even from the rst evolving steps of the framework and also treated by Oliver [10] through the introduction of the so-called auxiliary function '. The function 'is de ned as to allow for the existence of a new displacement vector ^u=u+'[juj]; (2) by letting'to have a single property: it will have the value of 1 on any node lying on the border of + and zero in those on the remaining border: 3 Basic'function requirements '=( 1x2 + 0x2 (3) This is the only formal mathematical requirement as stated in [10] apart from the basic fact that 'should be continuous. The most common form of 'taken from the literature is yet the simplest one to satisfy Eq. 3 for a given element geometry. E.g., for the linear tetrahedron, the simplest 'remains a composition of standard linear interpolation functions. In this case, 'also happens to be unique. Analogous constructions for'can be easily made for a 2D triangle and for 1D cases (Fig. 3). Figure 3: Example of the 'function for a constant stress triangle de ned by the nodes N1;N2;N3, with node 1 on the + domain and the others on . Note that for this case there is only one possible linear 'de nition If'complies with the previous de nition, then the new vector ^ucan be seen as a continuous fusion of the local displacement eld and the jump vector associated to the fracture, having border displacements that coincide with those of u. Thus, for the practical purposes of nodal displacements in a discretized element, ^ucontains all the correct information of border nodes and can be used for imposing boundary conditions in a direct and unequivocal fashion. Note that, on the other hand, the displacement eld within the element does not coincide with that of u.^ucan be seen as a special version of uwhere the discontinuity jump is continuously spread by'along the element domain. 'could be thus interpreted as a sort of regularisation function (see Fig. 2). With this de nition in mind, Eq. 1 is rewritten as u=^u+ (H') [juj]; (4) where the previous vector uhas been replaced by the newly de ned vector ^u, and the discontinuity base function now gets composed of both a Heaviside and '. From this point, authors generally stop referring to uand simply take ^uas the de facto standard solution vector for the reasons mentioned above, and Eq. 4 is taken as the reference for developing all derived mechanical elds. While the introduction of 'solves a fundamental consistency problem on boundary conditions imposition, it also responsible for other undesired e ects on general element kinematics and fracture mode coupling. A deeper discussion is presented on Section 4. 2.2. Kinematics at terminal separation conditions and meaning of [juj] Post-localisation behavior within an element is generally modeled by making use of one of two popular options: a continuous strong discontinuity approach (CSDA) or a discrete strong discontinuity approach (DSDA) [18, 19]. Regardless of the approach taken, usually there is an initial mechanical connection between the newly split bodies, which will harden and/or decay as the crack separation evolves depending on the governing post-localisation law. Eventually, there should be no in uencing forces from +to when the split bodies have suciently separated. This state of failure will be referred to as a terminal separation 4 condition . In such state, the traction happening at the fracture surface should be zero, and the internal fracture kinematics vector [ juj] will practically represent the rigid body motion of the +domain with respect to , depending on the de nition chosen for the vector [ juj]. [juj] in most of the current E-FEM literature has been simply de ned as a one to three components vector of constant values in space at the element scale describing a rigid body translation of +with respect to . However, as [juj] stands as an internal element de nition with no strong global restrictions, it could be arranged to describe more complex kinematic modes such as rigid body rotations, having non-constant component behavior through the failure surface  d. [juj] thus remains a free modeling choice in this sense. But, whatever this choice may be, its behavior should be consistent with the standard element kinematics (whether uor^u). An obvious example is that the magnitude of a constant [ juj] in the normal direction ^nto the discontinuity plane (crack normal separation), should be no greater than the magnitude of the relative nodal displacement in the element giving rise to this separation. This is shown in Fig. 4 for a 2D constant stress triangle (CST). A relative displacement dnis imposed on this element and, assuming post-localisation, a corresponding crack separation [ juj]ndevelops. By common sense, it should be expected that [ juj]ndn, where at terminal separation conditions we could probably have [juj]n'dndepending on the speci c model. Figure 4: Example of kinematic consistency that should hold between the fracture displacement vector [ juj] and the nodal displacement vector dfor a constant stress triangle. Considering this, it is important to remark the framework does notexplicitly guarantee any kinematic consistency between the internal variable [ juj] and a standard global variable ^uor any derived quantities such asdn. This is due to the same fact that [ juj] is an internal element de nition. This is specially true when reaching a terminal separation condition, where there's no further in uence of any post-localisation law derived through CSDA or DSDA that would imply any additional coercion to [ juj] in any way. In a terminal separation condition, there's no other governing law than that of the purely natural kinematics arising from the modeling choices made in the basic formulation of the strong discontinuity, such as '. In the end, it should be clear that the physical meaning of [ juj] should be carefully assessed with respect to the formulation choices made within the E-FEM framework. This will be further discussed in Section 4 and 5. Having all this in mind, the strain eld is extracted by taking the symmetrical gradient rs() in Eq. 4. This study will assume a constant vector [ juj] through  d, which implies that this strong discontinuity study will be limited to a rigid body translation of +with respect to . The strain eld can be thus 5 expressed as: "=rsu=rs^u+rs[(H') [juj]] )"=rs^u(r' [juj])s+ (rH [juj])s )"=rs^u(r' [juj])s |{z} ^"b+(^n [juj])s |{z} ^"u; (5) where the rst tensor product is identi ed as the bounded part of the strong discontinuity strain ^"band the second tensor product is the unbounded strain ^"u.is the Dirac delta, with a trigger position corresponding to the failure surface  d, just as the Heaviside described previously. 3. In-depth analysis of variational foundations The use of customized eld enhancement functions on the E-FEM framework requires a variational principle that grants sucient exibility to e ectively manage multiple independent mechanical elds when making element discretisation. After all, the E-FEM remains a mixed FE formulation. This is the main reason why the three- eld Hu-Washizu principle has been the preferred approach for almost all authors working on the framework. There are, however, many details usually omitted on the description of its appli- cation that merit formal mathematical justi cation. This will allow the reader to know why some numerical choices prevail in the E-FEM literature and also to better understand the measures taken for improving the formulation explained in later sections. We take the Hu-Washizu functional for three independent elds of displacement, strain and stress ( u;";) as the point of departure: IHW(u;";) =1 2Z ":C:"dVZ :"dV +Z :rudVZ fbudVZ @ tudA;(6) where the double bar notation speci es a mechanical eld on its 2nd order tensor format (the omission of it would mean that the eld is on the Voigt vector notation), Cis a fourth-order linear elasticity tensor, fbis a vector of body forces and tis a vector of traction forces as commonly found in the conventional FEM. To reach equilibrium, a stationary condition for this functional is required, and this implies taking the variation of it along with the variations associated to each of the displacement, strain and stress elds (IHW(u;";) = 0) , leading to three independent equations, already written using a Voigt vector notation: Hu-Washizu variational system for the E-FEM framework Z @utdVZ utfbdVZ @ uttdA= 0 (7a) Z et(@u")dV= 0 (7b) Z e"t(("))dV= 0 (7c) It's worth noticing the di erence made in Eq. 7c between ("), which is the stress eld calculated from the constitutive relations (in this case through an elastic relation (") =C"), and, which is the 6 independent stress eld to be found in the analysis. These elds are not necessarily equal, and the variational analysis will provide with a speci c relation between them. The Eq. Sys. 7 is the most commonly used form of the Hu-Washizu principle as seen on the E-FEM literature. The solution for the system, in this continuous format, inevitably returns the equilibrium of the body (already expressed on Eq. 7a), the strong kinematic compatibility and constitutive relations: ru"= 0 (8a) (")= 0 (8b) From this point, any further step requires stating a eld discretisation strategy, and exibility will be drawn from Eqs. 8a and 8b as required, having cases where these are not necessarily satis ed in a strong manner. 3.1. A word on the discretisation strategy When using the Hu-Washizu functional, it is important to remember that the three elds of displacement, strain and stress ( u;";) along with the other three elds of displacement variation, strain variation and stress variation ( u;";) are all completely independent. Therefore, each of the six elds will require a unique discretisation choice and its entirely up to each study to determine the most convenient eld discretisation setup, keeping in mind that the main goal of the framework is to represent the e ects of the strong discontinuity in the most physically consistent, yet numerically ecient way. Apart from this, the only formal restriction is to satisfy the Eq. Sys. 7 at all costs. It should be remembered that the discretisation strategy follows the outline of the SKON family of E-FEM approaches: the structure of the Hu-Washizu principle is exploited to render the formulation as physically consistent as possible, including both correct kinematics and statics simultaneously, resulting in a non-symmetric FE formulation. 3.1.1. Displacement eld discretisation For discretizing the displacement eld u, authors naturally take advantage of the facts discussed in Section 2.1, and take only ^uinto account: u=^u=Nd; (9) withNas a standard interpolation matrix and dthe standard nodal displacement vector. The variation of the displacement eld uis also discretized in the same way, taking the displacement vector variation d: u=Nd (10) Eq. 9 implies that the displacement eld to be found through this strategy does not exactly correspond to the initial de nition of uin the interior of the element, but will perfectly coincide with it on the nodes, and therefore the vector dcontains always the real displacement of the element through its nodes. This is considered as enough information to be extracted from the displacement eld as there's generally no interest in studying it in detail considering that the vector [ juj] already contains embedded crack displacement information. The advantages emerge when dealing with Eq. 7a, as the simpli cation of the displacement eld (especially the displacement variation eld) will allow to reach the classical standard FEM global solution form without any additional terms involved. Note that Eq. 9 does notimply Eq. 10 in any way: having the same interpolation matrix Nis totally a strategic choice that is convenient at this point. This present work will comply with it. 7 3.1.2. Strain eld discretisation Discretisation of the strain eld "is done following all the terms described in Eq. 5, giving this eld the role for modeling all strong discontinuity kinematics on the element, including bound and unbound terms: "=Bd+Gs[juj] =Bd+ (Gsb+Gsu) [juj]; (11) where Bis the standard matrix corresponding to the partial di erentiations of NandGsis the corresponding matrix for the strong discontinuity terms, including a bounded part Gsband an unbounded part Gsu. The speci c form of the matrix operators GsbandGsuis constructed by considering Eq. 5 alone. Their structure is xed and determined by the natural kinematics of the strong discontinuity model assumed. A strain eld following such idea within the E-FEM framework is referred to as a Kinematically Enhanced Strain (KES) . The choice for the strain variation eld stays analogous by considering the same standard strain term associated with standard displacement variations, but leaving a choice open for a matrix G swhich, for reasons explained later (and also a deliberate choice), will adopt the same bounded-unbounded composition G sb,G suand will be associated to the variation of the crack displacements [juj]: "=Bd+G s[juj] =Bd+ (G sb+G sn)[juj] (12) In this case, the speci c form of the matrix operators G sbandG suis not determined by strong discontinuity kinematics directly. It is rather derived from the demands of the variational principle itself through the satisfaction of Eq. 7, as it will be shown in the following sections. This means in principle that the components of this strain eld attend more to static considerations . A strain eld following such idea within the E-FEM framework is referred to as a Enhanced Assumed Strain (EAS) . 3.1.3. Stress eld discretisation Finally, the independent stress and stress variation elds are simply designated in a free-form fashion: =Ss (13a) =Ss; (13b) wheresis an independent stress vector with the state of stresses associated to each node of a given element (e.g., , a six-component vector for a 3D, constant eld). would be its corresponding variation. The interpolation matrices SandSare absolutely unrelated to any of the previously de ned entities in the other mechanical elds. 3.1.4. Calculated stress eld discretisation The stress eld coming from the constitutive relation (") will be calculated by taking only the bounded part of". The boundedness of (") despite of having unbounded terms within "has been discussed since the rst attempts of consolidating the E-FEM framework [15]. It remains an accepted practice in current theoretical developments. Considering this, this stress eld is expressed as: Calculated stress eld (") =C"b=C(Bd+Gsb[juj]) (14) 8 3.2. Application of the discretisation strategy This is probably the hardest and the most loosely described mathematical development section in the E-FEM framework while it determines most of the operational conditions of the numerical solution. The de nitions of the chosen discretisation strategy in the previous section are basically thrown into the varia- tional principle represented by Eq. 7 to study their mathematical implications. The intention in this section is to unfold key parts of the process that strongly regulate the structure of current formulations as well as of future proposals. 3.2.1. Basic orthogonality analysis: the role of S We will start by taking Eq. 7b, which brings the following expression after considering the de nitions of u,", the fact that [juj] is a constant and that sremains completely arbitrary: Z et(@u")dV=Z estSt(Bd[Bd+Gs[juj]])dV=stZ eStGsdV[juj] =0 )Z eStGsdV=0 (15) Eq. 15 demands an orthogonality between StandGsin a weak sense. Given that Gshas the speci c purpose of modelling the kinematics of the strong discontinuity, it shall not be modi ed to satisfy Eq. 15. Thus, Shas to comply. Note that this variational relation does not impose any form restrictions to critical functions, as Swill not give shape to any realmechanical eld. It could be stated that Sbecomes just a support interpolation matrix that allows Gsto exist. We can choose whatever kinematics we want for the strong discontinuity implied in Gs, as long as we can actually prove that this support interpolation matrix Ssatisfying Eq. 15 exists. At this point, the speci c form of 'comes into play, which is important to discuss because it determines the overall form of Gs.Gsin turn is composed by bounded and an unbounded parts GsbandGsu. Let us suppose now that we work on a linear tetrahedron, and thus with a linear '=a+bx+cy+dz, which is one of the most prominent choices made on the E-FEM literature. With this assumption and considering the structure of Eq. 5, expressions for matrix operators GsbandGsuon [juj] can be devised as: Gsb=2 6666664';x 0 0 0';y 0 0 0';z ';y';x 0 0';z';y ';z 0';x3 7777775=r';Gsu=2 6666664nx0 0 0ny0 0 0nz nynx0 0nzny nz0nx3 7777775=Hs (16) where'x;'y;'zare the partial derivatives of 'andnx;ny;nzare the components of the normal vector ^nto the crack plane. If we assume for a rst trial that Sremains constant, there will be only one single option to de ne it: the identity matrix I6(or a scalar multiple of it). Eq. 15 can be then worked using integral properties of the Dirac delta to show that: Z eStGsdV=Z eI6(r'+Hs)dV=Ver'+AHs; (17) whereVeis the element volume and Ais the area of the crack surface. It can be easily seen that, as r' andVdo not have any kind of relation to AandHs, Eq. 17 can never be asserted as zero and thus a constant Ssatisfying Eq. 15 does not exist. Note that a constant Swas chosen to be consistent with the natural order for a stress eld associated to a standard linear displacement eld on a linear tetrahedron. There is no option but to go for a non-conventional de nition for S. Fortunately, Hu-Washizu allows for this. 9 Using a linear Swould imply a linear interpolation matrix based on a series of four points ( x i;y i;z i), which would make a total of 12 free parameters to achieve the sought orthogonality. Let us de ne the four associated interpolation functions S i=a i+b ix+c iy+d izto these points. Eq. 15 can then be worked as: 0 @Z eStdV1 Ar'0 @Z AStdA1 AHs= 0; (18) leading to a linear system of the form: SVe i';j+SA inj= 0; i= 1;2;3;4j=x;y;z (19) whereSVe iare the interpolation functions S iintegrated through the element volume VeandSA iare these same functions but integrated over the crack surface A. As there are a total of 12 independent equations to ful ll and Shas 12 free parameters, a unique solution for Sexists. Note that this system will determine a very speci c set of nodal positions ( x i;y i;z i) that are absolutely unrelated to realelement nodal positions. The base of Sis not the standard element base and thus it will be entirely di erent of that coming from S, which has to be the real one. This means that, from this very point, the framework is changed from a classical Galerkin scheme to an asymmetric Petrov-Galerkin scheme , a fact that almost all authors in the current E-FEM literature ignore. For more complex choices of 'or even a non-constant [ juj], it can be shown by the aforementioned approach that a suciently elaborate Scan be chosen to satisfy Eq. 15. This last point is the main takeaway from the work on Eq. 7b. 3.2.2. The bridge between real and constitutive stress elds Eq 7c can be worked by using Eq. 12 and noting that dand[juj] remain independent, giving rise to two new equations: Z e"t(("))dV=Z e dtBt+[juj]tGt s ((")Ss)dV =dtZ eBt((")Ss)dV+[juj]tZ eGt s((")Ss)dV= 0 )Z Bt(")dV=Z BtSsdV (20a) )Z Gt s(")dV=Z Gt sSsdV (20b) Eq. 20a is the only equation that allows to establish a direct relation between the calculated stresses from the constitutive law (") and the interpolated realstress Ss, since the form of Bis known beforehand for a given standard element geometry. Eq. 20a is a weak equality involving Btas a weight factor. The form ofSalways remains a choice. For a constant real stress eld, Swill be a 66 identity matrix, whereas for a linear real stress interpolation, Swill be a dense 624 matrix, and so on. Based on this choice and the form ofB, Eq. 20a may or may not grant a strong equality between constitutive stresses and real interpolated stresses. In the case that Sis chosen to match the order of ("), Eq. 20a e ectively reduces to a strong equality. Otherwise, for the case we have a constant Band a constant S, a real stress vector scan still be calculated though a simple volume averaged integral regardless of the order of ("): 10 Volume-averaged real stress calculation s=1 VeZ (")dV (21) Note that, in any case, real stress calculations remain explicit and do not require for methods like a least-squares approach like in other mixed element variational frameworks [20, 1]. 3.2.3. The bridge between real and constitutive traction vectors Eq. 20b is treated di erently. As mentioned before, G sis assumed to be composed by bounded and unbounded terms as its counterpart Gs, and in this section it will be clari ed why. Usually, authors in the eldchoose at this point to enforce weak orthogonality conditions to each separate side of Eq. 20b since it becomes easier to devise the traction vector explicitly [14]. The present analysis will continue to do so. This creates once again two new relations: Z Gt sSsdV=Z Gt s(")dV= 0 (22) The main idea is to derive an expression for the realtraction vector Tbased on the realstress eld ( Ss), and then connecting it to the traction vector based on constitutive stress T", so thatTcan be explicitly calculated. This can be done by the developing the left and center side of Eq. 22, leading to: Z Gt sbSsdV+Z Gt suSsdV=Z Gt sb(")dV+Z Gt su(")dV (23) If we propose the speci c form of the unbounded G suas exactly the same as Gsu=Hs, we can use again the Dirac delta's properties to arrive to: Z Gt sbSsdV+Z dHt sSsdA=Z Gt sb(")dV+Z dHt s(")dA )Z dTdA=Z Gt sb[(")Ss]dV+Z dT"dA; (24) where it has been recognized that Ht sremains a stress projection operator over the crack surface having ^nas its normal, thus returning traction vectors at once operating over both stress elds. The main intent of assuming G suis to e ectively reveal the traction vectors TandT". Depending on the stress relation coming from Eq. 20a, Eq. 24 might be simpli ed. A strong equality between TandT"is possible only if we can ensure that the respective stress elds are constant and equal to each other. 3.2.4. EAS and static considerations - The patch test condition The zero implied on the left side of Eq. 22 allows to de ne G sunder the EAS approach. The arbitrary vectorscan be then taken out for reaching the following weak orthogonality condition: Z Gt sSdV= 0 (25) As the form Sis chosen based on real physics, it is clear that G smust be determined by it. The works on Eq. 24 have already required the unbounded part G suto be equal to H sin order to obtain the traction 11 vectors, so only the bounded part G sbcan be used to satisfy Eq. 25: Z Gt sbSdV+Ht sZ dSdA= 0; (26) where the nal form of G sbwill basically depend on the choice of S. At this stage, authors generally also seek to comply with the traditional FE patch test to ensure element density convergence for a given parent element by verifying that it is able to work with constant stress elds regardless of its underlying complexity [20, 1, 14]. This implies to additionally satisfy Eq. 26 for the case of a constant S, regardless of the form already chosen for it. This is the reason why most of the authors working on the E-FEM framework try to choose a constant Sfrom the beginning, so that only one single solution for Eq. 26 is required, satisfying the patch test for once and for all. For that case indeed, solution to Eq. 26 simpli es to the expression commonly found in E-FEM literature for the SKON, EAS approach [21, 22]: EAS G sboperator for constant real stress G sb=A VeHs; (27) for which G sbhas also been set as a constant matrix, as well. For other cases that are less restrictive, G sbhas to rise in complexity by simultaneously satisfying Eq. 26 for both a constant and a higher degree S, in a process analogous to that one followed during the analysis in Section 3.2.1. This leads to a linear system on the coecients for G sbsimilar to that of Eq. 19. Note that, no matter what the choice is, it inevitably leads to have Gsb6=G sb, again departing from a classical Galerkin scheme to a Petrov-Garlekin one. This is the only way to retain consistent kinematic representation while remaining variationally consistent under the light of Hu-Washizu. In the proposals made on the next sections, we will retain the choice of a constant real stress eld, along with a constant S and the validity of Eq. 27. 3.2.5. Final traction calculation If a constant G sbis retained, Eqs. 24 and 21 allow further simpli cations to get to a weak equality between traction vectors: T"=T=T; (28) where we have used the fact that a constant yields a constant T. Having cleared this link, the right hand of Eq. 22 nally leads a way to explicitly calculate a surface average of the traction vector T, for then applying any pertinent damage, softening or traction-separation laws. The calculation is simply stated as: T=1 AZ dT"dA=1 AZ Gt sb(")dV (29) From now on, traction vector calculations will be simply referring to T. 3.2.6. Internal-external force balance This part of the variational analysis has been left to the last as it doesn't bring any remarkable implica- tions. Eq. 7a can be stated as follows: Z BtdV=Z NtfbdV+Z NttdA; (30) 12 where Eq. 20a can be used literally as it is to vanish the real stress eld from the expression and leave all force balance calculations in terms of the calculated stress eld ("): Z Bt(")dV=Z NtfbdV+Z NttdA (31) Note that no strong equality between the stress elds is required to satisfy Eq. 31 in any way. With this last equation, the application of the complete variational principle for a wide range of cases typically considered on the E-FEM framework is nished. The present work will now take all the basics presented in Section 2 and 3 to describe the general inconsistency problems on commonly proposed versions of the formulation in a precise manner. 4. Current formulation approaches and associated pathologies The following analysis takes Eq. 29 as a departure point. After developing G sband("), the following form can be reached: General traction vector equation T=Te+M[juj]^n; (32) whereTeis a traction vector depending only on the standard element displacement vector d, and Mis a crack sti ness matrix. The vector [ juj]^nis the crack displacement vector expressed on the local frame ^n. Eq. 32, which is expressed entirely on the local crack base ( ^n;^t;^m), is by itself probably the most important yet one of the most overlooked variationally-based expressions in this strong discontinuity analysis framework. Indeed, it grants information about the in uence of the load-driven and the crack separation- driven parts on the overall evolution of the fracture traction T. The latter determines, along the chosen softening laws, the complete post-localisation mechanics of the element, including what happens in terminal separation conditions. Note that Eq. 32 has been derived merely from the discretisation strategy of the framework and the pure kinematics model proposed for the strong discontinuity itself. It has absolutely nothing to do with the application of any further traction-separation law schemes. The following subsections describe di erent ways to handle Eq. 32 depending on the general fracture modeling approach. 4.1. Single mode formulations Single mode formulations are those that consider a single kinematic mode associated to the fracture within an element, normally a rigid body normal separation or a parallel sliding distance (the latter assuming a uniform fracture plane  d). This implies the suppression of all other fracture kinematic modes. For illustrating the implications of this choice, Eq. 32 can be taken under terminal separation conditions: this is, assuming all averaged traction components are driven down to zero after the crack has suciently developed. We present the resulting equation as a system taking each of the local fracture plane directions: Tn=Ten+Mnn[juj]n+Mnt[juj]t+Mnm[juj]t= 0 Tt=Tet+Mtn[juj]n+Mtt[juj]t+Mtm[juj]t= 0 Tm=Tem+Mmn[juj]n+Mmt[juj]t+Mmm[juj]t= 0(33) As an example, considering a single normal separation mode formulation, the model will be entirely based on [juj]nand will automatically make [ juj]t= [juj]m= 0. For the system given by Eq. 33, this would imply: Tn=Ten+Mnn[juj]n= 0 Tt=Tet+Mtn[juj]n= 0 Tm=Tem+Mmn[juj]n= 0(34) 13 The updated Eq. Sys. 34 remains overconstrained and cannot be satis ed by a single value of [ juj]n. Authors choosing this approach will only take the traction-separation equation corresponding to the chosen fracture mode to be satis ed. In this case, only the rst equation in 34 depicting the general traction component Tnwill be driven to zero. A speci c [ juj]nvalue will ensure this condition. This [ juj]nwill evidently not satisfy the remaining traction separation relations in the system. As the components Tet;Temdepend on the arbitrary disposition of the nodal displacement vector d, general traction components Tt;Ttwill take totally non-physical values. This leads to the creation of undesired spurious internal forces that will produce severe stress locking. In order for single mode formulations to thrive during entire nonlinear FEM simulations, it is required todeliberately drive these crack surface traction components to zero during a numerical solution, along with any other stresses still associated to the crack. This can be done by simply applying an arbitrary decay explicitly to these or by setting them immediately to zero when localisation is detected on a given element [21]. The only conditions in which this kind of formulation will naturally drive all stresses down to zero is for a speci c combination of load (through d) and element geometry orientation with respect to  d(coecients Mnn;Mtn;Mmn). As it will be shown later in a more detailed analysis of Mfor a linear tetrahedron and a normal separation mode, these conditions correspond to purely axial strain on the ^ndirection, having an element face parallel to a planar  d. Because of the aforementioned reasons, single mode formulations do not ensure physical sense for crack internal variables. However, they have the evident advantage of their simplicity and numerical solution speed, as the system to solve for [ juj] becomes a single algebraic equation. As simple as it is, this kind of formulation has been successfully used for the modeling of complex fracture phenomena in heterogeneous, quasi-brittle materials [22, 23, 24]. It has granted satisfactory results for the global strength prediction and overall crack position of actual test specimens under tensile and compressive loads, even when considering triaxial con nement preloading conditions [25]. 4.2. Full crack translation formulations: the role of ' Full crack translation formulations consider all crack displacement components when articulating the traction-separation equation system. In terminal conditions, the full system depicted by Eq. 33 is to be solved for all variables [ juj]n;[juj]t;[juj]m. Regardless of the nature of the traction separation decay (exponential, polynomial, etc.), Eq. Sys. 33 remains strictly linear. A unique solution will always exist for [juj] given that Mis not singular. A detailed analysis of Mis thus relevant in this study. The crack sti ness matrix Mhas a direct dependence on G sband(") (Eq. 29).(") in turn depends directly on Gsb(Eq. 14), and thus on '(Eq. 16). As already discussed on Section 2.1, the nature of 'depends on the parent element, and authors typically take the simplest mathematical de nition of it satisfying the basic requirements in Eq. 3. The case of a linear tetrahedron parent element will be taken again as an illustrating reference. The simplest'is linear and easily constructed by stating linear interpolation functions taking the values of 1 or 0 depending on which nodes are on which side of the fracture surface. A 2D case is again depicted in Fig. 3. Note that it is not required to use the same function 'for all three displacement components [juj]n;[juj]t;[juj]m. Note that in this case there's only one possible de nition of a linear 'for this parent element, anyway. The standard deformation matrix Bcan be expressed as a block partition in a node-wise manner: B= B1B2B3B4 ; (35) 14 where each block Bican be further expressed as a function of the unit vector i= [ ix; iy; iz]tnormal to the opposing face to each node iin the tetrahedron (this is naturally coming from the basic FEM formulation of this parent element): Bi=Ai 3Ve2 6666664 ix 0 0 0 iy 0 0 0 iz iy ix 0 0 iz iy iz 0 ix3 7777775; (36) whereAiis the area of the opposite face to node i. Knowing that 'follows a similar structure as Bbut considering the selective condition in Eq. 3, Gsbcan be expressed as a superposition of the blocks de ned in Eq. 36: Bounded Gsboperator for a linear ' Gsb=NeX ipiBi; (37) wherepiadopts a value of 1 on the nodes lying on +, and zero otherwise. On the other hand, a constant eld will be assumed so that the EAS approach matrix G sbcan be obtained using Eq. 27. With Eq. 14 in mind, Eq. 29 can now be developed to get the following: Sti ness Matrix Mand load-driven traction Tecalculation 0 @1 VeZ HT sCBdV1 A |{z} Ked+0 @1 VeZ HT sCGsbdV1 A |{z} Km[juj] (38) T=RTKed+RTKmR[juj]^n=Te+M[juj]^n; (39) where Ris a rotation matrix allowing the passage from a global coordinate crack displacement vector [juj] to a local-based [ juj]^n. Working an explicit form for Myields nally: M=1 Ve2 4cn1^n cn2^t cn2^m cs^t cs^n 0 cs^m 0cs^n 3 5; (40) wheredenotes a dot product ( ^t =^tT ) and =AiP ifollows the superposition brought by Eq. 37. The constants cn1;cn2;csare a function of the linear elastic material model parameters E;. The reader can nd more details about the derivation of this compact expression for Mon the Appendix A. For now, let's assume a fracture plane that makes a 3-1 partition of the tetrahedron's nodes, having a single node on +(Fig. 5), so that the interpretation of becomes easier. For the case of the linear tetrahedron, it is geometrically impossible to have the orientation of the opposite face to the single node to be perpendicular to ^n, thus ^nwill be always di erent than zero. Given this fact and the structure given by Eq. 40, it can be safely assumed that Mwill not be singular. The system de ned in 33 will then have a unique solution for [ juj]n;[juj]t;[juj]m. It is relevant to note the presence of coupling fracture sti ness coecients between di erent directions (Mij;i6=j). This means that, even at terminal conditions (Eq. 33), a normal displacement [ juj]nof the 15 Figure 5: Tetrahedron element having a strong discontinuity with a basis ^n;^t;^mand a partition 3-1, node 1 located on + fractured piece of the element +will automatically induce a parallel rigid body motion [ juj]t;[juj]mas well. This does not make physical sense since at this stage +is already a completely independent body provisioned with free rigid body displacement modes. Another way to see this would be that the traction vector components at the fracture are strongly coupled through [ juj]n;[juj]t;[juj]m. Moving in a normal direction would mean to produce traction on the parallel direction. Again, this is not representative of a fractured body kinematics. Eq. 40 also states that all coecients Mijwill depend in general on the orientation of the parent element's geometry with respect to the fracture plane. This is, the Mijcoecients have an explicit mesh dependency . Thus, if we were to avoid any unwanted kinematic/traction coupling, the element geometry must have a speci c orientation with respect to the fracture plane. Otherwise, while a mathematical solution for the system given by Eq. 33 will be always assured, the overall physical sense of [ juj] will be once again compromised. Wells [9] had already identi ed the in uence of the function 'on the overall fracture kinematics and the induced strong coupling between traction components at the fracture plane. He proposed to x it by rede ning'in such a way that the Mcomponents having the dot products ^tand ^mwere driven to zero. This is, a 'capable of diagonalizing M. Illustrating the process in two dimensions, if an initial function'was de ned for a CST element on the xyplane having partial derivatives ';xand';y, the following rede nition for the partial derivatives was made: '0 ;x=n2 x';x+nxny';y (41a) '0 ;y=n2 y';y+nxny';x (41b) It is not necessary to check after 'directly, as only the partial derivatives of 'are needed during the numerical calculation. Indeed, such a rede nition of 'will avoid the unwanted coupling of traction components. However, this 'will no longer be able to satisfy the basic requirements of Eq. 3. Remember that this was the main reason for conceiving 'in the rst place. If this alternate de nition of 'is checked on the position x1shown in Fig. 6 where it is supposed to be zero, it can be shown that: '(x1)/j(^n 12)jj(x1^n)j; (42) wheredenotes a cross product, 12is the normal vector to the edge between nodes 1 and 2 and x1is the global position vector for node 1. This product of magnitudes nulli es whether the fracture line is aligned with the the normal of the opposed element edge or if it is also aligned with the position vector in question. The reader can nd more details on the Appendix B. 16 Figure 6: 2D CST element having a strong discontinuity with a vector ^nand a partition 2-1, node 3 located on +. The vector in this case corresponds to the normal to the surface de ned by nodes N2;N3;N4 . Again, a mesh dependency is introduced, albeit more silent and rather impacting the pursuit of global equilibrium. It is clear that a deep knowledge of the fundamentals of the E-FEM framework is required for any kind of improvement proposals. This is the reason why this work has been devoted on its rst half to describe its theoretical basis. 4.3. General crack kinematics formulations At terminal separation conditions, there should be no internal forces in the element associated to any rigid body motion of +with respect to . Otherwise, global strength assessment on a complete FEM model will never be realistic since these elemental contributions are out of control, regardless of their already damaged state. Liberalisation of the three rigid body translation modes for [ juj] ensures that it is possible to drive traction components Tn;Tm;Tndown to zero. Unfortunately, this does not ensure a complete release of element internal forces. Internal forces can be obtained extracting a part of the internal-external force balance on Eq. 31: fint=Z Bt(")dV=Z BtdV=VeBts; (43) where Eq. 20a and the fact that the present work stands for a constant stress eld ( =Ss=s) have been used. On the other hand, the kinematics proposed in Section 4.2 will only ensure the average traction vector Tacting upon  dis driven down to zero: T=1 AZ Gt sb(")dV=0 (44) For a constant stress parent element (constant G sb;B), it's obvious to see that Eq. 44 does not imply reaching zero on Eq. 43 since in general we have G sb/B. In this case, the whole real stress tensor on the element must be nulli ed so that all internal forces come to zero. Thanks to Eq. 44, we also know that some of the components of this stress tensor expressed on the local coordinate base ( ^n;^t;^m) correspond precisely to the traction vector components that have been already worked out in Section 4.2. There are thus four components left without any damage process: :2 4nn=Tnnt=Tt(((((nm=Tm tn=Tttttm (((((mn=Tmmtmm3 5 (45) These untouched components tt;tm;tt(only three since the tensor is symmetric) should be damaged in such a way that the framework retains mathematical consistency and physical meaningfulness. 17 In order to damage , Eq. 21 can be taken as a point of departure to relate a damaged state of stress to the constitutive stress ("). Then, Eq. 14 is used to see the implications on speci c terms within the formulation. At terminal conditions, nullifying will thus imply: Z dV=1 VeZ (")dV=1 VeZ C(Gsb[juj] +Bd)dV= 0 (46) 1 Ve0 @Z GsbdV1 A[juj] =1 Ve0 @Z BdV1 Ad (47) Gsb[juj] =Bd (48) At the end, Eq. 48 makes a relation between two average strain elds: one coming from standard dis- placements and the other coming from the strong discontinuity enrichment function. For the case of a linear tetrahedron, the average standard strain eld on the right side of Eq. 48 maps a R12vector to a R6subspace Sd. We also have B=B, asBis constant for this case. On the other hand, the left side enriched average strain maps a R3vector to a R6subspaceS[juj]. IfGsbhappens to be constant, then we would also have Gsb=Gsb. Despite the fact we can stretch more demands to the matrix Gsbthrough a richer de nition of', it can be shown that, in general, S[juj]kSd. Figure 7 illustrates this situation. It is impossible to nullify the strain subspace spanned by donly counting on a three rigid body displacement mode vector [juj]. Fracture kinematics hasto be enriched to at least 6 modes to do so. Thus, even full crack translation models might need to force a stress damaging mechanism such as in single mode formulations discussed in Section 4.1 to avoid stress locking. Figure 7: Generation of strain sub-spaces coming from dand [ juj]. The subspace generated by the fracture kinematics only represented through rigid body translations is not able to cover the subspace spanned by the nodal standard displacement vector. The idea of an extension for fracture kinematics variables on the E-FEM framework has been started by introducing the concept of a non-homogeneous (linear) crack displacement elds over a fracture surface in [26] in bi-dimensional elements. The concept was further consolidated by associating a set of fracture kinematic modes to this feature in [11, 12] and later applied in [8, 27, 28], also for two dimensional problems. The approach is to consider that each of the components of [ juj] are rather de ned using linear functions instead of constant values. Each of the parameters within these linear functions is associated to general kinematics for the split body representing the +domain. This can include rigid body translations, rigid body rotations and even simple axial strains in speci c directions. The centroid of the  dis generally taken as the reference position for the rigid body modes. In essence, fracture kinematics is enriched by providing + with improved modes that are able to represent more than just rigid body translations [ juj]n;[juj]t;[juj]m. For a 2D case, Fig. 8 illustrates how +kinematics on a CST element is enriched by considering two rigid body translations, a single rigid body rotation on the plane and a single extension (simple axial strain) on the parallel direction to  d. This is the type of enrichment studied by Armero and Linder [11]. For the sake 18 Figure 8: Enriched kinematics variables on the fracture local frame for a 2D CST element having a fracture line  dwith orientation vectors ^n;^tand nodesN1;N2lying on +. This kinematic variable set is composed by two rigid body translations [juj]n0;[juj]t0, one rigid body rotation and a single axial strain ton the direction parallel to  d. The centroid of  dis taken as the zero reference for rigid body motion description. of readability, all descriptions on the rest of this section will take this CST triangle example as reference, knowing that the proposal in later sections will make a 3D generalisation of this approach.The enrichment of fracture kinematics implies a rede nition for [ juj] that introduces a new intermediary matrix Jalong with a more general fracture kinematics vector : [juj] =J (49) =[juj]n0[juj]t0 tt(50) This enrichment has consequences on the mathematical framework developed on Sections 2 and 3, starting by the basic strain expression derived on Eq. 5, where the handling of gradients will change due to J. Additional interactions will also emerge between 'and this new matrix operator, yet giving rise to a di erent structure for the strain enrichment operators. Eqs. 47 and 48 will now be: 1 Ve0 @Z GsbJdV1 A=1 Ve0 @Z G0 sbdV1 A=1 Ve0 @Z BdV1 Ad (51) G0 sb=Bd; (52) where a new compound operator G0 sbalong with its averaged counterpart G0 sbhave been de ned. In order to build a new coherent structure for G0 sbunder this scheme, authors typically choose not to rework Eq. 5 and prefer to build G0 sbdirectly in a numerical, column by column fashion [11, 27]. Each of the columns of G0 sbrepresents a link between a fracture kinematic mode kand a speci c nodal displacement statedk: Kinematic mode linking, general relation G0 sbkk=Bdk (53) Here, the matrix G0 sbis given a physical meaning by linking kinematic states of the element standard nodes to each of the fracture kinematic modes. Since all other sections of the framework can be expressed as a function of G0 sb, this approach is enough to de ne the internal numerical solution process. 19 For example, suppose we would like to link the kinematic variable [ juj]n0to an actual set of nodal displacements d[juj]n0on the element. To do so, a real [ juj]n0is meant to represent a constant displacement of [juj]n0for nodesN1;N2lying on +in the ^ndirection. This speci c displacement state d[juj]n0will generate a standard strain through B=B. The negative of this strain has to be equal to the product of a speci c column of G0 sb(the rst one) and the kinematic mode [ juj]n0. The block structure of B=B1B2B3 can be used to our advantage: G0 sb[juj]n0[juj]n0=Bd[juj]n0 G0 sb[juj]n0[juj]n0=B1B2B32 4[juj]n0^n [juj]n0^n 03 5=NeX ipi(Bi^n) [juj]n0 G0 sb[juj]n0=NeX ipi(Bi^n) (54) The remaining columns G0 sb[juj]t0;G0 sb;G0 sbtcan be deducted following the same logic. This process allows to numerically de ne G0 sbcompletely based on basic parent element characteristics. One remark to this process is the nature of the linkdone between nodal displacements and kinematic variables (Eq. 53). A given mapping from the R12displacement space to the R6strain space is not unique , meaning that two or more di erent displacement states dkmight correspond to exactly the same standard strain vector. This means that, even if we create a link between a speci c fracture kinematic variable kto a desired displacement state through a given dk, the strain produced by dkcould be also produced by other entirely di erent d0 k, breaking the meaning of this link, and thus the physical meaning of k. This is specially true when trying to distinguish what happens on +or . There might be cases in which nodes lying on will contribute to a vector dcapable of exciting unexpected modes kthat were originally associated to the kinematics of +. Rigid body modes are re exive in the sense that they can be interpreted by xing a reference frame whether at or +, avoiding the risk of any mapping confusion. Simple axial strain modes are not, and these can be confused by this imperfect bijection. There is thus an inevitable kinematic coupling between +and that prevents the approach to be fully objective. There is some work that has been done, as an example, to consider a more general de nition for the factor pi: make it a continuous variable from 0 to 1 to control the coupling between +and in a more exible manner [29]. As this doesn't really solve the root problem, the authors of the present work have decided to keep the classical binary de nition for pi. Another issue with the approach, pointed out well by Armero and Linder [11], is the fact that it is assumed that the column decomposition of G0 sbyields a set of vectors that compose a linearly independent base for building a "space. For some node partitions on +; , speci cally when we have a single node on +, it is found that the base will not be linearly independent. This means that there is no way to know a unique solution for all kmodes for a given element kinematic state. Let's take the same CST triangle as an example, for a 1-2 node domain partition having node 1 on +. For a given displacement d1of this node, there are in nite ways to make a [ juj]n0;[juj]t0;;tdistribution in such a way that the resulting + motion returns the same d1displacement on node 1 (See Fig. 9). [11] propose a numerical workaround to this problem by making use of Lagrange multipliers to relax unconstrained parameters during the solution process. This still leaves the question of physical meaningfulness and legitimacy of one solution over any other. Finally, we already know the risk of forcing a de nition of any Gsb-related operators without making a direct assessment of '.None of the previous relations actually ensures satisfaction of any of the main requirements for 'in Eq. 3. [11] mention the possibility to integrate G0 sbto inquire about 'but this process is far from trivial given the multiple gradients implied in Eq. 5. A lot of information about the speci c 20 Figure 9: Illustration of linearly dependent kinematic mode decomposition for a CST having one single node on the +domain. Two di erent paths using di erent values for the kinematic modes are portrayed. Both paths (and endless others), lead to d1 in this case. 'giving rise to the given G0 sbhas been already lost at this point. The only evident fact would be that ' remains far from linear. The next section will take this development state as a departure point, and will continue to work on a 3D setting. 5. Formulation approach proposal for three dimensional problems This section makes a formulation proposal considering all the fundamentals from Sections 2 and 3 as well as the issues found during the evolution of di erent strong discontinuity modeling approaches described in section 4. At the end, this proposal is not meant to eradicate all the theoretical faults analysed so far. It is rather meant to make its best to reach a sound compromise between a theoretically sound E-FEM formulation but still operationally attractive. The key is to really retain mathematical robustness, ease of implementation and mesh independence as much as possible, which are the key characteristics that make the E-FEM framework attractive. 5.1. A consistent enrichment for ' In section 4.2 the role of the function 'has been discussed. 'allows for a consistent boundary condition imposition, but it also introduces a mesh dependency that generates an unwanted coupling between traction components at the fracture surface. This in turn compromises the physical meaning of fracture displacement variables [juj]n;[juj]t;[juj]m. [9] proposed a x, but is not completely e ective due to a violation of the mathematical framework. The theory described in Section 2 and 3 can be used to make a consistent 'proposal. The idea is to propose a'with enough exibility to both satisfy the requirements in Eq. 3 and the nulli cation of the out-of-diagonal terms of the Mmatrix in Eq. 40, as well as any other eventual requirements . Recalling the analysis made in Section 3.2.1, it is allowed to increase the complexity for the de nition for 'as long as it remains within the C0class of functions (continuous, derivative not necessarily continuous). It doesn't have to follow the same polynomial degree than the parent element. As an example, for a linear tetrahedron, a complete quadratic de nition for 'is perfectly valid: 21 Controlled/parametrised de nition for the 'function '=PT 2 = 0 1 2 3 12 23 13 11 22 33 P2=1x y z xy yz xz x2y2z2t(55) This increase in order for 'implies an increase in order for the element strain eld "through Gsb(Eq. 11), and thus for the calculated stress eld ("). If the realstress eldis chosen to follow the natural order in the parent element (a constant eld for a linear tetrahedron), the strong constitutive equality (Eq. 8b) will not be possible and stresses will have to be calculated through volume averages (Eq. 21). On the other hand, the advantage of having a discordance of (") with respect to a constant  eld is that absolutely nothing will change in the treatment of the EAS part of the framework (Section 3.2.4). This means that Eq. 27 can still be used for managing a constant expression for G sb. The integrand in Eq. 29 is not constant anymore. This means that the integration process for all traction calculations will have to consider the geometric parameters of  das well as the speci c geometries of the subvolumes +and . Eventually, after vanishing all variable terms of the base P2through de nite integration, an expression for the fracture sti ness matrix Mas a function of 'coecients will be attained. As per Eq. 40, we know that Mremains almost symmetric, where the respective out-of-diagonal term pairs (Mij;Mji) share the same internal products ( ^tT );(^mT ), and that there's already a pair ( Mtm;Mmt) which is zero. Thus, nullifying the out-of-diagonal terms in Mfor avoiding unwanted traction couplings implies the imposition of two additional linear equations on the coecients: Manual diagonalisation of the fracture sti ness matrix M ^tT ( ) = 0!Mnt( ) = 0; Mtn( ) = 0 (56a) ^mT ( ) = 0!Mmt( ) = 0; Mtm( ) = 0 (56b) Adding this to the basic 'requirements (Eq. 3, four linear equations on k') makes a total of 6 linear constraints. As a complete quadratic base P2allows for a total of 10 free parameters, it is completely feasible to design a 'function capable of satisfying BC imposition consistency at the same time that mechanical requirements for traction behavior. While the relations 56a and 56b ensure diagonalisation of the fracture sti ness matrix M, there's nothing that ensures that the magnitude of these diagonal components is physically meaningful. Each diagonal term Mjjindependently controls the fracture sti ness associated to each displacement [ juj]j. These displacements should be consistent with the element standard nodal displacement vector d. As expected, 'can also help on this errand by imposing additional conditions for the diagonal fracture sti ness terms. Note that, as per Eq. 40, the diagonal terms are notindependent, and only one single constraint would be needed to regulate them all. The case depicted in Fig. 4 is taken as a reference for devising such constraints. For a 2-1 node partition on a 2D triangle, it is easy to nd an exact relation between, for instance, the fracture normal separation [juj]nand the normal displacement  dnof a single isolated node partition on +: they should be equal . For other con gurations however, like a 2-2 node partition on a tetrahedron, it is not evident to assign a constraint on two independent nodes versus a single rigid body displacement: there are multiple ways to do 22 so and none of them would be rigorously correct. This work proposes to take the average separation  d between the groups of nodes on +and : d=d+d(57a) d=PpidiPpiP(1pi)diP(1pi); (57b) Taking the normal direction ^n, we require: [juj]n= dn=^nTd (58) At this point, it will be assumed that 'has been already arranged as to diagonalise Mthrough Eqs. 56a, 56b. With this, at terminal conditions, we can practically use the rst equation of the system in Eq. 33 to nally arrive to a linear constraint for Mnn: Ten+Mnn[juj]n= 0 ^nTKed+Mnn ^nTd = 0 Fracture sti ness attuning for terminal separation conditions Mnn( ) =Mnn1=^nTKed ^nT FtPpiIT 3(4)FT P(1pi) d(59) Proper block matrix algebra has been used to express this constraint in terms of the complete vector d, using some auxiliary matrix blocks: FT=p1I3p2I3p3I3p4I3 (60a) IT 3(4)=I3I3I3I3 ; (60b) having the 33 identity matrix I3as the basic block. Note that in Eq. 59, we are de ning a component of Min terms of a nodal displacement vector d, which is the load imposed to the element. As the solution for all coecients takes place only once when reaching localisation, the load vector taken into account is precisely that corresponding to the closest localisation state for the element. For the kinematics to remain perfectly consistent under the de nition of thisMnnsti ness, the load path would have to remain proportional . This is highly unlikely during a global fracture phenomenon, since there are di erent loading and relaxation stages within the same element that would makedto change direction in an aggressive fashion. Mnncan be updated as well in such cases. This will imply more than one solution for the coecients. Imposing the constraint on a di erent component Mtt;Mmmwill yield the same overall results given that a scaling factor of cs=cn1is used as per Eq. 40. As this represents only one additional linear constraint on the system for the coecients (making 7 constraints and 10 free parameters), once again it is still possible to build a proper 'function handling all requirements at once. The next discussion will continue to explore the exibility of 'under the light of fracture kinematics enrichment for 3D elements. 23 5.2. Fracture kinematics enrichment in 3D coexisting with ' In Section 4.3, an introduction has been made to a more general de nition of fracture kinematics. The approach has been partially illustrated for a 2D parent element following the work in [11], where an explicit 'is practically disregarded. This section will describe a detailed application for a 3D element, including the clear distinction and managing of the 'function. The process starts with the same de nition of [ juj] as a linear eld vector, but for three dimensions in the local frame: [juj] =2 4[juj]n [juj]t [juj]m3 5=2 4an+bn+cn+dn at+bt+ct+dt am+bm+cm+dm3 5; (61) where each of the parameters aj;bj;cj;djis to be related to di erent fracture kinematic modes k. It can be shown that such linear eld structure allows for the modelling of three rigid body translations [ juj]n0, [juj]t0, [juj]m0, three rigid body rotations n;t;mand three axial deformations (associated to +)n;t;m. Figure 10: Relations between the parameters ai;bi;ci;diand the rigid body rotation modes n;t;m. Note that two parameters may be in uenced by the same rotation mode. Each plot illustrates the e ect of having a small rotation in the axes corresponding ton;t;m, respectively. Figure 10 illustrates how to build the relations between all 3D rotation modes and the parameters aj;bj;cj;dj. The general idea is to make a careful geometric interpretation for the meaning of each pa- rameter by introducing a small variation while keeping all the remaining parameters at zero. This helps to build simple equations to relate these parameters to each fracture mode k. At the end, if small-angle approximations are allowed, the elds are nally expressed as: 2 4[juj]n [juj]t [juj]m3 5=2 4m+t+nn+ [juj]n0 ttn+m+ [juj]t0 n+mmt+ [juj]m03 5=J (62) 24 From this point, an expression for the intermediary matrix Jcan be readily made having a de nition for the vector: J=2 41 0 0 0   0 0 0 1 0000 0 0 10 0 0 3 52 6666666666664[juj]n0 [juj]t0 [juj]m0 n t m nn tt mm3 7777777777775(63) With this analysis, Eq. 5 can be worked for a general de nition for the strain eld "considering both  and'. As all structures having to do with the enriched kinematics and Jare in the local fracture surface frame (;; ), it makes sense to work all further developments on the local frame from now on. The same bounded-unbounded structure for "can be derived: G0 soperators, enriched modes formulation "=Bd+ [G0 sb+G0 su] (64a) G0 sb=HrJ('rJ)r'J (64b) G0 su=HsJ (64c) The term ('rJ) has been grouped because it is notthe product of a deformation gradient operator on Jand', but rather a compound operator. The Heaviside function Halso forces a distinction between a bounded operator de ned on +(G0+ sb) and on (G0 sb). In general, ( G0 sb) will be an explicit function of 'coecients . To start the kinematic linking process between fracture modes and nodal displacements followed in Section 4.3, the requirement stated on Eq. 52 can be worked on a subvolume-basis: G0 sb=1 VZ G0 sbdV=Bd; (65) where the averaged operators G0 sbwill be explicitly linear matrices on the coecients . The linking of each kinematic mode will establish a series of linear equations on each of the column elements of G0 sb. For instance, the rst mode [ juj]n0already studied in Section 4.3 for two dimensions, would require satisfaction of exactly the same form of Eq. 54. Fortunately, it can be shown that the rst column of G0 sbhas exactly the same zeros as Bi, so that in local base the following would be obtained: G0 sb[juj]n0=2 666666664G0 sb1[juj]n0( ) 0 0 G0 sb4[juj]n0( ) 0 G0 sb6[juj]n0( )3 777777775=NeX ipi2 6666664Bni 0 0 Bti 0 Bmi3 7777775; (66) where the rst index lin theG0 sbl[juj]n0scalars is just the row placement. By looking at this structure, it can be seen that linking this kinematic mode requires satisfaction of three linear equations on the coecients. 25 The speci c expressions for G0 sb1[juj]n0;G0 sb4[juj]n0;G0 sb6[juj]n0have to be retrieved by working 'andin Eq. 64b and then taking the volume averaged integral demanded in Eq. 65. Note that 'is to be expressed as a function of local coordinates ( ;; ), and subvolume coordinate description in the local frame will be also required for evaluating the volume-averaged integrals. The linking of some kinematic modes can be more demanding than others. Some will even repeat linear equations and these can be naturally omitted, e.g., the rst three columns of G0 sbwhich contain repetitions of the simple partial derivatives ';;';;';. This depends on the symmetry of the kinematics de ned and the model of the strong discontinuity overall. At the end, it is true that the demands for 'will do nothing but to increase, as the 'free parameters are the only way to satisfy these relations. As a result, a full quadratic base proposal P2might not be enough to ful ll all these requirements additionally from Eq. 3 and other constraints on the fracture sti ness matrix M: nulli cation of non-diagonal terms (Sys. 56) and terminal separation conditions (Eq. 59). While increasing the polynomial base order of 'might a quick and tempting option, there are also other options for improvement: •As mentioned already in Section 2.1, it is not required to have the same 'structure associated for all fracture displacement components [ juj]n;[juj]t;[juj]m. Indeed, three di erent functions 'n;'t;'m may coexist in the model: '2 4[juj]n [juj]t [juj]m3 5!2 4'n0 0 0't0 0 0'm3 52 4[juj]n [juj]t [juj]m3 5 (67) While this might apparently triple the amount of free parameters, some of the model symmetries that allowed for redundancy simpli cations will no longer emerge. As an example, the basic deformation gradient operator on '(the rst three columns of G0 sb) will now have all di erent terms with respect to the one developed with a single ': 2 6666664'; 0 0 0'; 0 0 0'; ';'; 0 0';'; '; 0';3 7777775!2 6666664'n; 0 0 0't; 0 0 0 'm;& 'n;'t; 0 0't;&'m; 'n;& 0'm;3 7777775(68) On the other hand, basic 'requirements in Eq. 3 will now require the triple of linear relations for boundary condition consistency: one set for each 'n;'t;'m. 'n(xi) =pi (69a) 't(xi) =pi (69b) 'm(xi) =pi (69c) While there is evidently a trade-o , there's still a gain in the e ective number of free parameters in the global'structure, without the need to modify the complexity of the algebraic base of it. •In general,'has only aC0continuity requirement. This means that while the function itself is required to be continuous through space, its derivatives are not. This allows for a piece-wise de nition for '. The most natural choice is to propose a rst function for +and then a second one for . Piece-wise'does not increase the number of equations for Eq. 3 requirements, but it still breaks some of the model symmetries. There will be also new linear equations to satisfy, which are associated to the basicC0continuity of 'at d: ' d='+ d; (70) which will depend on the nature of the base Pchosen for the structure of '. For aP2base,C0 continuity requires six linear relations. Again, the gain in free parameters is worth the trade-o . 26 By considering the combination of both options, 'can globally add up to 120 free parameters using a completeP3base. While this might seem just too much, the E-FEM model managed in [24] which features also a weak discontinuity to represent di erent material phases within the element, requires approximately 110 linear relations to reach full consistency as described in this work. This process is carried out only one time for each element that has reached localisation. This last model is the one de nitively used for later element simulations in this work. The reader can nd the detail on the mode linking expressions as well as general 'handling on Appendix C. 5.3. A comment on linear system handling Calculation of all 'parameters based on all kinematic considerations may represent a formidable linear system depending on the parent element needs. The more complex the overall de nition for 'is (Appendix C), the hardest is to ensure the well-posedness of this linear system on the coecients. This is basically the price to pay when building a more robust kinematic model depending on '. On one hand, the linear dependency issues arising for some conditions during the construction of the columns of G0 sb(discussed in Section 4.3) will still be present on the 3D proposal. This will render sometimes the linear system globally overconstrained (for the total number of equations) and locally underconstrained (impossible to nd an unique solution for some groups of parameters). On the other hand, for complex 'structures many symmetries in the overall kinematic model are lost. As an example, the structure of the Mmatrix deduced in Eq. 40 is not guaranteed for a piece-wise, direction dependent,P3de nition for '. This implies that more coecients have to be used to re-enforce some desired constraints. It is not known by the authors of this work if this leads to the introduction of additional linear relations that might contradict some others already present in the system, thus compromising well- posedness. To avoid the need of introducing complex numerical relaxation mechanisms during the solution of the coecients system, this work has used a Singular Value Decomposition (SVD) approach [30]. The SVD approach basically takes the system, whether under or overconstrained, and solves whether for the least- square-optimal solution for overconstraints or for the minimum value possible for free independent variables if underconstrained. The SVD decomposition is only made once the element reaches localisation. Once de nite values are obtained for the coecients, the G0 sbcan be numerically populated and these can be used for the rest of the load steps in the numerical analysis. 5.4. Further treatment of the traction separation law system The introduction of additional fracture variables poses questions on how to handle the traction-separation part of the framework. The procedure depends basically on how the virtual fracture displacement eld [juj] is discretized and the impact on the structure of Eq. 25. While the realfracture displacement eld vector [juj] e ectively requires a 9 variable interpolation matrix through 'andJ, the virtual counterpart [juj] does not have to share the same structure. [juj] can still be assumed as a three constant eld vector, or it can take the same detailed description as [ juj]. The Hu-Washizu framework allows this. In this case, Eq. 25 retains exactly the same structure having the real stress interpolation matrix Sand the EAS operator G sb. Despite the fact that now we have a G0 sbconsidering the related upgrades, the EAS de nition for G sbdoesn't have to consider any structures from the enriched fracture kinematics. The only requirement for the EAS de nition of G sbis just to have the same bounded-unbounded composition featuring a Dirac delta and a projection operator Hsas a minimum. This will allow the emergence and equality of the traction vector terms from Eq. 24 onwards. 27 Eq. 32 will return a new system with an enriched set of 9 variables in : T=Te+M; (71) whereTeremains exactly the same as per Eq. 38. Note that if the work is already done on the local frame, there is no need to use the rotation matrix Ranymore. The expression for base fracture sti ness matrix M will have an updated de nition: Crack sti ness matrix M, enriched modes formulation M=0 @1 VeZ HT sCG0 sbdV1 A=1 VeHT sCZ G0 sbdV M=1 VeHT sC V+G0+ sb+VG0 sb (72) Here, all matrices have been already assumed in the local frame to avoid further transformations, and the de nitions implied in Eq. 65 have been used. Note that matrix Mnow has dimensions 3 9. The original three traction separation laws associated to traction components Tn;Tt,Tmwill return an underdetermined system: Traction-separation system, enriched modes formulation Tn=Ten+X kMnkk=qn() (73a) Tt=Tet+X kMtkk=qt() (73b) Tm=Tem+X kMmkk=qm(); (73c) whereqn;qt;qmare decaying expressions controlling the overall magnitude of Tn;Tt;Tmuntil reaching fracture terminal conditions, where Tn;Tt;Tmshould be zero. To fully determine the fracture mechanics, six additional relations are needed. Three equations can be proposed to explicitly damage the ijterms left untouched by the original approach: tt;tm;mm(Eq. 45). This additional damaging process can be done through additional traction-separation equations with speci c decay behaviors. This brings a set of three additional equations to the system: Complementary damage equations tt=Ttt+X kMttkk=qtt() (74a) tm=Ttm+X kMtmkk=qtm() (74b) mm=Tmm+X kMmmkk=qmm(); (74c) 28 where the new load-driven terms Ttt;Ttm;Tmmare calculated using the same structure as Eq. 38, but now using stress projection operators Hst;Hsmin the remaining local directions ^t;^m, respectively: Hst=2 6666664tx0 0 0ty0 0 0tz tytx0 0tzty tz0tx3 7777775;in local frame : Hst=2 66666640 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 03 7777775(75) Hsm=2 6666664mx0 0 0my0 0 0mz mymx0 0mzmy mz0mx3 7777775;in local frame : Hsm=2 66666640 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 03 7777775; (76) and then using analogous expressions for new load-driven traction vectors Tet;Temfor nally projecting to ^t;^mand obtaining the desired components: Tet=Ketd=0 @1 VeZ HT stCBdV1 Ad=HT stCBd (77) Ttt=Tet^t; Ttm=Tet^m (78) Tem=Kemd=0 @1 VeZ HT smCBdV1 Ad=HT smCBd (79) Tmm=Tem^m; (80) The sti ness vectors Mtt;Mtm;Mmmare calculated the same way as the rows of the original Mmatrix, but using the newly de ned projection operators Hst;Hsm: Mt=1 VeHT stC V+G0+ sb+VG0 sb (81) Mtt=^tTMt;Mtm=^mTMt (82) Mm=1 VeHT smC V+G0+ sb+VG0 sb (83) Mmm=^mTMm (84) It's important to remark that introducing the system of Eq. 74 is not the same as enforcing an arti cial decay ontt;tm;ttas in Section 4.1. We're systematically involving all fracture modes kin both systems (Eqs. 73 and 74). Simultaneous solution of these systems will ensure kinematically and statically consistent values fornullifyingtt;tm;tt, which should not exist at all once the crack develops. The overall system in still remains underdetermined. The last remaining three equations are free to implement more e ects in fracture dynamics without having to resort to 'parameters exclusively. For example, three nal relations could be established to nish the uncoupling of rigid body translations within system in Eq. 73 by grouping and isolating the e ects of all other kinematic modes di erent than rigid body displacements: 29 System closing equations Mnnn+Mntt+Mnmm+Mnnn+Mntt+Mnmm= 0 (85a) Mtnn+Mttt+Mtmm+Mtnn+Mttt+Mtmm= 0 (85b) Mmnn+Mmtt+Mmmm+Mmnn+Mmtt+Mmmm= 0 (85c) This way, the system is completely closed and a complete solution for can be found at each load state dafter element localisation. The approach ensures that tt;tm;mmget progressively suppressed as the fracture develops. On the other hand, it might seem awkward to have a [juj] discretized very di erently with respect to [juj], but this already happens with other elds within the framework. 6. Elemental validations In this section, calculations with a single element have been done on many of the discussed formulations to illustrate the quality and pertinence of the results concerning fracture mechanics and element kinematics. Five models will be will be presented in a progressively complex fashion, going from a single mode approach up to the latest formulation proposal of Section 5. A graphical comparison will be made between all models, covering static behavior (evolution of the fracture traction and its related stress state) and kinematic behavior (consistency between internal kinematic modes and true element kinematics through d). This will allow to see how the framework gains robustness and physical meaningfulness overall as the formulations evolve. Many authors working on the 3D version of the framework perform element-level calculations with an idealized regular element geometry and fracture orientations. The load direction is also often well-aligned to element edges or faces. Fig 11 shows examples of such cases on a linear tetrahedron for tension and shear testing. As discussed in previous sections, these conditions will hide the framework de ciencies since these are often the conditions in which almost all formulation versions will work properly. These will be referred to as vanilla conditions. Obtaining meaningful results using vanilla conditions at element level does not guarantee successful numerical simulations with larger scale models. Figure 11: Schematics of typical vanilla elements taken for formulation testing at small scale. Left element for pure tension, right element for pure shear. It can be easily seen that even a simple large scale model such as a cubical material sample with a tetrahedral unstructured mesh in pure tension will generate far from vanilla conditions at local element level. If element fracture planes are generated through free and spontaneous element localisations, a complex 30 loading-unloading state will be observed with plane orientations having random partition types between nodes in all elements. In some elements, the combination of element geometry and the nodal displacement vectordwill be such that it will activate complex kinematic modes (Section 5.2), and a formulation not prepared for this will inevitably fail to return physical or even mathematically meaningful results. Arobust E-FEM formulation should be able to undertake such scenarios to be truly useful. From the view of the authors of this work, there's no point in moving forward with a formulation if it is not capable of passing a realistic element test. This is the reason why this section presents the results at element level. We have taken one of such non-vanilla element samples coming from a real large scale model, having a speci c geometry and fracture plane orientation. An approximate illustration of the chosen element (again, a linear tetrahedron) is shown in Fig. 12. It has a characteristic length of approximately 1mm. The reference frame is such that the normal axis to the fracture plane ^nis almost perfectly aligned with the global zdirection, but not to any of the tetrahedron's faces. The remaining parallel directions ^t;^mare not aligned with the x;ydirections. It has a 2-2 node partition with the volume of representing approximately 5% of total elemental volume while +having the rest. Figure 12: Two di erent views of a realistic element taken as the reference for the numerical simulations to test all formulations described in this work. The load displacement vector dhas been built in such a way that it helps to validate the kinematic consistency of the formulations. For this, load requirements are proposed as a function of the generalized fracture kinematic modes kdescribed in Section 5.2. A monotonic, nonlinear behavior is prescribed for each of these modes, reaching a nal target value. A load progression factor between 0 and 1 is used to calculate any intermediate value between the starting point (zero) and the target value. Figures 13a, 13b and 13c show the chosen behaviors for each k. For this case, a quadratic behavior has been assigned for rigid body translations, an inverted exponential behavior for rigid body rotations and a radical behavior for simple axial strains. The load demand overall remains a composition of strong rigid body displacements along with mild rotations and strains, having an emphasis on normal separation. Allkinematic modes are activated without exception. The load proposed ensures the element arrives to terminal separation conditions at target values. The load displacement vector dassociated to this particular evolution of the modes kcan be retrieved by knowing beforehand the nodes of the element corresponding to the +domain and the mode linking 31 equations discussed in Section 5.2. Once again, expressions have been explicitly detailed in the Appendix C. Starting from a load d0corresponding to element localisation, the e ective dprescription can be calculated as a superposition of all e ects: d( ) =d0+m;n;tX jd[juj]j( ) +m;n;tX jdj( ) +m;n;tX jdj( ); (86) whered[juj]jare the displacement vectors associated to rigid body translation demands, djfor rigid body rotations and djthe displacements related to simple axial strains. Since the load requirements are nally prescribed based precisely on a controlled local fracture kinematics behavior, a good E-FEM formulation should return a kinematics calculation as closest and most consistent as possible with the source fracture kinematics, now referred to as ref: the kinematic reference . As mentioned before, the element is assumed to have already gone through localisation with a given stress stateyat the fracture plane: 2 4ynnyntynm ytnyttytm ymnymtymm3 5!y=2 6666664ynn ynt ynm ytt ytm ymm3 7777775=2 66666649:0 0:0 0:0 2:5 0:138 1:4593 7777775MPa; whereynn;ytn;ymncorrespond to the components of the traction vector Tn;Tt;Tmat localisation. These are calculated through a speci c localisation criterion in the element, which in this case has been a Rankine criterion. The remaining components ytt;ytm;ymm are found by projecting the entire state of stresses to ^tand ^mexactly at this load level. The traction-separation equations consider an exponential decay law for all pertinent fracture trac- tion/stress components: qj=yjeynn GfI[juj]n0+ynt GfII[juj]t0+ynm GfII[juj]m0 ; (87) whereyjare each of the aforementioned initial fracture yield stress components. Note that the expo- nential argument is exclusively controlled by the fracture kinematics corresponding to rigid body trans- lations.GfI;GfIIare the associated fracture energies for fracture separation [ juj]n0and fracture sliding [juj]t0;[juj]m0, respectively. This is a particular proposal that will conveniently t all formulation types proposed so far, although richer proposals could be made for the generalized kinematics approaches. All calculations have been done using SageMath 8.7 [31] mathematical open software on a Jupyter notebook 5.7.6 [32] platform. This has included all symbolic integrations to handle all enhancement functions parameters and operators such as ',JandG0 sb. All numerical procedures such as SVD operations and linear system solutions have been treated in this platform, as well. Graphical visualisation of the realistic element in question (Fig. 12) was also originally done using these tools. In the following subsections, ve di erent formulation versions will be described in order of increasing complexity. At the end, a graphical comparison will be made on their kinematic and static behaviours: the evolution of the relevant fracture kinematic state variables and the fracture stress state variables over load progression starting from localisation. 32 0.0 0.2 0.4 0.6 0.8 1.0 Load factor β−15−10−50510152025Displacement ( µm)Prescribed rigid body displacements [|u|]n0- Target: 25 µm [|u|]t0- Target: -15 µm [|u|]m0- Target: 10 µm(a) Displacements: quadratic evolution 0.0 0.2 0.4 0.6 0.8 1.0 Load factor β−0.10−0.050.000.050.100.150.200.25rotation (mrad)Prescribed rigid body rotations θn- Target: 0.2 mrad θt- Target: -0.1 mrad θm- Target: 0.25 mrad (b) Rotations: inverse exponential evolution 0.0 0.2 0.4 0.6 0.8 1.0 Load factor β−0.2−0.10.00.10.20.3strain (10−3mm/mm)Prescribed simple axial strains /epsilon1n- Target: 1.0×10−4mm/mm /epsilon1t- Target: -2 .0×10−4mm/mm /epsilon1m- Target: 3.0×10−4mm/mm (c) Strains: radical evolution Figure 13: Proposed load evolution through a controlled behavior of generalized kinematic modes 6.1. Single mode formulation For a single mode approach, the calculation begins by considering rigid body translations [ juj]n, [juj]t, [juj]m, and having only one of them di erent than zero (Section 4.1). To follow the example, a normal separation calculation will be made, leaving thus [ juj]t= 0;[juj]m= 0. The'function remains linear and unique, so that Gsbadopts its simplest way possible as a constant matrix (Eq. 37). The approach can be summarized in Table 1. The approach considers only the rst equation of the traction-separation law system, and solves for [ juj]n. This solution implies nding the intersection between a straight line Ten+Mnn[juj]nand the nonlinear behavioryneynn GfI[juj]n. The solution can be expressed in a closed, explicit form using the main branch of the Lambert function W0[22]. After a solution for [ juj]nis obtained, the approach then proceeds to calculate the remaining statics of the fracture plane, mainly the three fracture traction components Tn;Tt;Tm. At last, the remaining components 33 Single Mode Formulation Kinematic modes : [juj]n Linear, known ':'= 0+ 01x+ 01y Fixed, constant G sb:Gsb=NeX ipiBi Traction-separation system :T=Te+M[juj]^n=q([juj]^n) Nullify [juj]t;[juj]m: [juj]t= 0;[juj]m= 0 Solve only for [juj]n:Ten+Mnn[juj]n=ynneynn GfI[juj]n Calculate Tn;Tt;Tm:Tn=Ten+Mnn[juj]n Tt=Tet+Mtn[juj]n Tm=Tem+Mmn[juj]n Calculate:=C"=C(Bd+Gsb^n[juj]n) Determine tt;tm;mm:tt=^tT^t; tm=^tT^m mm=^mT^m(88) Table 1: Single mode formulation approach de nition of the state of stresses tt;tm;mmare retrieved for reference by calculating the entire stress (Eq. 14) and its projection to the corresponding directions. 6.2. Full translation formulations The next step is to consider formulations that make complete liberalisation of all fracture displacement components [juj]n, [juj]t, [juj]m(Section 4.2). Two versions from this family are considered for elemental simulations, depending on their assumptions taken on the 'function. One considers a classical linear de nition for '(and therefore a constant Gsb), representing the direct extension of the single mode approach of the previous section. The other builds a quadratic 'that helps the diagonalisation of the Mmatrix as well as for kinematic consistency at terminal separation conditions. Regardless of the 'assumptions, these formulations require a solution for the three-variable nonlinear system: T=Te+M[juj] =yeynn GfI[juj]n+ynt GfII[juj]t+ynm GfII[juj]m ; (89) which still has a closed, analytical solution using condensation techniques and the same Lambert W function W0approach. The rst three-mode approach is summarized in Table 2. For this case, the calculation of the traction vector components Tn;Tt;Tmand the stress components tt;tm;mmis done exactly the same way as with the single mode approach using the solution values for [juj]n;[juj]t;[juj]m. The other three-mode approach will go forward and consider a more complex de nition for the 'function, namely a quadratic one, as discussed at the beginning of Section 5.1 supported by a complete P2quadratic polynomial base (Eq. 55). Ten free parameters ( coecients) are available with this 'proposal. These will be used to build a linear system satisfying speci c kinematic and static conditions: four of them to satisfy the basic boundary condition imposition requirements (Eq. 3), two for driving out-of-diagonal Mijfracture sti ness terms to 34 Three Modes Formulation, Fixed ' Kinematic modes : [juj]n;[juj]t;[juj]m Linear, known ':'= 0+ 01x+ 01y Fixed, constant G sb:Gsb=NeX ipiBi Traction-separation system & Solve for [juj]n;[juj]t;[juj]m:T=Te+M[juj]^n=q([juj]^n) Calculate Tn;Tt;Tm% Calculate:=C"=C(Bd+GsbR[juj]^n) Determine tt;tm;mm:tt=^tT^t; tm=^tT^m mm=^mT^m(90) Table 2: Three Modes Formulation approach de nition zero (Eq. 40) and a last one to achieve consistent terminal separation conditions according to Eq. 59. This makes seven linear constraints, which leaves three free parameters. An SVD approach is used (Section 5.3) to manage these conditions and to look for the coecient vector with the smallest norm possible. Once is solved for, the element internal solution process continues as usual, solving the three-variable system and calculating the fracture traction vector. The complete fracture stress has to be calculated using a volume average (Eq. 21) as Gsbis no longer constant in space. Please refer to Table 3. 6.3. Enriched kinematics formulation The last formulation type uses full enriched fracture kinematics as described in Sections 4.3 and 5.2 through a generalized mode vector : =[juj]n0[juj]t0[juj]m0ntmntmt; counting three rigid body displacements, three rigid body rotations and three simple axial strains. Once again, two di erent models from this family are included in the present analysis: one considering a xed de nition for 'and the other a more complex, parameterised one. For a xed 'approach, the idea is to numerically de ne the G0 sboperators by calculating each column corresponding to each of the modes kusing kinematic mode linking relations. The process was already discussed in Section 4.3 for 2D as originally conceived in [11]. The same idea is brought now to the current three-dimensional problem. The speci c expressions for each of the columns for G0 sbare explicitly detailed inAppendix C ("right side expressions" ), as well as all the rationale behind the kinematic mode linking in a 3D context. There is no need to look after the original shape for the 'function. This formulation is summarized in Table 4. Once the G0 sboperator is populated, the approach follows the ideas proposed in Section 5.4 for building the three equation system depicted in 71 for the traction vector components Tn;Tn;Tn. Three additional equations for damaging the fracture stress components tt;tt;ttare built as well following the structure in 74. The last required three equations are taken from 85, which help to uncouple fracture translations completely. Again, it is remarked that these last three auxiliary equations are proposed to close the system completely and to further help giving a physical sense to the kinematic modes k, using some of them as 35 Three Modes Formulation, Controlled ' Kinematic modes : [juj]n;[juj]t;[juj]m Quadratic, unknown single ':'=PT 2 Linear G sb, function of :Gsb=NeX ipir'( ;x;y;z ) Solution for coecients :'(xi; ) =pi; Mi6=j( ) = 0 Mnn( ) =Mnn1 Traction-separation system & Solve for [juj]n;[juj]t;[juj]m:T=Te+M[juj]^n=q([juj]^n) Calculate Tn;Tt;Tm% Calculate:=1 VeZ C[Bd+GsbR[juj]^n]dV Determine tt;tm;mm:tt=^tT^t; tm=^tT^m mm=^mT^m(91) Table 3: Three Modes formulation, Controlled 'approach de nition bu er variables. At the end, a 9 9 nonlinear system is obtained, regardless of the 'approach. Once the system is solved, all fracture statics variables ( Tn;Tt;Tm;tt;tm;mm) will be available directly. The controlled 'version of this model deals with a considerably more complex 'de nition compared to that used on the three-mode approach: a full cubic, piece-wise, triple 'function capable of delivering up to 120 free parameters. The G0 sboperators will be dependent on the coecients following the discussion of Section 5.2, and mode linking relations will serve to build linear restrictions on . The speci c structure of the complete linear system on for this proposal is detailed in the Appendix C. It makes a total of 117 linear constraints. An SVD solution approach is once again used to avoid under/over-constraint problems. Once the system is solved, the static solution process continues exactly the same as mentioned before. The approach is stated in Table 5. 6.4. Results and discussion A series of plots is presented next to make a direct comparison of the results obtained from the elemental simulations for the ve formulations types. The results are divided in two sections: static analysis and kine- matic analysis. The static analysis pertains the results of all fracture traction vector components Tn;Tt;Tm and the fracture stress state components tt;tm;mm. The kinematic analysis concerns the behavior of the pertinent kinematic modes (one for the single formulation, three for the full translation and nine for the completely enriched formulations) and their comparison to the kinematic reference already de ned at the beginning of this section. Static results The load pro le is conceived with enough crack separation and sliding for placing the element at terminal separation conditions. The +and the regions thus act as completely independent bodies, and no forces shall be exerted from one body to the other at all. This means that all fracture traction vector components Tn;Tt;Tmand all the remaining stress state components associated to the fracture interface tt;tm;mm should be zero at this point. The expected behavior for a good E-FEM formulation would be then to drive 36 Enriched Modes Formulation, Fixed ' Kinematic modes : [juj]n0;[juj]t0;[juj]m0 n;t;m;n;t;m Fixed (unknown shape) ': Not used explicitly Fixed columns for G0 sb:G0 sbkk=Bdk Traction-separation system :Tn=Ten+X Mnnkk=qnn() Tt=Tet+X Mntkk=qnt() Tm=Tem+X Mnmkk=qnm() Complementary damage :tt=Ttt+X Mttkk=qtt() tm=Ttm+X Mtmkk=qtm() mm=Tmm+X Mmmkk=qmm() Closing constraints :X Alkk=Bl; l= 1;2;3 Solve for Calculate Tn;Tt;Tm% Calculate tt;tm;mm(92) Table 4: Enriched Modes formulation, Fixed 'approach de nition all these static components to zero at the end of the load pro le. This is, the following plots evidence the damaging quality of the formulations. Figure 14 show the behavior of the traction vector components Tn,Tt,Tmas a function of the load factor . All formulations manage to damage Tncompletely: this is the most basic functionality of any E-FEM-based model. Generally speaking, the decaying rate depends on the fracture sti ness component Mnn, with the only exception of the three-mode, xed 'formulation. The value of Mnnwill di er slightly for formulations where normal separation consistency has been ensured through 'coecients (Eq. 59). The three-mode, xed 'approach will have an e ective normal fracture sti ness that will be the resulting e ect of combining Mnnas well as coupled sti ness terms Mnt,Mnm, which will produce for this case a very low e ective normal sti ness, so that the traction Tnis drove immediately down after starting the load application. ForTtandTm, all formulations managing at least three kinematic modes are capable of successfully damaging these parallel traction components. These already start at zero since their localisation states ynt;ynmhave begun at zero attending to the Rankine criterion assumed before. The task for the formu- lations is just then to keep these traction components stabilised at zero. The plots only point out the fact that a single mode formulation will never be able to so. Despite this fact, single mode calculations remain conveniently simple and these stress anomalies can be manually suppressed to continue with global nonlinear calculations in a large scale model. This is the approach followed by [22, 33, 34, 24]. Even with its partial loss of physical meaning at local level, its numerical robustness and fracture representation capabilities at large scale models have proven to be successful. Figure 15 show the behavior of the stress state components tt,tmandmm. Fracture stress components are only accessible when projecting the complete state to the ^tand ^mdirections. Only the enriched modes 37 Enriched Modes Formulation, Controlled ' Kinematic modes : [juj]n0;[juj]t0;[juj]m0 n;t;m;n;t;m Cubic, Piece-wise, triple ':' j=PT 3  j; j =n;t;m Quad. G0 sb, function of :G0 sb( ) =HrJ('rJ)r'J Solution for coecients :G0 sbk( )k=Bdk '(xi; ) =pi '+( ) d='( ) d Mi6=j( ) = 0; i;j =n;t;m Mnn( ) =Mnn1;Mtt( ) =Mtt1 Mmm( ) =Mmm1 Traction-separation system :Tn=Ten+X Mnnkk=qnn() Tt=Tet+X Mntkk=qnt() Tm=Tem+X Mnmkk=qnm() Complementary damage :tt=Ttt+X Mttii=qtt() tm=Ttm+X Mtmii=qtm() mm=Tmm+X Mmmii=qmm() Closing constraints :X Alkk=Bl; l= 1;2;3 Solve for Calculate Tn;Tt;Tm% Calculate tt;tm;mm(93) Table 5: Enriched Modes formulation, Controlled 'approach de nition formulations are designed to explicitly track these components and to damage them properly. They do so at di erent rates, again, due to 'implications on fracture sti ness. All remaining formulations simply make these stress components to grow without control, rendering the need for manual suppression necessary for avoiding stress locking problems. Based on this static results analysis, only enriched mode formulations exhibit proper damaging properties for this realistic element. Kinematic results For the kinematic consistency analysis, we already have the kinematic reference on which the load design was based. The following plots will show how the formulations manage to follow this reference for each kinematic mode k, knowing before-hand that the global approach taken for all formulations possesses mathematical ambiguities that will always prevent a perfect match (Section 4.3). Figure 16 show the analysis on rigid body displacements. The normal crack separation [ juj]n0plot is the one still having all ve formulations with correct predictions at once. The formulations having a controlled ' 38 0.0 0.2 0.4 0.6 0.8 1.0 Load factor β024681012Tn(MPa) Normal Traction ( ˆn) Single Mode Three modes, Fixed ϕ Three modes, Controlled ϕ Enriched kinematics, Fixed ϕ Enriched kinematics, Controlled ϕ 0.0 0.2 0.4 0.6 0.8 1.0 Load factor β−17.5−15.0−12.5−10.0−7.5−5.0−2.50.0Tt(MPa) Parallel Traction ( ˆt) Single Mode Three modes, Fixed ϕ Three modes, Controlled ϕ Enriched kinematics, Fixed ϕ Enriched kinematics, Controlled ϕ 0.0 0.2 0.4 0.6 0.8 1.0 Load factor β−6−5−4−3−2−101Tm(MPa) Parallel Traction (ˆm) Single Mode Three modes, Fixed ϕ Three modes, Controlled ϕ Enriched kinematics, Fixed ϕ Enriched kinematics, Controlled ϕFigure 14: Comparison of fracture traction components response for all ve formulations. will be naturally closer to the kinematic reference as the fracture sti ness has been designed to be consistent with normal displacement at terminal separation conditions. The other two having slightly di erent Mnn miss the target by 30%, but still retain the same behavior tendency. Again, all formulations were prepared, one way or another, to handle normal separation kinematics. The three-mode, xed 'again, has a very di erent e ective normal sti ness, which causes even the rst load step solution for [ juj]n0to be excessively large. The lack of kinematic consistency in the models starts to be truly appreciated in the analysis of crack sliding modes [juj]t0;[juj]m0. In this case there's only one single approach that gets the kinematic reference right: the enhanced kinematics with controlled '. This is because this approach is the only one having a de nition for 'complex enough allowing to have free parameters to help imposing terminal separation conditions on the fracture sliding sti ness terms Mtt;Mmm. It's important to note that, while the three- mode, controlled 'attempts to do the same, its limited 'structure does not brake the linear dependency observed on the diagonal of the Mmatrix. In this case, only one direction can be made to comply with terminal separation conditions, and it was chosen to be ^n. The three-mode xed 'approach once again, has mixed sti ness interactions which drive the e ective sliding sti ness to unrealistic levels in both directions. Formulations not being able to control sliding fracture sti ness will even predict a completely opposite 39 0.0 0.2 0.4 0.6 0.8 1.0 Load factor β0.00.51.01.52.02.53.03.5σtt(MPa) Fracture stress component σtt Single Mode Three modes, Fixed ϕ Three modes, Controlled ϕ Enriched kinematics, Fixed ϕ Enriched kinematics, Controlled ϕ 0.0 0.2 0.4 0.6 0.8 1.0 Load factor β0.0000.0250.0500.0750.1000.1250.1500.1750.200σtm(MPa) Fracture stress component σtm Single Mode Three modes, Fixed ϕ Three modes, Controlled ϕ Enriched kinematics, Fixed ϕ Enriched kinematics, Controlled ϕ 0.0 0.2 0.4 0.6 0.8 1.0 Load factor β0.000.250.500.751.001.251.501.752.00σmm(MPa) Fracture stress component σmm Single Mode Three modes, Fixed ϕ Three modes, Controlled ϕ Enriched kinematics, Fixed ϕ Enriched kinematics, Controlled ϕFigure 15: Comparison of the response of fracture stress state components tt;tm;mmfor all ve formulations. sliding direction to that demanded by the load vector d. This totally breaks the physical sense of these kinematic modes. It's interesting to note that, if the zero response from the single mode formulation was included, it would still be missing the kinematic target by less than the other formulations. In some conditions, the lack of response is indeed a good response, after all. Figure 17 show the analysis on rigid body rotations. From now on, only the last two formulations involving enriched modes are able to contend. The kinematic