General Consistency of Strong Discontinuity Kinematics in Embedded Finite Element Method (E-FEM) Formulations

Abstract

This paper performs an in-depth study of the theoretical basis behind the strong discontinuity methods to improve local fracture simulations using the Embedded Finite Element Method (E-FEM). The process starts from a review of the elemental enhancement functions found in current E-FEM literature, providing the reader a solid context of E-FEM fundamentals. A set of theoretical pathologies is then discussed, which prevent current frameworks from attaining full kinematic consistency and introduce unintended mesh dependencies. Based on this analysis, a new proposal of strong discontinuity enhancement functions is presented considering generalised fracture kinematics in a full tridimensional setting and a more robust definition of internal auxiliary functions. Element-level simulations are performed to compare the outputs within a group of selected E-FEM 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, demonstrating an increase in robustness that might drive the usefulness of E-FEM techniques for fracture simulations to a higher level.

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 fields. 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 efficient 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 performed using a traditional FEM solver framework, requiring no modifications to it unless the framework is deliberately enriched with a non-local analysis feature such as, for example, the use of a crack path field method to ensure continuity between elements [2].

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]. The latter have a more robust and deeper definition of their kinematics, attacking directly the definition 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 such as X-FEM or gained much ground in 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 from gaining a consensual solid stance as a reliable and accurate structural numerical simulation approach for material fractures. These have been recognised by many authors working in the field: 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 of how the basic definitions of strong discontinuity kinematics have been introduced so far on the E-FEM framework (Section 2 and Section 3) to identify the roots of its potential theoretical faults and related numerical issues (Section 4). Some authors provided their valuable insights on these challenges and proposed theoretical enhancements as a workaround [8,11,12,13]. This work attempts to consolidate these efforts 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 specific E-FEM strong discontinuity enhancement scheme (Section 5). Simulations at the element level are performed to compare the results to older enhancement proposals (Section 6), with encouraging conclusions (Section 7). For the reader’s convenience, Figure 1 describes the general layout and progression of the aforementioned discussions.

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 such as reclosure, local compression and explicit local friction, as well as the possibility to consider other kinds of discontinuities in a simultaneous fashion.

2. Analysis of Basic Internal Strong Discontinuity Kinematics

The starting point of every fracture representation within the E-FEM framework is to enrich the basic displacement field 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:

$$\mathit{u}=\overline{\mathit{u}}+{\mathcal{H}}_{{\mathrm{\Gamma }}_d}\left[\left|\mathit{u}\right|\right],$$

(1)

where the theoretical displacement of the body $\mathit{u}$ is set as the composition of a regular displacement $\overline{\mathit{u}}$ and a jump normally introduced by means of a Heaviside function having the location of the fault surface ${\mathrm{\Gamma }}_d$ as its main trigger. An internal fracture kinematics vector $\left[\left|\mathit{u}\right|\right]$ is used to represent the general motion of one of the fractured regions with respect to the other. Figure 2 and Figure 3 are very common illustrations used by authors in the field to visualise the discontinuous enrichment for the fracture of a body split into two regions ${\mathsf{\Omega }}^{+}$ and ${\mathsf{\Omega }}^{-}$.

The origin, shape and orientation of the surface ${\mathrm{\Gamma }}_d$ representing the fracture remain completely independent from this analysis, and will depend on the specific strain localisation criterion chosen for studying a given material. As the general philosophy of the E-FEM framework tends to aim for formulation simplicity, 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 performed on a general 3D domain unless specific remarks for 2D or 1D cases are made.

While Equation (1) already allows for the calculation of a valid strain field directly through a symmetric gradient operator ${\nabla }^s\left(•\right)=\frac{1}{2}\left[\nabla \left(•\right)+{\nabla }^t\left(•\right)\right]$, there is a couple of important basic remarks worth making at this point.

2.1. Kinematic Consistency of Boundary Condition Imposition

A discretised element whose kinematics are governed by Equation (1) does not allow a consistent imposition of boundary conditions on any traditional global FEM solution engine. This is because the displacements on the borders of the ${\mathsf{\Omega }}^{-}$ domain of the element (i.e., the nodes on it) are directly reachable through the vector $\overline{\mathit{u}}$, while the displacements on the borders of ${\mathsf{\Omega }}^{+}$ are rather a composition of $\overline{\mathit{u}}$ and $\left[\left|\mathit{u}\right|\right]$. We can use the vector $\overline{\mathit{u}}$ alone to impose boundary conditions on ${\mathsf{\Omega }}^{-}$, but not on ${\mathsf{\Omega }}^{+}$ since the actual displacement in this region remains a composition of two different vectors $\overline{\mathit{u}}+\left[\left|\mathit{u}\right|\right]$.

This issue has been well identified even from the first evolving steps of the framework and also treated by Oliver [10] through the introduction of the so-called auxiliary function $\phi$. The function $\phi$ is defined as to allow for the existence of a new displacement vector

$$\widehat{\mathit{u}}=\overline{\mathit{u}}+\phi \left[\left|\mathit{u}\right|\right],$$

(2)

by letting $\phi$ have a single property: it will have the value of 1 on any node lying on the border of ${\mathsf{\Omega }}^{+}$ and zero in those on the remaining ${\mathsf{\Omega }}^{-}$ border:

$$\phi =\left\{\begin{array}{cc}1\hfill & \mathit{x}\in {\mathsf{\Omega }}^{+}\hfill \\ 0\hfill & \mathit{x}\in {\mathsf{\Omega }}^{-}\hfill \end{array}\right.,i=1,..,N$$

(3)

where ${\mathit{x}}_i$ denotes a nodal position. This is the only formal mathematical requirement as stated in [10] apart from the basic fact that $\phi$ should be continuous. The most common form of $\phi$ taken from the literature is yet the simplest one to satisfy Equation (3) for a given element geometry. For example, for the linear tetrahedron, the simplest $\phi$ remains a composition of standard linear interpolation functions. In this case, $\phi$ also happens to be unique. Analogous constructions for $\phi$ can be easily made for a 2D triangle and for 1D cases (Figure 4).

If $\phi$ complies with the previous definition, then the new vector $\widehat{\mathit{u}}$ can be seen as a continuous fusion of the local displacement field and the jump vector associated with the fracture, having border displacements that coincide with those of $\mathit{u}$. Thus, for the practical purposes of nodal displacements in a discretised element, $\widehat{\mathit{u}}$ contains 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 field within the element does not coincide with that of $\mathit{u}$. $\widehat{\mathit{u}}$ can be seen as a special version of $\mathit{u}$ where the discontinuity jump is continuously fused by $\phi$ along the element domain. $\phi$ could be thus interpreted as a sort of regularisation function (see Figure 3).

With this definition in mind, Equation (1) is rewritten as

$$\mathit{u}=\widehat{\mathit{u}}+\left({\mathcal{H}}_{\mathrm{\Gamma }}-\phi \right)\left[\left|\mathit{u}\right|\right],$$

(4)

where the previous vector $\overline{\mathit{u}}$ has been replaced by the newly defined vector $\widehat{\mathit{u}}$, and the discontinuity base function now becomes composed of both a Heaviside and $\phi$. From this point, authors generally stop referring to $\overline{\mathit{u}}$ and simply take $\widehat{\mathit{u}}$ as the de facto standard solution vector for the reasons mentioned above, and Equation (4) is taken as the reference for developing all derived mechanical fields.

While the introduction of $\phi$ solves a fundamental consistency problem on boundary conditions imposition, it is also responsible for other undesired effects on general element kinematics and fracture mode coupling. A deeper discussion is presented in Section 4.

2.2. Kinematics at Terminal Separation Conditions and Meaning of $\left[\left|\mathit{u}\right|\right]$

Post-localisation behaviour within an element is generally modelled 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 influencing forces from ${\mathsf{\Omega }}^{+}$ to ${\mathsf{\Omega }}^{-}$ when the split bodies have sufficiently separated. This state of failure will be referred to as a terminal separation condition. In such state, the traction happening at the fracture surface should be zero, and the internal fracture kinematics vector $\left[\left|\mathit{u}\right|\right]$ will practically represent the rigid body motion of the ${\mathsf{\Omega }}^{+}$ domain with respect to ${\mathsf{\Omega }}^{-}$, depending on the definition chosen for the vector $\left[\left|\mathit{u}\right|\right]$.

$\left[\left|\mathit{u}\right|\right]$ in most of the current E-FEM literature has been simply defined as a one-to-three components vector of constant values in space at the element scale, describing a rigid body translation of ${\mathsf{\Omega }}^{+}$ with respect to ${\mathsf{\Omega }}^{-}$. However, as $\left[\left|\mathit{u}\right|\right]$ stands as an internal element definition with no strong global restrictions, it could be arranged to describe more complex kinematic modes such as rigid body rotations, having non-constant component behaviour through the failure surface ${\mathrm{\Gamma }}_d$.

$\left[\left|\mathit{u}\right|\right]$ thus remains a free modelling choice in this sense. However, whatever this choice may be, its behaviour should be consistent with the standard element kinematics (whether $\mathit{u}$ or $\widehat{\mathit{u}}$). An obvious example is that the magnitude of a constant $\left[\left|\mathit{u}\right|\right]$ in the normal direction $\widehat{\mathit{n}}$ to 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 Figure 5 for a 2D constant stress triangle (CST). A relative displacement $d_n$ is imposed on this element and, assuming post-localisation, a corresponding crack separation ${\left[\left|u\right|\right]}_n$ develops. By common sense, it should be expected that ${\left[\left|u\right|\right]}_n\le d_n$, where at terminal separation conditions, we could probably have ${\left[\left|u\right|\right]}_n\simeq d_n$ depending on the specific model.

Considering this, it is important to remark the framework does not explicitly guarantee any kinematic consistency between the internal variable $\left[\left|\mathit{u}\right|\right]$ and a standard global variable $\widehat{\mathit{u}}$ or any derived quantities such as $d_n$. This is due to the same fact that $\left[\left|\mathit{u}\right|\right]$ is an internal element definition. This is especially true when reaching a terminal separation condition, where there is no further influence of any post-localisation law derived through CSDA or DSDA that would imply any additional coercion to $\left[\left|\mathit{u}\right|\right]$ in any way. In a terminal separation condition, there is no other governing law than that of the purely natural kinematics arising from the modelling choices made in the basic formulation of the strong discontinuity, such as $\phi$.

In the end, it should be clear that the physical meaning of $\left[\left|\mathit{u}\right|\right]$ 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 Section 5.

Having all this in mind, the strain field is extracted by taking the symmetrical gradient ${\nabla }^s\left(•\right)$ in Equation (4). This study will assume a constant vector $\left[\left|\mathit{u}\right|\right]$ through ${\mathrm{\Gamma }}_d$, which implies that this strong discontinuity study will be limited to a rigid body translation of ${\mathsf{\Omega }}^{+}$ with respect to ${\mathsf{\Omega }}^{-}$. The strain field can be thus expressed as:

$$\begin{array}{c}\mathit{\epsilon }={\nabla }^s\mathit{u}={\nabla }^s\widehat{\mathit{u}}+{\nabla }^s\left[\left({\mathcal{H}}_{\mathrm{\Gamma }}-\phi \right)\left[\left|\mathit{u}\right|\right]\right]\\ \Rightarrow \mathit{\epsilon }={\nabla }^s\widehat{\mathit{u}}-{\left(\nabla \phi \otimes \left[\left|\mathit{u}\right|\right]\right)}^s+{\left(\nabla {\mathcal{H}}_{\mathrm{\Gamma }}\otimes \left[\left|\mathit{u}\right|\right]\right)}^s\\ \Rightarrow \mathit{\epsilon }={\nabla }^s\widehat{\mathit{u}}\underset{{\widehat{\mathit{\epsilon }}}_b}{\underbrace{-{\left(\nabla \phi \otimes \left[\left|\mathit{u}\right|\right]\right)}^s}}\underset{{\widehat{\mathit{\epsilon }}}_u}{\underbrace{+{\delta }_{\mathrm{\Gamma }}{\left(\widehat{\mathit{n}}\otimes \left[\left|\mathit{u}\right|\right]\right)}^s}},\end{array}$$

(5)

where the first tensor product is identified as the bounded part of the strong discontinuity strain ${\widehat{\mathit{\epsilon }}}_b$ and the second tensor product is the unbounded strain ${\widehat{\mathit{\epsilon }}}_u$. ${\delta }_{\mathrm{\Gamma }}$ is the Dirac delta, with a trigger position corresponding to the failure surface ${\mathrm{\Gamma }}_d$, just as the Heaviside was described previously.

3. In-Depth Analysis of Variational Foundations

The use of customised field enhancement functions on the E-FEM framework requires a variational principle that grants sufficient flexibility to effectively manage multiple independent mechanical fields when making element discretisation. After all, the E-FEM remains a mixed FE formulation. This is the main reason why the three-field Hu–Washizu principle has been the preferred approach for almost all authors working on the framework. There are, however, many details usually omitted in the description of its application that merit formal mathematical justification. 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 fields of displacement, strain and stress $\left(\mathit{u},\mathit{\epsilon },\mathit{\sigma }\right)$ as the point of departure:

$$\begin{array}{cc}\hfill & I_{HW}\left(\mathit{u},\mathit{\epsilon },\mathit{\sigma }\right)=\frac{1}{2}\underset{\mathsf{\Omega }}{\int }\overline{\overline{\mathit{\epsilon }}}:\mathcal{C}:\overline{\overline{\mathit{\epsilon }}}dV-\underset{\mathsf{\Omega }}{\int }\overline{\overline{\mathit{\sigma }}}:\overline{\overline{\mathit{\epsilon }}}dV\hfill \\ & +\underset{\mathsf{\Omega }}{\int }\overline{\overline{\mathit{\sigma }}}:\nabla \mathit{u}dV-\underset{\mathsf{\Omega }}{\int }{\mathit{f}}_b·\mathit{u}dV-\underset{\partial \mathsf{\Omega }}{\int }\mathit{t}·\mathit{u}dA,\hfill \end{array}$$

(6)

where the double bar notation $\overline{\overline{•}}$ specifies a mechanical field on its 2nd order tensor format (the omission of it would mean that the field is on the Voigt vector notation), $\mathcal{C}$ is a fourth-order linear elasticity tensor, ${\mathit{f}}_b$ is a vector of body forces and $\mathit{t}$ is 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 with each of the displacement, strain and stress fields $\left(\delta I_{HW}\left(\delta \mathit{u},\delta \mathit{\epsilon },\delta \mathit{\sigma }\right)=0\right)$, leading to three independent equations, already written using a Voigt vector notation:

$$\begin{array}{c}\hfill \underset{\mathsf{\Omega }}{\int }\partial \delta {\mathit{u}}^t\mathit{\sigma }dV-\underset{\mathsf{\Omega }}{\int }\delta {\mathit{u}}^t{\mathit{f}}_{b}dV-\underset{\partial \mathsf{\Omega }}{\int }\delta {\mathit{u}}^t\mathit{t}dA=0\end{array}$$

(7a)

$$\begin{array}{c}\hfill \underset{{\mathsf{\Omega }}_e}{\int }\delta {\mathit{\sigma }}^t\left(\partial \mathit{u}-\mathit{\epsilon }\right)dV=0\end{array}$$

(7b)

$$\begin{array}{c}\hfill \underset{{\mathsf{\Omega }}_e}{\int }\delta {\mathit{\epsilon }}^t\left(\mathit{\sigma }\left(\mathit{\epsilon }\right)-\mathit{\sigma }\right)dV=0\end{array}$$

(7c)

It is worth noticing the difference made in Equation (7c) between $\mathit{\sigma }\left(\mathit{\epsilon }\right)$, which is the stress field calculated from the constitutive relations (in this case through an elastic relation $\mathit{\sigma }\left(\mathit{\epsilon }\right)=\mathbf{C}\mathit{\epsilon }$), and $\mathit{\sigma }$, which is the independent stress field to be found in the analysis. These fields are not necessarily equal, and the variational analysis will provide a specific relation between them.

The system of Equations (7a)–(7c) is the most commonly used form of the Hu–Washizu principle as seen in the E-FEM literature. The solution for the system, in this continuous format, inevitably returns the equilibrium of the body (already expressed in Equation (7a)), the strong kinematic compatibility and constitutive relations:

$$\begin{array}{c}\hfill \nabla \mathit{u}-\overline{\overline{\mathit{\epsilon }}}=0\end{array}$$

(8a)

$$\begin{array}{c}\hfill \overline{\overline{\mathit{\sigma }}}\left(\mathit{\epsilon }\right)-\overline{\overline{\mathit{\sigma }}}=0\end{array}$$

(8b)

From this point, any further step requires stating a field discretisation strategy, and flexibility will be drawn from Equations (8a) and (8b) as required, having cases where these are not necessarily satisfied 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 fields of displacement, strain and stress $\left(\mathit{u},\mathit{\epsilon },\sigma \right)$ along with the other three fields of displacement variation, strain variation and stress variation $\left(\delta \mathit{u},\delta \mathit{\epsilon },\delta \sigma \right)$ are all completely independent. Therefore, each of the six fields will require a unique discretisation choice and its entirely up to each study to determine the most convenient field discretisation setup, keeping in mind that the main goal of the framework is to represent the effects of the strong discontinuity in the most physically consistent, yet numerically efficient way. Apart from this, the only formal restriction is to satisfy the system of Equations (7a)–(7c) 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 Field Discretisation

For discretising the displacement field $\mathit{u}$, authors naturally take advantage of the facts discussed in Section 2.1, and take only $\widehat{\mathit{u}}$ into account:

$$\mathit{u}=\widehat{\mathit{u}}=\mathbf{N}\mathit{d},$$

(9)

with $\mathbf{N}$ as a standard interpolation matrix and $\mathit{d}$ the standard nodal displacement vector. The variation of the displacement field $\delta \mathit{u}$ is also discretised in the same way, taking the displacement vector variation $\delta \mathit{d}$:

$$\delta \mathit{u}=\mathbf{N}\delta \mathit{d}$$

(10)

Equation (9) implies that the displacement field to be found through this strategy does not exactly correspond to the initial definition of $\mathit{u}$ in the interior of the element but will perfectly coincide with it on the nodes, and therefore the vector $\mathit{d}$ contains always the real displacement of the element through its nodes. This is considered as enough information to be extracted from the displacement field as there is generally no interest in studying it in detail considering that the vector $\left[\left|\mathit{u}\right|\right]$ already contains embedded crack displacement information. The advantages emerge when dealing with Equation (7a), as the simplification of the displacement field (especially the displacement variation field) will allow reaching the classical standard FEM global solution form without any additional terms involved. Note that Equation (9) does not imply Equation (10) in any way: having the same interpolation matrix $\mathbf{N}$ is totally a strategic choice that is convenient at this point. This present work will comply with it.

3.1.2. Strain Field Discretisation

Discretisation of the strain field $\epsilon$ is performed following all the terms described in Equation (5), giving this field the role for modelling all strong discontinuity kinematics on the element, including bound and unbound terms:

$$\mathit{\epsilon }=\mathbf{B}\mathit{d}+{\mathbf{G}}_s\left[\left|\mathit{u}\right|\right]=\mathbf{B}\mathit{d}+\left({\mathbf{G}}_{sb}+{\mathbf{G}}_{su}\right)\left[\left|\mathit{u}\right|\right],$$

(11)

where $\mathbf{B}$ is the standard matrix corresponding to the partial differentiations of $\mathbf{N}$ and ${\mathbf{G}}_s$ is the corresponding matrix for the strong discontinuity terms, including a bounded part ${\mathbf{G}}_{sb}$ and an unbounded part ${\mathbf{G}}_{su}$. The specific form of the matrix operators ${\mathbf{G}}_{sb}$ and ${\mathbf{G}}_{su}$ is constructed by considering Equation (5) alone. Their structure is fixed and determined by the natural kinematics of the strong discontinuity model assumed. A strain field following such idea within the E-FEM framework is referred to as a Kinematically Enhanced Strain (KES).

The choice for the strain variation field stays analogous by considering the same standard strain term associated with standard displacement variations, but leaving a choice open for a matrix ${\mathbf{G}}_s^{*}$ which, for reasons explained later (and also a deliberate choice), will adopt the same bounded-unbounded composition ${\mathbf{G}}_{sb}^{*}$, ${\mathbf{G}}_{su}^{*}$ and will be associated with the variation of the crack displacements $\delta \left[\left|\mathit{u}\right|\right]$:

$$\delta \mathit{\epsilon }=\mathbf{B}\delta \mathit{d}+{\mathbf{G}}_s^{*}\delta \left[\left|\mathit{u}\right|\right]=\mathbf{B}\delta \mathit{d}+\left({\mathbf{G}}_{sb}^{*}+{\mathbf{G}}_{sn}^{*}\right)\delta \left[\left|\mathit{u}\right|\right]$$

(12)

In this case, the specific form of the matrix operators ${\mathbf{G}}_{sb}^{*}$ and ${\mathbf{G}}_{su}^{*}$ is not determined by strong discontinuity kinematics directly. It is rather derived from the demands of the variational principle itself through the satisfaction of Equation (7), as it will be shown in the following sections. This means in principle that the components of this strain field attend more to static considerations. A strain field following such idea within the E-FEM framework is referred to as an Enhanced Assumed Strain (EAS).

3.1.3. Stress Field Discretisation

Finally, the independent stress $\mathit{\sigma }$ and stress variation fields $\delta \mathit{\sigma }$ are simply designated in a free-form fashion:

$$\begin{array}{cc}\hfill \mathit{\sigma }& =\mathbf{S}\mathit{s}\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill \delta \mathit{\sigma }& ={\mathbf{S}}^{*}\delta \mathit{s},\hfill \end{array}$$

(13b)

where $\mathit{s}$ is an independent stress vector with the state of stresses associated with each node of a given element (e.g., a six-component vector for a 3D, constant field). $\delta \mathit{\sigma }$ would be its corresponding variation. The interpolation matrices $\mathbf{S}$ and ${\mathbf{S}}^{*}$ are absolutely unrelated to any of the previously defined entities in the other mechanical fields.

3.1.4. Calculated Stress Field Discretisation

The stress field coming from the constitutive relation $\mathit{\sigma }\left(\mathit{\epsilon }\right)$ will be calculated by taking only the bounded part of $\epsilon$. The boundedness of $\mathit{\sigma }\left(\mathit{\epsilon }\right)$ despite having unbounded terms within $\mathit{\epsilon }$ has been discussed since the first attempts of consolidating the E-FEM framework [15]. It remains an accepted practice in current theoretical developments. Considering this, this stress field is expressed as:

$$\mathit{\sigma }\left(\mathit{\epsilon }\right)=\mathbf{C}{\mathit{\epsilon }}_b=\mathbf{C}\left(\mathbf{B}\mathit{d}+{\mathbf{G}}_{sb}\left[\left|\mathit{u}\right|\right]\right)$$

(14)

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 definitions of the chosen discretisation strategy in the previous section are basically thrown into the variational principle represented by Equation (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 ${\mathbf{S}}^{*}$

We will start by taking Equation (7b), which brings the following expression after considering the definitions of $\mathit{u}$, $\mathit{\epsilon }$, the fact that $\left[\left|\mathit{u}\right|\right]$ is a constant and that $\delta \mathit{s}$ remains completely arbitrary:

$$\begin{array}{c}\underset{{\mathsf{\Omega }}_e}{\int }\delta {\mathit{\sigma }}^t\left(\partial \mathit{u}-\mathit{\epsilon }\right)dV=\underset{{\mathsf{\Omega }}_e}{\int }\delta {\mathit{s}}^t{\mathit{S}}^{*t}\left(\mathit{B}\mathit{d}-\left[\mathit{B}\mathit{d}+{\mathbf{G}}_s\left[\left|\mathit{u}\right|\right]\right]\right)dV=\delta {\mathit{s}}^t\underset{{\mathsf{\Omega }}_e}{\int }{\mathit{S}}^{*t}{\mathbf{G}}_{s}dV\left[\left|\mathit{u}\right|\right]=\mathbf{0}\end{array}$$

$$\begin{array}{c}\Rightarrow \underset{{\mathsf{\Omega }}_e}{\int }{\mathit{S}}^{*t}{\mathbf{G}}_{s}dV=\mathbf{0}\end{array}$$

(15)

Equation (15) demands an orthogonality between ${\mathbf{S}}^{*t}$ and ${\mathbf{G}}_s$ in a weak sense. Given that ${\mathbf{G}}_s$ has the specific purpose of modelling the kinematics of the strong discontinuity, it shall not be modified to satisfy Equation (15). Thus, ${\mathbf{S}}^{*}$ has to comply. Note that this variational relation does not impose any form of restrictions to critical functions, as ${\mathbf{S}}^{*}$ will not give shape to any real mechanical field. It could be stated that ${\mathbf{S}}^{*}$ becomes just a support interpolation matrix that allows ${\mathbf{G}}_s$ to exist. We can choose whatever kinematics we want for the strong discontinuity implied in ${\mathbf{G}}_s$, as long as we can actually prove that this support interpolation matrix ${\mathbf{S}}^{*}$ satisfying Equation (15) exists.

At this point, the specific form of $\phi$ comes into play, which is important to discuss because it determines the overall form of ${\mathbf{G}}_s$. ${\mathbf{G}}_s$ in turn is composed of bounded and unbounded parts ${\mathbf{G}}_{sb}$ and ${\mathbf{G}}_{su}$. Let us suppose now that we work on a linear tetrahedron, and thus with a linear $\phi =a+bx+cy+dz$, which is one of the most prominent choices made in the E-FEM literature. With this assumption and considering the structure of Equation (5), expressions for matrix operators ${\mathbf{G}}_{sb}$ and ${\mathbf{G}}_{su}$ on $\left[\left|\mathit{u}\right|\right]$ can be devised as:

$${\mathbf{G}}_{sb}=-\left[\begin{array}{ccc}{\phi }_{,x}& 0& 0\\ 0& {\phi }_{,y}& 0\\ 0& 0& {\phi }_{,z}\\ {\phi }_{,y}& {\phi }_{,x}& 0\\ 0& {\phi }_{,z}& {\phi }_{,y}\\ {\phi }_{,z}& 0& {\phi }_{,x}\end{array}\right]=-{\nabla }^{\prime },{\mathbf{G}}_{su}=-{\delta }_{\mathrm{\Gamma }}\left[\begin{array}{ccc}{n}_x& 0& 0\\ 0& n_y& 0\\ 0& 0& n_z\\ n_y& n_x& 0\\ 0& n_z& n_y\\ n_z& 0& n_x\end{array}\right]=-{\delta }_{\mathrm{\Gamma }}{\mathbf{H}}_s$$

(16)

where ${\phi }_x,{\phi }_y,{\phi }_z$ are the partial derivatives of $\phi$ and $n_x,n_y,n_z$ are the components of the normal vector $\widehat{\mathit{n}}$ to the crack plane. If we assume for a first trial that ${\mathbf{S}}^{*}$ remains constant, there will be only one single option to define it: the identity matrix ${\mathbf{I}}_6$ (or a scalar multiple of it). Equation (15) can be then worked using integral properties of the Dirac delta to show that:

$$\underset{{\mathsf{\Omega }}_e}{\int }{\mathbf{S}}^{*t}{\mathbf{G}}_{s}dV=\underset{{\mathsf{\Omega }}_e}{\int }{\mathbf{I}}_6\left(-{\nabla }^{\prime }+{\delta }_{\mathrm{\Gamma }}{\mathbf{H}}_s\right)dV=V_e{\nabla }^{\prime }+A_{\mathrm{\Gamma }}{\mathbf{H}}_s,$$

(17)

where $V_e$ is the element volume and $A_{\mathrm{\Gamma }}$ is the area of the crack surface. It can be easily seen that, as ${\nabla }^{\prime }$ and V do not have any kind of relation to $A_{\mathrm{\Gamma }}$ and ${\mathbf{H}}_s$, Equation (17) can never be asserted as zero and thus a constant ${\mathbf{S}}^{*}$ satisfying Equation (15) does not exist. Note that a constant ${\mathbf{S}}^{*}$ was chosen to be consistent with the natural order for a stress field associated with a standard linear displacement field on a linear tetrahedron. There is no option but to go for a non-conventional definition for ${\mathbf{S}}^{*}$. Fortunately, Hu–Washizu allows for this.

Using a linear ${\mathbf{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 define the four associated interpolation functions $S_i^{*}=a_i^{*}+b_i^{*}x+c_i^{*}y+d_i^{*}z$ to these points. Equation (15) can then be worked as:

$$-\left(\underset{{\mathsf{\Omega }}_e}{\int }{\mathbf{S}}^{*t}dV\right){\mathbf{\nabla }}^{\prime }-\left(\underset{A_{\mathrm{\Gamma }}}{\int }{\mathbf{S}}^{*t}dA\right){\mathbf{H}}_s=0,$$

(18)

leading to a linear system of the form:

$$S_i^{*V_e}{\phi }_{,j}+S_i^{*A_{\mathrm{\Gamma }}}{n}_j=0,i=1,2,3,4j=x,y,z$$

(19)

where $S_i^{*V_e}$ are the interpolation functions $S_i^{*}$ integrated through the element volume $V_e$, and $S_i^{*A_{\mathrm{\Gamma }}}$ are these same functions but integrated over the crack surface $A_{\mathrm{\Gamma }}$. As there are a total of 12 independent equations to fulfil and ${\mathbf{S}}^{*}$ has 12 free parameters, a unique solution for ${\mathbf{S}}^{*}$ exists.

Note that this system will determine a very specific set of nodal positions $(x_i^{*},y_i^{*},z_i^{*})$ that are absolutely unrelated to real element nodal positions. The base of ${\mathbf{S}}^{*}$ is not the standard element base and thus it will be entirely different to that coming from $\mathbf{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 $\phi$ or even a non-constant $\left[\left|\mathit{u}\right|\right]$, it can be shown by the aforementioned approach that a sufficiently elaborate ${\mathbf{S}}^{*}$ can be chosen to satisfy Equation (15). This last point is the main takeaway from the work in Equation (7b).

3.2.2. The Bridge between Real and Constitutive Stress Fields

Equation (7c) can be worked by using Equation (12) and noting that $\delta \mathit{d}$ and $\delta \left[\left|\mathit{u}\right|\right]$ remain independent, giving rise to two new equations:

$$\begin{array}{c}\underset{{\mathsf{\Omega }}_e}{\int }\delta {\mathit{\epsilon }}^t\left(\mathit{\sigma }\left(\mathit{\epsilon }\right)-\mathit{\sigma }\right)dV=\underset{{\mathsf{\Omega }}_e}{\int }\left(\delta {\mathit{d}}^t{\mathbf{B}}^t+\delta {\left[\left|\mathit{u}\right|\right]}^t{\mathbf{G}}_s^{*t}\right)\left(\mathit{\sigma }\left(\mathit{\epsilon }\right)-\mathbf{S}\mathit{s}\right)dV\\ =\delta {\mathit{d}}^t\underset{{\mathsf{\Omega }}_e}{\int }{\mathbf{B}}^t\left(\mathit{\sigma }\left(\mathit{\epsilon }\right)-\mathbf{S}\mathit{s}\right)dV+\delta {\left[\left|\mathit{u}\right|\right]}^t\underset{{\mathsf{\Omega }}_e}{\int }{\mathbf{G}}_s^{*t}\left(\mathit{\sigma }\left(\mathit{\epsilon }\right)-\mathbf{S}\mathit{s}\right)dV=0\end{array}$$

$$\begin{array}{c}\Rightarrow \underset{\mathsf{\Omega }}{\int }{\mathbf{B}}^t\mathit{\sigma }\left(\mathit{\epsilon }\right)dV=\underset{\mathsf{\Omega }}{\int }{\mathbf{B}}^t\mathbf{S}\mathit{s}dV\end{array}$$

(20a)

$$\begin{array}{c}\Rightarrow \underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_s^{*t}\mathit{\sigma }\left(\mathit{\epsilon }\right)dV=\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_s^{*t}\mathbf{S}\mathit{s}dV\end{array}$$

(20b)

Equation (20a) is the only equation that allows to establish a direct relation between the calculated stresses from the constitutive law $\mathit{\sigma }\left(\mathit{\epsilon }\right)$ and the interpolated real stress $\mathbf{S}\mathit{s}$, since the form of $\mathbf{B}$ is known beforehand for a given standard element geometry. Equation (20a) is a weak equality involving ${\mathbf{B}}^t$ as a weight factor. The form of $\mathbf{S}$ always remains a choice. For a constant real stress field, $\mathbf{S}$ will be a $6\times 6$ identity matrix, whereas for a linear real stress interpolation, $\mathbf{S}$ will be a dense $6\times 24$ matrix, and so on. Based on this choice and the form of $\mathbf{B}$, Equation (20a) may or may not grant a strong equality between constitutive stresses and real interpolated stresses. In the case that $\mathbf{S}$ is chosen to match the order of $\mathit{\sigma }\left(\mathit{\epsilon }\right)$, Equation (20a) effectively reduces to a strong equality. Otherwise, for the case in which we have a constant $\mathbf{B}$ and a constant $\mathbf{S}$, a real stress vector $\mathit{s}$ can still be calculated through a simple volume averaged integral regardless of the order of $\mathit{\sigma }\left(\mathit{\epsilon }\right)$:

$$\mathit{s}=\frac{1}{V_e}\underset{\mathsf{\Omega }}{\int }\mathit{\sigma }\left(\mathit{\epsilon }\right)dV$$

(21)

Note that, in any case, real stress calculations remain explicit and do not require methods such as a least-squares approach as in other mixed element variational frameworks [1,20].

3.2.3. The Bridge between Real and Constitutive Traction Vectors

Equation (20b) is treated differently. As mentioned before, ${\mathbf{G}}_s^{*}$ is assumed to be composed of bounded and unbounded terms as its counterpart ${\mathbf{G}}_s$, and in this section, it will be clarified why. Usually, authors in the field choose at this point to enforce weak orthogonality conditions to each separate side of Equation (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:

$$\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_s^{*t}\mathbf{S}\mathit{s}dV=\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_s^{*t}\mathit{\sigma }\left(\mathit{\epsilon }\right)dV=0$$

(22)

The main idea is to derive an expression for the real traction vector $\mathit{T}$ based on the real stress field $\left(\mathbf{S}\mathit{s}\right)$, and then connecting it to the traction vector based on constitutive stress ${\mathit{T}}_{\mathit{\epsilon }}$, so that $\mathit{T}$ can be explicitly calculated. This can be achieved by developing the left and centre side of Equation (22), leading to:

$$\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{sb}^{*t}\mathbf{S}\mathit{s}dV+\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{su}^{*t}\mathbf{S}\mathit{s}dV=\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{sb}^{*t}\sigma \left(\mathit{\epsilon }\right)dV+\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{su}^{*t}\mathit{\sigma }\left(\mathit{\epsilon }\right)dV$$

(23)

If we propose the specific form of the unbounded ${\mathbf{G}}_{su}^{*}$ as exactly the same as ${\mathbf{G}}_{su}={\delta }_{\mathrm{\Gamma }}{\mathbf{H}}_s$, we can use again the Dirac delta’s properties to arrive at:

$$\begin{array}{c}\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{sb}^{*t}\mathbf{S}\mathit{s}dV+\underset{{\mathrm{\Gamma }}_d}{\int }{\mathbf{H}}_s^t\mathbf{S}\mathit{s}dA=\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{sb}^{*t}\mathit{\sigma }\left(\mathit{\epsilon }\right)dV+\underset{{\mathrm{\Gamma }}_d}{\int }{\mathbf{H}}_s^t\mathit{\sigma }\left(\mathit{\epsilon }\right)dA\end{array}$$

$$\begin{array}{c}\Rightarrow \underset{{\mathrm{\Gamma }}_d}{\int }\mathit{T}dA=\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{sb}^{*t}\left[\mathit{\sigma }\left(\mathit{\epsilon }\right)-\mathbf{S}\mathit{s}\right]dV+\underset{{\mathrm{\Gamma }}_d}{\int }{\mathit{T}}_{\mathit{\epsilon }}dA,\end{array}$$

(24)

where it has been recognised that ${\mathbf{H}}_s^t$ remains a stress projection operator over the crack surface having $\widehat{\mathit{n}}$ as its normal, thus returning traction vectors at once operating over both stress fields. The main intent of assuming ${\mathbf{G}}_{su}^{*}$ is to effectively reveal the traction vectors $\mathit{T}$ and ${\mathit{T}}_{\mathit{\epsilon }}$. Depending on the stress relation coming from Equation (20a), Equation (24) might be simplified. A strong equality between $\mathit{T}$ and ${\mathit{T}}_{\mathit{\epsilon }}$ is possible only if we can ensure that the respective stress fields 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 Equation (22) allows to define ${\mathbf{G}}_s^{*}$ under the EAS approach. The arbitrary vector $\mathit{s}$ can be then taken out for reaching the following weak orthogonality condition:

$$\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_s^{*t}\mathbf{S}dV=0$$

(25)

As the form $\mathbf{S}$ is chosen based on real physics, it is clear that ${\mathbf{G}}_s^{*}$ must be determined by it. The works on Equation (24) have already required the unbounded part ${\mathbf{G}}_{su}^{*}$ to be equal to ${\delta }_{\mathrm{\Gamma }}{\mathbf{H}}_s^{*}$ in order to obtain the traction vectors, so only the bounded part ${\mathbf{G}}_{sb}^{*}$ can be used to satisfy Equation (25):

$$\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{sb}^{*t}\mathbf{S}dV+{\mathbf{H}}_s^t\underset{{\mathrm{\Gamma }}_d}{\int }\mathbf{S}dA=0,$$

(26)

where the final form of ${\mathbf{G}}_{sb}^{*}$ will basically depend on the choice of $\mathbf{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 fields regardless of its underlying complexity [1,14,20]. This implies to additionally satisfy Equation (26) for the case of a constant $\mathbf{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 $\mathbf{S}$ from the beginning, so that only one single solution for Equation (26) is required, satisfying the patch test for once and for all. For that case indeed, the solution to Equation (26) simplifies to the expression commonly found in E-FEM literature for the SKON, EAS approach [21,22]:

$${\mathbf{G}}_{sb}^{*}=-\frac{A_{\mathrm{\Gamma }}}{V_e}{\mathbf{H}}_s,$$

(27)

for which ${\mathbf{G}}_{sb}^{*}$ has also been set as a constant matrix, as well. For other cases that are less restrictive, ${\mathbf{G}}_{sb}^{*}$ has to rise in complexity by simultaneously satisfying Equation (26) for both a constant and a higher degree $\mathbf{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 coefficients for ${\mathbf{G}}_{sb}^{*}$ similar to that of Equation (19).

Note that, no matter what the choice is, it inevitably leads to ${\mathbf{G}}_{sb}\ne {\mathbf{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 in the next sections, we will retain the choice of a constant real stress field, along with a constant $\mathbf{S}$ and the validity of Equation (27).

3.2.5. Final Traction Calculation

If a constant ${\mathbf{G}}_{sb}^{*}$ is retained, Equations (21) and (24) allow further simplifications to reach a weak equality between traction vectors:

$${\overline{\mathit{T}}}_{\epsilon }=\overline{\mathit{T}}=\mathit{T},$$

(28)

where we have used the fact that a constant $\sigma$ yields a constant $\mathit{T}$. Having cleared this link, the right hand of Equation (22) finally leads a way to explicitly calculate a surface average of the traction vector $\overline{\mathit{T}}$, for then applying any pertinent damage, softening or traction–separation laws. The calculation is simply stated as:

$$\mathit{T}=\frac{1}{A_{\mathrm{\Gamma }}}\underset{{\mathrm{\Gamma }}_d}{\int }{\mathit{T}}_{\mathit{\epsilon }}dA=-\frac{1}{A_{\mathrm{\Gamma }}}\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{sb}^{*t}\mathit{\sigma }\left(\epsilon \right)dV$$

(29)

From now on, traction vector calculations will be simply referring to $\mathit{T}$.

3.2.6. Internal–External Force Balance

This part of the variational analysis has been left to last as it does not bring any remarkable implications. Equation (7a) can be stated as follows:

$$\underset{\mathsf{\Omega }}{\int }{\mathbf{B}}^t\sigma dV=\underset{\mathsf{\Omega }}{\int }{\mathbf{N}}^t{\mathit{f}}_{b}dV+\underset{\mathsf{\Omega }}{\int }{\mathbf{N}}^t\mathit{t}dA,$$

(30)

where Equation (20a) can be used literally as it is to vanish the real stress field from the expression and leave all force balance calculations in terms of the calculated stress field $\sigma \left(\mathit{\epsilon }\right)$:

$$\underset{\mathsf{\Omega }}{\int }{\mathbf{B}}^t\sigma \left(\mathit{\epsilon }\right)dV=\underset{\mathsf{\Omega }}{\int }{\mathbf{N}}^t{\mathit{f}}_{b}dV+\underset{\mathsf{\Omega }}{\int }{\mathbf{N}}^t\mathit{t}dA$$

(31)

Note that no strong equality between the stress fields is required to satisfy Equation (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 finished. The present work will now take all the basics presented in Section 2 and Section 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 Equation (29) as a departure point. After developing ${\mathbf{G}}_{sb}^{*}$ and $\sigma \left(\epsilon \right)$, the following form can be reached:

$$\mathit{T}={\mathit{T}}_e+\mathbf{M}{\left[\left|\mathit{u}\right|\right]}_{\widehat{\mathit{n}}},$$

(32)

where ${\mathit{T}}_e$ is a traction vector depending only on the standard element displacement vector $\mathit{d}$, and $\mathbf{M}$ is a crack stiffness matrix. The vector ${\left[\left|\mathit{u}\right|\right]}_{\widehat{\mathit{n}}}$ is the crack displacement vector expressed on the local frame $\widehat{\mathit{n}}$.

Equation (32), which is expressed entirely on the local crack base ($\widehat{\mathit{n}},\widehat{\mathit{t}},\widehat{\mathit{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 influence of the load-driven and the crack separation-driven parts on the overall evolution of the fracture traction $\mathit{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 Equation (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 different ways to handle Equation (32) depending on the general fracture modelling approach.

4.1. Single Mode Formulations

Single mode formulations are those that consider a single kinematic mode associated with the fracture within an element, normally a rigid body normal separation or a parallel sliding distance (the latter assuming a uniform fracture plane ${\mathrm{\Gamma }}_d$). This implies the suppression of all other fracture kinematic modes. For illustrating the implications of this choice, Equation (32) can be taken under terminal separation conditions: that is, assuming all averaged traction components are driven down to zero after the crack has sufficiently developed. We present the resulting equation as a system taking each of the local fracture plane directions:

$$\begin{array}{cc}\hfill & T_n=T_{en}+M_{nn}{\left[\left|u\right|\right]}_n+M_{nt}{\left[\left|u\right|\right]}_t+M_{nm}{\left[\left|u\right|\right]}_t=0\hfill \\ \hfill & T_t=T_{et}+M_{tn}{\left[\left|u\right|\right]}_n+M_{tt}{\left[\left|u\right|\right]}_t+M_{tm}{\left[\left|u\right|\right]}_t=0\hfill \\ \hfill & T_m=T_{em}+M_{mn}{\left[\left|u\right|\right]}_n+M_{mt}{\left[\left|u\right|\right]}_t+M_{mm}{\left[\left|u\right|\right]}_t=0\hfill \end{array}$$

(33)

As an example, considering a single normal separation mode formulation, the model will be entirely based on ${\left[\left|u\right|\right]}_n$ and will automatically make ${\left[\left|u\right|\right]}_t={\left[\left|u\right|\right]}_m=0$. For the system given by Equation (33), this would imply:

$$\begin{array}{cc}\hfill & T_n=T_{en}+M_{nn}{\left[\left|u\right|\right]}_n=0\hfill \\ \hfill & T_t=T_{et}+M_{tn}{\left[\left|u\right|\right]}_n=0\hfill \\ \hfill & T_m=T_{em}+M_{mn}{\left[\left|u\right|\right]}_n=0\hfill \end{array}$$

(34)

The updated Equation (34) remains overconstrained and cannot be satisfied by a single value of ${\left[\left|u\right|\right]}_n$. Authors choosing this approach will only take the traction–separation equation corresponding to the chosen fracture mode to be satisfied. In this case, only the first equation in (34) depicting the general traction component $T_n$ will be driven to zero. A specific ${\left[\left|u\right|\right]}_n$ value will ensure this condition. This ${\left[\left|u\right|\right]}_n$ will evidently not satisfy the remaining traction–separation relations in the system. As the components $T_{et},T_{em}$ depend on the arbitrary disposition of the nodal displacement vector $\mathit{d}$, general traction components $T_t,T_t$ will 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 to deliberately drive these crack surface traction components to zero during a numerical solution, along with any other stresses still associated with the crack. This can be achieved 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 specific combination of load (through $\mathit{d}$) and element geometry orientation with respect to ${\mathrm{\Gamma }}_d$ (coefficients $M_{nn},M_{tn},M_{mn}$). As it will be shown later in a more detailed analysis of $\mathbf{M}$ for a linear tetrahedron and a normal separation mode, these conditions correspond to purely axial strain on the $\widehat{\mathit{n}}$ direction, having an element face parallel to a planar ${\mathrm{\Gamma }}_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 $\left[\left|\mathit{u}\right|\right]$ becomes a single algebraic equation. As simple as it is, this kind of formulation has been successfully used for the modelling 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 confinement preloading conditions [25].

4.2. Full Crack Translation Formulations: The Role of $\phi$

Full crack translation formulations consider all crack displacement components when articulating the traction–separation equation system. In terminal conditions, the full system depicted by Equation (33) is to be solved for all variables ${\left[\left|u\right|\right]}_n,{\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m$. Regardless of the nature of the traction–separation decay (exponential, polynomial, etc.), Equation (33) remains strictly linear. A unique solution will always exist for $\left[\left|\mathit{u}\right|\right]$ given that $\mathbf{M}$ is not singular. A detailed analysis of $\mathbf{M}$ is thus relevant in this study.

The crack stiffness matrix $\mathbf{M}$ has a direct dependence on ${\mathbf{G}}_{sb}^{*}$ and $\sigma \left(\mathit{\epsilon }\right)$ (Equation (29)). $\sigma \left(\mathit{\epsilon }\right)$ in turn depends directly on ${\mathbf{G}}_{sb}$ (Equation (14)), and thus on $\phi$ (Equation (16)). As already discussed in Section 2.1, the nature of $\phi$ depends on the parent element, and authors typically take the simplest mathematical definition of it satisfying the basic requirements in Equation (3).

The case of a linear tetrahedron, a parent element will be taken again as an illustrating reference. The simplest $\phi$ 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 Figure 4. Note that it is not required to use the same function $\phi$ for all three displacement components ${\left[\left|u\right|\right]}_n,{\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m$. Note that in this case there is only one possible definition of a linear $\phi$ for this parent element, anyway.

The standard deformation matrix $\mathbf{B}$ can be expressed as a block partition in a node-wise manner:

$$\mathbf{B}=\left[\begin{array}{cccc}{\mathbf{B}}_1& {\mathbf{B}}_2& {\mathbf{B}}_3& {\mathbf{B}}_4\end{array}\right],$$

(35)

where each block ${\mathbf{B}}_i$ can be further expressed as a function of the unit vector ${\mathit{\Psi }}_i={[{\Psi }_{ix},{\Psi }_{iy},{\Psi }_{iz}]}^t$ normal to the opposing face to each node i in the tetrahedron (this is naturally coming from the basic FEM formulation of this parent element):

$${\mathbf{B}}_i=\frac{A_i}{3V_e}\left[\begin{array}{ccc}{\Psi }_{ix}& 0& 0\\ 0& {\Psi }_{iy}& 0\\ 0& 0& {\Psi }_{iz}\\ {\Psi }_{iy}& {\Psi }_{ix}& 0\\ 0& {\Psi }_{iz}& {\Psi }_{iy}\\ {\Psi }_{iz}& 0& {\Psi }_{ix}\end{array}\right],$$

(36)

where $A_i$ is the area of the opposite face to node i. Knowing that $\phi$ follows a similar structure as $\mathbf{B}$ but considering the selective condition in Equation (3), ${\mathbf{G}}_{sb}$ can be expressed as a superposition of the blocks defined in Equation (36):

$${\mathbf{G}}_{sb}=-\sum _i^{N_e}{p}_i{\mathbf{B}}_i,$$

(37)

where $p_i$ adopts a value of 1 on the nodes lying on ${\mathsf{\Omega }}^{+}$, and zero otherwise.

On the other hand, a constant field $\sigma$ will be assumed so that the EAS approach matrix ${\mathbf{G}}_{sb}^{*}$ can be obtained using Equation (27). With Equation (14) in mind, Equation (29) can now be developed to obtain the following:

$$\begin{array}{c}\underset{{\mathbf{K}}_e}{\underbrace{\left(\frac{1}{V_e}\underset{\mathsf{\Omega }}{\int }{\mathbf{H}}_s^T\mathbf{C}\mathbf{B}dV\right)}}\mathit{d}+\underset{{\mathbf{K}}_m}{\underbrace{\left(\frac{1}{V_e}\underset{\mathsf{\Omega }}{\int }{\mathbf{H}}_s^T\mathbf{C}{\mathbf{G}}_{sb}dV\right)}}\left[\left|\mathit{u}\right|\right]\end{array}$$

(38)

$$\begin{array}{c}\mathit{T}={\mathbf{R}}^T{\mathbf{K}}_e\mathit{d}+{\mathbf{R}}^T{\mathbf{K}}_m\mathbf{R}{\left[\left|\mathit{u}\right|\right]}_{\widehat{\mathit{n}}}={\mathit{T}}_e+\mathbf{M}{\left[\left|\mathit{u}\right|\right]}_{\widehat{\mathit{n}}},\end{array}$$

(39)

where $\mathbf{R}$ is a rotation matrix allowing the passage from a global coordinate crack displacement vector $\left[\left|\mathit{u}\right|\right]$ to a local-based ${\left[\left|\mathit{u}\right|\right]}_{\widehat{\mathit{n}}}$. Working an explicit form for $\mathbf{M}$ yields finally:

$$\mathbf{M}=-\frac{1}{V_e}\left[\begin{array}{ccc}{c}_{n1}\widehat{\mathit{n}}·\mathit{\Psi }& c_{n2}\widehat{\mathit{t}}·\mathit{\Psi }& c_{n2}\widehat{\mathit{m}}·\mathit{\Psi }\\ c_s\widehat{\mathit{t}}·\mathit{\Psi }& c_s\widehat{\mathit{n}}·\mathit{\Psi }& 0\\ c_s\widehat{\mathit{m}}·\mathit{\Psi }& 0& c_s\widehat{\mathit{n}}·\mathit{\Psi }\end{array}\right],$$

(40)

where · denotes a dot product ($\widehat{\mathit{t}}·\mathit{\Psi }={\widehat{\mathit{t}}}^T\mathit{\Psi }$) and $\mathit{\Psi }=A_i\sum {\mathit{\Psi }}_i$ follows the superposition brought by Equation (37). The constants $c_{n1},c_{n2},c_s$ are a function of the linear elastic material model parameters $E,\nu$. The reader can find more details about the derivation of this compact expression for $\mathbf{M}$ in Appendix A.

For now, let us assume a fracture plane that makes a 3-1 partition of the tetrahedron’s nodes, having a single node on ${\mathsf{\Omega }}^{+}$ (Figure 6), so that the interpretation of $\mathit{\Psi }$ 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 $\widehat{\mathit{n}}$, thus $\mathit{\Psi }·\widehat{\mathit{n}}$ will be always different than zero. Given this fact and the structure given by Equation (40), it can be safely assumed that $\mathbf{M}$ will not be singular. The system defined in (33) will then have a unique solution for ${\left[\left|u\right|\right]}_n,{\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m$.

It is relevant to note the presence of coupling fracture stiffness coefficients between different directions ($M_{ij},i\ne j$). This means that, even at terminal conditions (Equation (33)), a normal displacement ${\left[\left|u\right|\right]}_n$ of the fractured piece of the element ${\mathsf{\Omega }}^{+}$ will automatically induce a parallel rigid body motion ${\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m$ as well. This does not make physical sense since at this stage ${\mathsf{\Omega }}^{+}$ 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 ${\left[\left|u\right|\right]}_n,{\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m$. Moving in a normal direction would mean to produce traction on the parallel direction. Again, this is not representative of fractured body kinematics.

Equation (40) also states that all coefficients $M_{ij}$ will depend in general on the orientation of the parent element’s geometry with respect to the fracture plane. That is, the $M_{ij}$ coefficients have an explicit mesh dependency. Thus, if we were to avoid any unwanted kinematic/traction coupling, the element geometry must have a specific orientation with respect to the fracture plane. Otherwise, while a mathematical solution for the system given by Equation (33) will always be assured, the overall physical sense of $\left[\left|\mathit{u}\right|\right]$ will be once again compromised. Wells [9] had already identified the influence of the function $\phi$ on the overall fracture kinematics and the induced strong coupling between traction components at the fracture plane. He proposed to fix it by redefining $\phi$ in such a way that the $\mathbf{M}$ components having the dot products $\mathit{\Psi }·\widehat{\mathit{t}}$ and $\mathit{\Psi }·\widehat{\mathit{m}}$ were driven to zero. That is, a $\phi$ capable of diagonalising $\mathbf{M}$. Illustrating the process in two dimensions, if an initial function $\phi$ was defined for a CST element on the $xy$ plane having partial derivatives ${\phi }_{,x}$ and ${\phi }_{,y}$, the following redefinition for the partial derivatives was made:

$$\begin{array}{c}{\phi }_{,x}^{\prime }=n_x^2{\phi }_{,x}+n_x{n}_y{\phi }_{,y}\end{array}$$

(41a)

$$\begin{array}{c}{\phi }_{,y}^{\prime }=n_y^2{\phi }_{,y}+n_x{n}_y{\phi }_{,x}\end{array}$$

(41b)

It is not necessary to check after $\phi$ directly, as only the partial derivatives of $\phi$ are needed during the numerical calculation.

Indeed, such a redefinition of $\phi$ will avoid the unwanted coupling of traction components. However, this $\phi$ will no longer be able to satisfy the basic requirements of Equation (3). Remember that this was the main reason for conceiving $\phi$ in the first place. If this alternate definition of $\phi$ is checked on the position ${\mathit{x}}_1$ shown in Figure 7 where it is supposed to be zero, it can be shown that:

$$\phi \left({\mathit{x}}_1\right)\propto \left|\left(\widehat{\mathit{n}}\times {\mathit{\Psi }}_{12}\right)\right|\left|\left({\mathit{x}}_1\times \widehat{\mathit{n}}\right)\right|,$$

(42)

where × denotes a cross product, ${\mathit{\Psi }}_{12}$ is the normal vector to the edge between nodes 1 and 2 and ${\mathit{x}}_1$ is the global position vector for node 1. This product of magnitudes nullifies whether the fracture line is aligned with the normal of the opposed element edge or if it is also aligned with the position vector in question. The reader can find more details in Appendix B.

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 in the first half to describing its theoretical basis.

4.3. General Crack Kinematics Formulations

At terminal separation conditions, there should be no internal forces in the element associated with any rigid body motion of ${\mathsf{\Omega }}^{+}$ with respect to ${\mathsf{\Omega }}^{-}$. Otherwise, the 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 $\left[\left|\mathit{u}\right|\right]$ ensures that it is possible to drive traction components $\overline{T_n},\overline{T_m},\overline{T_n}$ down 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 in Equation (31):

$$f_{int}=\underset{\mathsf{\Omega }}{\int }{\mathbf{B}}^t\mathit{\sigma }\left(\mathit{\epsilon }\right)dV=\underset{\mathsf{\Omega }}{\int }{\mathbf{B}}^t\mathit{\sigma }dV=V_e{\mathbf{B}}^t\mathit{s},$$

(43)

where Equation (20a) and the fact that the present work stands for a constant stress field ($\sigma =\mathbb{S}\mathit{s}=\mathit{s}$) have been used. On the other hand, the kinematics proposed in Section 4.2 will only ensure the average traction vector $\mathit{T}$ acting upon ${\mathrm{\Gamma }}_d$ is driven down to zero:

$$\mathit{T}=-\frac{1}{A_{\mathrm{\Gamma }}}\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{sb}^{*t}\mathit{\sigma }\left(\mathit{\epsilon }\right)dV=\mathbf{0}$$

(44)

For a constant stress parent element (constant ${\mathbf{G}}_{sb}^{*},\mathbf{B}$), it is obvious to see that Equation (44) does not imply reaching zero in Equation (43) since, in general, we have ${\mathbf{G}}_{sb}^{*}\propto /\mathbf{B}$. In this case, the whole real stress tensor on the element must be nullified so that all internal forces come to zero. Thanks to Equation (44), we also know that some of the components of this stress tensor expressed on the local coordinate base ($\widehat{\mathit{n}},\widehat{\mathit{t}},\widehat{\mathit{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:

$$\mathit{\sigma }:\left[\begin{array}{ccc}{\sigma }_{nn}=T_n\to 0& {\sigma }_{nt}=T_t\to 0& {\sigma }_{nm}=T_m\to 0\\ {\sigma }_{tn}=T_t\to 0& {\sigma }_{tt}& {\sigma }_{tm}\\ {\sigma }_{mn}=T_m\to 0& {\sigma }_{mt}& {\sigma }_{mm}\end{array}\right]$$

(45)

These untouched components ${\sigma }_{tt},{\sigma }_{tm},{\sigma }_{tt}$ (only three since the tensor is symmetric) should be damaged in such a way that the framework retains mathematical consistency and physical meaningfulness.

In order to damage $\sigma$, Equation (21) can be taken as a point of departure to relate a damaged state of stress $\sigma$ to the constitutive stress $\sigma \left(\mathit{\epsilon }\right)$. Then, Equation (14) is used to see the implications on specific terms within the formulation. At terminal conditions, nullifying $\sigma$ will thus imply:

$$\begin{array}{c}\underset{\mathsf{\Omega }}{\int }\mathit{\sigma }dV=\frac{1}{V_e}\underset{\mathsf{\Omega }}{\int }\mathit{\sigma }\left(\mathit{\epsilon }\right)dV=\frac{1}{V_e}\underset{\mathsf{\Omega }}{\int }\mathbf{C}\left({\mathbf{G}}_{sb}\left[\left|\mathit{u}\right|\right]+\mathbf{B}\mathit{d}\right)dV=0\end{array}$$

(46)

$$\begin{array}{c}\frac{1}{V_e}\left(\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{sb}dV\right)\left[\left|\mathit{u}\right|\right]=-\frac{1}{V_e}\left(\underset{\mathsf{\Omega }}{\int }\mathbf{B}dV\right)\mathit{d}\end{array}$$

(47)

$$\overline{{\mathbf{G}}_{sb}}\left[\left|\mathit{u}\right|\right]=-\overline{\mathbf{B}}\mathit{d}$$

(48)

In the end, Equation (48) makes a relation between two average strain fields: one coming from standard displacements and the other coming from the strong discontinuity enrichment function. For the case of a linear tetrahedron, the average standard strain field on the right side of Equation (48) maps a ${\mathbb{R}}^{12}$ vector to a ${\mathbb{R}}^6$ subspace ${\mathcal{S}}_{\mathit{d}}$. We also have $\overline{\mathbf{B}}=\mathbf{B}$, as $\mathbf{B}$ is constant for this case. On the other hand, the left side enriched average strain maps a ${\mathbb{R}}^3$ vector to a ${\mathbb{R}}^6$ subspace ${\mathcal{S}}_{\left[\left|\mathit{u}\right|\right]}$. If ${\mathbf{G}}_{sb}$ happens to be constant, then we would also have ${\mathbf{G}}_{sb}=\overline{{\mathbf{G}}_{sb}}$. Despite the fact that we can stretch more demands to the matrix ${\mathbf{G}}_{sb}$ through a richer definition of $\phi$, it can be shown that, in general, ${\mathcal{S}}_{\left[\left|\mathit{u}\right|\right]}⫆/{\mathcal{S}}_{\mathit{d}}$. Figure 8 illustrates this situation. It is impossible to nullify the strain subspace spanned by $\mathit{d}$ only counting on a three rigid body displacement mode vector $\left[\left|\mathit{u}\right|\right]$. Fracture kinematics have to 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.

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 fields 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 $\left[\left|\mathit{u}\right|\right]$ are rather defined using linear functions instead of constant values. Each of the parameters within these linear functions is associated with general kinematics for the split body representing the ${\mathsf{\Omega }}^{+}$ domain. This can include rigid body translations, rigid body rotations and even simple axial strains in specific directions. The centroid of the ${\mathrm{\Gamma }}_d$ is generally taken as the reference position for the rigid body modes. In essence, fracture kinematics is enriched by providing ${\mathsf{\Omega }}^{+}$ with improved modes that are able to represent more than just rigid body translations ${\left[\left|u\right|\right]}_n,{\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m$. For a 2D case, Figure 9 illustrates how ${\mathsf{\Omega }}^{+}$ 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 ${\mathrm{\Gamma }}_d$. This is the type of enrichment studied by Armero and Linder [11]. For the sake 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 redefinition for $\left[\left|\mathit{u}\right|\right]$ that introduces a new intermediary matrix $\mathbf{J}$ along with a more general fracture kinematics vector $\xi$:

$$\begin{array}{c}\left[\left|\mathit{u}\right|\right]=\mathbf{J}\xi \end{array}$$

(49)

$$\begin{array}{c}\mathit{\xi }={\left[\begin{array}{cccc}{\left[\left|u\right|\right]}_{n0}& {\left[\left|u\right|\right]}_{t0}& \theta & {ϵ}_t\end{array}\right]}^t\end{array}$$

(50)

This enrichment has consequences on the mathematical framework developed in Section 2 and Section 3, starting by the basic strain expression derived in Equation (5), where the handling of gradients will change due to $\mathbf{J}$. Additional interactions will also emerge between $\phi$ and this new matrix operator, yet giving rise to a different structure for the strain enrichment operators. Equations (47) and (48) will now be:

$$\begin{array}{c}\frac{1}{V_e}\left(\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{sb}\mathbf{J}dV\right)\xi =\frac{1}{V_e}\left(\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{sb}^{\prime }dV\right)\xi =-\frac{1}{V_e}\left(\underset{\mathsf{\Omega }}{\int }\mathbf{B}dV\right)\mathit{d}\end{array}$$

(51)

$$\begin{array}{c}\overline{{\mathbf{G}}_{sb}^{\prime }}\xi =-\overline{\mathbf{B}}\mathit{d},\end{array}$$

(52)

where a new compound operator ${\mathbf{G}}_{sb}^{\prime }$ along with its averaged counterpart $\overline{{\mathbf{G}}_{sb}^{\prime }}$ have been defined.

In order to build a new coherent structure for ${\mathbf{G}}_{sb}^{\prime }$ under this scheme, authors typically choose not to rework Equation (5) and prefer to build $\overline{{\mathbf{G}}_{sb}^{\prime }}$ directly in a numerical, column by column fashion [11,27]. Each of the columns of $\overline{{\mathbf{G}}_{sb}^{\prime }}$ represent a link between a fracture kinematic mode ${\xi }_k$ and a specific nodal displacement state ${\mathit{d}}_k$:

$$\overline{{\mathbf{G}}_{s{b}_k}^{\prime }}{\xi }_k=-\mathbf{B}{\mathit{d}}_k$$

(53)

Here, the matrix $\overline{{\mathbf{G}}_{sb}^{\prime }}$ is given 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 $\overline{{\mathbf{G}}_{sb}^{\prime }}$, this approach is enough to define the internal numerical solution process.

For example, suppose we would like to link the kinematic variable ${\left[\left|u\right|\right]}_{n0}$ to an actual set of nodal displacements ${\mathit{d}}_{{\left[\left|u\right|\right]}_{n0}}$ on the element. To do so, a real ${\left[\left|u\right|\right]}_{n0}$ is meant to represent a constant displacement of ${\left[\left|u\right|\right]}_{n0}$ for nodes $N_1,N_2$ lying on ${\mathsf{\Omega }}^{+}$ in the $\widehat{\mathit{n}}$ direction. This specific displacement state ${\mathit{d}}_{{\left[\left|u\right|\right]}_{n0}}$ will generate a standard strain through $\mathbf{B}=\overline{\mathbf{B}}$. The negative of this strain has to be equal to the product of a specific column of $\overline{{\mathbf{G}}_{sb}^{\prime }}$ (the first one) and the kinematic mode ${\left[\left|u\right|\right]}_{n0}$. The block structure of $\mathbf{B}=\left[\begin{array}{ccc}{\mathbf{B}}_1& {\mathbf{B}}_2& {\mathbf{B}}_3\end{array}\right]$ can be used to our advantage:

$$\begin{array}{c}{\overline{{\mathbf{G}}_{sb}^{\prime }}}_{{\left[\left|u\right|\right]}_{n0}}{\left[\left|u\right|\right]}_{n0}=-\mathbf{B}{\mathit{d}}_{{\left[\left|u\right|\right]}_{n0}}\\ {\overline{{\mathbf{G}}_{sb}^{\prime }}}_{{\left[\left|u\right|\right]}_{n0}}{\left[\left|u\right|\right]}_{n0}=-\left[\begin{array}{ccc}{\mathbf{B}}_1& {\mathbf{B}}_2& {\mathbf{B}}_3\end{array}\right]\left[\begin{array}{c}{\left[\left|u\right|\right]}_{n0}\widehat{\mathit{n}}\\ {\left[\left|u\right|\right]}_{n0}\widehat{\mathit{n}}\\ 0\end{array}\right]=-\sum _i^{N_e}{p}_i\left({\mathbf{B}}_i\widehat{\mathit{n}}\right){\left[\left|u\right|\right]}_{n0}\end{array}$$

$$\begin{array}{c}{\overline{{\mathbf{G}}_{sb}^{\prime }}}_{{\left[\left|u\right|\right]}_{n0}}=-\sum _i^{N_e}{p}_i\left({\mathbf{B}}_i\widehat{\mathit{n}}\right)\end{array}$$

(54)

The remaining columns ${\overline{{\mathbf{G}}_{sb}^{\prime }}}_{{\left[\left|u\right|\right]}_{t0}},{\overline{{\mathbf{G}}_{sb}^{\prime }}}_{\theta },{\overline{{\mathbf{G}}_{sb}^{\prime }}}_{{ϵ}_t}$ can be deducted following the same logic. This process allows to numerically define $\overline{{\mathbf{G}}_{sb}^{\prime }}$ completely based on basic parent element characteristics.

One remark to this process is the nature of the link made between nodal displacements and kinematic variables (Equation (53)). A given mapping from the ${\mathbb{R}}^{12}$ displacement space to the ${\mathbb{R}}^6$ strain space is not unique, meaning that two or more different displacement states ${\mathit{d}}_k$ might correspond to exactly the same standard strain vector. This means that, even if we create a link between a specific fracture kinematic variable ${\xi }_k$ to a desired displacement state through a given ${\mathit{d}}_k$, the strain produced by ${\mathit{d}}_k$ could also be produced by other entirely different ${\mathit{d}}_k^{\prime }$, breaking the meaning of this link, and thus the physical meaning of ${\xi }_k$.

This is especially true when trying to distinguish what happens on ${\mathsf{\Omega }}^{+}$ or ${\mathsf{\Omega }}^{-}$. There might be cases in which nodes lying on ${\mathsf{\Omega }}^{-}$ will contribute to a vector $\mathit{d}$ capable of exciting unexpected modes ${\xi }_k$ that were originally associated with the kinematics of ${\mathsf{\Omega }}^{+}$. Rigid body modes are reflexive in the sense that they can be interpreted by fixing a reference frame whether at ${\mathsf{\Omega }}^{-}$ or ${\mathsf{\Omega }}^{+}$, 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 ${\mathsf{\Omega }}^{+}$ and ${\mathsf{\Omega }}^{-}$ that prevents the approach from being fully objective. There is some work that has been undertaken, as an example, to consider a more general definition for the factor $p_i$: make it a continuous variable from 0 to 1 to control the coupling between ${\mathsf{\Omega }}^{+}$ and ${\mathsf{\Omega }}^{-}$ in a more flexible manner [29]. As this does not really solve the root problem, the authors of the present work have decided to keep the classical binary definition for $p_i$.

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 $\overline{{\mathbf{G}}_{sb}^{\prime }}$ yields a set of vectors that compose a linearly independent base for building a $\mathit{\epsilon }$ space. For some node partitions on ${\mathsf{\Omega }}^{+},{\mathsf{\Omega }}^{-}$, specifically when we have a single node on ${\mathsf{\Omega }}^{+}$, 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 ${\xi }_k$ modes for a given element kinematic state. Let us take the same CST triangle as an example, for a 1-2 node domain partition having node 1 on ${\mathsf{\Omega }}^{+}$. For a given displacement ${\mathit{d}}_1$ of this node, there are infinite ways to make a ${\left[\left|u\right|\right]}_{n0},{\left[\left|u\right|\right]}_{t0},\theta ,{ϵ}_t$ distribution in such a way that the resulting ${\mathsf{\Omega }}^{+}$ motion returns the same ${\mathit{d}}_1$ displacement on node 1 (See Figure 10). Linder and Armero [11] proposed 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 definition of any ${\mathbf{G}}_{sb}$-related operators without making a direct assessment of $\phi$. None of the previous relations actually ensures satisfaction of any of the main requirements for $\phi$ in Equation (3). Linder and Armero [11] mention the possibility of integrating $\overline{{\mathbf{G}}_{sb}^{\prime }}$ to inquire about $\phi$, but this process is far from trivial given the multiple gradients implied in Equation (5). A lot of information about the specific $\phi$ giving rise to the given $\overline{{\mathbf{G}}_{sb}^{\prime }}$ has been already lost at this point. The only evident fact would be that $\phi$ 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 Section 2 and Section 3 as well as the issues found during the evolution of different strong discontinuity modelling 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 do its best to reach a sound compromise between a theoretically sound E-FEM formulation but still remain 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 $\phi$

In Section 4.2, the role of the function $\phi$ has been discussed. $\phi$ 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 ${\left[\left|u\right|\right]}_n,{\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m$. Wells [9] proposed a fix, but it is not completely effective due to a violation of the mathematical framework.

The theory described in Section 2 and Section 3 can be used to make a consistent $\phi$ proposal. The idea is to propose a $\phi$ with enough flexibility to both satisfy the requirements in Equation (3) and the nullification of the out-of-diagonal terms of the $\mathbf{M}$ matrix in Equation (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 definition for $\phi$ as long as it remains within the $C^0$ class of functions (continuous, derivative not necessarily continuous). It does not have to follow the same polynomial degree as the parent element. As an example, for a linear tetrahedron, a complete quadratic definition for $\phi$ is perfectly valid:

$$\begin{array}{c}\phi ={\mathit{P}}_2^T\alpha \\ \mathit{\alpha }=\left[\begin{array}{cccccccccc}{\alpha }_0& {\alpha }_1& {\alpha }_2& {\alpha }_3& {\alpha }_{12}& {\alpha }_{23}& {\alpha }_{13}& {\alpha }_{11}& {\alpha }_{22}& {\alpha }_{33}\end{array}\right]\\ {\mathit{P}}_2={\left[\begin{array}{cccccccccc}1& x& y& z& xy& yz& xz& x^2& y^2& z^2\end{array}\right]}^t\end{array}$$

(55)

This increase in order for $\phi$ implies an increase in order for the element strain field $\mathit{\epsilon }$ through ${\mathbf{G}}_{sb}$ (Equation (11)), and thus for the calculated stress field $\sigma \left(\mathit{\epsilon }\right)$. If the real stress field $\sigma$ is chosen to follow the natural order in the parent element (a constant field for a linear tetrahedron), the strong constitutive equality (Equation (8b)) will not be possible and stresses $\sigma$ will have to be calculated through volume averages (Equation (21)).

On the other hand, the advantage of having a discordance of $\sigma \left(\mathit{\epsilon }\right)$ with respect to a constant $\sigma$ field is that absolutely nothing will change in the treatment of the EAS part of the framework (Section 3.2.4). This means that Equation (27) can still be used for managing a constant expression for ${\mathbf{G}}_{sb}^{*}$.

The integrand in Equation (29) is not constant anymore. This means that the integration process for all traction calculations will have to consider the geometric parameters of ${\mathrm{\Gamma }}_d$ as well as the specific geometries of the subvolumes ${\mathsf{\Omega }}^{+}$ and ${\mathsf{\Omega }}^{-}$. Eventually, after vanishing all variable terms of the base ${\mathit{P}}_2$ through definite integration, an expression for the fracture stiffness matrix $\mathbf{M}$ as a function of $\phi$ coefficients $\alpha$ will be attained. As per Equation (40), we know that $\mathbf{M}$ remains almost symmetric, where the respective out-of-diagonal term pairs $(M_{ij},M_{ji})$ share the same internal products $\left({\widehat{\mathit{t}}}^T\mathit{\Psi }\right),\left({\widehat{\mathit{m}}}^T\mathit{\Psi }\right)$, and that there is already a pair $(M_{tm},M_{mt})$ that is zero. Thus, nullifying the out-of-diagonal terms in $\mathbf{M}$ for avoiding unwanted traction couplings implies the imposition of two additional linear equations on the $\alpha$ coefficients:

Adding this to the basic $\phi$ requirements (Equation (3), four linear equations on ${\mathit{k}}_{\phi }$) makes a total of 6 linear constraints. As a complete quadratic base ${\mathit{P}}_2$ allows for a total of 10 free parameters, it is completely feasible to design a $\phi$ function capable of satisfying BC imposition consistency and at the same time the mechanical requirements for traction behaviour.

While the relations ((56a) and (56b)) ensure diagonalisation of the fracture stiffness matrix $\mathbf{M}$, there is nothing that ensures that the magnitude of these diagonal components is physically meaningful. Each diagonal term $M_{jj}$ independently controls the fracture stiffness associated with each displacement ${\left[\left|u\right|\right]}_j$. These displacements should be consistent with the element standard nodal displacement vector $\mathit{d}$. As expected, $\phi$ can also help on this errand by imposing additional conditions for the diagonal fracture stiffness terms. Note that, as per Equation (40), the diagonal terms are not independent, and only one single constraint would be needed to regulate them all.

The case depicted in Figure 5 is taken as a reference for devising such constraints. For a 2-1 node partition on a 2D triangle, it is easy to find an exact relation between, for instance, the fracture normal separation ${\left[\left|u\right|\right]}_n$ and the normal displacement $\Delta d_n$ of a single isolated node partition on ${\mathsf{\Omega }}^{+}$: they should be equal. For other configurations, however, such as 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 so and none of them would be rigorously correct. This work proposes to take the average separation $\Delta \overline{\mathit{d}}$ between the groups of nodes on ${\mathsf{\Omega }}^{+}$ and ${\mathsf{\Omega }}^{-}$:

$$\begin{array}{c}\Delta \overline{\mathit{d}}=\overline{{\mathit{d}}^{+}}-\overline{{\mathit{d}}^{-}}\end{array}$$

(57a)

$$\begin{array}{c}\Delta \overline{\mathit{d}}=\frac{\sum p_i{\mathit{d}}_i}{\sum p_i}-\frac{\sum \left(1-p_i\right){\mathit{d}}_i}{\sum \left(1-p_i\right)},\end{array}$$

(57b)

Taking the normal direction $\widehat{\mathit{n}}$, we require:

$${\left[\left|u\right|\right]}_n=\Delta d_n={\widehat{\mathit{n}}}^T\Delta \overline{\mathit{d}}$$

(58)

At this point, it will be assumed that $\phi$ has been already arranged as to diagonalise $\mathbf{M}$ through Equations (56a) and (56b). With this, at terminal conditions, we can practically use the first equation of the system in Equation (33) to finally arrive to a linear constraint for $M_{nn}$:

$$\begin{array}{c}{T}_{en}+M_{nn}{\left[\left|u\right|\right]}_n=0\\ {\widehat{\mathit{n}}}^T{\mathbf{K}}_e\mathit{d}+M_{nn}\left({\widehat{\mathit{n}}}^T\Delta \overline{\mathit{d}}\right)=0\end{array}$$

$$M_{nn}\left(\mathit{\alpha }\right)=M_{n{n}_{\infty }}=-\frac{{\widehat{\mathit{n}}}^T{\mathbf{K}}_e\mathit{d}}{{\widehat{\mathit{n}}}^T\left[\frac{{\mathbf{F}}^t}{\sum p_i}-\frac{{\mathbf{I}}_{3\left(4\right)}^T-{\mathbf{F}}^T}{\sum \left(1-p_i\right)}\right]\mathit{d}}$$

(59)

A proper block matrix algebra has been used to express this constraint in terms of the complete vector $\mathit{d}$, using some auxiliary matrix blocks:

$$\begin{array}{c}{\mathbf{F}}^T=\left[\begin{array}{cccc}{p}_1{\mathbf{I}}_3& p_2{\mathbf{I}}_3& p_3{\mathbf{I}}_3& p_4{\mathbf{I}}_3\end{array}\right]\end{array}$$

(60a)

$$\begin{array}{c}{\mathbf{I}}_{3\left(4\right)}^T=\left[\begin{array}{cccc}{\mathbf{I}}_3& {\mathbf{I}}_3& {\mathbf{I}}_3& {\mathbf{I}}_3\end{array}\right],\end{array}$$

(60b)

having the $3\times 3$ identity matrix ${\mathbf{I}}_3$ as the basic block.

Note that in Equation (59), we are defining a component of $\mathbf{M}$ in terms of a nodal displacement vector $\mathit{d}$, which is the load imposed to the element. As the solution for all $\alpha$ coefficients 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 definition of this $M_{nn}$ stiffness, the load path would have to remain proportional. This is highly unlikely during a global fracture phenomenon, since there are different loading and relaxation stages within the same element that would make $\mathit{d}$ change direction in an aggressive fashion. $M_{nn}$ can be updated as well in such cases. This will imply more than one solution for the $\alpha$ coefficients. Imposing the constraint on a different component $M_{tt},M_{mm}$ will yield the same overall results given that a scaling factor of $c_s/c_{n1}$ is used as per Equation (40). As this represents only one additional linear constraint on the system for the $\alpha$ coefficients (making 7 constraints and 10 free parameters), once again it is still possible to build a proper $\phi$ function handling all requirements at once.

The next discussion will continue to explore the flexibility of $\phi$ under the light of fracture kinematics enrichment for 3D elements.

5.2. Fracture Kinematics Enrichment in 3D Coexisting with $\phi$

In Section 4.3, an introduction has been made to a more general definition of fracture kinematics. The approach has been partially illustrated for a 2D parent element following the work in [11], where an explicit $\phi$ is practically disregarded. This section will describe a detailed application for a 3D element, including the clear distinction and managing of the $\phi$ function.

The process starts with the same definition of $\left[\left|\mathit{u}\right|\right]$ as a linear field vector, but for three dimensions in the local frame:

$$\left[\left|\mathit{u}\right|\right]=\left[\begin{array}{c}{\left[\left|u\right|\right]}_n\\ {\left[\left|u\right|\right]}_t\\ {\left[\left|u\right|\right]}_m\end{array}\right]=\left[\begin{array}{c}{a}_n\eta +b_n\zeta +c_n\xi +d_n\\ a_t\eta +b_t\zeta +c_t\xi +d_t\\ a_m\eta +b_m\zeta +c_m\xi +d_m\end{array}\right],$$

(61)

where each of the parameters $a_j,b_j,c_j,d_j$ is to be related to different fracture kinematic modes ${\xi }_k$. It can be shown that such linear field structure allows for the modelling of three rigid body translations ${\left[\left|u\right|\right]}_{n_0}$, ${\left[\left|u\right|\right]}_{t_0}$, ${\left[\left|u\right|\right]}_{m_0}$, three rigid body rotations ${\theta }_n,{\theta }_t,{\theta }_m$ and three axial deformations (associated with ${\mathsf{\Omega }}^{+}$) ${ϵ}_n,{ϵ}_t,{ϵ}_m$.

Figure 11 illustrates how to build the relations between all 3D rotation modes and the parameters $a_j,b_j,c_j,d_j$. The general idea is to make a careful geometric interpretation for the meaning of each parameter 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 ${\xi }_k$. At the end, if small-angle approximations are allowed, the fields are finally expressed as:

$$\left[\begin{array}{c}{\left[\left|u\right|\right]}_n\\ {\left[\left|u\right|\right]}_t\\ {\left[\left|u\right|\right]}_m\end{array}\right]=\left[\begin{array}{c}-{\theta }_m\eta +{\theta }_t\zeta +{ϵ}_n\xi +{\left[\left|u\right|\right]}_{n_0}\\ {ϵ}_t\eta -{\theta }_n\zeta +{\theta }_m\xi +{\left[\left|u\right|\right]}_{t_0}\\ {\theta }_n\eta +{ϵ}_m\zeta -{\theta }_t\xi +{\left[\left|u\right|\right]}_{m_0}\end{array}\right]=\mathbf{J}\xi$$

(62)

From this point, an expression for the intermediary matrix $\mathbf{J}$ can be readily made having a definition for the vector $\mathit{\xi }$:

$$\mathbf{J}\mathit{\xi }=\left[\begin{array}{ccccccccc}1& 0& 0& 0& \zeta & -\eta & \xi & 0& 0\\ 0& 1& 0& -\zeta & 0& \xi & 0& \eta & 0\\ 0& 0& 1& \eta & -\xi & 0& 0& 0& \zeta \end{array}\right]\left[\begin{array}{c}{\left[\left|u\right|\right]}_{n_0}\\ {\left[\left|u\right|\right]}_{t_0}\\ {\left[\left|u\right|\right]}_{m_0}\\ {\theta }_n\\ {\theta }_t\\ {\theta }_m\\ {ϵ}_n\\ {ϵ}_t\\ {ϵ}_m\end{array}\right]$$

(63)

With this analysis, Equation (5) can be worked for a general definition for the strain field $\epsilon$ considering both $\xi$ and $\phi$. As all structures having to do with the enriched kinematics and $\mathbf{J}$ are in the local fracture surface frame $(\xi ,\eta ,\zeta )$, it makes sense to work all further developments on the local frame from now on. The same bounded–unbounded structure for $\mathit{\epsilon }$ can be derived:

$$\begin{array}{c}\mathit{\epsilon }=\mathbf{B}\mathit{d}+\left[{\mathbf{G}}_{sb}^{\prime }+{\mathbf{G}}_{su}^{\prime }\right]\xi \end{array}$$

(64a)

$$\begin{array}{c}{\mathbf{G}}_{sb}^{\prime }={\mathcal{H}}_{\mathrm{\Gamma }}\nabla \mathbf{J}-\left(\phi \nabla \mathbf{J}\right)-\nabla \phi \mathbf{J}\end{array}$$

(64b)

$$\begin{array}{c}{\mathbf{G}}_{su}^{\prime }={\delta }_{\mathrm{\Gamma }}{\mathbf{H}}_s\mathbf{J}\end{array}$$

(64c)

The term $\left(\phi \nabla \mathbf{J}\right)$ has been grouped because it is not the product of a deformation gradient operator on $\mathbf{J}$ and $\phi$, but rather a compound operator. The Heaviside function ${\mathcal{H}}_{\mathrm{\Gamma }}$ also forces a distinction between a bounded operator defined on ${\mathsf{\Omega }}^{+}$ (${\mathbf{G}}_{sb}^{\prime +}$) and on ${\mathsf{\Omega }}^{-}$ (${\mathbf{G}}_{sb}^{\prime -}$).

In general, (${\mathbf{G}}_{sb}^{\prime \pm }$) will be an explicit function of $\phi$ coefficients $\alpha$. To start the kinematic linking process between fracture modes and nodal displacements followed in Section 4.3, the requirement stated in Equation (52) can be worked on a subvolume-basis:

$$\overline{{\mathbf{G}}_{sb}^{\prime \pm }}\xi =\frac{1}{V^{\pm }}\underset{{\mathsf{\Omega }}^{\pm }}{\int }{\mathbf{G}}_{sb}^{\prime \pm }dV\xi =-\overline{\mathbf{B}}\mathit{d},$$

(65)

where the averaged operators $\overline{{\mathbf{G}}_{sb}^{\prime \pm }}$ will be explicitly linear matrices on the coefficients $\alpha$.

The linking of each kinematic mode will establish a series of linear equations on each of the column elements of $\overline{{\mathbf{G}}_{sb}^{\prime \pm }}$. For instance, the first mode ${\left[\left|u\right|\right]}_{n_0}$ already studied in Section 4.3 for two dimensions, would require satisfaction of exactly the same form of Equation (54). Fortunately, it can be shown that the first column of $\overline{{\mathbf{G}}_{sb}^{\prime \pm }}$ has exactly the same zeros as ${\mathbf{B}}_i$ so that in the local base, the following would be obtained:

$${\overline{{\mathbf{G}}_{sb}^{\prime \pm }}}_{{\left[\left|u\right|\right]}_{n0}}=\left[\begin{array}{c}{\overline{G_{sb}^{\prime \pm }}}_{1{\left[\left|u\right|\right]}_{n0}}\left(\alpha \right)\\ 0\\ 0\\ {\overline{G_{sb}^{\prime \pm }}}_{4{\left[\left|u\right|\right]}_{n0}}\left(\alpha \right)\\ 0\\ {\overline{G_{sb}^{\prime \pm }}}_{6{\left[\left|u\right|\right]}_{n0}}\left(\alpha \right)\end{array}\right]=-\sum _i^{N_e}{p}_i\left[\begin{array}{c}{B}_{ni}\\ 0\\ 0\\ B_{ti}\\ 0\\ B_{mi}\end{array}\right],$$

(66)

where the first index l in the ${\overline{G_{sb}^{\prime \pm }}}_{l{\left[\left|u\right|\right]}_{n0}}$ scalars 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 $\alpha$ coefficients. The specific expressions for ${\overline{G_{sb}^{\prime \pm }}}_{1{\left[\left|u\right|\right]}_{n0}},{\overline{G_{sb}^{\prime \pm }}}_{4{\left[\left|u\right|\right]}_{n0}},{\overline{G_{sb}^{\prime \pm }}}_{6{\left[\left|u\right|\right]}_{n0}}$ have to be retrieved by working $\phi$ and $\xi$ in Equation (64b) and then taking the volume averaged integral demanded in Equation (65). Note that $\phi$ is to be expressed as a function of local coordinates $(\xi ,\eta ,\zeta )$, and the subvolume coordinate description in the local frame will also be 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 first three columns of ${\mathbf{G}}_{sb}^{\prime \pm }$, which contain repetitions of the simple partial derivatives ${\phi }_{,\xi },{\phi }_{,\eta },{\phi }_{,\zeta }$. This depends on the symmetry of the kinematics defined and the model of the strong discontinuity overall. At the end, it is true that the demands for $\phi$ will do nothing but increase, as the $\phi$ free parameters are the only way to satisfy these relations. As a result, a full quadratic base proposal ${\mathit{P}}_2$ might not be enough to fulfil all these requirements additionally from Equation (3) and other constraints on the fracture stiffness matrix $\mathbf{M}$: nullification of non-diagonal terms (Sys. (56)) and terminal separation conditions (Equation (59)). While increasing the polynomial base order of $\phi$ might be 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 $\phi$ structure associated for all fracture displacement components ${\left[\left|u\right|\right]}_n,{\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m$. Indeed, three different functions ${\phi }_n,{\phi }_t,{\phi }_m$ may coexist in the model:

  $$\phi \left[\begin{array}{c}{\left[\left|u\right|\right]}_n\\ {\left[\left|u\right|\right]}_t\\ {\left[\left|u\right|\right]}_m\end{array}\right]\to \left[\begin{array}{ccc}{\phi }_n& 0& 0\\ 0& {\phi }_t& 0\\ 0& 0& {\phi }_m\end{array}\right]\left[\begin{array}{c}{\left[\left|u\right|\right]}_n\\ {\left[\left|u\right|\right]}_t\\ {\left[\left|u\right|\right]}_m\end{array}\right]$$

  (67)
  While this might apparently triple the amount of free parameters, some of the model symmetries that are allowed for redundancy simplifications will no longer emerge. As an example, the basic deformation gradient operator on $\phi$ (the first three columns of ${\mathbf{G}}_{sb}^{\prime \pm }$) will now have all different terms with respect to the one developed with a single $\phi$:

  $$\left[\begin{array}{ccc}{\phi }_{,\xi }& 0& 0\\ 0& {\phi }_{,\eta }& 0\\ 0& 0& {\phi }_{,\zeta }\\ {\phi }_{,\eta }& {\phi }_{,\xi }& 0\\ 0& {\phi }_{,\zeta }& {\phi }_{,\eta }\\ {\phi }_{,\zeta }& 0& {\phi }_{,\xi }\end{array}\right]\to \left[\begin{array}{ccc}{\phi }_{n,\xi }& 0& 0\\ 0& {\phi }_{t,\eta }& 0\\ 0& 0& {\phi }_{m,\varsigma }\\ {\phi }_{n,\eta }& {\phi }_{t,\xi }& 0\\ 0& {\phi }_{t,\varsigma }& {\phi }_{m,\eta }\\ {\phi }_{n,\varsigma }& 0& {\phi }_{m,\xi }\end{array}\right]$$

  (68)
  On the other hand, basic $\phi$ requirements in Equation (3) will now require the triple of linear relations for boundary condition consistency: one set for each ${\phi }_n,{\phi }_t,{\phi }_m$.

  $$\begin{array}{cc}\hfill {\phi }_n\left({\mathit{x}}_i\right)& =p_i\hfill \end{array}$$

  (69a)

  $$\begin{array}{cc}\hfill {\phi }_t\left({\mathit{x}}_i\right)& =p_i\hfill \end{array}$$

  (69b)

  $$\begin{array}{cc}\hfill {\phi }_m\left({\mathit{x}}_i\right)& =p_i\hfill \end{array}$$

  (69c)
  While there is evidently a trade-off, there is still a gain in the effective number of free parameters in the global $\phi$ structure, without the need to modify the complexity of the algebraic base of it.
• In general, $\phi$ has only a $C^0$ continuity 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 definition for $\phi$. The most natural choice is to propose a first function for ${\mathsf{\Omega }}^{+}$ and then a second one for ${\mathsf{\Omega }}^{-}$. Piece-wise $\phi$ does not increase the number of equations for Equation (3) requirements, but it still breaks some of the model symmetries. There will also be new linear equations to satisfy, which are associated with the basic $C^0$ continuity of $\phi$ at ${\mathrm{\Gamma }}_d$:

  $${\left.{\phi }^{-}\right|}_{{\mathrm{\Gamma }}_d}={\left.{\phi }^{+}\right|}_{{\mathrm{\Gamma }}_d},$$

  (70)

  which will depend on the nature of the base $\mathit{P}$ chosen for the structure of $\phi$. For a ${\mathit{P}}_2$ base, $C^0$ continuity requires six linear relations. Again, the gain in free parameters is worth the trade-off.

By considering the combination of both options, $\phi$ can globally add up to 120 free parameters using a complete ${\mathit{P}}_3$ base. While this might seem just too much, the E-FEM model managed in [24], which also features a weak discontinuity to represent different 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 definitively used for later element simulations in this work. The reader can find the detail on the mode linking expressions as well as general $\phi$ handling on Appendix C.

5.3. A Comment on Linear System Handling

Calculation of all $\phi$ parameters based on all kinematic considerations may represent a formidable linear system depending on the parent element needs. The more complex the overall definition for $\phi$ is (Appendix C), the harder it is to ensure the well-posedness of this linear system on the $\alpha$ coefficients. This is basically the price to pay when building a more robust kinematic model depending on $\phi$.

On the one hand, the linear dependency issues arising for some conditions during the construction of the columns of $\overline{{\mathbf{G}}_{sb}^{\prime \pm }}$ (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 find a unique solution for some groups of parameters).

On the other hand, for complex $\phi$ structures, many symmetries in the overall kinematic model are lost. As an example, the structure of the $\mathbf{M}$ matrix deduced in Equation (40) is not guaranteed for a piece-wise, direction dependent, ${\mathit{P}}_3$ definition for $\phi$. This implies that more $\alpha$ coefficients 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 $\alpha$ coefficients 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 if overconstrained or for the minimum value possible for free independent variables if underconstrained. The SVD decomposition is only made once the element reaches localisation. Once definite values are obtained for the $\alpha$ coefficients, the $\overline{{\mathbf{G}}_{sb}^{\prime \pm }}$ can 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 field $\delta \left[\left|\mathit{u}\right|\right]$ is discretised and the impact on the structure of Equation (25). While the $\mathit{real}$ fracture displacement field vector $\left[\left|\mathit{u}\right|\right]$ effectively requires a 9 variable interpolation matrix through $\phi$ and $\mathbf{J}$, the virtual counterpart $\delta \left[\left|\mathit{u}\right|\right]$does not have to share the same structure. $\delta \left[\left|\mathit{u}\right|\right]$ can still be assumed as a three constant field vector, or it can take the same detailed description as $\left[\left|\mathit{u}\right|\right]$. The Hu–Washizu framework allows this.

In this case, Equation (25) retains exactly the same structure having the real stress interpolation matrix $\mathbf{S}$ and the EAS operator ${\mathbf{G}}_{sb}^{*}$. Despite the fact that now we have a ${\mathbf{G}}_{sb}^{\prime \pm }$ considering the $\xi$ related upgrades, the EAS definition for ${\mathbf{G}}_{sb}^{*}$ does not have to consider any structures from the enriched fracture kinematics. The only requirement for the EAS definition of ${\mathbf{G}}_{sb}^{*}$ is just to have the same bounded-unbounded composition featuring a Dirac delta and a projection operator ${\mathbf{H}}_s$ as a minimum. This will allow the emergence and equality of the traction vector terms from Equation (24) onwards.

Equation (32) will return a new system with an enriched set of 9 variables in $\xi$:

$$\mathit{T}={\mathit{T}}_e+\mathbf{M}\xi ,$$

(71)

where ${\mathit{T}}_e$ remains exactly the same as per Equation (38). Note that if the work has already been performed on the local frame, there is no need to use the rotation matrix $\mathbf{R}$ anymore. The expression for base fracture stiffness matrix $\mathbf{M}$ will have an updated definition:

$$\begin{array}{c}\mathbf{M}=\left(\frac{1}{V_e}\underset{\mathsf{\Omega }}{\int }{\mathbf{H}}_s^T\mathbf{C}{\mathbf{G}}_{sb}^{\prime }dV\right)\xi =\frac{1}{V_e}{\mathbf{H}}_s^T\mathbf{C}\underset{\mathsf{\Omega }}{\int }{\mathbf{G}}_{sb}^{\prime \pm }dV\end{array}$$

$$\begin{array}{c}\mathbf{M}=\frac{1}{V_e}{\mathbf{H}}_s^T\mathbf{C}\left(V^{+}\overline{{\mathbf{G}}_{sb}^{\prime +}}+V^{-}\overline{{\mathbf{G}}_{sb}^{\prime -}}\right)\end{array}$$

(72)

Here, all matrices have been already assumed in the local frame to avoid further transformations, and the definitions implied in Equation (65) have been used. Note that matrix $\mathbf{M}$ now has dimensions $3\times 9$.

The original three traction–separation laws associated with traction components $T_n,T_t$, $T_m$ will return an underdetermined system:

$$\begin{array}{c}{T}_n=T_{en}+\sum _k{M}_{nk}{\xi }_k=q_n\left(\mathit{\xi }\right)\end{array}$$

(73a)

$$\begin{array}{c}{T}_t=T_{et}+\sum _k{M}_{tk}{\xi }_k=q_t\left(\mathit{\xi }\right)\end{array}$$

(73b)

$$\begin{array}{c}{T}_m=T_{em}+\sum _k{M}_{mk}{\xi }_k=q_m\left(\mathit{\xi }\right),\end{array}$$

(73c)

where $q_n,q_t,q_m$ are decaying expressions controlling the overall magnitude of $T_n,T_t,T_m$ until reaching fracture terminal conditions, where $T_n,T_t,T_m$ should be zero.

To fully determine the fracture mechanics, six additional relations are needed. Three equations can be proposed to explicitly damage the ${\sigma }_{ij}$ terms left untouched by the original approach: ${\sigma }_{tt},{\sigma }_{tm},{\sigma }_{mm}$ (Equation (45)). This additional damaging process can be performed through additional traction–separation equations with specific decay behaviours. This brings a set of three additional equations to the system:

$$\begin{array}{c}{\sigma }_{tt}=T_{tt}+\sum _k{M}_{t_{tk}}{\xi }_k=q_{tt}\left(\mathit{\xi }\right)\end{array}$$

(74a)

$$\begin{array}{c}{\sigma }_{tm}=T_{tm}+\sum _k{M}_{t_{mk}}{\xi }_k=q_{tm}\left(\mathit{\xi }\right)\end{array}$$

(74b)

$$\begin{array}{c}{\sigma }_{mm}=T_{mm}+\sum _k{M}_{m_{mk}}{\xi }_k=q_{mm}\left(\mathit{\xi }\right),\end{array}$$

(74c)

where the new load-driven terms $T_{tt},T_{tm},T_{mm}$ are calculated using the same structure as Equation (38), but now using stress projection operators ${\mathbf{H}}_{st},{\mathbf{H}}_{sm}$ in the remaining local directions $\widehat{\mathit{t}},\widehat{\mathit{m}}$, respectively:

$${\mathbf{H}}_{st}=\left[\begin{array}{ccc}{t}_x& 0& 0\\ 0& t_y& 0\\ 0& 0& t_z\\ t_y& t_x& 0\\ 0& t_z& t_y\\ t_z& 0& t_x\end{array}\right],\mathrm{in}\mathrm{local}\mathrm{frame}:{\mathbf{H}}_{st}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\\ 1& 0& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right]$$

(75)

$$\begin{array}{c}{\mathbf{H}}_{sm}=\left[\begin{array}{ccc}{m}_x& 0& 0\\ 0& m_y& 0\\ 0& 0& m_z\\ m_y& m_x& 0\\ 0& m_z& m_y\\ m_z& 0& m_x\end{array}\right],\mathrm{in}\mathrm{local}\mathrm{frame}:{\mathbf{H}}_{sm}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 1\\ 0& 0& 0\\ 0& 1& 0\\ 1& 0& 0\end{array}\right],\end{array}$$

(76)

and then using analogous expressions for new load-driven traction vectors ${\mathit{T}}_{et},{\mathit{T}}_{em}$ for finally projecting to $\widehat{\mathit{t}},\widehat{\mathit{m}}$ and obtaining the desired components:

$$\begin{array}{c}{\mathit{T}}_{et}={\mathbf{K}}_{et}\mathit{d}=\left(\frac{1}{V_e}\underset{\mathsf{\Omega }}{\int }{\mathbf{H}}_{st}^T\mathbf{C}\mathbf{B}dV\right)\mathit{d}={\mathbf{H}}_{st}^T\mathbf{C}\mathbf{B}\mathit{d}\end{array}$$

(77)

$$\begin{array}{c}{T}_{tt}={\mathit{T}}_{et}·\widehat{\mathit{t}},T_{tm}={\mathit{T}}_{et}·\widehat{\mathit{m}}\end{array}$$

(78)

$$\begin{array}{c}{\mathit{T}}_{em}={\mathbf{K}}_{em}\mathit{d}=\left(\frac{1}{V_e}\underset{\mathsf{\Omega }}{\int }{\mathbf{H}}_{sm}^T\mathbf{C}\mathbf{B}dV\right)\mathit{d}={\mathbf{H}}_{sm}^T\mathbf{C}\mathbf{B}\mathit{d}\end{array}$$

(79)

$$\begin{array}{c}{T}_{mm}={\mathit{T}}_{em}·\widehat{\mathit{m}},\end{array}$$

(80)

The stiffness vectors ${\mathit{M}}_{tt},{\mathit{M}}_{tm},{\mathit{M}}_{mm}$ are calculated in the same way as the rows of the original $\mathbf{M}$ matrix, but using the newly defined projection operators ${\mathbf{H}}_{st},{\mathbf{H}}_{sm}$:

$$\begin{array}{c}{\mathbf{M}}_t=\frac{1}{V_e}{\mathbf{H}}_{st}^T\mathbf{C}\left(V^{+}\overline{{\mathbf{G}}_{sb}^{\prime +}}+V^{-}\overline{{\mathbf{G}}_{sb}^{\prime -}}\right)\end{array}$$

(81)

$$\begin{array}{c}{\mathit{M}}_{tt}={\widehat{\mathit{t}}}^T{\mathbf{M}}_t,{\mathit{M}}_{tm}={\widehat{\mathit{m}}}^T{\mathbf{M}}_t\end{array}$$

(82)

$$\begin{array}{c}{\mathbf{M}}_m=\frac{1}{V_e}{\mathbf{H}}_{sm}^T\mathbf{C}\left(V^{+}\overline{{\mathbf{G}}_{sb}^{\prime +}}+V^{-}\overline{{\mathbf{G}}_{sb}^{\prime -}}\right)\end{array}$$

(83)

$$\begin{array}{c}{\mathit{M}}_{mm}={\widehat{\mathit{m}}}^T{\mathbf{M}}_m\end{array}$$

(84)

It is important to remark that introducing the system of Equation (74) is not the same as enforcing an artificial decay on ${\sigma }_{tt},{\sigma }_{tm},{\sigma }_{tt}$ as in Section 4.1. We are systematically involving all fracture modes ${\xi }_k$ in both systems (Equations (73) and (74)). Simultaneous solutions of these systems will ensure kinematically and statically consistent values for $\xi$ nullifying ${\sigma }_{tt},{\sigma }_{tm},{\sigma }_{tt}$, which should not exist at all once the crack develops.

The overall system in $\xi$ still remains underdetermined. The last remaining three equations are free to implement more effects in fracture dynamics without having to resort to $\phi$ parameters exclusively. For example, three final relations could be established to finish the uncoupling of rigid body translations within the system in Equation (73) by grouping and isolating the effects of all other kinematic modes different than rigid body displacements:

$$\begin{array}{c}{M}_{n{\theta }_n}{\theta }_n+M_{n{\theta }_t}{\theta }_t+M_{n{\theta }_m}{\theta }_m+M_{n{ϵ}_n}{ϵ}_n+M_{n{ϵ}_t}{ϵ}_t+M_{n{ϵ}_m}{ϵ}_m=0\end{array}$$

(85a)

$$\begin{array}{c}{M}_{t{\theta }_n}{\theta }_n+M_{t{\theta }_t}{\theta }_t+M_{t{\theta }_m}{\theta }_m+M_{t{ϵ}_n}{ϵ}_n+M_{t{ϵ}_t}{ϵ}_t+M_{t{ϵ}_m}{ϵ}_m=0\end{array}$$

(85b)

$$\begin{array}{c}{M}_{m{\theta }_n}{\theta }_n+M_{m{\theta }_t}{\theta }_t+M_{m{\theta }_m}{\theta }_m+M_{m{ϵ}_n}{ϵ}_n+M_{m{ϵ}_t}{ϵ}_t+M_{m{ϵ}_m}{ϵ}_m=0\end{array}$$

(85c)

This way, the system is completely closed and a complete solution for $\xi$ can be found at each load state $\mathit{d}$ after element localisation. The approach ensures that ${\sigma }_{tt},{\sigma }_{tm},{\sigma }_{mm}$ becomes progressively suppressed as the fracture develops. On the other hand, it might seem awkward to have a $\delta \left[\left|\mathit{u}\right|\right]$ discretised very differently with respect to $\left[\left|\mathit{u}\right|\right]$, but this already happens with other fields within the framework.

6. Elemental Validations

In this section, calculations with a single element have been performed 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 behaviour (evolution of the fracture traction and its related stress state) and kinematic behaviour (consistency between internal kinematic modes and true element kinematics through $\mathit{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 idealised regular element geometry and fracture orientations. The load direction is also often well-aligned to element edges or faces. Figure 12 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 deficiencies 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 the element level does not guarantee successful numerical simulations with larger scale models.

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 the local element level. If element fracture planes are generated through free and spontaneous element localisations, a complex 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 vector $\mathit{d}$ will 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.

A robust E-FEM formulation should be able to undertake such scenarios to be truly useful. From the view of the authors of this work, there is 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 the element level. We have taken one of such non-vanilla element samples coming from a real large scale model, having a specific geometry and fracture plane orientation. An approximate illustration of the chosen element (again, a linear tetrahedron) is shown in Figure 13. It has a characteristic length of approximately 1 mm. The reference frame is such that the normal axis to the fracture plane $\widehat{\mathit{n}}$ is almost perfectly aligned with the global z direction, but not to any of the tetrahedron’s faces. The remaining parallel directions $\widehat{\mathit{t}},\widehat{\mathit{m}}$ are not aligned with the $x,y$ directions. It has a 2-2 node partition with the volume of ${\mathsf{\Omega }}^{-}$ representing approximately 5% of total elemental volume while ${\mathsf{\Omega }}^{+}$ has the rest.

The load displacement vector $\mathit{d}$ has 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 generalised fracture kinematic modes ${\xi }_k$ described in Section 5.2. A monotonic, nonlinear behaviour is prescribed for each of these modes, reaching a final target value. A load progression factor $\beta$ between 0 and 1 is used to calculate any intermediate value between the starting point (zero) and the target value. Figure 14a–c show the chosen behaviours for each ${\xi }_k$. For this case, a quadratic behaviour has been assigned for rigid body translations, an inverted exponential behaviour for rigid body rotations and a radical behaviour 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. All kinematic modes are activated without exception. The load proposed ensures the element arrives at the terminal separation conditions at target values.

The load displacement vector $\mathit{d}$ associated with this particular evolution of the modes ${\xi }_k$ can be retrieved by knowing beforehand the nodes of the element corresponding to the ${\mathsf{\Omega }}^{+}$ domain and the mode linking equations discussed in Section 5.2. Once again, expressions have been explicitly detailed in Appendix C. Starting from a load ${\mathit{d}}_0$ corresponding to element localisation, the effective $\mathit{d}$ prescription can be calculated as a superposition of all effects:

$$\mathit{d}\left(\beta \right)={\mathit{d}}_0+\sum _j^{m,n,t}{\mathit{d}}_{{\left[\left|u\right|\right]}_j}\left(\beta \right)+\sum _j^{m,n,t}{\mathit{d}}_{{\theta }_j}\left(\beta \right)+\sum _j^{m,n,t}{\mathit{d}}_{{ϵ}_j}\left(\beta \right),$$

(86)

where ${\mathit{d}}_{{\left[\left|u\right|\right]}_j}$ are the displacement vectors associated with rigid body translation demands, ${\mathit{d}}_{{\theta }_j}$ for rigid body rotations and ${\mathit{d}}_{{ϵ}_j}$ the displacements related to simple axial strains.

Since the load requirements are finally prescribed based precisely on a controlled local fracture kinematics behaviour, a good E-FEM formulation should return a kinematics calculation $\xi$ as close to and most consistent as possible with the source fracture kinematics, now referred to as ${\xi }_{ref}$: the kinematic reference.

As mentioned before, the element is assumed to have already gone through localisation with a given stress state ${\sigma }_y$ at the fracture plane:

$$\left[\begin{array}{ccc}{\sigma }_{y_{nn}}& {\sigma }_{y_{nt}}& {\sigma }_{y_{nm}}\\ {\sigma }_{y_{tn}}& {\sigma }_{y_{tt}}& {\sigma }_{y_{tm}}\\ {\sigma }_{y_{mn}}& {\sigma }_{y_{mt}}& {\sigma }_{y_{mm}}\end{array}\right]\to {\sigma }_y=\left[\begin{array}{c}{\sigma }_{y_{nn}}\\ {\sigma }_{y_{nt}}\\ {\sigma }_{y_{nm}}\\ {\sigma }_{y_{tt}}\\ {\sigma }_{y_{tm}}\\ {\sigma }_{y_{mm}}\end{array}\right]=\left[\begin{array}{c}9.0\\ 0.0\\ 0.0\\ 2.5\\ 0.138\\ 1.459\end{array}\right]\mathrm{MPa},$$

where ${\sigma }_{y_{nn}},{\sigma }_{y_{tn}},{\sigma }_{y_{mn}}$ correspond to the components of the traction vector $T_n,T_t,T_m$ at localisation. These are calculated through a specific localisation criterion in the element, which in this case has been a Rankine criterion. The remaining components ${\sigma }_{ytt},{\sigma }_{ytm},{\sigma }_{ymm}$ are found by projecting the entire state of stresses to $\widehat{\mathit{t}}$ and $\widehat{\mathit{m}}$ exactly at this load level.

The traction–separation equations consider an exponential decay law for all pertinent fracture traction/stress components:

$$q_j={\sigma }_{y_j}{e}^{-\left(\frac{{\sigma }_{y_{nn}}}{G_{fI}}{\left[\left|u\right|\right]}_{n_0}+\frac{{\sigma }_{y_{nt}}}{G_{fII}}{\left[\left|u\right|\right]}_{t_0}+\frac{{\sigma }_{y_{nm}}}{G_{fII}}{\left[\left|u\right|\right]}_{m_0}\right)},$$

(87)

where ${\sigma }_{y_j}$ are each of the aforementioned initial fracture yield stress components. Note that the exponential argument is exclusively controlled by the fracture kinematics corresponding to rigid body translations. $G_{fI},G_{fII}$ are the associated fracture energies for fracture separation ${\left[\left|u\right|\right]}_{n_0}$ and fracture sliding ${\left[\left|u\right|\right]}_{t_0},{\left[\left|u\right|\right]}_{m_0}$, respectively. This is a particular proposal that will conveniently fit all formulation types proposed so far, although richer proposals could be made for the generalised kinematics approaches.

All calculations have been made 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 $\phi$, $\mathbf{J}$ and $\overline{{\mathbf{G}}_{sb}^{\prime }}$. 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 (Figure 13) was also originally performed using these tools.

In the following subsections, five different 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.

6.1. Single Mode Formulation

For a single mode approach, the calculation begins by considering rigid body translations ${\left[\left|u\right|\right]}_n$, ${\left[\left|u\right|\right]}_t$, ${\left[\left|u\right|\right]}_m$ and having only one of them different than zero (Section 4.1). To follow the example, a normal separation calculation will be made, leaving thus ${\left[\left|u\right|\right]}_t=0,{\left[\left|u\right|\right]}_m=0$. The $\phi$ function remains linear and unique, so that ${\mathbf{G}}_{sb}$ adopts its simplest way possible as a constant matrix (Equation (37)). The approach can be summarised in Equation (89).

The approach considers only the first equation of the traction–separation law system, and solves for ${\left[\left|u\right|\right]}_n$. This solution implies finding the intersection between a straight line $T_{en}+M_{nn}{\left[\left|u\right|\right]}_n$ and the nonlinear behaviour ${\sigma }_{yn}{e}^{-\frac{{\sigma }_{ynn}}{G_{fI}}{\left[\left|u\right|\right]}_n}$. The solution can be expressed in a closed, explicit form using the main branch of the Lambert function $W_0$ [22].

After a solution for ${\left[\left|u\right|\right]}_n$ is obtained, the approach then proceeds to calculate the remaining statics of the fracture plane, mainly the three fracture traction components $T_n,T_t,T_m$. At last, the remaining components of the state of stresses ${\sigma }_{tt},{\sigma }_{tm},{\sigma }_{mm}$ are retrieved for reference by calculating the entire stress $\sigma$ (Equation (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 ${\left[\left|u\right|\right]}_n$, ${\left[\left|u\right|\right]}_t$, ${\left[\left|u\right|\right]}_m$ (Section 4.2). Two versions from this family are considered for elemental simulations, depending on their assumptions made on the $\phi$ function. One considers a classical linear definition for $\phi$ (and therefore a constant ${\mathbf{G}}_{sb}$), representing the direct extension of the single mode approach of the previous section. The other builds a quadratic $\phi$ that helps the diagonalisation of the $\mathbf{M}$ matrix as well as for kinematic consistency at terminal separation conditions. Regardless of the $\phi$ assumptions, these formulations require a solution for the three-variable nonlinear system:

$$\mathit{T}={\mathit{T}}_e+\mathbf{M}\left[\left|\mathit{u}\right|\right]={\mathit{\sigma }}_y{e}^{-\left(\frac{{\sigma }_{y_{nn}}}{G_{fI}}{\left[\left|u\right|\right]}_n+\frac{{\sigma }_{y_{nt}}}{G_{fII}}{\left[\left|u\right|\right]}_t+\frac{{\sigma }_{y_{nm}}}{G_{fII}}{\left[\left|u\right|\right]}_m\right)},$$

(88)

which still has a closed, analytical solution using condensation techniques and the same Lambert W function $W_0$ approach. The first three-mode approach is summarised in Equation (90).

For this case, the calculation of the traction vector components $T_n,T_t,T_m$ and the stress components ${\sigma }_{tt},{\sigma }_{tm},{\sigma }_{mm}$ is performed exactly the same way as with the single mode approach using the solution values for ${\left[\left|u\right|\right]}_n,{\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m$.

The other three-mode approach will go forward and consider a more complex definition for the $\phi$ function, namely a quadratic one, as discussed at the beginning of Section 5.1 supported by a complete ${\mathit{P}}_2$ quadratic polynomial base (Equation (55)).

Ten free parameters ($\alpha$ coefficients) are available with this $\phi$ proposal. These will be used to build a linear system satisfying specific kinematic and static conditions: four of them to satisfy the basic boundary condition imposition requirements (Equation (3)), two for driving out-of-diagonal $M_{ij}$ fracture stiffness terms to zero (Equation (40)) and a last one to achieve consistent terminal separation conditions according to Equation (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 $\alpha$ coefficient vector with the smallest norm possible.

$$\begin{array}{c}\underset{¯}{\mathbf{Single\; Mode\; Formulation}}\\ \begin{array}{cccccc}\hfill \mathbf{Kinematic}\mathbf{modes}:& \hfill & & \hfill & & {\left[\left|u\right|\right]}_n\hfill \\ \hfill \mathbf{Linear},\mathbf{known}\phi :& \hfill & & \hfill & & \phi ={\alpha }_0+{\alpha }_{01}x+{\alpha }_{01}y\hfill \\ \hfill \mathbf{Fixed},\mathbf{constant}{\mathbf{G}}_{sb}:& \hfill & & \hfill & & {\mathbf{G}}_{sb}=-\sum _i^{N_e}{p}_i{\mathbf{B}}_i\hfill \\ \hfill \mathbf{Traction}–\mathbf{separation}\mathbf{system}:& \hfill & & \hfill & & \mathit{T}={\mathit{T}}_e+\mathbf{M}{\left[\left|\mathit{u}\right|\right]}_{\widehat{\mathit{n}}}=\mathbf{q}\left({\left[\left|\mathit{u}\right|\right]}_{\widehat{\mathit{n}}}\right)\hfill \\ \hfill \mathbf{Nullify}{\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m:& \hfill & & \hfill & & {\left[\left|u\right|\right]}_t=0,{\left[\left|u\right|\right]}_m=0\hfill \\ \hfill \mathbf{Solve}\mathbf{only}\mathbf{for}{\left[\left|u\right|\right]}_n:& \hfill & & \hfill & & T_{en}+M_{nn}{\left[\left|u\right|\right]}_n={\sigma }_{y_{nn}}{e}^{-\frac{{\sigma }_{y_{nn}}}{G_{fI}}{\left[\left|u\right|\right]}_n}\hfill \\ \hfill \mathbf{Calculate}{T}_n,T_t,T_m:& \hfill & & \hfill & & T_n=T_{en}+M_{nn}{\left[\left|u\right|\right]}_n\hfill \\ \hfill & \hfill & & \hfill & & T_t=T_{et}+M_{tn}{\left[\left|u\right|\right]}_n\hfill \\ \hfill & \hfill & & \hfill & & T_m=T_{em}+M_{mn}{\left[\left|u\right|\right]}_n\hfill \\ \hfill \mathbf{Calculate}\mathit{\sigma }:& \hfill & & \hfill & & \mathit{\sigma }=\mathbf{C}\mathit{\epsilon }=\mathbf{C}\left(\mathbf{B}\mathit{d}+{\mathbf{G}}_{sb}\widehat{\mathit{n}}{\left[\left|u\right|\right]}_n\right)\hfill \\ \hfill \mathbf{Determine}{\sigma }_{tt},{\sigma }_{tm},{\sigma }_{mm}:& \hfill & & \hfill & & {\sigma }_{tt}={\widehat{\mathit{t}}}^T\overline{\overline{\sigma }}\widehat{\mathit{t}},{\sigma }_{tm}={\widehat{\mathit{t}}}^T\overline{\overline{\sigma }}\widehat{\mathit{m}}\hfill \\ \hfill & \hfill & & \hfill & & {\sigma }_{mm}={\widehat{\mathit{m}}}^T\overline{\overline{\sigma }}\widehat{\mathit{m}}\hfill \end{array}\end{array}$$

(89)

$$\begin{array}{c}\underset{¯}{\mathbf{Three\; modes\; formulation\; approach\; definition}}\\ \begin{array}{cccccc}\hfill \mathbf{Kinematic}\mathbf{modes}:& \hfill & & \hfill & & {\left[\left|u\right|\right]}_n,{\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m\hfill \\ \hfill \mathbf{Linear},\mathbf{known}\phi :& \hfill & & \hfill & & \phi ={\alpha }_0+{\alpha }_{01}x+{\alpha }_{01}y\hfill \\ \hfill \mathbf{Fixed},\mathbf{constant}{\mathbf{G}}_{sb}:& \hfill & & \hfill & & {\mathbf{G}}_{sb}=-\sum _i^{N_e}{p}_i{\mathbf{B}}_i\hfill \\ \hfill \mathbf{Traction}-\mathbf{separation}\mathbf{system}\searrow & \hfill & & \hfill & & \\ \hfill \mathbf{Solve}\mathbf{for}{\left[\left|u\right|\right]}_n,{\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m:& \hfill & & \hfill & & \mathit{T}={\mathit{T}}_e+\mathbf{M}{\left[\left|\mathit{u}\right|\right]}_{\widehat{\mathit{n}}}=\mathbf{q}\left({\left[\left|\mathit{u}\right|\right]}_{\widehat{\mathit{n}}}\right)\hfill \\ \hfill \mathbf{Calculate}{T}_n,T_t,T_m\nearrow & \hfill & & \hfill & & \\ \hfill \mathbf{Calculate}\mathit{\sigma }:& \hfill & & \hfill & & \mathit{\sigma }=\mathbf{C}\mathit{\epsilon }=\mathbf{C}\left(\mathbf{B}\mathit{d}+{\mathbf{G}}_{sb}\mathbf{R}{\left[\left|\mathit{u}\right|\right]}_{\widehat{\mathit{n}}}\right)\hfill \\ \hfill \mathbf{Determine}{\sigma }_{tt},{\sigma }_{tm},{\sigma }_{mm}:& \hfill & & \hfill & & {\sigma }_{tt}={\widehat{\mathit{t}}}^T\overline{\overline{\sigma }}\widehat{\mathit{t}},{\sigma }_{tm}={\widehat{\mathit{t}}}^T\overline{\overline{\sigma }}\widehat{\mathit{m}}\hfill \\ \hfill & \hfill & & \hfill & & {\sigma }_{mm}={\widehat{\mathit{m}}}^T\overline{\overline{\sigma }}\widehat{\mathit{m}}\hfill \end{array}\end{array}$$

(90)

Once $\mathit{\alpha }$ is solved, the element internal solution process continues as usual, solving the three-variable system and calculating the fracture traction vector. The complete fracture stress $\mathit{\sigma }$ has to be calculated using a volume average (Equation (21)) as ${\mathbf{G}}_{sb}$ is no longer constant in space. Please refer to Equation (91).

$$\begin{array}{c}\underset{¯}{\mathbf{Three\; modes\; formulation,\; controlled}\phi }\\ \begin{array}{cccccc}\hfill \mathbf{Kinematic}\mathbf{modes}:& \hfill & & \hfill & & {\left[\left|u\right|\right]}_n,{\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m\hfill \\ \hfill \mathbf{Quadratic},\mathbf{unknown}\mathbf{single}\phi :& \hfill & & \hfill & & \phi ={\mathit{P}}_2^T\mathit{\alpha }\hfill \\ \hfill \mathbf{Linear}{\mathbf{G}}_{sb},\mathbf{function}\mathbf{of}\mathit{\alpha }:& \hfill & & \hfill & & {\mathbf{G}}_{sb}=-\sum _i^{N_e}{p}_i\nabla \mathit{\phi }\left(\mathit{\alpha },x,y,z\right)\hfill \\ \hfill \mathbf{Solution}\mathbf{for}\mathit{\alpha }\mathbf{coefficients}:& \hfill & & \hfill & & \phi ({\mathit{x}}_i,\mathit{\alpha })=p_i,M_{i\ne j}\left(\mathit{\alpha }\right)=0\hfill \\ \hfill & \hfill & & \hfill & & M_{nn}\left(\mathit{\alpha }\right)=M_{n{n}_{\infty }}\hfill \\ \hfill \mathbf{Traction}–\mathbf{separation}\mathbf{system}\searrow & \hfill & & \hfill & & \\ \hfill \mathbf{Solve}\mathbf{for}{\left[\left|u\right|\right]}_n,{\left[\left|u\right|\right]}_t,{\left[\left|u\right|\right]}_m:& \hfill & & \hfill & & \mathit{T}={\mathit{T}}_e+\mathbf{M}{\left[\left|\mathit{u}\right|\right]}_{\widehat{\mathit{n}}}=\mathit{q}\left({\left[\left|\mathit{u}\right|\right]}_{\widehat{\mathit{n}}}\right)\hfill \\ \hfill \mathbf{Calculate}{T}_n,T_t,T_m\nearrow & \hfill & & \hfill & & \\ \hfill \mathbf{Calculate}\mathit{\sigma }:& \hfill & & \hfill & & \mathit{\sigma }=\frac{1}{V_e}\mathrm{\int }\mathbf{C}\left[\mathbf{B}\mathit{d}+{\mathbf{G}}_{sb}\mathbf{R}{\left[\left|u\right|\right]}_{\widehat{\mathit{n}}}\right]dV\hfill \\ \hfill \mathbf{Determine}{\sigma }_{tt},{\sigma }_{tm},{\sigma }_{mm}:& \hfill & & \hfill & & {\sigma }_{tt}={\widehat{\mathit{t}}}^T\overline{\overline{\sigma }}\widehat{\mathit{t}},{\sigma }_{tm}={\widehat{\mathit{t}}}^T\overline{\overline{\sigma }}\widehat{\mathit{m}}\hfill \\ \hfill & \hfill & & \hfill & & {\sigma }_{mm}={\widehat{\mathit{m}}}^T\overline{\overline{\sigma }}\widehat{\mathit{m}}\hfill \end{array}\end{array}$$

(91)

6.3. Enriched Kinematics Formulation

The last formulation type uses full enriched fracture kinematics as described in Section 4.3 and Section 5.2 through a generalised mode vector $\xi$:

$$\xi ={\left[\begin{array}{ccccccccc}{\left[\left|u\right|\right]}_{n0}& {\left[\left|u\right|\right]}_{t0}& {\left[\left|u\right|\right]}_{m0}& {\theta }_n& {\theta }_t& {\theta }_m& {ϵ}_n& {ϵ}_t& {ϵ}_m\end{array}\right]}^t,$$

counting three rigid body displacements, three rigid body rotations and three simple axial strains. Once again, two different models from this family are included in the present analysis: one considering a fixed definition for $\phi$ and the other a more complex, parameterised one.

For a fixed $\phi$ approach, the idea is to numerically define the $\overline{{\mathbf{G}}_{sb}^{\prime \pm }}$ operators by calculating each column corresponding to each of the modes ${\xi }_k$ using 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 specific expressions for each of the columns for $\overline{{\mathbf{G}}_{sb}^{\prime \pm }}$ are explicitly detailed in Appendix C, 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 $\phi$ function. This formulation is summarised in Equation (92).

$$\begin{array}{c}\underset{¯}{\mathbf{Enriched\; Modes\; Formulation,\; Fixed}\phi }\\ \begin{array}{cccccc}\hfill \mathbf{Kinematic}\mathbf{modes}:& \hfill & & \hfill & & {\left[\left|u\right|\right]}_{n_0},{\left[\left|u\right|\right]}_{t_0},{\left[\left|u\right|\right]}_{m_0}\hfill \\ \hfill & \hfill & & \hfill & & {\theta }_n,{\theta }_t,{\theta }_m;{ϵ}_n,{ϵ}_t,{ϵ}_m\hfill \\ \hfill \mathbf{Fixed}(\mathbf{unknown}\mathbf{shape})\phi :& \hfill & & \hfill & & \mathrm{Not}\mathrm{used}\mathrm{explicitly}\hfill \\ \hfill \mathbf{Fixed}\mathbf{columns}\mathbf{for}\overline{{\mathbf{G}}_{sb}^{\prime \pm }}:& \hfill & & \hfill & & {\overline{{\mathbf{G}}_{sb}^{\prime \pm }}}_k{\xi }_k=-\mathbf{B}{\mathit{d}}_k\hfill \\ \hfill \mathbf{Traction}–\mathbf{separation}\mathbf{system}:& \hfill & & \hfill & & T_n=T_{en}+\sum M_{n_{nk}}{\xi }_k=q_{nn}\left(\mathit{\xi }\right)\hfill \\ \hfill & \hfill & & \hfill & & T_t=T_{et}+\sum M_{n_{tk}}{\xi }_k=q_{nt}\left(\mathit{\xi }\right)\hfill \\ \hfill & \hfill & & \hfill & & T_m=T_{em}+\sum M_{n_{mk}}{\xi }_k=q_{nm}\left(\mathit{\xi }\right)\hfill \\ \hfill \mathbf{Complementary}\mathbf{damage}:& \hfill & & \hfill & & {\sigma }_{tt}=T_{tt}+\sum M_{t_{tk}}{\xi }_k=q_{tt}\left(\mathit{\xi }\right)\hfill \\ \hfill & \hfill & & \hfill & & {\sigma }_{tm}=T_{tm}+\sum M_{t_{mk}}{\xi }_k=q_{tm}\left(\mathit{\xi }\right)\hfill \\ \hfill & \hfill & & \hfill & & {\sigma }_{mm}=T_{mm}+\sum M_{m_{mk}}{\xi }_k=q_{mm}\left(\mathit{\xi }\right)\hfill \\ \hfill \mathbf{Closing}\mathbf{constraints}:& \hfill & & \hfill & & \sum A_{lk}{\xi }_k=B_l,l=1,2,3\hfill \\ \hfill \mathbf{Solve}\mathbf{for}\mathit{\xi }& \hfill & & \hfill & & \searrow \hfill \\ \hfill \mathbf{Calculate}{T}_n,T_t,T_m& \hfill & & \hfill & & \to \mathrm{From}\mathrm{the}\mathrm{set}\mathrm{of}9\mathrm{equations}\mathrm{just}\mathrm{above}\hfill \\ \hfill \mathrm{Calculate}{\sigma }_{tt},{\sigma }_{tm},{\sigma }_{mm}& \hfill & & \hfill & & \nearrow \hfill \end{array}\end{array}$$

(92)

Once the $\overline{{\mathbf{G}}_{sb}^{\prime \pm }}$ operator is populated, the approach follows the ideas proposed in Section 5.4 for building the three equation system depicted in Equation (71) for the traction vector components $T_n,T_n,T_n$. Three additional equations for damaging the fracture stress components ${\sigma }_{tt},{\sigma }_{tt},{\sigma }_{tt}$ are built as well following the structure in Equation (74). The last required three equations are taken from Equation (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 ${\xi }_k$, using some of them as buffer variables. At the end, a $9\times 9$ nonlinear system is obtained, regardless of the $\phi$ approach. Once the system is solved, all fracture statics variables $(T_n,T_t,T_m,{\sigma }_{tt},{\sigma }_{tm},{\sigma }_{mm})$ will be available directly.

The controlled $\phi$ version of this model deals with a considerably more complex $\phi$ definition compared to that used on the three-mode approach: a full cubic, piece-wise, triple $\phi$ function capable of delivering up to 120 free parameters. The ${\mathbf{G}}_{sb}^{\prime \pm }$ operators will be dependent on the $\mathit{\alpha }$ coefficients following the discussion of Section 5.2, and mode linking relations will serve to build linear restrictions on $\mathit{\alpha }$. The specific structure of the complete linear system on $\mathit{\alpha }$ for this proposal is detailed in 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 $\mathit{\alpha }$ system is solved, the static solution process continues exactly the same as mentioned before. The approach is stated in Equation (93).

$$\begin{array}{c}\underset{¯}{\mathbf{Enriched\; Modes\; Formulation,\; Controlled}\phi }\\ \begin{array}{cccccc}\hfill \mathbf{Kinematic}\mathbf{modes}:& \hfill & & \hfill & & {\left[\left|u\right|\right]}_{n_0},{\left[\left|u\right|\right]}_{t_0},{\left[\left|u\right|\right]}_{m_0}\hfill \\ \hfill & \hfill & & \hfill & & {\theta }_n,{\theta }_t,{\theta }_m;{ϵ}_n,{ϵ}_t,{ϵ}_m\hfill \\ \hfill \mathbf{Cubic},\mathbf{Piece-wise},\mathbf{triple}\phi :& \hfill & & \hfill & & {\phi }_j^{\pm }={\mathit{P}}_3^T{\mathit{\alpha }}_j^{\pm },j=n,t,m\hfill \\ \hfill \mathbf{Quad}.{\mathbf{G}}_{sb}^{\prime \pm },\mathbf{function}\mathbf{of}\mathit{\alpha }:& \hfill & & \hfill & & {\mathbf{G}}_{sb}^{\prime \pm }\left(\mathit{\alpha }\right)={\mathcal{H}}_{\mathrm{\Gamma }}\nabla \mathbf{J}-\left(\mathit{\phi }\nabla \mathbf{J}\right)-\nabla \mathit{\phi }\mathbf{J}\hfill \\ \hfill \mathbf{Solution}\mathbf{for}\mathit{\alpha }\mathbf{coefficients}:& \hfill & & \hfill & & {\overline{{\mathbf{G}}_{sb}^{\prime \pm }}}_k\left(\mathit{\alpha }\right){\xi }_k=-\mathbf{B}{\mathit{d}}_k\hfill \\ \hfill & \hfill & & \hfill & & \phi ({\mathit{x}}_i,\mathit{\alpha })=p_i\hfill \\ \hfill & \hfill & & \hfill & & {\left.{\phi }^{+}\left(\mathit{\alpha }\right)\right|}_{{\mathrm{\Gamma }}_d}={\left.{\phi }^{-}\left(\mathit{\alpha }\right)\right|}_{{\mathrm{\Gamma }}_d}\hfill \\ \hfill & \hfill & & \hfill & & M_{i\ne j}\left(\mathit{\alpha }\right)=0,i,j=n,t,m\hfill \\ \hfill & \hfill & & \hfill & & M_{nn}\left(\mathit{\alpha }\right)=M_{n{n}_{\infty }},M_{tt}\left(\mathit{\alpha }\right)=M_{t{t}_{\infty }}\hfill \\ \hfill & \hfill & & \hfill & & M_{mm}\left(\mathit{\alpha }\right)=M_{m{m}_{\infty }}\hfill \\ \hfill \mathbf{Traction}-\mathbf{separation}\mathbf{system}:& \hfill & & \hfill & & T_n=T_{en}+\sum M_{n_{nk}}{\xi }_k=q_{nn}\left(\mathit{\xi }\right)\hfill \\ \hfill & \hfill & & \hfill & & T_t=T_{et}+\sum M_{n_{tk}}{\xi }_k=q_{nt}\left(\xi \right)\hfill \\ \hfill & \hfill & & \hfill & & T_m=T_{em}+\sum M_{n_{mk}}{\xi }_k=q_{nm}\left(\mathit{\xi }\right)\hfill \\ \hfill \mathbf{Complementary}\mathbf{damage}:& \hfill & & \hfill & & {\sigma }_{tt}=T_{tt}+\sum M_{t{t}_i}{\xi }_i=q_{tt}\left(\mathit{\xi }\right)\hfill \\ \hfill & \hfill & & \hfill & & {\sigma }_{tm}=T_{tm}+\sum M_{t{m}_i}{\xi }_i=q_{tm}\left(\mathit{\xi }\right)\hfill \\ \hfill & \hfill & & \hfill & & {\sigma }_{mm}=T_{mm}+\sum M_{m{m}_i}{\xi }_i=q_{mm}\left(\mathit{\xi }\right)\hfill \\ \hfill \mathbf{Closing}\mathbf{constraints}:& \hfill & & \hfill & & \sum A_{lk}{\xi }_k=B_l,l=1,2,3\hfill \\ \hfill \mathbf{Solve}\mathbf{for}\mathit{\xi }& \hfill & & \hfill & & \searrow \hfill \\ \hfill \mathbf{Calculate}{T}_n,T_t,T_m& \hfill & & \hfill & & \to \mathrm{From}\mathrm{the}\mathrm{set}\mathrm{of}9\mathrm{equations}\mathrm{just}\mathrm{above}\hfill \\ \hfill \mathbf{Calculate}{\sigma }_{tt},{\sigma }_{tm},{\sigma }_{mm}& \hfill & & \hfill & & \nearrow \hfill \end{array}\end{array}$$

(93)

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 five formulations types. The results are divided into two sections: static analysis and kinematic analysis. The static analysis pertains the results of all fracture traction vector components $T_n,T_t,T_m$ and the fracture stress state components ${\sigma }_{tt},{\sigma }_{tm},{\sigma }_{mm}$. The kinematic analysis concerns the behaviour 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 defined at the beginning of this section.

6.4.1. Static Results

The load profile is conceived with enough crack separation and sliding for placing the element at terminal separation conditions. The ${\mathsf{\Omega }}^{+}$ and the ${\mathsf{\Omega }}^{-}$ 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 $T_n,T_t,T_m$ and all the remaining stress state components associated with the fracture interface ${\sigma }_{tt},{\sigma }_{tm},{\sigma }_{mm}$ should be zero at this point. The expected behaviour for a good E-FEM formulation would be then to drive all these static components to zero at the end of the load profile. That is, the following plots evidence the damaging quality of the formulations. Figure 15 shows the behaviour of the traction vector components $T_n$, $T_t$, $T_m$ as a function of the load factor $\beta$. All formulations manage to damage $T_n$ completely: this is the most basic functionality of any E-FEM-based model. Generally speaking, the decaying rate depends on the fracture stiffness component $M_{nn}$, with the only exception of the three-mode, fixed $\phi$ formulation. The value of $M_{nn}$ will differ slightly for formulations where normal separation consistency has been ensured through $\phi$ coefficients (Equation (59)). The three-mode, fixed $\phi$ approach will have an effective normal fracture stiffness that will be the resulting effect of combining $M_{nn}$ as well as coupled stiffness terms $M_{nt}$, $M_{nm}$, which will produce for this case a very low effective normal stiffness, so that the traction $T_n$ is driven down immediately after starting the load application.

For $T_t$ and $T_m$, 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 ${\sigma }_{ynt},{\sigma }_{ynm}$ have begun at zero attending to the Rankine criterion assumed before. The task for the formulations is just 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 do 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,24,33,34]. Even with its partial loss of physical meaning at the local level, its numerical robustness and fracture representation capabilities at large scale models have proven to be successful.

Figure 16 shows the behaviour of the stress state components ${\sigma }_{tt}$, ${\sigma }_{tm}$ and ${\sigma }_{mm}$. Fracture stress components are only accessible when projecting the complete state $\sigma$ to the $\widehat{\mathit{t}}$ and $\widehat{\mathit{m}}$ directions. Only the enriched mode formulations are designed to explicitly track these components and to damage them properly. They do so at different rates, again, due to $\phi$ implications on fracture stiffness. All remaining formulations simply make these stress components 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.

6.4.2. 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 ${\xi }_k$, knowing beforehand that the global approach taken for all formulations possesses mathematical ambiguities that will always prevent a perfect match (Section 4.3).

Figure 17 shows the analysis of rigid body displacements. The normal crack separation ${\left[\left|u\right|\right]}_{n_0}$ plot is the one still having all five formulations with correct predictions at once. The formulations having a controlled $\phi$ will be naturally closer to the kinematic reference as the fracture stiffness has been designed to be consistent with normal displacement at terminal separation conditions. The other two having slightly different $M_{nn}$ miss the target by 30% but still retain the same behaviour tendency. Again, all formulations were prepared, one way or another, to handle normal separation kinematics. The three-mode, fixed $\phi$ again, has a very different effective normal stiffness, which causes even the first load step solution for ${\left[\left|u\right|\right]}_{n_0}$ to be excessively large.

The lack of kinematic consistency in the models starts to be truly appreciated in the analysis of crack sliding modes ${\left[\left|u\right|\right]}_{t_0},{\left[\left|u\right|\right]}_{m_0}$. In this case, there is only one single approach that achieves the correct the kinematic reference: the enhanced kinematics with controlled $\phi$. This is because this approach is the only one with a definition for $\phi$ complex enough to allow free parameters to help impose terminal separation conditions on the fracture sliding stiffness terms $M_{tt},M_{mm}$. It is important to note that, while the three-mode, controlled $\phi$ attempts to do the same, its limited $\phi$ structure does not brake the linear dependency observed on the diagonal of the $\mathbf{M}$ matrix. In this case, only one direction can be made to comply with terminal separation conditions, and it was chosen to be $\widehat{\mathit{n}}$. The three-mode fixed $\phi$ approach, once again, has mixed stiffness interactions, which drive the effective sliding stiffness to unrealistic levels in both directions.

Formulations not being able to control sliding fracture stiffness will even predict a completely opposite sliding direction to that demanded by the load vector $\mathit{d}$. This totally breaks the physical sense of these kinematic modes. It is 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 18 shows the analysis on rigid body rotations. From now on, only the last two formulations involving enriched modes are able to contend. The kinematic reference remains low compared to the displacement demands, so that the task for the formulations is basically to keep the response as low as possible. This is related to the fact that rotation modes, as well as axial strain modes, are mainly conceived as buffer variables that help kinematic consistency efforts to focus on fracture displacement modes.

In all cases, both formulations follow the same rotation tendency, having the controlled $\phi$ version exhibit a considerably lower response. The torsion ${\theta }_n$ presents the most notable deviations, with almost ${83}^{\circ }$ for the case of a fixed $\phi$ and $8^{\circ }$ for a controlled one.

Finally, the same trend is observed with the simple axial strains, as shown in Figure 19. The fixed $\phi$ approach presents excessively high values for lateral strains. No formulation is able to contain the strain response at zero, but the results suggest the controlled $\phi$ formulation fairs better at controlling this buffer kinematic mode.

To better explain the tendency of results in both Figure 18 and Figure 19, it is worth noting that the exponential decay proposed (Equation (87)) for damaging all traction vector and stress components related to the local fracture only involves the rigid body translation fracture kinematic modes ${\left[\left|u\right|\right]}_{n_0},{\left[\left|u\right|\right]}_{t_0},{\left[\left|u\right|\right]}_{m_0}$. This implies that, when the local fracture develops, the normal separation and parallel sliding modes are the ones dominating the development and the final onset of the element segmentation. Thus, it is also implied that rotation and simple axial strain modes $({\theta }_n,{\theta }_t,{\theta }_m,{ϵ}_n,{ϵ}_t,{ϵ}_m)$ will not have the same degree of influence on the overall fracture process and that the formulation as it is intends not to involve them during this part of the solution for all local element variables. This is a very particular choice made by the authors in this work, and a different choice could be made.

Such choice will sacrifice some of the capability to achieve kinematic consistency for rigid body rotations and simple axial strains up to a certain degree. In the case of the fixed and controlled $\phi$ formulations portrayed in Figure 18 and Figure 19, it is seen that both have the natural tendency to diverge from low values as the external load is increased. However, the most complex formulation is able to deal better with such partial loss of robustness and is rendered acceptable by the authors of this work. Thus, this choice allows for a reasonable compromise between obtaining near perfect kinematic consistency for an entire 3D E-FEM framework and keeping it mathematically manageable as to make a large-scale implementation process as efficient and less error-prone as possible. As the analysis in previous cases, a single mode formulation would technically obtain the correct kinematic consistency, as having these modes by default as non-existent implies a zero value for them and therefore a much closer match with respect to the kinematic reference in Figure 18 and Figure 19. However, it should be kept in mind that single mode formulations come with the inconvenience of inherent static inconsistency, and they will not function properly without a deliberate element erosion process. The latest generation of E-FEM formulations avoids such undesired traits.

7. Conclusions

An extensive study has been presented on the E-FEM strong discontinuity analysis framework, stating a solid basis of its fundamentals and discussing the most common recent formulation approaches. This base has helped to analyse general flaws observed in current formulations and to propose reasonable options for evolving the framework. Many authors had already hinted at some improvements. This work has consolidated this knowledge, added some enhancements and conformed a generalised working proposal in three-dimensional models.

A realistic test element has been taken as a reference to compare the performance of all described formulations. Simulations at this level have been performed considering a prescription of fracture kinematic modes just right after localisation, producing a nodal displacement profile. This load has been designed to make the element reach terminal separation conditions. A comparison of all responses allowed to study relevant behaviours on fracture kinematics and statics.

As more complexity is introduced in the formulations, an increase in consistency is perceived in both the kinematic and static outputs of the element. The most elaborate proposal (fully enriched kinematics with explicit $\phi$ parameters) is the only one to report a relatively clean and sound fracture kinematics vector in its entirety. However, its mathematical complexity and computational cost would raise questions on its practicality and ease of implementation. For the case developed in this study, this last formulation requires solving (only once) a linear system of 120 equations to finish the definition of the structure of the enhancement functions. The previous formulation following in degree of complexity (fully enriched kinematics, direct $\phi$) does not require solving such system at all. Still, it is able to satisfactorily remove all relevant fracture stresses, avoiding all numerical locking nuisances observed with other formulations.

Is it really worth adding this much complexity to the formulations just for the sake of obtaining cleaner fracture statics and kinematics? If it is desired to use the fracture kinematics information for more advanced modelling of internal fracture phenomena such as frictional sliding, fault reclosure and the effects of compression, then having these outputs correct becomes important. These complex calculations will rely on having a more realistic state of how the fracture moves. On the other hand, if the modelling goal is just to have a general idea of how a global fracture pattern starts to develop in a model, an approximate path and an estimation of a global strength, the fixed $\phi$ approaches will suffice.

In all cases, while none of the formulations are perfectly robust and mechanically correct, it remains satisfactory to see that the E-FEM framework can be taken this far for the modelling of strong discontinuities, gaining a lot in physical meaningfulness and still retaining the charms that normally attract many computational mechanics researchers.