Linear smoothed extended finite element method

TL;DR

This paper introduces a linear smoothing extension to the XFEM that improves the accuracy of stress intensity factors in fracture mechanics problems while maintaining similar convergence rates and computational efficiency.

Contribution

The paper develops and evaluates a linear smoothing approach for XFEM, enhancing accuracy in numerical integration involving singularities and non-polynomial functions.

Findings

Linear smoothed XFEM yields more accurate stress intensity factors.

All methods show similar convergence rates.

Linear smoothing maintains computational efficiency.

Abstract

The extended finite element method (XFEM) was introduced in 1999 to treat problems involving discontinuities with no or minimal remeshing through appropriate enrichment functions. This enables elements to be split by a discontinuity, strong or weak and hence requires the integration of discontinuous functions or functions with discontinuous derivatives over elementary volumes. Moreover, in the case of open surfaces and singularities, special, usually non-polynomial functions must also be integrated. A variety of approaches have been proposed to facilitate these special types of numerical integration, which have been shown to have a large impact on the accuracy and the convergence of the numerical solution. The smoothed extended finite element method (SmXFEM) [1], for example, makes numerical integration elegant and simple by transforming volume integrals into surface integrals. However, it was reported in [1, 2] that the strain smoothing is inaccurate when non-polynomial functions are in the basis. This is due to the constant smoothing function used over the smoothing domains which destroys the effect of the singularity. In this paper, we investigate the benefits of a recently developed Linear smoothing procedure [3] which provides better approximation to higher order polynomial fields in the basis. Some benchmark problems in the context of linear elastic fracture mechanics (LEFM) are solved to compare the standard XFEM, the constant-smoothed XFEM (Sm-XFEM) and the linear-smoothed XFEM (LSm-XFEM). We observe that the convergence rates of all three methods are the same. The stress intensity factors (SIFs) computed through the proposed LSm-XFEM are however more accurate than that obtained through Sm-XFEM. To conclude, compared to the conventional XFEM, the same order of accuracy is achieved at a relatively low computational effort.

Tables

Table 1. Table 1: Number of sub-cells used to compute the stiffness matrix for the constant smoothed XFEM (Sm-XFEM) and the linear smoothed XFEM (LSm-XFEM). In case of Sm-XFEM, the smoothing function is chosen as f ​ ( 𝐱 ) = 1 𝑓 𝐱 1 f(\mathbf{x})=1 , whilst in case of LSm-XFEM, a complete set of polynomials is chosen. For example, f ​ ( 𝐱 ) = { 1 ​ x ​ y } 𝑓 𝐱 1 𝑥 𝑦 f(\mathbf{x})=\{1\ x\ y\} for two dimensions and f ​ ( 𝐱 ) = { 1 ​ x ​ y ​ z } 𝑓 𝐱 1 𝑥 𝑦 𝑧 f(\mathbf{x})=\{1\ x\ y\ z\} for three dimensions as smoothing function.

	Type of element	Sm-XFEM	LSm-XFEM
two dimensions	Standard elements	4	1
	Tip enriched elements	5	5
	Split enriched elements	8	8
three dimensions	Standard elements	6	1
	Tip enriched elements	24	24
	Split enriched elements	12	12

Equations

$$\nabla \cdot \mathbf{\sigma }+\mathbf{b}=0\text{in }\Omega$$

$$\mathbf{\sigma }\cdot \mathbf{n}=0\text{on }{\Gamma }_c$$

$$\mathbf{\sigma }\cdot \mathbf{n}=\hat{\mathbf{t}}\text{on }{\Gamma }_t$$

$$\mathbf{u}=\bar{\mathbf{u}}\text{on }{\Gamma }_u$$

$$\mathbf{\epsilon }={\nabla }_s\mathbf{u}$$

$$\underset{\Omega }{\int }\mathbf{\epsilon }(\mathbf{u}):\mathbf{C}:\mathbf{\epsilon }(\mathbf{v})\mathrm{d}\Omega =\underset{{\Gamma }_t}{\int }\mathbf{tv}\mathrm{d}\Gamma$$

$$\mathcal{U}$$

$$\mathcal{V}$$

$$\forall {\mathrm{v}}^h\in {\mathcal{V}}^{h}a({\mathrm{u}}^h,{\mathrm{v}}^h)=\ell ({\mathrm{v}}^h),$$

$$\mathrm{K}{\mathrm{u}}^h$$

$$\mathrm{K}=\sum _h{\mathrm{K}}^h$$

$$\mathrm{f}=\sum _h{\mathrm{f}}^h$$

$${\mathrm{u}}^h(\mathbf{x})=\sum _{I\in {\mathcal{N}}^{std}}{\mathbf{N}}_I(\mathrm{x}){\mathbf{u}}_I+\sum _{J\in {\mathcal{N}}^{hev}}{\mathbf{N}}_J(\mathrm{x})H(\mathbf{x}){\mathbf{a}}_J+\sum _{K\in {\mathcal{N}}^{tip}}{\mathbf{N}}_K(\mathrm{x})\left(\sum _{m=1}^4{F}_m(r,\theta ){\mathbf{b}}_K^m\right)$$

$$F_m(r,\theta )$$

$$\mathbf{K}{\mathrm{u}}^h=\mathbf{F}$$

$${\mathbf{K}}^e=\left[\begin{array}{ccc}{\mathbf{K}}^{uu}& {\mathbf{K}}^{ua}& {\mathbf{K}}^{ub}\\ {\mathbf{K}}^{au}& {\mathbf{K}}^{aa}& {\mathbf{K}}^{ab}\\ {\mathbf{K}}^{bu}& {\mathbf{K}}^{ba}& {\mathbf{K}}^{bb}\end{array}\right]$$

$${\tilde{\epsilon }}_{ij}^h({\mathbf{x}}_c)=\underset{{\Omega }^h}{\int }{\epsilon }_{ij}^h(\mathbf{x})\Phi (\mathbf{x}-{\mathbf{x}}_c)\mathrm{d}\Omega$$

$${\tilde{N}}_{I,i}^h(\mathbf{x})=\underset{{\Omega }^h}{\int }{N}_{I,i}^h(\mathbf{x})\Phi (\mathbf{x})\mathrm{d}\Omega$$

$$\underset{{\Omega }^h}{\int }{N}_{I,i}^h(\mathbf{x})\Phi (\mathbf{x})\mathrm{d}\Omega =\underset{{\Gamma }^h}{\int }{N}_I^h(\mathbf{x})\Phi (x)\mathbf{n}(\mathbf{x})\mathrm{d}\Gamma$$

$$\underset{{\Omega }^h}{\int }{N}_{I,i}^h(\mathbf{x})\Phi (\mathbf{x})\mathrm{d}\Omega =\underset{{\Gamma }^h}{\int }{N}_I^h(\mathbf{x})\Phi (x)\mathbf{n}(\mathbf{x})\mathrm{d}\Gamma -\underset{{\Omega }^h}{\int }{N}_I^h(\mathbf{x}){\Phi }^{\prime }(\mathbf{x})\mathrm{d}\Omega$$

$${\tilde{N}}_{I,i}^h({\mathbf{x}}_c)=\frac{1}{A_c}\underset{{\Gamma }^h}{\int }{N}_I^h(\mathbf{x})\mathbf{n}(\mathbf{x})\mathrm{d}\Gamma$$

$${\tilde{\epsilon }}^h({\mathbf{x}}_c)={\tilde{B}}_c({\mathbf{x}}_c)\mathbf{q}$$

$${\tilde{B}}_c=\left[{\tilde{\mathbf{B}}}_{c1}{\tilde{\mathbf{B}}}_{c2}\cdots {\tilde{\mathbf{B}}}_{nc}\right]$$

$${\tilde{N}}_{I,i}^h({\mathbf{x}}_c)=\frac{1}{V_c}\sum _{b=1}^{nb}\left(\begin{array}{ccc}{N}_I^h({\mathbf{x}}_b^G){\mathbf{n}}_{\mathbf{x}}& 0& 0\\ 0& N_I^h({\mathbf{x}}_b^G){\mathbf{n}}_{\mathbf{y}}& 0\\ 0& 0& N_I^h({\mathbf{x}}_b^G){\mathbf{n}}_{\mathbf{z}}\\ N_I^h({\mathbf{x}}_b^G){\mathbf{n}}_{\mathbf{y}}& N_I^h({\mathbf{x}}_b^G){\mathbf{n}}_{\mathbf{x}}& 0\\ 0& N_I^h({\mathbf{x}}_b^G){\mathbf{n}}_{\mathbf{z}}& N_I^h({\mathbf{x}}_b^G){\mathbf{n}}_{\mathbf{y}}\\ N_I^h({\mathbf{x}}_b^G){\mathbf{n}}_{\mathbf{z}}& 0& N_I^h({\mathbf{x}}_b^G){\mathbf{n}}_{\mathbf{x}}\end{array}\right)A_b^c$$

$$\underset{{\Omega }^h}{\int }{\mathbf{N}}_{I,x}^h(\mathbf{x})\mathrm{d}\Omega =\underset{{\Gamma }^h}{\int }{\mathbf{N}}_I^h(\mathbf{x})n_x\mathrm{d}\Gamma$$

$$\underset{{\Omega }^h}{\int }{\mathbf{N}}_{I,x}^h(\mathbf{x})x\mathrm{d}\Omega =\underset{{\Gamma }^h}{\int }{\mathbf{N}}_I^h(\mathbf{x})x{n}_x\mathrm{d}\Gamma -\underset{{\Omega }^h}{\int }{\mathbf{N}}_I^h(\mathbf{x})\mathrm{d}\Omega$$

$$\underset{{\Omega }^h}{\int }{\mathbf{N}}_{I,x}^h(\mathbf{x})y\mathrm{d}\Omega =\underset{{\Gamma }^h}{\int }{\mathbf{N}}_I^h(\mathbf{x})y{n}_x\mathrm{d}\Gamma$$

$$\underset{{\Omega }^h}{\int }{\mathbf{N}}_{I,x}^h(\mathbf{x})z\mathrm{d}\Omega =\underset{{\Gamma }^h}{\int }{\mathbf{N}}_I^h(\mathbf{x})z{n}_x\mathrm{d}\Gamma$$

$$\underset{{\Omega }^h}{\int }{\mathbf{N}}_{I,y}^h(\mathbf{x})\mathrm{d}\Omega =\underset{{\Gamma }^h}{\int }{\mathbf{N}}_I^h(\mathbf{x})n_y\mathrm{d}\Gamma$$

$$\underset{{\Omega }^h}{\int }{\mathbf{N}}_{I,y}^h(\mathbf{x})x\mathrm{d}\Omega =\underset{{\Gamma }^h}{\int }{\mathbf{N}}_I^h(\mathbf{x})x{n}_y\mathrm{d}\Gamma$$

$$\underset{{\Omega }^h}{\int }{\mathbf{N}}_{I,y}^h(\mathbf{x})y\mathrm{d}\Omega =\underset{{\Gamma }^h}{\int }{\mathbf{N}}_I^h(\mathbf{x})y{n}_y\mathrm{d}\Gamma -\underset{{\Omega }^h}{\int }{\mathbf{N}}_I^h(\mathbf{x})\mathrm{d}\Omega$$

$$\underset{{\Omega }^h}{\int }{\mathbf{N}}_{I,y}^h(\mathbf{x})z\mathrm{d}\Omega =\underset{{\Gamma }^h}{\int }{\mathbf{N}}_I^h(\mathbf{x})z{n}_y\mathrm{d}\Gamma$$

$$\underset{{\Omega }^h}{\int }{\mathbf{N}}_{I,z}^h(\mathbf{x})\mathrm{d}\Omega =\underset{{\Gamma }^h}{\int }{\mathbf{N}}_I^h(\mathbf{x})n_z\mathrm{d}\Gamma$$

Full text

††journal: arXiv

1 Introduction

The major difficulty associated with solving problems involving evolving discontinuities is the meshing and re-meshing required as the discontinuities evolve in time. This difficulty is exacerbated when singularities are also present, as is the case in crack growth simulations. These difficulties are somewhat alleviated by the introduction of enrichment functions to represent the discontinuities and the singularities at the patch level, in finite elements or meshfree methods. A first approach to treat discontinuities without an explicit meshing was proposed as early as 1995 in [4]. A much more versatile approach was presented a few years later in the form of the extended finite element method (XFEM) [5] [6] by exploiting the partition of unity property identified by Melenk and Babuška [7]. Partition of unity enrichment for problems with discontinuous solutions is now widely used both in academia and in industrial practice [8, 9, 10] and is known under various names, including the generalized finite element method (GFEM) [11, 12] and the extended finite element method (XFEM). The approach has also been widely used in the form of enriched meshfree methods [13].

Another problem associated with partition of unity methods involving non-polynomial basis functions is to integrate the resulting fields accurately. These enriched methods, also carry along the element mapping involved in building the system matrices. The regularity and positive definiteness of the isoparametric mapping poses a number of restrictions on the allowable shapes of the finite elements: for example, the element should be convex. Meshfree methods also have to face such problems associated with the regularity, distortion and clustering in the point cloud. Under large distortions, meshfree methods face numerical instabilities and low accuracy [14]. Nodal integration also leads to instabilities in cases where higher order shape functions are used. This is due to the fact that in the meshfree methods each node would be associated with a support domain. And the shape functions derivatives would be integrated in this support domain. This implies that each integration domain would be associated with only one integration point (i.e the node). Hence when only one integration is point is considered for higher order functions (other than constant strain) the integration would be similar to the inadequate reduced integration which in turn causes instabilities.

To alleviate these instabilities, the strain smoothing concept was introduced for meshfree methods [15]. The basic idea of strain smoothing is to transform numerical integration over volumes to integration over surfaces, thereby removing instabilities due to the evaluation of the shape function derivatives at the nodes. This approach was later extended to finite element methods by Liu et al [16]. The resulting method was coined the Smoothed finite element method (SFEM), was discussed in a number of papers [17, 2, 18, 19, 20, 21] and applied to a wide variety of problems. The SFEM reduces the mesh sensitivity to some extent by avoiding the necessity of evaluating the Jacobian. Since the derivatives are not needed, the iso-parametric mapping is also avoided.

The SFEM involves computation of a smoothed strain from the standard compatible strain field. The smoothed strain is evaluated as a spatial average of the standard strain field over smoothing subcells which cover the domain andthat can be chosen independently from the mesh structure. These smoothing cells are typically constructed from the mesh following different approaches. This gives rise to a number of methods including cell-based SFEM (CS-FEM) [16, 22, 21, 17], node-based SFEM (NS-FEM) [19], edge-based SFEM (ESFEM) [18], face-based SFEM (FS-FEM).

The method was later extended to solve problems with field discontinuities, both strong and weak, by Bordas et al [1]. This was achieved by extending strain smoothing to the partition of unity framework [23, 7]. Though the smoothed FEM did make the integration procedure more elegant, it was also reported in [2] that the error levels are higher due to the inaccurate approximation of the near tip singular fields. Similar errors were also encountered in higher order elements and polygons [24]. It is noteworthy that similar difficulties are also present in meshfree methods, as addressed in [25] by employing the Discrete Divergence Consistency (DDC) requirement by establishing a compatibility relation between the shape function and its derivatives. This was recently extended for the FEM in [3] and named: Linear smoothing (LS) scheme. It was also reported that the linear smoothing scheme provides an improved accuracy compared to the standard constant-based smoothing of the SFEM.

The present paper aims at investigating how the use of higher-order smoothing, in particular linear smoothing, resolves the lack of accuracy observed when constant smoothing is employed with non-polynomial bases functions. The paper is organized as follows. After presenting the governing equations for elasto-statics, a brief overview of the extended finite element method is given in Section 2. Section 3 presents the linear smoothing technique. A few standard benchmark problems in linear elastic fracture mechanics, solved by using the developed method are presented and the accuracy, convergence and the efficiency of the proposed method are discussed in Section 4, followed by concluding remarks in the last section.

2 Theoretical formulation

2.1 Governing equations for elastostatics

Consider a linear elastic body as shown in Figure 1, with a discontinuity. Let the domain be $\Omega\in\mathbb{R}^{d}$ bounded by $\Gamma$. In this case the boundary has three parts namely $\Gamma_{u}$, where the displacement boundary conditions are applied, $\Gamma_{t}$, where the tractive boundary conditions are applied and $\Gamma_{c}$, which is the traction free surface representing the discontinuity, such that $\Gamma=\Gamma_{u}\cup\Gamma_{t}\cup\Gamma_{c}$ and $\Gamma_{u}\cap\Gamma_{t}=\emptyset$.

The governing equation to be solved is

[TABLE]

The boundary conditions are as follows

[TABLE]

where, is the gradient operator, $\bm{\sigma}$ is the Cauchy stress tensor, $\mathbf{b}$ is the body force per unit volume, $\bf{n}$ is the unit outward normal and $\mathbf{t}$ is the applied tractive stress. For a body undergoing small displacements and strains, the strain displacement equation reads as

[TABLE]

where, $\nabla_{s}$ is the symmetric part of the gradient operator. By substituting the constitutive relations and the strain-displacement relations the weak form of the above Equation (1) becomes Equation(4) in the absence of the body forces: find $\bm{\mathrm{u}}\in\mathscr{U}$ such that

[TABLE]

where, $\bf{u}$ and $\bf{v}$ are the trial and the test function spaces, respectively. For the aforementioned problem, the function spaces includes functions that are discontinuous across $\Gamma_{c}$.

[TABLE]

where the space $\mathcal{W}(\Omega)$ includes linear displacement fields. The domain is partitioned into elements $\Omega^{h}$, and on using shape functions $\phi_{a}$ that span at least the linear space, we substitute vector-valued trial and test functions $\bm{\mathrm{u}}^{h}=\sum_{a}N_{a}\bm{\mathrm{u}}_{a}$ and $\bm{\mathrm{v}}^{h}=\sum_{b}N_{b}\bm{\mathrm{v}}_{b}$, respectively, into Equation (4) and apply a standard Galerkin procedure to obtain the discrete weak form: find $\bm{\mathrm{u}}^{h}\in\mathscr{U}^{h}$ such that

[TABLE]

which leads to the following system of linear equations:

[TABLE]

where $\bm{\mathrm{K}}$ is the assembled stiffness matrix, $\bm{\mathrm{f}}$ the assembled nodal force vector, $\bm{\mathrm{u}}^{h}$ the assembled vector of nodal displacements, $\bm{\mathrm{N}}$ is the matrix of shape functions, $\bm{\mathrm{C}}$ is the constitutive matrix for an isotropic linear elastic material, and $\bm{\mathrm{B}}=\bm{\nabla}_{s}\bm{\mathrm{N}}$ is the strain-displacement matrix that is computed using the derivatives of the shape functions.

2.2 eXtended Finite Element Method

With the regular FEM, the mesh topology has to conform to the discontinuous surface. The introduction of the XFEM has alleviated these difficulties by allowing the discontinuities to be independent of the underlying discretization. Within the framework of the eXtended Finite Element Method (XFEM), the trial functions take the following form:

[TABLE]

where $I$ is the set of all the nodes in the system, $J$ is the set of nodes which are completely cut by the crack, $K$ is the set of nodes which contain the crack tips as shown in Figure 2. $\mathbf{N}_{I}(\mathbf{x})$ are the standard shape functions associated with the standard DOF $\mathbf{u}_{I}$, $H(\mathbf{x})$ is the Heaviside function associated with the enriched DOF, $\mathbf{a}_{J}$ and $F_{m}(r,\theta)$ are the tip enrichment functions associated with the DOF, $\mathbf{b}_{K}^{m}$ that span the near tip asymptotic fields:

[TABLE]

Following the Galerkin procedure, this modification to the trial function space leads to an enlarged problem to solve:

[TABLE]

where

[TABLE]

where, the superscript $uu$ refers to standard FEM components, $aa$ refers to the Heaviside enrichment terms, $bb$ refers to the asymptotic enrichment terms and other terms can be defined similarly. The augmentation of non-polynomial functions to the trial function space, makes the numerical integration non-trivial. This has been of particular interest among research community, for example, equivalent polynomial approach by Ventura [26] and Ventura et al., [27], conformal mapping [28], Duffy transformation [29], generalized Gaussian quadrature [30], strain smoothing technique [2], exponentially convergent mapping [31],polar mapping [32] and very recently by using Euler’s homogeneous function theorem and Stoke’s theorem [33]. In [2], the strain smoothing technique was combined with the XFEM, coined as smoothed XFEM (Sm-XFEM) to integrate over elements intersected with discontinuous surface. The main advantages of the Sm-XFEM are that no subdivision of the split elements is required and that the derivatives of the shape functions (including the enrichment functions) are not required. However, it was observed that the error level was greater when compared to the conventional XFEM, whilst yielding similar convergence rates.

3 Linear smoothing in the XFEM

The strain smoothing was introduced in [15] for the meshfree methods, which was later extended to the FEM by Liu and co-workers [16]. The basic idea is to compute a strain field, referred to as ‘smoothed’ strain field by evaluating the weighted average of the standard (or compatible) strain field. The support domain of the associated material point can be chosen based on surrounding cells, nodes or edges. In this paper, we restrict our discussion only to the cell based strain smoothed FEM. Within this framework, at a point $x_{c}$ in element $\Omega^{h}$ the smoothed strain is given below

[TABLE]

In terms of the standard element shape function derivatives, $N^{h}_{I,i}(\mathbf{x})$, the smoothed derivatives are given by

[TABLE]

where, $\Phi$ is the smoothing function and $i=x,y,z$. By invoking the Divergence theorem Eq. 2 can be written as

[TABLE]

This form of the strain has the following advantages

Domain integration is reduced to a boundary integration along the smoothing cells

• 2.

Does not require the derivatives of the shape functions and hence does not need the Jacobian

• 3.

Does not need isoparametric mapping there by giving a leverage on the distortion level of the mesh

The choice of the smoothing function and the integration order used, decide the accuracy of the smoothed strain field. If a constant smoothing function is used, the method is termed the SFEM. It was shown in [2, 24] that the gradient becomes inaccurate for non-polynomials and higher order polynomial functions. Same issue was also faced with in the context of meshfree approximations and in [25] this inaccuracy was addressed by introducing an additional domain integral term which ensures consistency between the shape functions and their derivatives. This modified equation was termed the Divergence Consistency (DC). It was also shown that such consistency requirement is implicitly satisfied, if linear field is used. It can be seen that Equation(15) would reduce to Equation(14), if $\Phi$ is a constant.

[TABLE]

where, $\Gamma^{h}$ is the contour of the smoothing cell. Here the domain integral term vanishes as the smoothing function is constant over the domain. Since we assumed linear displacement functions, the strain would be a constant and a unique value can be computed using a single equation. Hence requiring just one interior Gauß points. This can be written as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where, $n_{c}$ is the number of sub-cells in an element, $V_{c}$ is the volume of the sub-cell, $n_{b}$ is the number of boundary surfaces of the sub-cell, $A_{b}^{c}$ and $x_{b}^{G}$ are the area and Gauß  point of the boundary surface $b$. The smoothing technique has been very efficient for polyhedral elements since the polyhedrons can be divided into number of sub-cells (usually tetrahedrons) and the stiffness matrix is summed up over each sub-cell. It can be seen in Equation(16d) that the derivatives of the shape functions are not needed in order to evaluate the strains. Hence, the computation of Jacobian is avoided. This also avoids the associated isoparametric mapping. The stiffness matrix is evaluated as in the regular finite element method by replacing the terms in the strain gradient matrix with the terms evaluated above and summing it up over the sub-cells. The constant smoothing technique when applied to elements other than Constant Strain elements (3-noded triangles and 4-noded tetrahedrons) yields results which are bounded by the results of reduced integration procedure (smoothing over one sub-cell) and full integration procedure (smoothing over infinite number of sub-cells). The method is hence not ‘variationally consistent’ for any number of sub-cells other than 1 and $\infty$ [17], whereas the linear smoothing procedure is variationally consistent. The constant smoothing and linear smoothing schemes differ in the choice of the smoothing function. In the linear smoothing scheme the basis function used is $f=\left[1\ x\ y\ z\ xy\ yz\ zx\ xyz\right]^{T}$ in case of hexahedral subcells and $f=\left[1\ x\ y\ z\right]^{T}$if tetrahedral sub-cells are used. Figure 3 shows one possible division of hexahedral elements into tetrahedral elements (also referred as subcells in the literature) for the purpose of numerical integration. The number of terms in the basis function should be consistent with the number of Gauß  points to obtain a unique solution. Since a linear basis function is being used the domain integral term which results as a consequence of the DC does not vanish and hence it has to be computed by using the appropriate order of Gaussian integration. In the case of tetrahedral sub-cells the system of equations for a linear basis would be

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

for $\mathbf{N}^{h}_{I,x}(\mathbf{x})$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

for $\mathbf{N}^{h}_{I,y}(\mathbf{x})$.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

for $\mathbf{N}^{h}_{I,z}(\mathbf{x})$. Here $\mathbf{N}_{I}$ represents the shape function associated with the $I^{th}$ node of the parent element. It is independent of the sub-cell. The location of the gauss points for the boundary integration and domain integration in a tetrahedral sub-cell are shown in Figure 4. Let the natural coordinates of the $m^{th}$ interior gauss point of a sub-cell be $\mathbf{p}_{m}=(x_{m},y_{m},z_{m})$ and its associated gauss weight be $w_{m}$; coordinates of the $k^{th}$ boundary of the sub-cell be $\mathbf{c}^{k}_{g}=(x^{k}_{g},y^{k}_{g},z^{k}_{g})$ and the associated weights be $v^{k}_{g}$ . The unit outward normal associated with the $g^{th}$ gauss point of the $k^{th}$ boundary of the sub-cell is denoted by $\mathbf{n}^{k}=(n^{k}_{x},n^{k}_{y},n^{k}_{z})$. The smoothed derivatives are computed numerically as follows

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The smoothed derivative of the $i^{th}$ shape function evaluated at the four interior Gauß points of a tetrahedral sub-cell is given by

[TABLE]

The same procedure is to be repeated for all the shape functions of the parent element. For the $m^{th}$ interior gauss point of a sub-cell of a $n$ sided polygon the smoothed strain displacement matrix is given by

[TABLE]

[TABLE]

For the displacement approximation given by Equation (8), the compatible strain field is given by:

[TABLE]

where $\mathbf{q}=\{\mathbf{u}\ \mathbf{a}\ \mathbf{b}\}$ is the vector of degrees of freedom, $\mathbf{B}_{\rm fem}$, $\mathbf{B}_{\rm hev}$ and $\mathbf{B}_{\rm tip}$ contains the strain displacement terms corresponding to the regular finite element part, Heaviside enriched part and the tip enriched part. The components of the compatible strain field are:

[TABLE]

where,

[TABLE]

The smoothed counterpart of the above compatible strain field can be obtained by following the procedure outlined earlier. The elements that are intersected by the discontinuous surface is divided into tetrahedra and a linear smoothing basis, $f(\mathbf{x})=\left\{1\ x\ y\ z\right\}$ is chosen to evaluate the smoothed strain.

Remark 1

In case of two dimensions, the subcell is a triangle and the smoothing procedure can be derived from the linear basis

[TABLE]

with derivative

[TABLE]

4 Numerical Examples

In this section, the accuracy and the convergence properties of the proposed formulation is numerically studied within the framework of linear elastic fracture mechanics (LEFM) in both two and three dimensions. The domain is discretized with four noded quadrilateral and eight noded hexahedral elements in two and three dimensions, respectively. The numerical results from the present formulation is compared with the conventional XFEM and the SmXFEM [2]. The following convention is adopted to compute the stiffness matrix within the framework of the smoothing technique:

For the conventional XFEM, the elements that are intersected by the discontinuous surface is triangulated and a higher order triangular quadrature scheme is adopted. And for the standard elements, 2$\times$2 Gauß  quadrature rule is adopted. To estimate the error and to study the convergence properties, the $L^{2}$ norm and the $H^{1}$ semi-norm is used.

4.1 Infinite plate with center crack under far-field tension

Consider an infinite plate with a centre crack subjected to far field tension under plane strain condition has been considered. Consider an infinite plate as shown in Figure 5 A small section ABCD has been solved. The effect of the far-field stress has been accounted by prescribing equivalent displacements as given by following closed form solution Equation (32) in polar coordinates centered at the crack tip.

[TABLE]

where $K_{I}$ $=$ $\sigma$ $\sqrt{\pi a}$, the mode I stress intensity factor, $\nu$ is the Poisson’s ratio, $E$ is the Young’s modulus. The dimension has been chosen as 10 x 10 mm. ‘a’ is chosen as 100 mm. The convergence of the relative error in the displacement ($L^{2}$ norm) and the stress intensity factor is shown in Figure 6. It can be seen that in general all the three methods yields a rate of convergence of 1 in the $L^{2}$ norm and 0.5 in the $H^{1}$ semi-norm. For a given dof, the conventional XFEM yields slightly accurate results than the Sm-XFEM or the LSm-XFEM but the errors are within the same order. Moreover it is noted that in the XFEM, 13 integration points per triangle (for the tip element) is used when compared to three integration points in case of LSm-XFEM and one integration point in case of Sm-XFEM. The sub-optimal rate of convergence is due to the fact that we are employing topological enrichment scheme as opposed to geometric enrichment.

4.2 Finite Plate with an edge crack subject to tensile and shear stresses

Next, consider a finite element with an edge crack subjected to tensile and shear stresses as shown in Figure 7. A consistent system of units is used for the analysis.

Plate subjected to tensile stress

In this case, the dimension of the plate is 1 x 2 units. The Youngs’ modulus, E and Poisson’s ratio, $\nu$ are taken as 1000 and 0.3 respectively. A state of plane strain condition is assumed. The crack width is taken as 0.5 units. The obtained SIF are compared with the reference empirical solution[46]:

[TABLE]

where, $f(\alpha)=1.12-0.231\alpha+10.55\alpha^{2}-21.72\alpha^{3}+30.39\alpha^{4}$, $\alpha=a/W$ is the crack width ratio, $a$ is the half-crack width and $w$ is the plate width. The convergence of the relative error in the stress intensity factor is shown in Figure 8. It can be seen that the all the three methods converge at almost identical rates ($\approx$ 0.5). The results of LS scheme are better than the CS scheme and are almost equal with the conventional XFEM.

Plate subjected to shear stresses

In this case, the dimensions of the plate are taken as $W=$ 7 units and $L=$ 16 units. The plate is subjected to shear stress on the top edges, while the displacements are constrained at the bottom edge. The crack width is taken as 3.5 units. The Youngs’ modulus, $E$ and Poisson’s ratio, $\nu$ are taken as 3$\times$ 107 and 0.25 respectively. Plane strain condition is assumed. The reference SIF is taken from [46], which is $K_{I}$ = 34 units, $K_{II}$ = 4.55 units. The convergence of the $K_{I}$ and the $K_{II}$ are presented in Figure 9. It is again seen that all the three methods have similar convergence rates. The LS scheme is also more accurate than the CS scheme with a very minor additional computational expense. It is again recalled that the additional integration points still require only the shape function values which can be obtained by linear interpolation along the boundary. The error can be attributed to the inadequate approximation space in the local crack tip region, i.e, the asymptotic fields are approximated by a linear field.

4.3 Plate with an inclined center crack subject to tensile stresses

Next, to illustrate the efficacy of the formulation SIFs in case of mixed-mode loading conditions, consider an inclined center crack subjected to far field tension (see Figure 10). The dimensions of the plate are taken as 20 $\times$ 20. The crack width, $2a$ is chosen as 2 units. A far field uniform tensile stress, 1$\times$104 units is applied with Young’s modulus, $E=$ 1$\times$107 and Poisson’s ratio, $\nu=$ 0.3. The accuracy of the numerically computed SIFs are compared with analytical SIFs given by:

[TABLE]

where $\beta$ is the inclination of the crack measured anti-clockwise from the $x-$ axis. Based on a progressive refinement, it was observed that a structured mesh of 100 $\times$ 100 quadrilateral mesh is adequate. The influence of the crack angle and different modelling approaches, viz., XFEM, Sm-XFEM, LSm-XFEM on the SIFs are shown in Figure 11. It can be seen that the results from the proposed approach are accurate and comparable with the conventional XFEM and slightly more accurate than the Sm-XFEM.

4.4 Finite plate with a through-thickness edge-crack subject to tensile stresses

As a last example, the linear smoothing technique is extended to three dimensional domain with a through-the-thickness edge crack subjected to uniform tensile stress as shown in Figure 12 with dimensions $W/a=$ 1 and $H/W=$ 3. The displacement at the bottom face is constrained in all directions and a uniform tensile stress $\sigma=$ 1$\times$104 is applied on the top face. The material properties are: Young’s modulus $E=$ 1$\times$107 and Poisson’s ratio $\nu=$ 0.3. The domain is discretized with structured eight noded hexahedral elements and the normalised SIF from [47] is taken as the reference solution.

The smoothed strain field over a standard element is computed without any further sub-divisions and with $f(\bm{\mathrm{x}})=\left[1\ x\ y\ z\ xy\ yz\ zx\ xyz\right]$ as a smoothing function. For the elements that are intersected by the discontinuous surface, the element is sub-divided into tetrahedra and $f(\bm{\mathrm{x}})=\left[1\ x\ y\ z\right]$ is chosen as the smoothing function. In case of the LSm-XFEM a total of 96 Gauß  points are used in case of tip enriched elements where as 300 Gauß  points are used in the conventional XFEM. In the case of Heaviside enriched elements 48 Gauß  points are used in case of LSm-XFEM and 60 Gauß  points are used in case of the conventional XFEM. The convergence of the relative error in the normalised stress intensity factor is shown in Figure 13. It can be seen that the LSm-XFEM is more accurate than the Sm-XFEM and is in good agreement with the conventional XFEM.

5 Conclusions

In this paper the Linear smoothing (Second order smoothing) was discussed and a method to couple it with the extended finite element method was presented. The developed method was used to solve problems with discontinuities and singularities in both two and three dimensions. The method also involves a rational integration procedure employing the Green’s theorem. The performance of the linear smoothing scheme for enriched approximation space was studied. Through numerical examples it was shown that the Linear smoothing scheme is more accurate than its constant counterpart. The linear smoothing scheme leads to almost identical results to the standard extended finite element method, but it requires fewer quadrature points, viz., approximately one-fourth to what is required with the conventional XFEM.

The constant smoothing and the linear smoothing technique is extended to three dimensions for the first time. Although the presented example in three dimensions is for straight crack, it can be easily extended to other crack profiles. The superior accuracy of the linear smoothing technique is also obtained in the three dimensional case. These results are attributed to the superior approximation properties of the linear smoothing compared to the constant strain smoothing, which is immediately apparent for problems involving complex, non-polynomial, enrichment functions. The remaining, incompressible, error level is attributed to the inadequate approximation space in the smoothed strain, i.e. to the inability of a linear smoothed strain to approximate the singular strains provided by the enriched approximations. Future, ongoing work includes the enrichment of the smoothing space with suitable enrichment functions in order to investigate any additional accuracy improvements as well as the introduction of the approach in recently developed stable extended finite element schemes [48, 49].