On the role of enrichment and statical admissibility of recovered fields in
a-posteriori error estimation for enriched finite element methods

Octavio A. González-Estrada1, Juan José Ródenas2, Stéphane P.A. Bordas1, Marc
Duflot3, Pierre Kerfriden1, Eugenio Giner2

1Institute of Mechanics and Advanced Materials. Cardiff School of Engineering,
Cardiff University, The Parade, Cardiff CF24 3AA Wales, UK.
2Centro de Investigación de Tecnología de Vehículos(CITV),
Universitad Politècnica de València, E-46022-Valencia, Spain.
3CENAERO, Rue des Frères Wright 29, B-6041 Gosselies, Belgium

###### Abstract

Purpose – This paper aims at assessing the effect of (1) the statical
admissibility of the recovered solution; (2) the ability of the recovered
solution to represent the singular solution; on the accuracy, local and global
effectivity of recovery-based error estimators for enriched finite element
methods (e.g. the extended finite element method, XFEM).
Design/methodology/approach – We study the performance of two recovery
techniques. The first is a recently developed superconvergent patch recovery
procedure with equilibration and enrichment (SPR-CX). The second is known as
the extended moving least squares recovery (XMLS), which enriches the
recovered solutions but does not enforce equilibrium constraints. Both are
extended recovery techniques as the polynomial basis used in the recovery
process is enriched with singular terms for a better description of the
singular nature of the solution.
Findings – Numerical results comparing the convergence and the effectivity
index of both techniques with those obtained without the enrichment
enhancement clearly show the need for the use of extended recovery techniques
in Zienkiewicz-Zhu type error estimators for this class of problems. The
results also reveal significant improvements in the effectivities yielded by
statically admissible recovered solutions.
Originality/value – This work shows that both extended recovery procedures and
statical admissibility are key to an accurate assessment of the quality of
enriched finite element approximations.

Keywords extended finite element method; error estimation; linear elastic
fracture mechanics; statical admissibility; extended recovery
Paper Type Research Paper

##  1 Introduction

Engineering structures, in particular in aerospace engineering, are intended
to operate with flawless components, especially for safety critical parts.
However, there is always a possibility that cracking will occur during
operation, risking catastrophic failure and associated casualties.

The mission of Damage Tolerance Assessment (DTA) is to assess the influence of
these defects, cracks and damage on the ability of a structure to perform
safely and reliably during its service life. An important goal of DTA is to
estimate the fatigue life of a structure, i.e. the time during which it
remains safe given a pre-existing flaw.

Damage tolerance assessment relies on the ability to accurately predict crack
paths and growth rates in complex structures. Since the simulation of three-
dimensional crack growth is either not supported by commercial software, or
requires significant effort and time for the analysts and is generally not
coupled to robust error indicators, reliable, quality-controlled, industrial
damage tolerance assessment is still a major challenge in engineering
practice.

The extended finite element method (XFEM) (Moës et al., 1999) is now one of
many successful numerical methods to solve fracture mechanics problems. The
particular advantage of the XFEM, which relies on the partition of unity (PU)
property (Melenk and Babuška, 1996) of finite element shape functions, is its
ability to model cracks without the mesh conforming to their geometry. This
allows the crack to split the background mesh arbitrarily, thereby leading to
significantly increased freedom in simulating crack growth.

This feat is achieved by adding new degrees of freedom to (i) describe the
discontinuity of the displacement field across the crack faces, within a given
element, and (ii) reproduce the asymptotic fields around the crack tip. Thanks
to the advances made in the XFEM during the last years, the method is now
considered to be a robust and accurate means of analysing fracture problems,
has been implemented in commercial codes (ABAQUS, 2009, CENAERO, 2011) and is
industrially in use for damage tolerance assessment of complex three
dimensional structures (Bordas and Moran, 2006, Wyart et al., 2007, Duflot,
2007).

While XFEM allows to model cracks without (re)meshing the crack faces as they
evolve and yields exceptionally accurate stress intensity factors (SIF) for 2D
problems, Bordas and Moran (2006), Wyart et al. (2007), Duflot (2007) show
that for realistic 3D structures, a “very fine” mesh is required to accurately
capture the complex three-dimensional stress field and obtain satisfactory
stress intensity factors (SIFs), the main drivers of linear elastic crack
propagation.

Consequently, based on the mesh used for the stress analysis, a new mesh
offering sufficient refinement throughout the whole potential path of the
crack must be constructed _a priori_ , i.e. before the crack path is known.
Practically, this is done by running preliminary analyses on coarse meshes to
obtain an approximative crack path and heuristically refining the mesh around
this path to increase accuracy.

Typically, this refinement does not rely on sound error measures, thus the
heuristically chosen mesh is in general inadequately suited and can cause
large inaccuracies in the crack growth path, especially around holes, which
can lead to non-conservative estimates of the safe life of the structure.

Thus, it is clear that although XFEM simplifies the treatment of cracks
compared to the standard FEM by lifting the burden of a geometry conforming
mesh, it still requires iterations and associated user intervention and it
employs heuristics which are detrimental to the robustness and accuracy of the
simulation process.

It would be desirable to minimise the changes to the mesh topology, thus user
intervention, while ensuring the stress fields and SIFs are accurately
computed at each crack growth step. This paper is one more step in attempting
to control the discretization error committed by enriched FE approximations,
decrease human intervention in damage tolerance assessment of complex
industrial structures and enhance confidence in the results by providing
enriched FEMs with sound error estimators which will guarantee a predetermined
accuracy level and suppress recourse to manual iterations and heuristics.

The error assessment procedures used in finite element analysis are well known
and can be usually classified into different families (Ainsworth and Oden,
2000): residual based error estimators, recovery based error estimators, dual
techniques, etc. Residual based error estimators were substantially improved
with the introduction of the _residual equilibration_ by Ainsworth and Oden
(2000). These error estimators have a strong mathematical basis and have been
frequently used to obtain lower and upper bounds of the error (Ainsworth and
Oden, 2000, Díez et al., 2004). Recovery based error estimates were first
introduced by Zienkiewicz and Zhu (1987) and are often preferred by
practitioners because they are robust and simple to use and provide an
enhanced solution. Further improvements were made with the introduction of new
recovery processes such as the superconvergent patch recovery (SPR) technique
proposed by Zienkiewicz and Zhu (1992a, b). Dual techniques based on the
evaluation of two different fields, one compatible in displacements and
another equilibrated in stresses, have also been used by Pereira et al. (1999)
to obtain bounds of the error. Herein we are going to focus on recovery based
techniques which follow the ideas of the Zienkiewicz-Zhu (ZZ) error estimator
proposed by Zienkiewicz and Zhu (1987).

The literature on error estimation techniques for mesh based partition of
unity methods is still scarce. One of the first steps in that direction was
made in the context of the Generalized Finite Element Method (GFEM) by
Strouboulis et al. (2001). The authors proposed a recovery-based error
estimator which provides good results for h-adapted meshes. In a later work
two new a posteriori error estimators for GFEM were presented (Strouboulis et
al., 2006). The first one was based on patch residual indicators and provided
an accurate theoretical upper bound estimate, but its computed version
severely underestimated the exact error. The second one was an error estimator
based on a recovered displacement field and its performance was closely
related to the quality of the GFEM solution. This recovery technique
constructs a recovered solution on patches using a basis enriched with
handbook functions.

In order to obtain a recovered stress field that improves the accuracy of the
stresses obtained by XFEM, Xiao and Karihaloo (2006) proposed a moving least
squares fitting adapted to the XFEM framework which considers the use of
statically admissible basis functions. Nevertheless, the recovered stress
field was not used to obtain an error indicator.

Pannachet et al. (2009) worked on error estimation for mesh adaptivity in XFEM
for cohesive crack problems. Two error estimates were used, one based on the
error in the energy norm and another that considers the error in a local
quantity of interest. The error estimation was based on solving a series of
local problems with prescribed homogeneous boundary conditions.

Panetier et al. (2010) presented an extension to enriched approximations of
the constitutive relation error (CRE) technique already available to evaluate
error bounds in FEM (Ladevèze and Pelle, 2005). This procedure has been used
to obtain local error bounds on quantities of interest for XFEM problems.

Bordas and Duflot (2007) and Bordas et al. (2008) proposed a recovery based
error estimator for 2D and 3D XFEM approximations for cracks known as the
extended moving least squares (XMLS). This method intrinsically enriched an
MLS formulation to include information about the singular fields near the
crack tip. Additionally, it used a diffraction method to introduce the
discontinuity in the recovered field. This error estimator provided accurate
results with effectivity indices close to unity (optimal value) for 2D and 3D
fracture mechanics problems. Later, Duflot and Bordas (2008) proposed a global
derivative recovery formulation extended to XFEM problems. The recovered
solution was sought in a space spanned by the near tip strain fields obtained
from differentiating the Westergaard asymptotic expansion. Although the
results provided by this technique were not as accurate as those in Bordas and
Duflot (2007), Bordas et al. (2008), they were deemed by the authors to
require less computational power.

Ródenas et al. (2008) presented a modification of the superconvergent patch
recovery (SPR) technique tailored to the XFEM framework called
$\textrm{SPR}_{\textrm{XFEM}}$. This technique was based on three key
ingredients: (1) the use of a _singular_ +_smooth_ stress field decomposition
procedure around the crack tip similar to that described by Ródenas et al.
(2006) for FEM; (2) direct calculation of recovered stresses at integration
points using the partition of unity, and (3) use of different stress
interpolation polynomials at each side of the crack when a crack intersects a
patch. In order to obtain an equilibrated smooth recovered stress field, a
simplified version of the SPR-C technique presented by Ródenas et al. (2007)
was used. This simplified SPR-C imposed the fulfilment of the boundary
equilibrium equation at boundary nodes but did not impose the satisfaction of
the internal equilibrium equations.

The numerical results presented by Bordas and Duflot (2007), Bordas et al.
(2008) and Ródenas et al. (2008) showed promising accuracy for both the XMLS
and the $\textrm{SPR}_{\textrm{XFEM}}$ techniques. It is apparent in those
papers that $\textrm{SPR}_{\textrm{XFEM}}$ led to effectivity indices
remarkably close to unity. Yet, these papers do not shed any significant light
on the respective roles played by two important ingredients in those error
estimators, namely:

  1. 1.

The statistical admissibility of the recovered solution;

  2. 2.

The enrichment of the recovered solution.

Moreover, the differences in the test cases analysed and in the quality
measures considered in each of the papers makes it difficult to objectively
compare the merits of both methods.

The aim of this paper is to assess the role of statistical admissibility and
enrichment of the recovered solution in recovery based error estimation of
enriched finite element approximation for linear elastic fracture. To do so,
we perform a systematic study of the results obtained when considering the
different features of two error estimation techniques: XMLS and SPR-CX. SPR-CX
is an enhanced version of the $\textrm{SPR}_{\textrm{XFEM}}$ presented by
Ródenas et al. (2008) for 2D problems. Advantages and disadvantages of each
method are also provided.

The outline of the paper is as follows. In Section 2, the XFEM is briefly
presented. Section 3 deals with error estimation and quality assessment of the
solution. In Sections 4 and 5 we introduce the error estimators used in XFEM
approximations: the SPR-CX and the XMLS, respectively. Some numerical examples
analysing both techniques and the effect of the enrichment functions in the
recovery process are presented in Section 6. Finally some concluding remarks
are provided in Section 7.

##  2 Reference problem and XFEM solution

Let us consider a 2D linear elastic fracture mechanics (LEFM) problem on a
bounded domain $\Omega\subset\mathbb{R}^{2}$. The unknown displacement field
$\bm{\mathrm{u}}$ is the solution of the boundary value problem

$\displaystyle\nabla\cdot\boldsymbol{\upsigma}(\bm{\mathrm{u}})+\bm{\mathrm{b}}$
$\displaystyle=\mathbf{0}$ $\displaystyle\textrm{in }\Omega$ (1)
$\displaystyle\boldsymbol{\upsigma}(\bm{\mathrm{u}})\cdot\bm{\mathrm{n}}$
$\displaystyle=\bm{\mathrm{t}}$ $\displaystyle\textrm{on }\Gamma_{N}$ (2)
$\displaystyle\boldsymbol{\upsigma}(\bm{\mathrm{u}})\cdot\bm{\mathrm{n}}$
$\displaystyle=\mathbf{0}$ $\displaystyle\textrm{on }\Gamma_{C}$ (3)
$\displaystyle\bm{\mathrm{u}}$ $\displaystyle=\bar{\bm{\mathrm{u}}}$
$\displaystyle\textrm{on }\Gamma_{D}$ (4)

where $\Gamma_{N}$ and $\Gamma_{D}$ are the Neumann and Dirichlet boundaries
and $\Gamma_{C}$ represents the crack faces, with
$\partial\Omega=\Gamma_{N}\cup\Gamma_{D}\cup\Gamma_{C}$ and
$\Gamma_{N}\cap\Gamma_{D}\cap\Gamma_{C}=\varnothing$. $\bm{\mathrm{b}}$ are
the body forces per unit volume, $\bm{\mathrm{t}}$ the tractions applied on
$\Gamma_{N}$ (being $\bm{\mathrm{n}}$ the normal vector to the boundary) and
$\bar{\bm{\mathrm{u}}}$ the prescribed displacements on $\Gamma_{D}$. The weak
form of the problem reads: Find $\bm{\mathrm{u}}\in V$ such that:

$\forall\bm{\mathrm{v}}\in V\qquad
a(\bm{\mathrm{u}},\bm{\mathrm{v}})=l(\bm{\mathrm{v}}),$ (5)

where $V$ is the standard test space for elasticity problems such that
$V=\\{\bm{\mathrm{v}}\;|\;\bm{\mathrm{v}}\in
H^{1}(\Omega),\bm{\mathrm{v}}|_{\Gamma_{D}}(\bm{\mathrm{x}})=\mathbf{0},\bm{\mathrm{v}}\;\textrm{discontinuous
on }\Gamma_{C}\\}$, and

$\displaystyle a(\bm{\mathrm{u}},\bm{\mathrm{v}})$
$\displaystyle:=\int_{\Omega}\boldsymbol{\upsigma}(\bm{\mathrm{u}}):\bm{\mathrm{\epsilon}}(\bm{\mathrm{v}})d\Omega=\int_{\Omega}\boldsymbol{\upsigma}(\bm{\mathrm{u}}):\boldsymbol{\mathsf{S}}:\boldsymbol{\upsigma}(\bm{\mathrm{v}})d\Omega$
(6) $\displaystyle l(\bm{\mathrm{v}})$
$\displaystyle:=\int_{\Omega}\bm{\mathrm{b}}\cdot\bm{\mathrm{v}}d\Omega+\int_{\Gamma_{N}}\bm{\mathrm{t}}\cdot\bm{\mathrm{v}}d\Gamma,$
(7)

where $\boldsymbol{\mathsf{S}}$ is the compliance tensor,
$\boldsymbol{\upsigma}$ and $\bm{\mathrm{\epsilon}}$ represent the stress and
strain operators.

LEFM problems are denoted by the singularity present at the crack tip. The
following expressions represent the first term of the asymptotic expansion
which describes the displacements and stresses for combined loading modes I
and II in 2D. These expressions can be found in the literature (Szabó and
Babuška, 1991, Ródenas et al., 2008) and are reproduced here for completeness:

$\displaystyle u_{1}(r,\phi)$ $\displaystyle=\frac{K_{\rm
I}}{2\mu}\sqrt{\frac{r}{2\pi}}\cos\frac{\phi}{2}\left(\kappa-\cos\phi\right)+\frac{K_{\rm
II}}{2\mu}\sqrt{\frac{r}{2\pi}}\sin\frac{\phi}{2}\left(2+\kappa+\cos\phi\right)$
(8) $\displaystyle u_{2}(r,\phi)$ $\displaystyle=\frac{K_{\rm
I}}{2\mu}\sqrt{\frac{r}{2\pi}}\sin\frac{\phi}{2}\left(\kappa-\cos\phi\right)+\frac{K_{\rm
II}}{2\mu}\sqrt{\frac{r}{2\pi}}\cos\frac{\phi}{2}\left(2-\kappa-\cos\phi\right)$

$\displaystyle\sigma_{11}(r,\phi)$ $\displaystyle=\frac{K_{\rm I}}{\sqrt{2\pi
r}}\cos\frac{\phi}{2}\left(1-\sin\frac{\phi}{2}\sin\frac{3\phi}{2}\right)-\frac{K_{\rm
II}}{\sqrt{2\pi
r}}\sin\frac{\phi}{2}\left(2+\cos\frac{\phi}{2}\cos\frac{3\phi}{2}\right)$ (9)
$\displaystyle\sigma_{22}(r,\phi)$ $\displaystyle=\frac{K_{\rm I}}{\sqrt{2\pi
r}}\cos\frac{\phi}{2}\left(1+\sin\frac{\phi}{2}\sin\frac{3\phi}{2}\right)+\frac{K_{\rm
II}}{\sqrt{2\pi r}}\sin\frac{\phi}{2}\cos\frac{\phi}{2}\cos\frac{3\phi}{2}$
$\displaystyle\sigma_{12}(r,\phi)$ $\displaystyle=\frac{K_{\rm I}}{\sqrt{2\pi
r}}\sin\frac{\phi}{2}\cos\frac{\phi}{2}\cos\frac{3\phi}{2}+\frac{K_{\rm
II}}{\sqrt{2\pi
r}}\cos\frac{\phi}{2}\left(1-\sin\frac{\phi}{2}\sin\frac{3\phi}{2}\right)$

where $r$ and $\phi$ are the crack tip polar coordinates, $K_{\rm I}$ and
$K_{\rm II}$ are the stress intensity factors for modes I and II, $\mu$ is the
shear modulus, and $\kappa$ the Kolosov’s constant, defined in terms of the
parameters of material $E$ (Young’s modulus) and $\upsilon$ (Poisson’s ratio),
according to the expressions:

$\mu=\frac{E}{2\left(1+\upsilon\right)},\qquad\kappa=\left\\{\begin{array}[]{l
c r}{\displaystyle 3-4\upsilon\qquad\textrm{plane strain}}\\\ \vskip 6.0pt
plus 2.0pt minus
2.0pt\cr{\displaystyle\frac{3-\upsilon}{1+\upsilon}\,\qquad\textrm{plane
stress}}\end{array}\right.$

This type of problem is difficult to model using a standard FEM approximation
as the mesh needs to explicitly conform to the crack geometry. With the XFEM
the discontinuity of the displacement field along the crack faces is
introduced by adding degrees of freedom to the nodes of the elements
intersected by the crack. This tackles the problem of adjusting the mesh to
the geometry of the crack (Moës et al., 1999, Stolarska et al., 2001).
Additionally, to describe the solution around the crack tip, the numerical
model introduces a basis that spans the near tip asymptotic fields. The
following expression is generally used to interpolate the displacements at a
point of coordinates $\bm{\mathrm{x}}$ accounting for the presence of a crack
tip in a 2D XFEM approximation:

$\bm{\mathrm{u}}_{h}(\bm{\mathrm{x}})=\sum_{i\in\mathcal{I}}N_{i}(\bm{\mathrm{x}})\textbf{a}_{i}+\sum_{j\in\mathcal{J}}N_{j}(\bm{\mathrm{x}})H(\bm{\mathrm{x}})\textbf{b}_{j}+\sum_{m\in\mathcal{M}}N_{m}(\bm{\mathrm{x}})\left(\sum_{\ell=1}^{4}F_{\ell}(\bm{\mathrm{x}})\textbf{c}_{m}^{\ell}\right)$
(10)

where $N_{i}$ are the shape functions associated with node $i$,
$\textbf{a}_{i}$ represent the conventional nodal degrees of freedom,
$\textbf{b}_{j}$ are the coefficients associated with the discontinuous
enrichment functions, and $\bm{\mathrm{c}}_{m}$ those associated with the
functions spanning the asymptotic field. In the above equation, $\mathcal{I}$
is the set of all the nodes in the mesh, $\mathcal{M}$ is the subset of nodes
enriched with crack tip functions, and $\mathcal{J}$ is the subset of nodes
enriched with the discontinuous enrichment (see Figure 1). In (10), the
Heaviside function $H$, with unitary modulus and a change of sign on the crack
face, describes the displacement discontinuity if the finite element is
intersected by the crack. The $F_{\ell}$ are the set of branch functions used
to represent the asymptotic expansion of the displacement field around the
crack tip seen in (8). For the 2D case, the following functions are used
(Belytschko and Black, 1999):

$\left\\{F_{\ell}\left(r,\phi\right)\right\\}\equiv\sqrt{r}\left\\{\sin\frac{\phi}{2},\cos\frac{\phi}{2},\sin\frac{\phi}{2}\sin\phi,\cos\frac{\phi}{2}\sin\phi\right\\}$
(11)

Figure 1: Classification of nodes in XFEM. Fixed enrichment area of radius
$r_{e}$

The main features of the XFEM implementation considered to evaluate the
numerical results is described in detail in Ródenas et al. (2008) and can be
summarized as follows:

  * •

Use of bilinear quadrilaterals.

  * •

Decomposition of elements intersected by the crack into integration subdomains
that do not contain the crack. Alternatives which do not required this
subdivision are proposed by Ventura (2006), Natarajan et al. (2010).

  * •

Use of a quasi-polar integration with a $5\times 5$ quadrature rule in
triangular subdomains for elements containing the crack tip.

  * •

No correction for blending elements. Methods to address blending errors are
proposed by Chessa et al. (2003), Gracie et al. (2008), Fries (2008), Tarancón
et al. (2009).

###  2.1 Evaluation of stress intensity factors

The stress intensity factors (SIFs) in LEFM represent the amplitude of the
singular stress fields and are key quantities of interest to simulate crack
growth in LEFM. Several post–processing methods, following local or global
(energy) approaches, are commonly used to extract SIFs (Banks-Sills, 1991) or
to calculate the energy release rate $G$. Energy or global methods are
considered to be the most accurate and efficient methods (Banks-Sills, 1991,
Li et al., 1985). Global methods based on the equivalent domain integral (EDI
methods) are specially well-suited for FEM and XFEM analyses as they are easy
to implement and can use information far from the singularity. In this paper,
the interaction integral as described in Shih and Asaro (1988), Yau et al.
(1980) has been used to extract the SIFs. This technique provides $K_{\rm I}$
and $K_{\rm II}$ for problems under mixed mode loading conditions using
auxiliary fields. Details on the implementation of the interaction integral
can be found for example in Moës et al. (1999), Ródenas et al. (2008), Giner
et al. (2005).

##  3 Error estimation in the energy norm

The approximate nature of the FEM and XFEM approximations implies a
discretization error which can be quantified using the error in the energy
norm for the solution
$\lVert\bm{\mathrm{e}}\rVert=\lVert\bm{\mathrm{u}}-\bm{\mathrm{u}}^{h}\rVert$.
To obtain an estimate of the discretization error
$\lVert\bm{\mathrm{e}}_{es}\rVert$, in the context of elasticity problems
solved using the FEM, the expression for the ZZ estimator is defined as
(Zienkiewicz and Zhu, 1987) (in matrix form):

$\lVert\bm{\mathrm{e}}\rVert^{2}\approx\lVert\bm{\mathrm{e}}_{es}\rVert^{2}=\int_{\Omega}\left(\bm{\mathrm{\sigma}}^{*}-\bm{\mathrm{\sigma}}^{h}\right)^{T}\bm{\mathrm{D}}^{-1}\left(\bm{\mathrm{\sigma}}^{*}-\bm{\mathrm{\sigma}}^{h}\right)d\Omega$
(12)

or alternatively for the strains

$\lVert\bm{\mathrm{e}}_{es}\rVert^{2}=\int_{\Omega}\left(\bm{\mathrm{\varepsilon}}^{*}-\bm{\mathrm{\varepsilon}}^{h}\right)^{T}\bm{\mathrm{D}}\left(\bm{\mathrm{\varepsilon}}^{*}-\bm{\mathrm{\varepsilon}}^{h}\right)d\Omega,$

where the domain $\Omega$ refers to the full domain of the problem or a local
subdomain (element), $\bm{\mathrm{\sigma}}^{h}$ represents the stress field
provided by the FEM, $\bm{\mathrm{\sigma}}^{*}$ is the recovered stress field,
which is a better approximation to the exact solution than
$\bm{\mathrm{\sigma}}^{h}$ and $\bm{\mathrm{D}}$ is the elasticity matrix of
the constitutive relation
$\bm{\mathrm{\sigma}}=\bm{\mathrm{D}}\bm{\mathrm{\varepsilon}}$.

The recovered stress field $\bm{\mathrm{\sigma}}^{*}$ is usually interpolated
in each element using the shape functions $\bm{\mathrm{N}}$ of the underlying
FE approximation and the values of the recovered stress field calculated at
the nodes $\bar{\bm{\mathrm{\sigma}}}^{*}$

$\bm{\mathrm{\sigma}}^{*}(\bm{\mathrm{x}})=\sum_{i=1}^{n_{e}}N_{i}(\bm{\mathrm{x}})\bar{\bm{\mathrm{\sigma}}}^{*}_{i}(\bm{\mathrm{x}}_{i}),$
(13)

where $n_{e}$ is the number of nodes in the element under consideration and
$\bar{\bm{\mathrm{\sigma}}}^{*}_{i}(\bm{\mathrm{x}}_{i})$ are the stresses
provided by the least squares technique at node $i$. The components of
$\bar{\bm{\mathrm{\sigma}}}^{*}_{i}$ are obtained using a polynomial
expansion, $\bar{\sigma}_{i,j}=\bm{\mathrm{p}}\bm{\mathrm{a}}$ (with
$j=xx,yy,xy$), defined over a set of contiguous elements connected to node $i$
called patch, where $\bm{\mathrm{p}}$ is the polynomial basis and
$\bm{\mathrm{a}}$ are the unknown coefficients.

The ZZ error estimator is asymptotically exact if the recovered solution used
in the error estimation converges at a higher rate than the finite element
solution (Zienkiewicz and Zhu, 1992b).

In this paper, we are interested in the role played by statical admissibility
and enrichment of the recovered solution in estimating the error committed by
XFEM. To do so, we study the performance of two recovery-based error
estimators, which exhibit different features, that have been recently
developed for XFEM:

  * •

The SPR-CX derived from the error estimator developed by Ródenas et al. (2008)
and summarized in Section 4 (two other versions of the SPR-CX technique have
been also considered);

  * •

The XMLS proposed by Bordas and Duflot (2007) and summarized in Section 5.

Both estimators provide a recovered stress field in order to evaluate the
estimated error in the energy norm by means of the expression shown in (12).

##  4 SPR-CX error estimator

The SPR-CX error estimator is an enhancement of the error estimator first
introduced by Ródenas et al. (2008), which incorporates the ideas proposed in
Ródenas et al. (2007) to guarantee the exact satisfaction of the equilibrium
locally on patches. In Ródenas et al. (2008) a set of key ideas are proposed
to modify the standard SPR by Zienkiewicz and Zhu (1992a), allowing its use
for singular problems.

The recovered stresses $\bm{\mathrm{\sigma}}^{*}$ are directly evaluated at an
integration point $\bm{\mathrm{x}}$ through the use of a partition of unity
procedure, properly weighting the stress interpolation polynomials obtained
from the different patches formed at the vertex nodes of the element
containing $\bm{\mathrm{x}}$:

$\bm{\mathrm{\sigma}}^{*}(\bm{\mathrm{x}})=\sum_{i=1}^{n_{v}}N_{i}(\bm{\mathrm{x}})\bm{\mathrm{\sigma}}^{*}_{i}(\bm{\mathrm{x}}),$
(14)

where $N_{i}$ are the shape functions associated to the vertex nodes $n_{v}$.

One major modification is the introduction of a splitting procedure to perform
the recovery. For singular problems the exact stress field
$\bm{\mathrm{\sigma}}$ is decomposed into two stress fields, a smooth field
$\bm{\mathrm{\sigma}}_{smo}$ and a singular field
$\bm{\mathrm{\sigma}}_{sing}$:

$\bm{\mathrm{\sigma}}=\bm{\mathrm{\sigma}}_{smo}+\bm{\mathrm{\sigma}}_{sing}.$
(15)

Then, the recovered field $\bm{\mathrm{\sigma}}^{*}$ required to compute the
error estimate given in (12) can be expressed as the contribution of two
recovered stress fields, one smooth $\bm{\mathrm{\sigma}}^{*}_{smo}$ and one
singular $\bm{\mathrm{\sigma}}^{*}_{sing}$:

$\bm{\mathrm{\sigma}}^{*}=\bm{\mathrm{\sigma}}^{*}_{smo}+\bm{\mathrm{\sigma}}^{*}_{sing}.$
(16)

For the recovery of the singular part, the expressions in (9) which describe
the asymptotic fields near the crack tip are used. To evaluate
$\bm{\mathrm{\sigma}}^{*}_{sing}$ from (9) we first obtain estimated values of
the stress intensity factors $K_{\rm I}$ and $K_{\rm II}$ using the
interaction integral as indicated in Section 2.1. The recovered part
$\bm{\mathrm{\sigma}}^{*}_{sing}$ is an equilibrated field as it satisfies the
internal equilibrium equations.

Once the field $\bm{\mathrm{\sigma}}^{*}_{sing}$ has been evaluated, an FE
approximation to the smooth part $\bm{\mathrm{\sigma}}^{h}_{smo}$ can be
obtained subtracting $\bm{\mathrm{\sigma}}^{*}_{sing}$ from the raw FE field:

$\bm{\mathrm{\sigma}}^{h}_{smo}=\bm{\mathrm{\sigma}}^{h}-\bm{\mathrm{\sigma}}^{*}_{sing}.$
(17)

Then, the field $\bm{\mathrm{\sigma}}^{*}_{smo}$ is evaluated applying an
SPR-C recovery procedure over the field $\bm{\mathrm{\sigma}}^{h}_{smo}$.

For patches intersected by the crack, the recovery technique uses different
stress interpolation polynomials on each side of the crack. This way it can
represent the discontinuity of the recovered stress field along the crack
faces, which is not the case for SPR that smoothes out the discontinuity (see,
e.g. Bordas and Duflot (2007)).

In order to obtain an equilibrated recovered stress field
$\bm{\mathrm{\sigma}}^{*}_{smo}$, the SPR-CX enforces the fulfilment of the
equilibrium equations locally on each patch. The constraint equations are
introduced via Lagrange multipliers into the linear system used to solve for
the coefficients of the polynomial expansion of the recovered stresses on each
patch. These include the satisfaction of the:

  * •

Internal equilibrium equations.

  * •

Boundary equilibrium equation: A point collocation approach is used to impose
the satisfaction of a second order approximation to the tractions along the
Neumann boundary.

  * •

Compatibility equation: This additional constraint is also imposed to further
increase the accuracy of the recovered stress field.

To evaluate the recovered field, quadratic polynomials have been used in the
patches along the boundary and crack faces, and linear polynomials for the
remaining patches. As more information about the solution is available along
the boundary, polynomials one degree higher are useful to improve the quality
of the recovered stress field.

The enforcement of equilibrium equations provides an equilibrated recovered
stress field locally on patches. However, the process used to obtain a
continuous field $\bm{\mathrm{\sigma}}^{*}$ shown in (14) introduces a small
lack of equilibrium as explained in Ródenas et al. (2010a). The reader is
referred to Ródenas et al. (2010a), Díez et al. (2007) for details.

##  5 XMLS error estimator

In the XMLS the solution is recovered through the use of the _moving least
squares_ (MLS) technique, developed by mathematicians to build and fit
surfaces. The XMLS technique extends the work of Tabbara et al. (1994) for FEM
to enriched approximations. The general idea of the XMLS is to use the
displacement solution provided by XFEM to obtain a recovered strain field
(Bordas and Duflot, 2007, Bordas et al., 2008). The smoothed strains are
recovered from the derivative of the MLS-smoothed XFEM displacement field:

$\displaystyle\bm{\mathrm{u}}^{*}(\bm{\mathrm{x}})$
$\displaystyle=\sum_{i\in\mathcal{N}_{x}}\psi_{i}(\bm{\mathrm{x}})\bm{\mathrm{u}}_{i}^{h}$
(18) $\displaystyle\bm{\mathrm{\varepsilon}}^{*}(\bm{\mathrm{x}})$
$\displaystyle=\sum_{i\in\mathcal{N}_{x}}\nabla_{s}(\psi_{i}(\bm{\mathrm{x}})\bm{\mathrm{u}}_{i}^{h}),$
(19)

where the $\bm{\mathrm{u}}_{i}^{h}$ are the raw nodal XFEM displacements and
$\psi_{i}(\bm{\mathrm{x}})$ are the MLS shape function values associated with
a node $i$ at a point $\bm{\mathrm{x}}$. $\mathcal{N}_{x}$ is the set of
$n_{x}$ nodes in the domain of influence of point $\bm{\mathrm{x}}$,
$\nabla_{s}$ is the symmetric gradient operator, $\bm{\mathrm{u}}^{*}$ and
$\bm{\mathrm{\varepsilon}}^{*}$ are the recovered displacement and strain
fields respectively. At each point $\bm{\mathrm{x}}_{i}$ the MLS shape
functions $\psi_{i}$ are evaluated using weighting functions $\omega_{i}$ and
an enriched basis $\bm{\mathrm{p}}(\bm{\mathrm{x}}_{i})$. The total $n_{x}$
non-zero MLS shape functions at point $\bm{\mathrm{x}}$ are evaluated as:

$(\psi_{i}(\bm{\mathrm{x}}))_{1\leq i\leq
n_{x}}=(\bm{\mathrm{A}}^{-1}(\bm{\mathrm{x}})\bm{\mathrm{p}}(\bm{\mathrm{x}}))^{T}\bm{\mathrm{p}}(\bm{\mathrm{x}}_{i})\omega_{i}(\bm{\mathrm{x}})$
(20)

where
$\bm{\mathrm{A}}(\bm{\mathrm{x}})=\sum_{i=1}^{n_{x}}\omega_{i}(\bm{\mathrm{x}})\bm{\mathrm{p}}(\bm{\mathrm{x}}_{i})\bm{\mathrm{p}}(\bm{\mathrm{x}}_{i})^{T}$
is a matrix to be inverted at every point $\bm{\mathrm{x}}$ (see Bordas and
Duflot (2007), Bordas et al. (2008) for further details).

For each supporting point $\bm{\mathrm{x}}_{i}$, the weighting function
$\omega_{i}$ is defined such that:

$\omega_{i}(s)=f_{4}(s)=\begin{cases}1-6s^{2}+8s^{3}-3s^{4}&\text{if
}\left|s\right|\leq 1\\\ 0&\text{if }\left|s\right|>1\end{cases}$ (21)

where $s$ is the normalized distance between the supporting point
$\bm{\mathrm{x}}_{i}$ and a point $\bm{\mathrm{x}}$ in the computational
domain. In order to describe the discontinuity, the distance $s$ in the weight
function defined for each supporting point is modified using the _diffraction
criterion_ (Belytschko et al., 1996). The basic idea of this criterion is
depicted in Figure 2. The weight function is continuous except across the
crack faces since the points at the other side of the crack are not considered
as part of the support. Near the crack tip the weight of a node $i$ over a
point of coordinates $\bm{\mathrm{x}}$ diminishes as the crack hides the
point. When the point $\bm{\mathrm{x}}$ is hidden by the crack the following
expression is used:

$s=\frac{\left\|\bm{\mathrm{x}}-\bm{\mathrm{x}}_{C}\right\|+\left\|\bm{\mathrm{x}}_{C}-\bm{\mathrm{x}}_{i}\right\|}{d_{i}}$
(22)

where $d_{i}$ is the radius of the support.

Figure 2: Diffraction criteria to introduce the discontinuity in the XMLS
approximation.

The MLS shape functions can reproduce any function in their basis. The basis
$\bm{\mathrm{p}}$ used is a linear basis enriched with the functions that
describe the first order asymptotic expansion at the crack tip as indicated in
(11):

$\bm{\mathrm{p}}=\left[1,x,y,\left[F_{1}(r,\phi),F_{2}(r,\phi),F_{3}(r,\phi),F_{4}(r,\phi)\right]\right]$
(23)

Note that although the enriched basis can reproduce the singular behaviour of
the solution around the crack tip, the resulting recovered field not
necessarily satisfies the equilibrium equations.

##  6 Numerical results

In this section, numerical experiments are performed to verify the behaviour
of both XFEM recovery based error estimators considered in this paper.

Babuška et al. (1994a, b, 1997) proposed a robustness patch test for quality
assessment of error estimators. However, the use of this test is not within
the scope of this paper and furthermore, to the authors’ knowledge, it has not
been used in the context of XFEM. Therefore, the more traditional approach of
using benchmark problems is considered here to analyse the response of the
different error estimators.

The accuracy of the error estimators is evaluated both locally and globally.
This evaluation has been based on the effectivity of the error in the energy
norm, which is quantified using the _effectivity index_ $\theta$ defined as:

$\theta=\frac{\lVert\textbf{e}_{es}\rVert}{\lVert\textbf{e}\rVert}.$ (24)

When the enhanced or recovered solution is close to the analytical solution
the effectivity approaches the theoretical value of 1, which indicates that it
is a good error estimator, i.e. the approximate error is close to the exact
error.

To assess the quality of the estimator at a local level, the local effectivity
$D$, inspired on the _robustness index_ found in Babuška et al. (1994a), is
used. For each element $e$, $D$ represents the variation of the effectivity
index in this element, $\theta^{e}$, with respect to the theoretical value
(the error estimator can be considered to be of good quality if it yields $D$
values close to zero). $D$ is defined according to the following expression,
where superscripts e indicate the element $e$:

$\begin{array}[]{c}{\displaystyle D=\theta^{e}-1\qquad{\rm
if}\qquad\theta^{e}\geq 1}\\\ {\displaystyle
D=1-\frac{1}{\theta^{e}}\qquad{\rm
if}\qquad\theta^{e}<1}\end{array}\qquad\qquad{\rm
with}\qquad\theta^{e}=\frac{\left\|\textbf{e}_{es}^{e}\right\|}{\left\|\textbf{e}^{e}\right\|}.$
(25)

To evaluate the overall quality of the error estimator we use the global
effectivity index $\theta$, the mean value $m\left(|D|\right)$ and the
standard deviation $\sigma\left(D\right)$ of the local effectivity. A good
quality error estimator yields values of $\theta$ close to one and values of
$m\left(|D|\right)$ and $\sigma\left(D\right)$ close to zero.

The techniques can be used in practical applications. However, in order to
properly compare their performance we have used an academic problem with exact
solution. In the analysis we solve the Westergaard problem (Gdoutos, 1993) as
it is one of the few problems in LEFM under mixed mode with an analytical
solution. In Giner et al. (2005), Ródenas et al. (2008) we can find explicit
expressions for the stress fields in terms of the spatial coordinates. In the
next subsection we show a description of the Westergaard problem and XFEM
model, taken from Ródenas et al. (2008) and reproduced here for completeness.

###  6.1 Westergaard problem and XFEM model

The Westergaard problem corresponds to an infinite plate loaded at infinity
with biaxial tractions $\sigma_{x\infty}=\sigma_{y\infty}=\sigma_{\infty}$ and
shear traction $\tau_{\infty}$, presenting a crack of length $2a$ as shown in
Figure 3. Combining the externally applied loads we can obtain different
loading conditions: pure mode I, II or mixed mode.

Figure 3: Westergaard problem. Infinite plate with a crack of length $2a$
under uniform tractions $\sigma_{\infty}$ (biaxial) and $\tau_{\infty}$.
Finite portion of the domain $\Omega_{0}$, modelled with FE.

In the numerical model only a finite portion of the domain ($a=1$ and $b=4$ in
Figure 3) is considered. The projection of the stress distribution
corresponding to the analytical Westergaard solution for modes I and II, given
by the expressions below, is applied to its boundary:

$\begin{array}[]{r@{\hspace{1ex}}c@{\hspace{1ex}}l}{\sigma_{x}^{I}}(x,y)\hskip
4.30554pt&{=}\hfil\hskip
4.30554pt&{\displaystyle\frac{\sigma_{\infty}}{\sqrt{\left|t\right|}}\bigg{[}\left(x\cos\frac{\phi}{2}-y\sin\frac{\phi}{2}\right)+y\frac{a^{2}}{\left|t\right|^{2}}\left(m\sin\frac{\phi}{2}-n\cos\frac{\phi}{2}\right)\bigg{]}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr{\sigma_{y}^{I}}(x,y)\hskip
4.30554pt&{=}\hfil\hskip
4.30554pt&{\displaystyle\frac{\sigma_{\infty}}{\sqrt{\left|t\right|}}\bigg{[}\left(x\cos\frac{\phi}{2}-y\sin\frac{\phi}{2}\right)-y\frac{a^{2}}{\left|t\right|^{2}}\left(m\sin\frac{\phi}{2}-n\cos\frac{\phi}{2}\right)\bigg{]}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr{\tau_{xy}^{I}}(x,y)\hskip
4.30554pt&{=}\hfil\hskip 4.30554pt&{\displaystyle
y\frac{a^{2}\sigma_{\infty}}{\left|t\right|^{2}\sqrt{\left|t\right|}}\left(m\cos\frac{\phi}{2}+n\sin\frac{\phi}{2}\right)}\end{array}$
(26)

$\begin{array}[]{r@{\hspace{1ex}}c@{\hspace{1ex}}l}{\sigma_{x}^{II}}(x,y)\hskip
4.30554pt&{=}\hfil\hskip
4.30554pt&{\displaystyle\frac{\tau_{\infty}}{\sqrt{\left|t\right|}}\bigg{[}2\left(y\cos\frac{\phi}{2}+x\sin\frac{\phi}{2}\right)-y\frac{a^{2}}{\left|t\right|^{2}}\left(m\cos\frac{\phi}{2}+n\sin\frac{\phi}{2}\right)\bigg{]}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr{\sigma_{y}^{II}}(x,y)\hskip
4.30554pt&{=}\hfil\hskip 4.30554pt&{\displaystyle
y\frac{a^{2}\tau_{\infty}}{\left|t\right|^{2}\sqrt{\left|t\right|}}\left(m\cos\frac{\phi}{2}+n\sin\frac{\phi}{2}\right)}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr{\tau_{xy}^{II}}(x,y)\hskip
4.30554pt&{=}\hfil\hskip
4.30554pt&{\displaystyle\frac{\tau_{\infty}}{\sqrt{\left|t\right|}}\bigg{[}\left(x\cos\frac{\phi}{2}-y\sin\frac{\phi}{2}\right)+y\frac{a^{2}}{\left|t\right|^{2}}\left(m\sin\frac{\phi}{2}-n\cos\frac{\phi}{2}\right)\bigg{]}}\end{array}$
(27)

where the stress fields are expressed as a function of $x$ and $y$, with
origin at the centre of the crack. The parameters $t$, $m$, $n$ and $\phi$ are
defined as

$\begin{split}t&=(x+iy)^{2}-a^{2}=(x^{2}-y^{2}-a^{2})+i(2xy)=m+in\\\
m&=\textrm{Re}(t)=\textrm{Re}(z^{2}-a^{2})=x^{2}-y^{2}-a^{2}\\\
n&=\textrm{Im}(t)=(z^{2}-a^{2})=2xy\\\
\phi&=\textrm{Arg}(\bar{t})=\textrm{Arg}(m-in)\qquad\textrm{with
}\phi\in\left[-\pi,\pi\right],\;i^{2}=-1\end{split}$ (28)

For the problem analysed, the exact value of the SIF is given as

$K_{{\rm I},ex}=\sigma\sqrt{\pi a}\qquad\qquad K_{{\rm II},ex}=\tau\sqrt{\pi
a}$ (29)

Three different loading configurations corresponding to the _pure mode I_
($\sigma_{\infty}=100,\;\tau_{\infty}=0$), _pure mode II_
($\sigma_{\infty}=0,\;\tau_{\infty}=100$), and _mixed mode_
($\sigma_{\infty}=30,\;\tau_{\infty}=90$) cases of the Westergaard problem are
considered. The geometric models and boundary conditions are shown in Figure
4.

Figure 4: Model for an infinite plate with a crack subjected to biaxial
tractions $\sigma_{\infty},\;\tau_{\infty}$ in the infinite.

Bilinear elements are used in the models, with a _singular_ +_smooth_
decomposition area of radius $\rho=0.5$ equal to the radius $r_{e}$ used for
the fixed enrichment area (note that the splitting radius $\rho$ should be
greater or equal to the enrichment radius $r_{e}$, Ródenas et al. (2008)). The
radius for the Plateau function used to extract the SIF is $r_{q}=0.9$.
Young’s modulus is $E=10^{7}$ and Poisson’s ratio $\upsilon=0.333$. For the
numerical integration of standard elements we use a $2\times 2$ Gaussian
quadrature rule. For split elements we use 7 Gauss points in each triangular
subdomain, and a $5\times 5$ quasipolar integration (Béchet et al., 2005) in
the subdomains of the element containing the crack tip.

Regarding global error estimation, the evolution of global parameters in
sequences of uniformly refined structured (Figure 5) and unstructured (Figure
6) meshes is studied. In the first case, the mesh sequence is defined so that
the crack tip always coincided with a node. For a more general scenario, this
condition has not been applied for the unstructured meshes.

Figure 5: Sequence of structured meshes. Figure 6: Sequence of unstructured
meshes.

###  6.2 Mode I and structured meshes

The first set of results presented in this study are for the Westergaard
problem under mode I load conditions. Figure 7 shows the values for the local
effectivity index $D$ for the first mesh in the sequence analysed. In the
figure, the size of the enrichment area is denoted by a circle. It can be
observed that far from the enrichment area both techniques yield similar
results, however, the XMLS technique exhibits a higher overestimation of the
error close to the singularity.

Figure 7: Mode I, structured mesh 1. Local effectivity index $D$ (ideal value
$D=0$). The circles denote the enriched zones around the crack tips. The
superiority of the statically admissible recovery (SPR-CX) compared to the
standard XMLS is clear.

The same behaviour is also observed for more refined meshes in the sequence.
Figure 8 shows a zoom in the enriched area of a finer mesh where, as before, a
higher overestimation of the error can be observed for the results obtained
with the XMLS error estimator. In this example and in order to study the
evolution of the accuracy of the recovered stress field when considering
different features in the recovery process, we consider two additional
versions:

  * •

The first considers a recovery procedure which enforces both internal and
boundary equilibrium, but does not include the singular+smooth splitting
technique (SPR-C). This approach is very similar in form to conventional SPR-
based recovery techniques widely used in FEM. It is well-known that this type
of recovery process will produce unreliable results when the stress field
contains a singularity because the polynomial representation of the recovered
stresses is not able to describe the singular field (e.g. Bordas and Duflot
(2007)).

  * •

The second version of the error estimator performs the splitting but does not
equilibrate the recovered field (SPR-X). As previously commented,
$\bm{\mathrm{\sigma}}^{*}_{sing}$ is an equilibrated field but
$\bm{\mathrm{\sigma}}^{*}_{smo}$ is not equilibrated. The aim of this second
version is to assess the influence of enforcing the equilibrium constraints.

Table 1 summarizes the main features of the different recovery procedures we
considered.

| Singular | Equilibrated sin- | Locally equili-
---|---|---|---
| functions | gular functions | brated field
SPR-X | yes | yes | no
SPR-C | no | not applicable | yes
SPR-CX | yes | yes | yes
XMLS | yes | no | no
Table 1: Comparison of features for the different recovery procedures Figure
8: Mode I, structured mesh 4. Local effectivity index $D$ (ideal value $D=0$).
The circles denote the enriched zones around the crack tips. The results show
the need for the inclusion of the near-tip fields in the recovery process
(SPR-X is far superior to SPR-C). It is also clear that enforcing statical
admissibility (SPR-CX) greatly improves the accuracy along the crack faces and
the Neumann boundaries. Also notice that XMLS leads to slightly improved
effectivities compared to SPR-CX around the crack faces, between the enriched
zone and the left boundary. SPR-C is clearly not able to reproduce the
singular fields, which is shown by large values of $D$ inside the enriched
region.

The performance of SPR-X in Figure 8 is poor, especially along the Neumann
boundaries, where the equilibrium equations are not enforced, as they are in
the SPR-CX. Particularly interesting are the results for the SPR-C technique,
where considerable overestimation of the exact error is observed _in the whole
enriched region_. This is due to the inability of polynomial functions to
reproduce the near-tip fields, even relatively far from the crack tip.

The above compared the relative performance of the methods locally. It is also
useful to have a measure of the global accuracy of the different error
estimators, Figure 9(left) shows the convergence of the estimated error in the
energy norm $\lVert\textbf{e}_{es}\rVert$ (see equation (12)) evaluated with
the proposed techniques, the convergence of the exact error
$\lVert\textbf{e}\rVert$111The exact error measures the error of the raw XFEM
compared to the exact solution. The theoretical convergence rate is
$\mathcal{O}(h^{p})$, $h$ being the element size and $p$ the polynomial degree
(1 for linear elements,…), or $\mathcal{O}(-\text{dof}^{p/2})$ if we consider
the number of degrees of freedom. is shown for comparison. It can be observed
that SPR-CX leads to the best approximation to the exact error in the energy
norm. The best approach is that which captures as closely as possible the
exact error. Non-singular formulations such as SPR-C greatly overestimate the
error, especially close to the singular point. Only SPR-C suffers from the
characteristic suboptimal convergence rate (factor 1/2, associated with the
strength of the singularity)

In order to evaluate the quality of the recovered field
$\bm{\mathrm{\sigma}}^{*}$ obtained with each of the recovery techniques, the
convergence in energy norm of the approximate error on the recovered field
$\lVert\textbf{e}^{*}\rVert=\lVert\bm{\mathrm{u}}-\bm{\mathrm{u}}^{*}\rVert$
is compared for the error estimators in Figure 9(right). It can be seen that
the XMLS, the SPR-CX and the SPR-X provide convergence rates higher than the
convergence rate for the exact error $\lVert\textbf{e}\rVert$, which is an
indicator of the asymptotic exactness of the error estimators based on these
recovery techniques (Zienkiewicz and Zhu, 1992b).

Figure 9: Mode I, structured meshes: Convergence of the estimated error in the
energy norm $\|\mathbf{e}_{es}\|$ (left). All methods which include the near-
tip fields do converge with the optimal convergence rate of 1/2. The best
method is SPR-CX. XMLS performs worse than SPR-X for coarse meshes, but it
becomes equivalent to SPR-CX for fine meshes whilst the SPR-X error increases
slightly. Convergence of the error for the recovered field
$\|\mathbf{e}_{es}^{*}\|$ (right). For the recovered field, SPR-CX is the best
method closely followed by SPR-X. XMLS errors are half an order of magnitude
larger than SPR-CX/SPR-X, but superior to the non-enriched recovery method
SPR-C. For all methods, except SPR-C, the error on the recovered solution
converges faster than the error on the raw solution, which shows that the
estimators are asymptotically exact.

Figure 10 shows the evolution of the parameters $\theta$, $m\left(|D|\right)$
y $\sigma\left(D\right)$ with respect to the number of degrees of freedom for
the error estimators. It can be verified that although the SPR-C technique
includes the satisfaction of the equilibrium equations, the effectivity of the
error estimator does not converge to the theoretical value ($\theta=1$). This
is due to the absence of the near-tip fields in the recovered solution. The
influence of introducing singular functions in the recovery process can be
observed in the results provided by SPR-X. In this case, the convergence is
obtained both locally and globally. This shows the importance of introducing a
description of the singular field in the recovery process, i.e. of using
_extended recovery techniques_ in XFEM.

Figure 10: Global indicators $\theta$, $m\left(|D|\right)$ and
$\sigma\left(D\right)$ for mode I and structured meshes.

Regarding the enforcement of the equilibrium equations, we can see that an
equilibrated formulation (SPR-CX) leads to better effectivities compared to
other non-equilibrated configurations (SPR-X, XMLS), which is an indication of
the advantages associated with equilibrated recoveries in this context.

Still in Figure 10, the results for the SPR-CX and XMLS show that the values
of the global effectivity index are close to (and less than) unity for the
SPR-CX estimator, whilst the XMLS is a lot less effective and shows
effectivities larger than 1.

The next step in our analysis is the verification of the convergence of the
mean value and standard deviation of the local effectivity index $D$. From a
practical perspective, one would want the local effectivity to be equally
good, on average, everywhere inside the domain; one would also want this
property to improve with mesh refinement. In other words, the average local
effectivity $m$ and its standard deviation $\sigma$ should decrease with mesh
refinement: $m$ and $\sigma$ measure the average behaviour of the method and
how far the results deviate from the mean, i.e. how spatially consistent they
are. The idea is to spot areas where there may be compensation between
overestimated and underestimated areas that could produce an apparently
’accurate’ global estimation. It can be confirmed that, although the curves
for $m(|D|)$ and $\sigma(D)$ for both estimators tend to zero with mesh
refinement, the SPR-CX is clearly superior to the XMLS. These results can be
explained by the fact that the SPR-CX enforces the fulfilment of the
equilibrium equations in patches, and evaluates the singular part of the
recovered field using the equilibrated expression that represents the first
term of the asymptotic expansion, whereas the XMLS enrichment uses only a set
of singular functions that would be able to reproduce this term but are not
necessarily equilibrated.

Although the results for the SPR-CX proved to be more accurate, the XMLS-type
techniques will prove useful to obtain error bounds. There is an increasing
interest in evaluating upper and lower bounds of the error for XFEM
approximations. Some work has been already done to obtain upper bounds using
recovery techniques as indicated in Ródenas et al. (2010a), where the authors
showed that an upper bound can be obtained if the recovered field is
continuous and equilibrated. However, they demonstrated that due to the use of
the conjoint polynomials process to enforce the continuity of
$\bm{\mathrm{\sigma}}^{*}$, some residuals in the equilibrium equations
appear, and the recovered field is only nearly equilibrated. Then, correction
terms have to be evaluated to obtain the upper bound. An equilibrated version
of the XMLS could provide a recovered stress field which would be continuous
and equilibrated, thus facilitating the evaluation of the upper bound. Initial
results for this class of techniques can be found in Ródenas et al. (2009,
2010b).

In this first example we have shown the effect of using the SPR-C and SPR-X
recovery techniques in the error estimation. Considering that these two
techniques are special cases of the SPR-CX with inferior results, for further
examples we will focus only on the SPR-CX and XMLS techniques.

###  6.3 Mode II and structured meshes

Figure 11 presents the results considering mode II loading conditions for the
local effectivity index $D$ on the fourth mesh of the sequence. Similarly to
the results for mode I, the XMLS estimator presents a higher overestimation of
the error near the enrichment area. This same behaviour is observed for the
whole set of structured meshes analysed under pure mode II.

Figure 11: Mode II, structured mesh 4. Local effectivity index $D$ (ideal
value $D=0$). SPR-CX leads to much better local effectivities than XMLS. It is
also remarkable that the effectivities are clearly worse in the enriched
region (circle) for both the XMLS and the SPR-CX, and that this effect is more
pronounced in the former method. Note the slightly worse results obtained by
SPR-CX around the crack faces between the boundary of the enriched region and
the left boundary, as in mode I.

As for the mode I case, the evolution of global accuracy parameters in mode II
exhibits the same behaviour seen for mode I loading conditions. Figure 12
shows the results for the convergence of the error in the energy norm. Figure
13 shows the evolution of global indicators $\theta$, $m\left(|D|\right)$ and
$\sigma\left(D\right)$. Again, SPR-CX provides the best results for the
Westergaard problem under pure mode II, using structured meshes.

Figure 12: Mode II, structured meshes: Convergence of the estimated error in
the energy norm $\|\mathbf{e}_{es}\|$ (left). The convergence rates are very
similar to those obtained in the mode I case, and, SPR-CX is still the most
accurate method, i.e. that for which the estimated error is closest to the
exact error. Convergence of the error in the recovered field
$\|\mathbf{e}_{es}^{*}\|$ (right). It can be noticed that the SPR-CX error
still converges faster than the exact error, but with a lower convergence rate
(0.59 versus 0.73) than in the mode I case. There is no such difference for
the XMLS results, which converge at practically the same rate (0.71 versus
0.72). The error level difference between SPR-CX and XMLS is about half an
order of magnitude, as in the mode I case.  Figure 13: Global indicators
$\theta$, $m\left(|D|\right)$ and $\sigma\left(D\right)$ for mode II and
structured meshes. The results are qualitatively and quantitatively similar to
those obtained in mode I.

###  6.4 Mixed mode and unstructured meshes

Considering a more general problem, Figure 14 shows the local effectivity
index $D$ for the fourth mesh in a sequence of unstructured meshes, having a
number of degrees of freedom (dof) similar to the mesh represented previously
for load modes I and II. The XMLS recovered field presents higher
overestimation of the error around the crack tip which corroborates the
results found in previous load cases.

Figure 14: Mixed mode, unstructured mesh 4. Local effectivity index $D$ (ideal
value $D=0$). Note the improved results compared to the structured case with
$D$ values ranging from -1.5 to 1.5 as opposed to -4 to 4.

Figures 15 and 16 represent the evolution of global parameters for
unstructured meshes and mixed mode. Once more, the best results for the error
estimates are obtained using the SPR-CX technique.

Figure 15:  Mixed mode, unstructured meshes: Convergence of the estimated
error in the energy norm $\|\mathbf{e}_{es}\|$ (left). Convergence of the
error for the recovered field $\|\mathbf{e}_{es}^{*}\|$ (right). The results
are almost identical to the structured mesh case, except for the faster
convergence obtained for SPR-CX in the unstructured compared to the structured
case (0.77 versus 0.59), keeping in mind that optimal convergence rates are
only formally obtained for structured meshes.  Figure 16: Global indicators
$\theta$, $m\left(|D|\right)$ and $\sigma\left(D\right)$ for mixed mode and
unstructured meshes. Notice that the XMLS performs more closely to SPR-CX
using those indicators, for unstructured than for structured meshes,
especially when measuring the mean value of the effectivities
$m\left(|D|\right)$.

##  7 Conclusions

The aim of this paper was to assess the accuracy gains provided by

  1. 1.

Including relevant enrichment functions in the recovery process;

  2. 2.

Enforcing statical admissibility of the recovered solution.

We focused on two recovery-based error estimators already available for LEFM
problems using the XFEM. The first technique called SPR-CX is an enhancement
of the SPR-based error estimator presented by Ródenas et al. (2008), where the
stress field is split into two parts (singular and smooth) and equilibrium
equations are enforced locally on patches. The second technique is the XMLS
proposed by Bordas and Duflot (2007), Bordas et al. (2008) which enriches the
basis of MLS shape functions, and uses a diffraction criterion, in order to
capture the discontinuity along the crack faces and the singularity at the
crack tip.

To analyse the behaviour of the ZZ error estimator using both techniques and
to assess the quality of the recovered stress field, we have evaluated the
effectivity index considering problems with an exact solution. Convergence of
the estimated error in the energy norm and other local error indicators are
also evaluated. To analyse the influence of the special features introduced in
the recovery process we have also considered two additional configurations:
SPR-C and SPR-X.

The results indicate that both techniques, SPR-CX and XMLS, provide error
estimates that converge to the exact value and can be considered as
asymptotically exact. Both techniques could be effectively used to estimate
the error in XFEM approximations, while other conventional recovery procedures
which do not include the enrichment functions in the recovery process have
proved not to converge to the exact error. This shows (albeit only
numerically) the need for the use of extended recovery techniques for accuracy
assessment in the XFEM context.

Better results are systematically obtained when using the SPR-CX to recover
the stress field, specially in the areas close to the singular point. For all
the different test cases analysed, the XMLS produced higher values of the
effectivity index in the enriched area, where the SPR-CX proved to be more
accurate. This can be ascribed to the fact that the SPR-CX technique recovers
the singular part of the solution using the known equilibrated exact
expressions for the asymptotic fields around the crack tip and enforces the
fulfilment of the equilibrium equations on patches. Further work currently in
progress will include the development of an equilibrated XMLS formulation,
which could provide a continuous and globally equilibrated recovered stress
field which can then be used to obtain upper bounds of the error in energy
norm.

The aim of our project is to tackle practical engineering problems such as
those presented in Bordas and Moran (2006), Wyart et al. (2007) where the
accuracy of the stress intensity factor is the target. The superiority of SPR-
CX may then be particularly relevant. To minimise the error on the stress
intensity factors we will target this error directly, through goal-oriented
error estimates. This will be reported in a forthcoming publication. However,
although the SPR-CX results are superior in general, an enhanced version of
the XMLS technique presented in this paper where we enforce equilibrium
conditions could result useful to evaluate upper error bounds when considering
goal-oriented error estimates, as it directly produces a continuous
equilibrated recovered stress field.

Our next step will be the comparison of available error estimators in three
dimensional settings in terms of accuracy versus computational cost to
minimise the error on the crack path and damage tolerance of the structure.

##  8 Acknowledgements

Stéphane Bordas would like to thank the support of the Royal Academy of
Engineering and of the Leverhulme Trust for his Senior Research Fellowship
entitled “Towards the next generation surgical simulators” as well as the
support of EPSRC under grant EP/G042705/1 Increased Reliability for
Industrially Relevant Automatic Crack Growth Simulation with the eXtended
Finite Element Method.
This work has been carried out within the framework of the research projects
DPI 2007-66773-C02-01, DPI2010-20542 and DPI2010-20990 of the Ministerio de
Ciencia e Innovación (Spain). Funding from FEDER, Universitat Politècnica de
València and Generalitat Valenciana is also acknowledged.

## References

  * ABAQUS (2009) ABAQUS (2009).  ABAQUS/standard user’s manual, v.6.9.  Hibbitt, Karlsson & Sorensen, Inc., Pawtucket, Rhode Island.
  * Ainsworth and Oden (2000) Ainsworth, M. and Oden, J. T. (2000).  _A posteriori_ Error Estimation in Finite Element Analysis.  John Wiley & Sons, Chichester.
  * Babuška et al. (1994a) Babuška, I., Strouboulis, T., and Upadhyay, C. S. (1994a).  A model study of the quality of _a posteriori_ error estimators for linear elliptic problems. Error estimation in the interior of patchwise uniform grids of triangles.  Computer Methods in Applied Mechanics and Engineering, 114(3-4):307–378.
  * Babuška et al. (1997) Babuška, I., Strouboulis, T., and Upadhyay, C. S. (1997).  A model study of the quality of _a posteriori_ error estimators for finite element solutions of linear elliptic problems, with particular reference to the behaviour near the boundary.  International Journal for Numerical Methods in Engineering, 40(14):2521–2577.
  * Babuška et al. (1994b) Babuška, I., Strouboulis, T., Upadhyay, C. S., Gangaraj, J., and Copps, K. (1994b).  Validation of _a posteriori_ error estimators by numerical approach.  International Journal for Numerical Methods in Engineering, 37(7):1073–1123.
  * Banks-Sills (1991) Banks-Sills, L. (1991).  Application of the finite element method to linear elastic fracture mechanics.  Applied Mechanics Review, 44(10):447–461.
  * Béchet et al. (2005) Béchet, E., Minnebo, H., Moës, N., and Burgardt, B. (2005).  Improved implementation and robustness study of the X-FEM method for stress analysis around cracks.  International Journal for Numerical Methods in Engineering, 64(8):1033–1056.
  * Belytschko and Black (1999) Belytschko, T. and Black, T. (1999).  Elastic crack growth in finite elements with minimal remeshing.  International Journal for Numerical Methods in Engineering, 45(5):601–620.
  * Belytschko et al. (1996) Belytschko, T., Krongauz, Y., Fleming, M., Organ, D., and Liu, W. K. (1996).  Smoothing and accelerated computations in the element free Galerkin method.  Journal of Computational and Applied Mathematics, 74(1-2):111–126.
  * Bordas and Duflot (2007) Bordas, S. and Duflot, M. (2007).  Derivative recovery and _a posteriori_ error estimate for extended finite elements.  Computer Methods in Applied Mechanics and Engineering, 196(35-36):3381–3399.
  * Bordas et al. (2008) Bordas, S., Duflot, M., and Le, P. (2008).  A simple error estimator for extended finite elements.  Communications in Numerical Methods in Engineering, 24(11):961–971.
  * Bordas and Moran (2006) Bordas, S. and Moran, B. (2006).  Enriched finite elements and level sets for damage tolerance assessment of complex structures.  Engineering Fracture Mechanics, 73(9):1176–1201.
  * CENAERO (2011) CENAERO (2011).  Morfeo.  CENAERO, Gosselies, Belgium.
  * Chessa et al. (2003) Chessa, J., Wang, H., and Belytschko, T. (2003).  On the construction of blending elements for local partition of unity enriched finite elements.  International Journal of Numerical Methods, 57(7):1015–1038.
  * Duflot (2007) Duflot, M. (2007).  A study of the representation of cracks with level sets.  International Journal for Numerical Methods in Engineering, 70(11):1261–1302.
  * Duflot and Bordas (2008) Duflot, M. and Bordas, S. (2008).  _A posteriori_ error estimation for extended finite element by an extended global recovery.  International Journal for Numerical Methods in Engineering, 76(8):1123–1138.
  * Díez et al. (2004) Díez, P., Parés, N., and Huerta, A. (2004).  Accurate upper and lower error bounds by solving flux-free local problems in stars.  Revue européenne des éléments finis, 13(5-6-7):497.
  * Díez et al. (2007) Díez, P., Ródenas, J. J., and Zienkiewicz, O. C. (2007).  Equilibrated patch recovery error estimates: simple and accurate upper bounds of the error.  International Journal for Numerical Methods in Engineering, 69(10):2075–2098.
  * Fries (2008) Fries, T. (2008).  A corrected XFEM approximation without problems in blending elements.  International Journal for Numerical Methods in Engineering, 75(5):503–532.
  * Gdoutos (1993) Gdoutos, E. E. (1993).  Fracture Mechanics: and Introduction. Solid Mechanics and its Applications.  Kluwer Academic Publishers, Dordrecht.
  * Giner et al. (2005) Giner, E., Fuenmayor, F. J., Baeza, L., and Tarancón, J. E. (2005).  Error estimation for the finite element evaluation of ${G}_{\text{I}}$ and ${G}_{\text{II}}$ in mixed-mode linear elastic fracture mechanics.  Finite Elements in Analysis and Design, 41(11-12):1079–1104.
  * Gracie et al. (2008) Gracie, R., Wang, H., and Belytschko, T. (2008).  Blending in the extended finite element method by discontinuous galerkin and assumed strain methods.  International Journal for Numerical Methods in Engineering, 74(11):1645–1669.
  * Ladevèze and Pelle (2005) Ladevèze, P. and Pelle, J. P. (2005).  Mastering calculations in linear and nonlinear mechanics.  Mechanical engineering series. Springer Science.
  * Li et al. (1985) Li, F. Z., Shih, C. F., and Needleman, A. (1985).  A comparison of methods for calculating energy release rates.  Engineering Fracture Mechanics, 21(2):405–421.
  * Melenk and Babuška (1996) Melenk, J. M. and Babuška, I. (1996).  The partition of unity finite element method: Basic theory and applications.  Computer Methods in Applied Mechanics and Engineering, 139(1-4):289–314.
  * Moës et al. (1999) Moës, N., Dolbow, J., and Belytschko, T. (1999).  A finite element method for crack growth without remeshing.  International Journal for Numerical Methods in Engineering, 46(1):131–150.
  * Natarajan et al. (2010) Natarajan, S., Mahapatra, D. R., and Bordas, S. P. A. (2010).  Integrating strong and weak discontinuities without integration subcells and example applications in an XFEM/GFEM framework.  International Journal for Numerical Methods in Engineering, 83(3):269–294.
  * Panetier et al. (2010) Panetier, J., Ladevèze, P., and Chamoin, L. (2010).  Strict and effective bounds in goal-oriented error estimation applied to fracture mechanics problems solved with XFEM.  International Journal for Numerical Methods in Engineering, 81(6):671–700.
  * Pannachet et al. (2009) Pannachet, T., Sluys, L. J., and Askes, H. (2009).  Error estimation and adaptivity for discontinuous failure.  International Journal for Numerical Methods in Engineering, 78(5):528–563.
  * Pereira et al. (1999) Pereira, O. J. B. A., de Almeida, J. P. M., and Maunder, E. A. W. (1999).  Adaptive methods for hybrid equilibrium finite element models.  Computer Methods in Applied Mechanics and Engineering, 176(1-4):19–39.
  * Ródenas et al. (2006) Ródenas, J. J., Giner, E., Tarancón, J. E., and González, O. A. (2006).  A recovery error estimator for singular problems using _singular+smooth_ field splitting.  In Topping, B. H. V., Montero, G., and Montenegro, R., editors, Fifth International Conference on Engineering Computational Technology, Stirling, Scotland. Civil-Comp Press.
  * Ródenas et al. (2010a) Ródenas, J. J., González-Estrada, O. A., Díez, P., and Fuenmayor, F. J. (2010a).  Accurate recovery-based upper error bounds for the extended finite element framework.  Computer Methods in Applied Mechanics and Engineering, 199(37-40):2607–2621.
  * Ródenas et al. (2009) Ródenas, J. J., González-Estrada, O. A., Fuenmayor, F. J., and Chinesta, F. (2009).  Accurate upper bound of the error in the finite element method using an equilibrated moving least squares recovery technique.  In Bouillard, P. and Díez, P., editors, ADMOS 2009. IV International Conference on Adaptive Modeling and Simulation, Barcelona, Spain. CINME.
  * Ródenas et al. (2010b) Ródenas, J. J., González-Estrada, O. A., Fuenmayor, F. J., and Chinesta, F. (2010b).  Upper bounds of the error in X-FEM based on a moving least squares (MLS) recovery technique.  In Khalili, N., Valliappan, S., Li, Q., and Russell, A., editors, 9th. World Congress on Computational Mechanics (WCCM9). 4th. Asian Pacific Congress on Computational Methods (APCOM2010). Centre for Infrastructure Engineering and Safety.
  * Ródenas et al. (2008) Ródenas, J. J., González-Estrada, O. A., Tarancón, J. E., and Fuenmayor, F. J. (2008).  A recovery-type error estimator for the extended finite element method based on singular+smooth stress field splitting.  International Journal for Numerical Methods in Engineering, 76(4):545–571.
  * Ródenas et al. (2007) Ródenas, J. J., Tur, M., Fuenmayor, F. J., and Vercher, A. (2007).  Improvement of the superconvergent patch recovery technique by the use of constraint equations: the SPR-C technique.  International Journal for Numerical Methods in Engineering, 70(6):705–727.
  * Shih and Asaro (1988) Shih, C. and Asaro, R. (1988).  Elastic-plastic analysis of cracks on bimaterial interfaces: Part I \- small scale yielding.  Journal of Applied Mechanics, 8:537–545.
  * Stolarska et al. (2001) Stolarska, M., Chopp, D. L., Moës, N., and Belytschko, T. (2001).  Modelling crack growth by level sets and the extended finite element method.  International Journal for Numerical Methods in Engineering, 51(8):943–960.
  * Strouboulis et al. (2001) Strouboulis, T., Copps, K., and Babuška, I. (2001).  The generalized finite element method.  Computer Methods in Applied Mechanics and Engineering, 190(32-33):4081–4193.
  * Strouboulis et al. (2006) Strouboulis, T., Zhang, L., Wang, D., and Babuška, I. (2006).  A posteriori error estimation for generalized finite element methods.  Computer Methods in Applied Mechanics and Engineering, 195(9-12):852–879.
  * Szabó and Babuška (1991) Szabó, B. A. and Babuška, I. (1991).  Finite Element Analysis.  John Wiley & Sons, New York.
  * Tabbara et al. (1994) Tabbara, M., Blacker, T., and Belytschko, T. (1994).  Finite element derivative recovery by moving least square interpolants.  Computer Methods in Applied Mechanics and Engineering, 117(1-2):211–223.
  * Tarancón et al. (2009) Tarancón, J. E., Vercher, A., Giner, E., and Fuenmayor, F. J. (2009).  Enhanced blending elements for XFEM applied to linear elastic fracture mechanics.  International Journal for Numerical Methods in Engineering, 77(1):126–148.
  * Ventura (2006) Ventura, G. (2006).  On the elimination of quadrature subcells for discontinuous functions in the eXtended Finite-Element Method.  International Journal for Numerical Methods in Engineering, 66(5):761–795.
  * Wyart et al. (2007) Wyart, E., Coulon, D., Duflot, M., Pardoen, T., Remacle, J. ., and Lani, F. (2007).  A substructured FE-shell/XFE-3D method for crack analysis in thin-walled structures.  International Journal for Numerical Methods in Engineering, 72(7):757–779.
  * Xiao and Karihaloo (2006) Xiao, Q. Z. and Karihaloo, B. L. (2006).  Improving the accuracy of XFEM crack tip fields using higher order quadrature and statically admissible stress recovery.  International Journal for Numerical Methods in Engineering, 66(9):1378–1410.
  * Yau et al. (1980) Yau, J., Wang, S., and Corten, H. (1980).  A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity.  Journal of Applied Mechanics, 47(2):335–341.
  * Zienkiewicz and Zhu (1987) Zienkiewicz, O. C. and Zhu, J. Z. (1987).  A simple error estimation and adaptive procedure for practical engineering analysis.  International Journal for Numerical Methods in Engineering, 24:337–357.
  * Zienkiewicz and Zhu (1992a) Zienkiewicz, O. C. and Zhu, J. Z. (1992a).  The superconvergent patch recovery and _a posteriori_ error estimates. Part I: The recovery technique.  International Journal for Numerical Methods in Engineering, 33(7):1331–1364.
  * Zienkiewicz and Zhu (1992b) Zienkiewicz, O. C. and Zhu, J. Z. (1992b).  The superconvergent patch recovery and _a posteriori_ error estimates. Part II: Error estimates and adaptivity.  International Journal for Numerical Methods in Engineering, 33(7):1365–1382.

Corresponding author
Octavio Andrés González-Estrada can be contacted at<EMAIL_ADDRESS>