This paper analyses in detail the use of the Embedded Finite Element Method (E-FEM) to simulate local material heterogeneities. The work starts by a making a short review on the evolution of weak discontinuity models within the E-FEM framework to discuss how they account for the presence of multiple materials within a single element structure. A theoretical basis is introduced through some mathematical weak discontinuity de nitions and the Hu-Washizu variational principle, for then establishing a set of requirements for retaining variational and kinematic consistency for any weak discontinuity enhancement proposal. From a general de nition of a displacement enhancement eld, two particular enhancement functions are derived by considering di erent consistency requirements: one which has been typically used in previous works and other which truly possesses variational consistency. A discussion is held on enhancement stability properties and the impact to global nite element solution processes. In the end, numerical simulations are made to assess the performance of each of these enhancements on the task of modelling a classical bi-material layered 3D tension problem. The nal discussion evaluates both model performance and ease of implementation. Keywords: E-FEM modelling, embedded weak discontinuity, enhanced nite elements, incompatible modes, local heterogeneity modelling, nite element enhancement functions, mixed nite element 1. Introduction The numerical modeling of the heterogeneous nature of some materials is important for studying and predicting complex features of their physical behavior, including their mechanical response and resistance under speci c conditions. Typical numerical analysis techniques such as the nite element or nite di erence methods generally approach this problem considering homogeneous base domains. This requires a homogenisation procedure for a determined representative patch of the heterogeneous material that is big enough to exhibit the material's behavior as if it were a whole continuum but still small enough to be able to reveal its heterogeneous structure. Such is the basis of representative volume element (RVE) approaches [1, 2]. At some point, every approach following this line will require a realistic modeling of a limited domain on the scale in which the heterogeneities of a given material can be geometrically described in an accurate way. For such multi-scale simulation processes, a classical FEM approach will require an adapted mesh (or at least a speci c treatment to it) for the small scale to consistently capture the geometrical distribution of material heterogeneities in such domain. While sophisticated meshing adaptation techniques for hetero- geneous objects are still an active subject of study [3, 4], the approach remains computationally expensive and mathematically complicated, depending always on the arbitrary shapes of the di erent material phases present on the heterogeneous structure. This is specially true when a study requires the analysis of a large amount of heterogeneity distribution samples for a meaningful statistical treatment, such as in the execution of Monte Carlo methods that require repetitive samplings for the homogenisation process [5]. 1 Alternative approaches for the numerical modelling of heterogeneous material domains have emerged, such as Voronoi cell techniques [6], discrete elements for granular rocks [7] or reinforced concrete [8], or the advanced nite element methods [9, 10, 11]. Some of these approaches will focus more on the material interfaces, like the Voronoi cells that make use of mixed 1-D nite elements to represent the presence of di erent material domains and the strength of the mechanical connection between them. Typically, it is the nite element methods in two or three dimensions that will grant a more meaningful representation of the state of stresses in continuous material models, since they attempt a more direct and accurate geometrical description of material heterogeneities. In regards to the latter, the meshfree or non-adaptive mesh approaches have been introduced as an attractive solution to this modelling problem. To this e ect, applications of many advanced FEM techniques can be found, such as the Generalized Finite Element Method (G-FEM) [11], the Extended Finite Element Method (X-FEM) [12], the Base Force Element Method (B-FEM) [13] or the Embedded Finite Element Method (E-FEM) [14], the last one being precisely the scope of study in this work. It is recognized by the authors of this study that some approaches such as the X-FEM, in the state as they are at the moment of writing, have a broader capacity for the representation of complex multiphasic material distributions within a nite element domain than that of the E-FEM element technology. The interest in this work, however, will remain focused on the E-FEM framework. The authors consider worth the e ort to keep the advantages of this approach, such as the simplicity and hermeticity of the mathematical enhancements that leave the global nite element solution process practically untouched. A brief analysis on the evolution of the use of the weak discontinuity model on the E-FEM approach will be presented to the reader in Section 2. A detailed analysis of the theoretical foundations for the de nition of weak discontinuity enhancements for the modelling of heterogeneous materials will be introduced in Section 3. This will help to establish a set of basic consistency requirements for de ning these enhancements under the light of the Hu-Washizu variational framework. In Section 4, two particular weak discontinuity enhancement functions will be derived based on this consistency analysis and the amount of requirements chosen to be satis ed. The rst of them is the one typically managed in the reference E-FEM works, and a second one is proposed to purposefully maximize variational and kinematic consistency. Finally, basic comparative numerical simulations are made between the weak discontinuity enhancement proposals in Section 5. The analytical solution for a basic heterogeneous model will serve as a reference to assess the performance of the enhancements, as well as any further comments on particular behaviours. A concluding discussion will follow in Section 6, considering all theoretical and practical aspects of the developments presented in this work. 2. The role of weak discontinuity enhancements on the E-FEM framework The use of weak discontinuity enhancements started actually as one of the rst embedded nite element approaches for the modelling of shear instability bands, with some pioneering studies paving the way for consolidating the E-FEM approach as a whole [15, 16, 17]. The main idea was to model a shear band through the use of two parallel strain discontinuity lines that would cross a non-adapted mesh, typically having a uniform geometry. The elements having sub-domains enclosed by the shear band would have di erent constitutive properties to represent the local instability happening inside. It introduces a jump on the strain eld, which translates into a sudden change in slope for its corresponding displacement eld without breaking its continuity (this is thus the reason of calling it a weak discontinuity enhancement ). Figure 1 shows an example for a constant stress triangle (CST) element. The localisation band model represented a nite and continuous region of a fracture process. The weak discontinuity approach for this kind of material failure modelling got diversi ed afterwards with the introduction of regularization processes [18, 19] to avoid scale dependencies, especially concerning 2 the problem of setting an arbitrary shear band thickness. Eventually, attention was diverted towards the strong discontinuity enhancement equipped with a discrete post-localisation law as the method of choice for the modelling of internal element fracture on the E-FEM framework [20, 21, 22, 23]. The tendency of this evolution was practically to drive the shear band model thickness to zero while maintaining mathematical and physical coherence on the formulation. This translates the discontinuity to the displacement eld directly (thus the reason of naming it a strong discontinuity ). The strong discontinuity enhancement was indeed proven to be a more pragmatic and robust way to avoid mesh dependencies as possible, granting more objectivity to the approach. Nonetheless, the application of embedded weak discontinuity enhancements for the modelling of shear bands still gathers some interest in recent works, such as ductile material failure simulations under dynamic conditions [24, 25]. Figure 1: Basic schematic of weak discontinuity enhancements used to represent a shear band within a triangular 2D element having a local frame ^n,^t. The shear band possesses a set of di erent (damaged) mechanical properties Elwhile the rest of the element retains its original elastic behvior E. Note the introduction of a thickness dcharacterising the shear band. The authors of this study consider that the application of weak discontinuity enhancements to model material heterogeneities really started with the works on fracture simulations of cementitious materials on the mesoscale [26, 27, 28]. These developments began by modeling a single weak discontinuity on 1-D beam elements to represent the presence of two di erent linear elastic sti ness domains coexisting on the same element. The perspective was di erent to that of Voronoi cell constructions [6] in the sense that no regular inclusion recognition had to be made on a material matrix to assign one beam per interface. A totally random, unstructured 3-D mesh was built with beam elements and a heterogeneity distribution in space was just projected directly onto these. Some elements would fall entirely on the domain of one material phase or other, while others would be found in a region where there was an interface between materials. Is is those elements that were enriched with a weak discontinuity enhancement function. The model was also equipped with a strong discontinuity enhancement at the same location as the weak discontinuity to represent eventual failure and separation of the domains. In this sense, the work was also innovating from the perspective of integrating both discontinuity enhancements for entirely di erent roles. While this model allowed an explicit use of the weak discontinuity to nally model mesh-independent hetero- geneities, no objective state of stresses was described in the domains as no spatially accurate representation of the continuum is possible by only making use of 1-D beam elements. It was until the work of Roubin [14] that this application of the weak discontinuity model was devised for 3-D elements, inspired on the works in [29, 27]. The main idea was to establish a piece-wise displacement eld enrichment that, once being processed through the application of a symmetrical gradient operator rsym, it would comply with the Maxwell interface strain compatibility conditions [15]. The model counts with one, two or three internal variables characterising the strain jump between materials depending on the dimension of the problem. As his predecessors, Roubin also appended a strong discontinuity enhancement to integrate a fracture model, but only considering a single fracture kinematic mode: normal separation. This 3-D development was later taken as a base in [30, 31, 32] to perform simulations for heterogeneous rocks and cementitious materials in a similar fashion. A variety of fracture phenomena was explored, such as 3 plane sliding, crack reclosure and multi-scale analyses, among other developments. Further applications of these ideas can be found in the domain of poromechanics and electromechanics [33, 34]. The use of the weak discontinuity on this format acquired yet more relevance with the work of Stamati et al. [31], where image processing techniques and X-Ray tomography made possible to project realistic heterogeneity distributions coming from actual samples used for experimental campaigns, reaching a new level of predictability and model validation procedures. This application of the weak discontinuity model for 3-D geometries as seen in the work of Roubin [14] has been taken as the point of departure for the present study. In the next section, the theoretical basis behind it will be scrutinised in detail. 3. Theoretical foundations and consistency analysis The basic construction of a weak discontinuity for the modelling of material heterogeneities starts with the assumption that a heterogeneous displacement eld u, referred from now on to as the physical displacement, can be expressed as the composition of an average, homogenized base eld uand a eld enhancement eu carrying the mathematical weak discontinuity: u=u+eu (1) Figure 2 illustrates a typical dual material partition for a 4-node tetrahedral element in domains +,  with a boundary @ and having a plane  das an interface. The base work in [14, 35] considers linear elastic properties for each domain such as Young moduli E+,Eand Poisson ratios +,. A local coordinate system ( ^n,^t,^m) de nes the orientation of the material interface, having ^nas the unit vector normal to  d. Figure 2: Basic schematic of a weak discontinuity in 3-D for the modelling of material heterogeneities within a tetrahedral element The homogeneous base eld uis determined entirely by the displacement of the boundary nodes and the natural interpolation functions of the element. The de nition of the eld euis determined by internal variables keeping in mind that its corresponding strain function should introduce the strain jump associated with the change of material domains. The strain elds, as second order tensors, are obtained through a symmetric gradient operator rs() =1 2h r()T+r()i : "=rsu=rsu+rseu="+e" (2a) "+="+e"+;x2 +(2b) "="+e";x2 (2c) where a distinction has been done between the strain elds on the +and domains at each side of  d. Note that the base eld "remains invariant by the de nition of u. 4 3.1. Enhancement requirements and variational analysis To retain kinematic and variational consistency, the weak discontinuity model has to comply with certain requirements through both displacement and strain elds. The approach in this study will be to determine the possible function space for the enhanced displacement eld by introducing and applying these constraints, also noting the set of constraints e ectively considered in the work of [14] that shapes the most typical choice for it in those works. From now on, the analysis will take place on the local reference frame ( ^n,^t,^m) unless stated otherwise. Its coordinate variables will be denoted as ;; . The most basic constraint pertains the physical displacement eld u: it shall not lose continuity through the material interface. Given that the base eld uis already continuous by de nition, this implies that that enhanced displacement eualso has to be continuous: Enhancement requirement 1: displacement continuity eu+ d=eu d(3) The next involves an analysis of the strain eld and the de nition of the strain discontinuity jump. A strain discontinuity jump  e"is de ned the di erence of strain elds e"+ande"at the material interface d, resulting also in a second order tensor: e"=e"+ de" d=2 4e"nne"nte"nm sym e"tte"tm sym sym  e"mm3 5 (4) The components of the strain jump  e"are not obliged to respect full continuity as their parent dis- placement eld, but must still comply with Maxwell strain compatibility conditions [15] to be coherent with it. For this, the projection of  e"on the normal direction ^nwill be allowed to be free of constraints, while all other unrelated components of the tensor will be driven down to zero. Recall that in a local coordinate setting, the projections can be easily obtained by just extracting the line-column corresponding to a given direction within the tensor. Thus: e"^n= e"nne"nte"nmT6=0 e"tt= e"tm= e"mm= 0(5) This leaves the strain jump tensor  e"with only three active components  e"nn, e"nt, e"nm. These will be rede ned as [ "]n, ["]tand ["]m, respectively. These are indeed the internal variables that de ne the weak discontinuity model. With this, a new enhancement requirement is de ned: Enhancement requirement 2: strain jump e"=e"+ de" d=2 4["]n["]t["]m ["]t0 0 ["]m0 03 5 (6) The requirements to follow need considerations coming from the variational analysis. As in [36, 14, 30] and in most of other works on the E-FEM framework, the Hu-Washizu variational principle is chosen due 5 to its exibility to handle element eld enhancements through the independence of displacement, strain and stress elds: Z @utdVZ utfbdVZ @ uttdA= 0 (7a) Z et(@u")dV= 0 (7b) Z e"t(("))dV= 0 (7c) where the real elds have been denoted as ( u;";) and the eld variations (or virtual elds) as (u;";), having a boundary traction vector tand body forces fb. It's important to note that the real stress eld , in general, is di erent from the stress coming from constitutive law calculations ("). The same can be said from the real strain eld "and@u, where@is a eld gradient operator. A Voigt vector notation has been already implied for all elds involved in the presentation of Eqs. 7a-7c. All elds are independent from each other, in the sense that they do not have to necessarily follow direct gradient relations such is in Eq. 2a. This is, the Hu-Washizu variational principle allows for exible eld discretization strategies. However, the elds should retain enough physical meaningfulness to be able to correctly model the phenomenon in question. Authors working on this framework generally choose a discretisation strategy as to render the model as manageable and ecient as possible sacri cing the minimal amount of mechanical representation quality. This choice also considers the ease of an integration process with other models that might have a similar eld discretisation approach (such as a strong discontinuity model). The displacement and displacement variation elds u,uare commonly discretized taking only the standard displacement eld u: u=u=Nd (8a) u=Nd (8b) withNas a standard interpolation matrix and dthe standard nodal displacement vector. dis the cor- responding variation. This strategy means that only the eld uis used for describing node positions and imposing boundary conditions. In such case, it should be clear that, in order to ensure that dretains the correct nodal information, the eld ushould have the same value as uat the boundaries @ of all the element (i.e. the nodes on it): uj@ =uj@ (9) Given that we already have a de nition as stated in Eq. 1, this implies: uj@ =uj@ +euj@ )euj@ = 0 (10) This analysis de nes the next constraint for the weak discontinuity model: Enhancement requirement 3: null value on borders eujx=xi= 0; i = 1;2;:::;N (11) wherexiare nodal positions and Nis the number of nodes of the element. This requirement stands as the most overlooked in the current literature of this family of formulations applied to the modelling of material 6 heterogeneities. It is also the one that will make a signi cant di erence in the enhancement function shape with respect to the one used in typical heterogeneous E-FEM studies. Going forward with the discretisation strategy, the domain dependent strain eld "and its variation "conserve all kinematics description terms as stated in Eqs. 2a-2b. Their enhanced sections ( e"ande", respectively), which depend on the internal variables [ "]n, ["]t, ["]m, are stated through the de nition of aweak discontinuity vector [j"j] =["]n["]t["]mTand its variation [j"j]. This discretisation strat- egy also allows the opportunity to use di erent interpolation matrices G w;G wfor the real and variation enhancements, respectively: "=( Bd+G+ w[j"j]x2 + Bd+G w[j"j]x2 (12a) "=( Bd+G+ w[j"j]x2 + Bd+G w[j"j]x2 (12b) Note that, until now, no speci c form for G w;G whas been still assigned. In the original work of Roubin [14], it is actually assumed that G w=G w. The reason for this choice will be explained shortly. The stress eld and its variation are just interpolated using single independent stress vectors sand sthrough the use of interpolation matrices SandS, respectively: =Ss (13a) =Ss; (13b) The de nition for the stress eld coming from the constitutive law (") is based on the assumption that each of the material domains possesses its own linear elastic constitutive law considering separate second order linear elastic constitutive tensors C+andC. These linear operators act upon di erent regions of therealstrain eld". The de nition for this stress eld is thus devised as: (") =( C+(Bd+G+ w[j"j])x2 + C(Bd+G w[j"j])x2 (14) It should be emphasized that while the constitutive stress (") is by default a domain-dependent de nition, therealstress eldis not necessarily obliged to follow the same characteristics. Having set a discretisation strategy, the variational analysis takes place by using all previous equations to develop the Hu-Washizus system through Eqs. 7a-7c. For instance, after working Eq. 7b and considering that independence of [ j"j] as well as the fact that the variation vector sis to remain arbitrary, it follows that: stZ eStG wdV[j"j] =0)Z eStG wdV=0 (15) Eq. 15 demands orthogonality between the virtual stress interpolation matrix Sand the enhancement matrix operators G won each respective subvolume V+;V. Since the matrix Sdoes not represent any kind constraint on any real eld in the framework, most authors just take this relation as granted assuming that for any convenient choice of shape for the realG wmatrix operators, there will always exist a de nition forScomplex enough to satisfy Eq. 15. Indeed this is case for usual enhancement function proposals, but the demonstration will not be done on the present work. It will be only stressed that, since Swill be de ned following exclusively the structure of a chosen G w, it is clear that virtual stresses and real stresses will be interpolated in an asymmetrical fashion. Therefore, the framework already departs from this point from a classical Bubnov-Galerkin variational approach to a Petrov-Galerkin one. 7 Next, the works with Eq. 7c lead to two independent equations thanks again to the arbitrary nature of dand[j"j]: dtZ eBt((")Ss)dV+[j"j]tZ eGt w((")Ss)dV= 0 (16a) !dtZ eBt((")Ss)dV= 0!Z Bt(")dV=Z BtSsdV (16b) ![j"j]tZ eGt w((")Ss)dV= 0!Z Gt w(")dV=Z Gt wSsdV (16c) Eq. 16b is the only relation in the framework that establishes an explicit connection between real and constitutive stresses. It stands as a weak equality involving a known standard interpolation gradient matrix B. In the case that the choices in the weak discontinuity model render the shapes ( e.g., polynomial order) of(") andto be the same, then a strong equality between the elds will hold. This also implies that the real stresseswill have to be domain dependent as ("). For other cases, notably for an assumed constant real stress eld, a volume-averaged expression can be used for its calculation: s=1 VeZ (")dV; (17) where the shape, order or domain distribution of (") is irrelevant. For working Eq. 16c, authors commonly assume that both left and right sides can be set to zero, creating orthogonality constraints. On the right side, given that the realstress interpolation matrix Sis xed based onrealnode data, the orthogonality de nes the structure of Gt w. In general we will have: Z eGT wSsdV=Z +G+T wSdVs++Z GT wSdVs= 0 (18) Depending on the chosen form for, the de nition for the GT woperators will change attending to Eq. 18. Additionally, there is the requirement of passing a nite element patch test [29, 37]. The patch test states that it should be possible to satisfy all variational equations for a constant real stress eld regardless of the form chosen for this eld in the rst place. In other words, Eq. 18 must be satis ed for both a case in which the interpolation matrix Sis the identity matrix and the case in which Scorresponds to the shape functions truly desired to describe the stress eld. With this in mind, it will be required to have s+=s=sfor validating the case of some uniform and arbitrary-valued sover all the element domain e. This yields another requirement for the weak discontinuity model enhancement, even if this part of the discussion focuses on its virtual de nition: Enhancement requirement 4: Patch Test Z +G+T wdV+Z GT wdV=0 (19) For instance, the works in [14, 35] take the simplest approach for a tetrahedron element, which is to assume a constant real stress eld. In this case, the patch test requirement actually obliges the entire real stress eld to be constant and unique through all the element. Thus, it will have to be always calculated using a volume average (Eq. 17). On the other hand, this lowers the burden of also having to satisfy Eq. 18 for a higher order S. This leaves Eq. 19 as the only requirement for G w, which allows the possibility to make constant matrix de nitions, yielding: V+G+ w+VG w=0 (20) 8 Setting the left side of Eq. 16c to zero brings the nal relation for calculating the weak discontinuity internal variables [ "]n, ["]t, ["]mas a function of a displacement input d. This is, if a de nite shape for all weak discontinuity enhancement matrix operators G wandG whave been established at this point. Eq. 14 can be used to develop: Z +G+T wC+ Bd+G+ w[j"j] dV+Z GT wC Bd+G w[j"j] dV=0 (21a) Kwbd+Kww[j"j] =0 (21b) Kwb=Z +G+T wC+BdV+Z GT wCBdV (21c) Kww=Z +G+T wC+G+ wdV+Z GT wCG wdV (21d) [j"j] =K1 wwKwbd; (21e) where speci c enhancement sti ness matrices Kwb;Kwwhave been de ned. The direct equality between both sides of equation Eq. 16c brings additional information on some internal constraints for the stress elds that will emerge implicitly without impacting further enhancement de nitions on the weak discontinuity model overall. The implications are interesting to know but are not operationally useful. For instance, if constant matrix de nitions are done for the virtual G woperators, one can use Eq. 20 to express everything a function of the G+ woperator, isolate and nally reach: 1 V+Z ++(")dV=1 VZ (")dV (22) If enhancements are de ned in such a way that the constitutive stress elds are constant on each domain, then Eq. 22 practically dictates a strong equality and uniformity of the constitutive stress elds in all the element. At the same time, the strain elds are assured to be di erent on each domain by de nition. While this emerging equality might appear counter-intuitive to the reader, this would resemble a very simplistic situation such as in a couple of springs connected in series, where the strain is di erent but the force through them is exactly the same. Indeed, Eq. 22 has been observed to hold for a large number of elements during actual numerical simulations based on the formulation of Roubin [14] for large models managed in the works of Stamati et al. [31]. 4. Weak discontinuity enhancement proposals Now that all relevant constraints for de ning weak enhancement functions have been introduced, a particularisation of the model will take place, deriving two di erent enhancement eld functions considering slightly di erent ways of achieving the satisfaction of consistency requirements. For the sake of simplicity and coherence with the background literature, a linear tetrahedron will be set as the base element from now on. For now, it will be assumed that a constant stress eld is sought. 4.1. Typical enhancement analysis Authors managing the modelling approach in [14, 35] decide to make the weak discontinuity enhancement completely symmetrical by letting G w=G w. This automatically renders the weak discontinuity model variationally symmetric at the expense of removing the exibility of having a virtual enhancement with di erent characteristics. At the same time, only requirements 1, 2 and 4 (Eqs. 3, 6, 19) are explicitly imposed to this unique enhancement function. Instead of requirement 3 (Eq. 11), a general zero reference for the enhancement is set at the interface plane  d. This last imposition is absolutely unrelated to any considerations on variational consistency. This line of approach also chooses the simplest de nition for the model: a linear eld euand therefore constant operators G w. 9 Considering these restrictions, it will be demonstrated that the possible function space for eureduces to a unique expression. Let the following linear de nitions for the piece-wise enhanced displacement eld be: eu+=a++b++c++d+ (23a) eu=a+b+c+d (23b) where each vector has components contributing to each local direction ( ^n;^t;^m),e.g.,eu+= eu+ neu+ teu+ mT. The goal is to particularize the vectors a,b,c,das a function of basic element data and the weak discontinuity variables [ j"j]n, [j"j]tand [j"j]m. In local coordinates, it is not hard to see that the interface plane dis simply described by the equation = 0. This eases the application of requirement 1: eu+ =0=eu =0(24a) a++c++d+=a+c+d (24b) )a+=a=a;c+=c=c;d+=d=d (24c) where thewill be omitted from now on in all coecients a;c;d. The strain jump requirement (Eq. 6) requires the calculation of the symmetric gradients on each part of the enhanced displacement eld: e"+e"=2 4["]n["]t["]m ["]t0 0 ["]m0 03 5 (25a) e"=2 664u n;1 2 u n;+u t; 1 2 u n;+u m; sym u t;1 2 u t;+u m; sym sym u m;3 775=2 4b n1 2 cn+b t1 2(dn+b m) symct1 2(dt+cm) sym sym dm3 5 (25b) e"+e"=2 4b+ nb n1 2 b+ tb t1 2(b+ mb m) 1 2 b+ tb t 0 0 1 2(b+ mb m) 0 03 5 (25c) )[j"j] =2 4["]n ["]t ["]m3 5=2 4b+ nb n 1 2 b+ tb t 1 2(b+ mb m)3 5 (25d) Here, the results of Eq. 24c have been used. Note that for this reason, the zeros required in the strain jump matrix are produced naturally. The application of the patch test (Eq. 19) becomes easier if both sides of Eq. 20 are multiplied (con- tracted) by the weak discontinuity variable vector [ j"j] to recover enhanced strain elds, but in vector format: V+G+ w[j"j]|{z} e"++VG w[j"j]|{z} e"=0[j"j] =0 (26a) V+2 6666664b+ n ct dm 1 2 cn+b+ t 1 2(dt+cm) 1 2(dn+b+ m)3 7777775=V2 6666664b n ct dm 1 2 cn+b t 1 2(dt+cm) 1 2(dn+b m)3 7777775(26b) )ct=dm= 0 (26c) 10 Afterwards, the original idea for this formulation coming from Markovic [29] involves the imposition of a requirement that is not mandatory for variational consistency. The enhanced displacement euis prescribed with a value of zero through all the interface plane ( = 0). While Markovic did never reveal any particular reasons for this decision in his research, it is later found in this work through further mathematical analysis (Sections 4.3, 4.4, 4.5) that this constraint provides some operational bene ts to the framework. Indeed, the imposition of a zero enhanced displacement reference at  dis an aggressive constraint that will simplify the enhancement function to a great extent: u+ =0=u =0= 0 (27a) 2 4an at am3 5+2 4cn 0 cm3 5+2 4dn dt 03 5=2 40 0 03 5 (27b) )2 4an at am3 5=2 40 0 03 5;2 4cn 0 cm3 5=2 40 0 03 5;2 4dn dt 03 5=2 40 0 03 5 (27c) It should not be forgotten that this condition is not a replacement of requirement 3 (Eq. 11). With these results and the expressions coming from the application of the patch test requirement (Eq. 26b), one can calculate the bcoecients directly as a function of weak discontinuity internal variables: b+ n=V V["]n; b+ t= 2V V["]t; b+ m= 2V V["]m; (28a) b n=V+ V["]n; b t=2V+ V["]t; b m=2V+ V["]m; (28b) In the end, the enhancement function reduces to just a set of bcoecients multiplying the coordinate , which is the normal distance from the interface plane  d: 2 4eu n eu t eu m3 5=2 4b n b t b m3 5=V V2 4["]n 2 ["]t 2 ["]m3 5 (29a) )eu=V V ["]n^n+ 2 ["]t^t+ 2 ["]m^m (29b) If this expression is reverted to global coordinates, the original weak discontinuity enhancement version as used in Roubin [14] is recovered: Weak discontinuity enhancement: typical de nition eu= ^n(xxd) ["]n^n+ 2 ["]t^t+ 2 ["]m^m (30)  =( +=V Vx2 + =V+ Vx2 (31) Here the coordinate has been expressed as a projection of a distance with respect to to the  dplane and a domain-dependent scalar  containing the volume ratios has been de ned. It's important to note that the intent of this work to derive this already-known enhancement eld shape is to demonstrate that, instead of starting with a seemingly arbitrary de nition [29, 14], its nal form comes rather from the application of a de nite set of constraints. This sheds light on the theoretical basis on which this family of enhancements is built upon. 11 From Eq. 30, general expressions for the G woperators can be found on the global reference frame by applying the symmetrical gradient operator in global coordinates. Einstein index notation is useful to reach the following typical strain eld tensor expression in terms of symmetric dyadic products: e"= h ["]n(n ^n)s+ 2 ["]m(n ^m)s+ 2 ["]t n ^tsi ; (32) where factors of 2 have been accounted for in the shear strain-related terms to ease the use of 3D constitutive laws. This expression is converted into a Voigt format to nally obtain: e"=G w[j"j] =  Hw[j"j];Hw=2 6666664n2 x nxmx nxtx n2 y nymy nyty n2 z nzmz nztz nxny+nynxnxmy+nymxnxty+nytx nzny+nynznzmy+nymznzty+nytz nznx+nxnznzmx+nxmznztx+nxtz3 7777775(33) Note that the domain-dependent term in all these de nitions stands as a single scalar  taking the form of domain volume ratios. 4.2. Consistent enhancement analysis For the case of a more variationally consistent enhancement eld considering requirement 3 (Eq. 11), it is more practical to start the analysis by expressing the enhancement as a piece-wise de nition of two linear elds using classical linear interpolation functions: eu+=eu+ 1 1+eu+ 2 2+eu+ 3 3+eu+ 4 4 (34a) eu=eu 1 1+eu 2 2+eu 3 3+eu 4 4 (34b) where the interpolation functions associated each node iof a base linear tetrahedron have been de ned as i=ai+bi+ci+di, where all coecients ai;bi;ci;diare known. This time, the goal of the model particularisation is to nd the value of all the nodal enhanced displacements eu ias a function of nodal coordinate information and the weak discontinuity internal variables. An auxiliary variable piis de ned as a position indicator between the +and the domains as follow: pi=( 1xi2 + 0xi2 i=f1;2;3;4g (35) Using this variable, a general mixed velocity variable euiis de ned: eui= (1pi)eu+ i+pieu i; i =f1;2;3;4g (36) Having all that, we can start by the application of requirement 3 (Eq. 11) in a very straight-forward fashion by just nullifying some nodal enhanced displacements in their corresponding domains. This leads to: eu+=NX i(1pi)eu+ i i (37a) eu=NX ipieu i i (37b) This also implies that the mixed variable euiactually captures the set of all non-zeroeu ivariables to be solved for in this process. 12 Next, displacement continuity (requirement 1) is applied: eu+ =0=eu =0(38a) NX i(1pi)eu+ i(ai+ci+di) =NX ipieu i(ai+ci+di) (38b) where the mixed variable euican be used to reach the following: NX i(12pi)aieui= 0NX i(12pi)cieui= 0;NX i(12pi)dieui= 0 (39) The strain jump requirement 2 follows exactly the same process followed previously in Eqs. 25a-25d to nd expressions relating enhanced nodal displacements to weak discontinuity internal variables: NX ibieuin= ["]nNX ibieuit= 2 ["]t;NX ibieuim= 2 ["]m(40) whereeuin;euit;euimare the three components of each non-zero euion the local directions associated to each nodei. At this point, it's worth stopping to make a variable summary on the linear system being currently built for all non-zero enhanced nodal displacements eui. All non-zero variables corresponding to each direction ^n;^t;^mcan be grouped in single vectors eun,eut,eumas follows: eun=2 664eu1n eu2n eu3n eu4n3 775;eut=2 664eu1t eu2t eu3t eu4t3 775;eum=2 664eu1m eu2m eu3m eu4m3 775(41) The system can then be summarized using block matrix de nitions: 2 4Ce0404 04Ce04 0404Ce3 52 4eun eut eum3 5=2 4["]n;e ["]m;e ["]t;e3 5; (42) where: Ce=2 664(12p1)a1(12p2)a2(12p3)a3(12p4)a4 (12p1)b1(12p2)b2(12p3)b3(12p4)b4 (12p1)c1(12p2)c2(12p3)c3(12p4)c4 (12p1)d1(12p2)d2(12p3)d3(12p4)d43 775;04=2 6640 0 0 0 0 0 0 0 0 0 0 0 0 0 0 03 775(43a) ["]n;e=2 6640 ["]n 0 03 775;["]t;e= 22 6640 ["]t 0 03 775;["]m;e= 22 6640 ["]m 0 03 775(43b) The system already counts with 12 enhanced nodal displacement variables and 12 equations, which render it closed with a unique solution if the coecient matrix in Eq. 42 is not singular. Further application of requirement 4 (the patch test) does not add any new variables to the system. Therefore, it can be stated that for a linear de nition of the enhanced eld eu, it is not possible to make a variationally symmetric de nition for the Gwmatrix operators having a unique base enhanced displacement eld. To be fully consistent while keeping linear de nitions, the framework requires to make Gw6=G w, taking requirement 4 (Eq. 19) as 13 the guidance to de ne G wand the other requirements to de ne Gwin a separate way. The typical weak discontinuity enhanced eld (Eq. 30.) thus cannot be, by de nition, variationally consistent. The analysis on this section will continue to particularize the enhanced eld with the system proposed in Eq. 42. As the system is block-diagonal, a compact-closed solution is found: eun=2 664C1 1;2 C1 2;2 C1 3;2 C1 4;23 775["]n;eut= 22 664C1 1;2 C1 2;2 C1 3;2 C1 4;23 775["]t;eum= 22 664C1 1;2 C1 2;2 C1 3;2 C1 4;23 775["]m; (44) where theC1 i;2coecients come from the second column of the inverse of the Cematrix. The particularized enhanced eld can then be expressed as: eu+=NeX i(1pi)C1 i;2 i2 4["]n 2 ["]t 2 ["]m3 5;eu=NeX ipiC1 i;2 i2 4["]n 2 ["]t 2 ["]m3 5 (45) Finally, the consistent weak enhancement eld can still be written in the familiar format: Weak discontinuity enhancement: consistent de nition eu=  ["]n^n+ 2 ["]t^m+ 2 ["]m^t (46)  =( +=PNe i(1pi)C1 i;2 ix2 + =PNe ipiC1 i;2 ix2 (47) It's important to note that, while  stays as a constant in the typical enhancement model, it becomes a variable parameter on the consistent model depending on nodal coordinates embedded in the interpolation functions i. TheG woperators can be devised again by making use of the symmetric gradient operator. Taking, for instance, the enhanced strain eld on the +domain one can reach an analogous tensor expression to that of the typical model: e"+=NeX i(1pi)C1 i;2h ["]n(ei ^n)s+ 2 ["]t(ei ^m)s+ 2 ["]t ei ^tsi (48a) ei=bicidiT(48b) Here, the vector eicoming from local interpolation function coecients, has to be transformed (rotated) to global coordinates as needed. Again, expressing in a Voigt vector eld format: e"+=G+ w[j"j] =NeX i(1pi)C1 i;2HW;i[j"j] (49a) HW;i=2 6666664binxbimxbitx cinycimycity dinzdimzditz biny+cinxbimy+cimxbity+citx diny+cinzdimy+cimzdity+citz dinx+binzdimx+bimzditx+bitz3 7777775(49b) 14 where the coecients bi;ci;diare already taken from a rotated eivector in this global de nition. G wfollows in a similar way: e"=G w[j"j] =piC1 i;2HW;i[j"j] (50) The virtual operators G wonly have the goal of complying with the patch test (requirement 4). If the simplest, constant de nition for them is adopted, satisfaction of Eq. 20 allows an in nite amount of choices. G wcan even take the form of the typical enhancement (Eq. 33), which has been already built to satisfy Eq. 20. 4.3. A discussion on enhancement stability properties The stability of the weak discontinuity model can be assessed by observing that in Eq. 21e, the cal- culation of the weak discontinuity variables [ "]t;["]t;["]mas a function of nodal standard displacements d depends on the inverse of a Kwwsti ness matrix. The behavior of Kww, depending on the form of the enhancement operators, will determine if the formulation becomes unstable under certain conditions. No other mathematical stability sources are identi ed. For a more direct analysis, it is convenient to work in the local frame ( ^n;^t;^m). It will be assumed that the consistent enhancement will use the same operators G was the typical enhancement, while retaining the real operators G was devised in Section 4.2. If this is the case, both formulations will share the same Hwoperator for virtual elds, which in local coordinates reduces to: HT w=2 41 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 13 5 (51) For a typical weak discontinuity enhancement, Eq. 21d can be developed to get to the following: Kww=V+V V2HT w C+V+CV+ Hw; (52) where linear elastic constitutive matrices in three dimensions can be assumed for the corresponding materials on +and as: C+=2 6666664c+ 1c+ 2c+ 20 0 0 c+ 2c+ 1c+ 20 0 0 c+ 2c+ 2c+ 10 0 0 0 0 0 c+ s0 0 0 0 0 0 c+ s0 0 0 0 0 0 c+ s3 7777775;C=2 6666664c 1c 2c 20 0 0 c 2c 1c 20 0 0 c 2c 2c 10 0 0 0 0 0 c s0 0 0 0 0 0 c s0 0 0 0 0 0 c s3 7777775(53) for some real, positive constants c 1;c 2;c s. Taking these de nitions, the Kwwmatrix associated to the typical enhancement can be calculated in a diagonal, compact expression: Kww=2 4V+c+ 1+Vc 1 0 0 0V+c+ s+Vc s 0 0 0 V+c+ s+Vc s3 5 (54) As the constants c 1;c sstay real and positive as well as the subvolumes V+;V, Eq. 54 reveals that the typical weak discontinuity enhancement turns out to be unconditionally stable , no matter what the orientation of the interface, the subvolume partition and the material elasticity characteristics are. This is one of the most prominent reasons why, despite not being variationally consistent by de nition, authors in [14, 30, 31] have been keen to keep it on their heterogeneity modelling approaches, as it eases the implementation process and the numerical solution control. 15 Working with the consistent enhancement, Eq. 21d returns the following expression: Kww=V+V VHT wNeX iC1 i;2 (1pi)C+piC Hw;i; (55) where the node-dependant matrix Hwican be reduced to: HT wi=2 4bi0 0ci0di 0ci0bidi0 0 0di0cibi3 5 (56) Taking the same constitutive matrix de nitions in Eq. 53, the Kwwfor this case can be devised as: Kww=V+V VNeX iC1 i;22 4bik1icik2idik2i ciksibiksi 0 diksi 0biksi3 5 (57) where thekjiparameters are de ned as follows: k1i= (12pi)c(12pi) 1 (58a) k1i= (12pi)c(12pi) 2 (58b) k1i= (12pi)c(12pi) s (58c) The structure of Eq. 57 is considerably more complex than that of Eq. 54, where in the former we can appre- ciate a more extensive participation of elemental parameters such as interpolation function coecients and functions that are domain dependent (interface location/orientation dependent) within nodal summations. No unconditional stability can be readily assured in Eq. 57. While the speci c subspace of parameters that drive Kwwunstable for the consistent enhancement will not be calculated in a rigorous fashion on this work, it is not hard to see that there's a condition that will intuitively introduce mathematical ambiguity problems: when the interface plane crosses exactly or very close to one or more element nodes, where it makes a sudden change of behaviour. If a given nite element model manages a non-structured mesh with random orientations while also having random material interfaces, instability or near-instability conditions will be certain to happen in some elements with a suciently large mesh. Implementation e orts have to consider this fact. 4.4. Sti ness matrices and impact the to global solution process One of the most attractive features of the E-FEM framework is the ability to limit the work with eld enhancements and all their associated variables within internal element calculation routines. All sti ness calculations derived from the special internal calculations associated to the discontinuities can be condensed and integrated to the standard elemental sti ness matrix. This way, no formal degrees of freedom are added to the nite element global solution process, so that the methods, routines and the solution platform shall remain untouched. However, the classical sti ness matrix properties normally identi ed in standard nite elements may change depending on the structure of the embedded enhancement. For the weak discontinuity enhancements presented in this work, the numerical solution process and the construction of an equivalent sti ness matrix can be started by taking Eq. 7a (elemental force balance) and Eq. 21b (main relation between nodal displacements and weak discontinuity variables) to build the following linear system: Kbbd+Kbw[j"j] =fe ext (59a) Kwbd+Kww[j"j] =0 (59b) 16 Here, Eq. 7a has been integrated using the previously de ned operators depending on the enhancement version to de ne sti ness matrices KbbandKbw. Condensation of the system then takes place by reducing [ j"j] from Eq. 59b and substituting on Eq. 59a, giving rise to a de nition of an equivalent elemental sti ness matrix Kscmultiplying the standard normal displacements vector d: Kscd=fe ext (60a) Ksc=KbbKbwK1 wwKwb (60b) With this, a global sti ness matrix assembly process may then be performed by taking the corresponding matrices Kscassociated to each element on a given model. For the case of a typical weak discontinuity enhancement, it is not hard to see that: Kbw=V+V VBT C+C Hw=V+V VHT w C+C BT =KT wb (61) so that: KbwK1 wwKwbT=KT wbKT wwKT bw=KbwK1 wwKwb (62) Given that Kbbis already symmetric, it can be concluded that the condensation process will always return asymmetric K sc. On the other hand, with the consistent enhancement, this is not the case: ( BTNeX iC1 i;2 (1pi)V++piV CHw;i)T 6=V+V VHT w C+C B (63a) Kbw6=KT wb (63b) This is, the consistent enhancement will, in general , return an asymmetrical sti ness matrix. This will introduce the need to use asymmetric solvers during a global numerical solution, with all computational and implementation implications that come along. This could be already perceived from the variational formulation itself, where the typical enhancement always managed the same operators for both virtual and real elds while the consistent enhancement didn't. 4.5. Variational consistency errors in the typical enhancement Now that it is known that the typical enhancement in general will not comply with basic requirement 3 (Eq. 11), it is relevant to discuss the conditions under which the formulation might produce large variational errors and the ones in which these errors will be kept within a reasonable range. For a case of a 1-D element, like in the applications done by Benkemoun [38] or Melnyk [26], there can be only one node on each side of the material interface, and its orientation will always be normal to the line de ning the body of the element. Under these conditions, it can be shown that the only di erence between a typical and a consistent enhancement is only a constant o set  euon the displacement function (refer to Fig. 3). Given that the slopes coincide and the operators of the formulation are based on eld derivatives, it can be concluded that the typical formulation e ectively complies with requirement 3 and thus also remains fully variationally consistent. For 2D and 3D elements in general, this is not the case. The slope of the typical enhancement will always be aligned to the orientation of the interface  d. If more than one node is present on one of the domains, the enhancement will not be able to return the same eld value on all nodes simultaneously, no matter what o set is given to the eld. The only condition in which this might happen is when the nodes within 17 (a) Original enhanced eld (b) Enhanced eld with an o set  eu Figure 3: Graphical analysis on 1D typical weak discontinuity enhancements. Nulli cation of the eld at element nodes can be achieved by adding an o set  eu a domain are all located on the same coordinate, which would mean having element geometry aligned to the material interface. As the respective nodal coordinates start to divert, the typical eld will miss to opportunity to nullify the values at the boundaries. Figure 4 illustrates this for a 2D constant stress triangle (CST). The o set  eucan be arranged to minimize this variational error by making the eld to roughly pass through zero at an averageposition of all nodes on a given domain. (a) Typical enhanced eld over a CST (b) Consistent enhanced eld over a CST Figure 4: Graphical analysis on 2D typical weak discontinuity enhancements taking a constant stress triangle (CST) as the base element, compared to a consistent enhancement. For the typical enhancement, it is assumed that an o set  eunhas been given to the eld to minimize the error when trying to nullify eld values on the CST nodes. Based on this rationale, the typical enhancement can be perceived as an average estimation of a fully consistent enhancement that will be closer to it under certain mesh geometry conditions. It can be expected that for an unstructured mesh with good aspect ratios, this estimation will grant reasonable results compared to a fully consistent approach. For a heavily distorted mesh having a very large disparity on coordinates within a domain, the variational errors induced will certainly get larger. 5. Numerical simulations In this section, a numerical comparison takes place between both weak discontinuity enhancements to discuss their performance with respect to an analytical solution to a given problem. 5.1. Numerical model description Numerical simulations have been done on a simple cube model made up of two material layers separated by a planar interface. The interface is parallel to two of the cube faces. If interface concentrations are 18 neglected, the total axial reaction associated to a normal displacement on one of the faces can be calculated by means of the classical theory of mechanics of materials, representing the system as two springs in series accounting for the axial sti ness of each layer. Linear elastic behavior is assumed for both materials, characterized by Young Moduli E+;Eand Poisson ratios +;. Figure 5 illustrates this simple mechanical system. Figure 5: Basic description of a bi-material layered prismatic body with total height Land cross section A, treated as a two-spring mechanical system in series having constants k+andk. The sti ness partition depends of the position hof the interface plane. The idea of the present study is to compare how each of the weak discontinuity enhancements can model this ideal bi-material layout by comparing to the classical analytical solution. The interface plane location hwill be varied taking regular steps from having a zero position at one of the cube faces until reaching the opposite side of the cube. This will represent situations in which the cube starts completely homogeneous with one of the material phases and gradually becomes entirely lled with the other material phase. The cube will feature an unstructured mesh, totally independent from the planar interface. The interface will cross a certain amount of elements on random edges and positions, and these elements will be enhanced with one of the weak discontinuity eld functions studied in previous sections. Special care has been taken with the mesh density: the size of the elements should be small enough to generate enough enhanced elements, but these special elements should cover a signi cant amount of cube volume in order for them to really have a meaningful contribution to the global response of the model. Otherwise, if the elements are too small, we might e ectively get a large number of enhanced elements near the interface, but also a much larger amount of normal elements having completely one material phase or the other, and thus the global sti ness response of the numerical sample will be dominated by the standard nite elements instead, which are not the object of this study. The mesh nally selected for this study is shown in Figure 6, also highlighting the number of enhanced elements resulting from having the planar interface at 30% height from the designed bottom position. It's pertinent to mention that the quality of the mesh has been kept rather high, with no aspect ratios going beyond 3. The cube has dimensions of 10 mm 10 mm10 mm. The material properties chosen for these simulations have been those normally associated to a simpli ed concrete mixture: a material phase of mortar (Young modulusE= 14000 MPa, Poisson ratio = 0:2) and an aggregate material ( E+= 70000 MPa, += 0:2). No other material properties are needed since all simulations will be made in static conditions. Boundary conditions have been set as to retain an ideal axial prism problem as much as possible without signi cant e ects of near-interface eld concentrations. The load has been imposed as a uniform displacement of 0:015 mm on the free upper face. Figure 7 illustrates the details of the model. FEAP (Finite Element Analysis Program) (REF) has been used as the nite element numerical solution platform to implement both enhancements described in previous sections and for simulating the problem mentioned beforehand. 21 static-implicit simulations have been performed considering each weak disconti- nuity enhancement approach having 21 uniformly separated positions for the interface plane, going from the lowerzface of the cube ( h= 0) all the way to the opposite face ( h= 10 mm) taking steps of 0.5 mm. The 19 (a) General view of entire model mesh. Interface in white outline. (b) Mesh view isolating all enriched elements. Interface in white outline. Figure 6: General mesh description and a view of all enriched elements crossed by the material interface plane. 92 elements and 43 nodes are used for the entire model. Enriched elements make up approximately for 50% of total volume. solution is strictly linear elastic, where the only solver-speci c di erence between each enhancement case has been the use of symmetric and asymmetric sti ness matrix handling routines, which are used only one time at the beginning of each analysis. Direct linear equation system solvers are used in either case. Figure 7: Description of model details, as simulated in the FEAP program. 5.2. Results and discussion Two di erent kinds of results have been considered for the current discussion in this work. One is the total vertical force reaction associated to the imposed displacement load for each case of interface plane position. The analytical calculation can be easily done attending to the representation in Figure 5 and nding the total vertical reaction through an equivalent sti ness keq: Fz=uzkeq (64a) keq=k+k k++k(64b) k=SE L(64c) 20 The other output of interest is the average strain eld value on each side of the interface. Analytically, these values are easily obtainable by just making: "= E=Feq SE(65) Numerical simulation results coming from both enhancement types for the vertical force reaction can be appreciated in Figure 8. The rst and last points in this plot represent the cases in which homogeneous material distributions are given for one or the other material phases (sti er case and more compliant case, respectively), in which all numerical and analytical models naturally coincide. In the rst part of the curves where the sti er material is predominant, both weak discontinuity enhancements coincide for a while having a consistent error with respect to the analytical curve. At some point after having at least 20% volume fraction of the ( ) material region, the consistent enhancement error drops almost entirely, closely sticking to the analytical curve until the end of the graph. The typical enhancement maintains a smooth behavior with a consistent error, which also fades at the end when the ( ) material dominates completely. A maximum error of approximately 14% is observed with the typical enhancement. The sudden variations on the consistent enhancement can be explained by its natural stability conditions depending on node positions relative to the cutting plane interface, already discussed in Section 4.3. 0 2 4 6 8 10 Interface height position (mm)200040006000800010000Vertical reaction (N)Vertical Reaction Comparison Analytical solution Typical enhancement Consistent enhancement Figure 8: Vertical reaction on the entire lower face of the cube model for both enhancement types and the analytical solution. The displacement load imposed to this simple model will ideally produce a piece-wise, constant strain eld. The numerical approaches should be able to produce these constant strain regions taking the contri- bution of all elements on each side of the interface, aiming to have the least dispersion as possible. For this, theaverage strain eld and its dispersion (standard deviation) have been calculated on each side of the cube model interface for both enhancement approaches, and results have been compared to the analytical model. The analytical results, of course, do not show any kind of eld dispersion as they only exhibit a unique strain value. Figure 9 shows this comparison for the case of the (+) material side. Again, at homogeneous conditions all models coincide. Both enhancements start diverging at the beginning and the consistent enhancement quickly catches up the analytical behavior with a mild error. The typical enhancement, once again, has a smoother curve keeping a sustained error. The last data point at h= 10 mm is not shown since the (+) region ceases to exist. The strain dispersion, represented as a translucent cloud around each enhancement curve, seems de nitely more controlled for the consistent enhancement, albeit with apparently more erratic uctuations coming from the stability nature of this formulation. 21 0 2 4 6 8 10 Interface height position (mm)0.00040.00060.00080.00100.00120.00140.0016Strain (mm/mm)Average strain comparison in the + domain Analytical solution Typical enhancement Consistent enhancement Strain dispersion: typical Strain dispersion: consistentFigure 9: Average strain and corresponding dispersion for both weak discontinuity enhancements compared to analytical calculations on the (+) (sti er material) domain. Figure 10 shows the analysis on the remaining material region. In this case, it is the rst point at zero that is missing since there is no ( ) material. The average strain again favours the consistent enhancement that keeps a lower error through all conditions. The dispersion of the consistent model is also remarkably low, with the exception of one point at 30% volume for the ( ) material, where an outlier data point occurs. After model inspections, it is indeed found that at this height many nodes lie very close to the interface with separations as low as 0.03 mm, which seemingly compromise the stability of this enhancement as discussed in section 4.3. The impact of this outlier is not noticeable for the case of the analysis on the (+) domain, since ath= 3 mm there are considerably more homogeneous elements made up of the (+) material, which help to smooth these statistics. 6. Conclusions A detailed analysis has been made on the use of weak discontinuity approaches within the E-FEM framework to model material heterogeneities. Kinematic and variational foundations have been stated to identify a set of consistency requirements for the general construction of weak discontinuity displacement enhancement elds. Based on the consideration of these requirements, two eld proposals have been derived: one which has been already used in the works of [14, 35, 31] and other which is fully consistent with all requirements. Simulations have been made on a simple bi-material system subject to an axial load to assess the performance of the enhancements. It has been seen that, in general, the consistent enhancement has a better performance than the typical enhancement at representing the expected response, but is more prone to uncontrolled uctuations. It can be also argued, however, that under high quality mesh conditions and moderate di erence between the elastic properties of material phases (a ratio of 5 explored in this work), the typical eld will have a reasonable behavior that will allow for sound estimations at local and global level mechanical outputs. The formulation by itself naturally presents unconditional stability and retains the symmetry of global sti ness matrices. The internal calculations required to particularise the function parameters are also considerably simpler than in the case of the consistent enhancement version. For these reasons, the authors of this study recommend its use whenever possible, especially when dealing with large 22 0 2 4 6 8 10 Interface height position (mm)0.0020.0030.0040.0050.006Strain (mm/mm)Average strain comparison in the Ω−domain Analytical solution Typical enhancement Consistent enhancement Strain dispersion: typical Strain dispersion: consistentFigure 10: Average strain and corresponding dispersion for both weak discontinuity enhancements compared to analytical calculations on the ( ) (more compliant) domain. and complex numerical models where solution times and stability are crucial for the success of the numerical analysis project. The use of the consistent weak discontinuity enhancement may nd a better niche on problems where high accuracy is required for eld shape calculations and where the use of an asymmetric solver poses no problem for a given FEM numerical solution platform. The authors in this work have found that the stability nuisances on the consistent formulation can actually be avoided at a large extent by just locally displacing the interface location sightly away from its compromising position, in such a way that the elemental volume distribution is just mildly changed without signi cantly a ecting the physical representation aimed by the element in question. The performance of the models has not been tested for badly shaped elements, where it is expected to have a worse behavior for the typical enhancement. Both enhancements also tend to have a slightly worse performance for locally describing sti er materials within a bi-material set. It would also be relevant to study material sets with more prominent di erences in elastic properties to assess the enhancements' behavior. It can also be shown that a higher polynomial degree proposal for the weak discontinuity enhancement will be able to ful ll allconsistency requirements at once . The study of its particular shape and properties could be an interesting topic for future work. This would allow an enhancement function both as robust as the typical and as accurate as the consistent one. In either case, the authors of this work nish by stating that the use of weak discontinuity enhancements for the representation of material heterogeneities remains a reliable and ecient numerical method for approaching the problem of multiple material phases representation featuring non-adapted meshes.