text:
\begin{document}

\begin{frontmatter}

\title{Numerical evaluation of stress intensity factors and T-stress for interfacial cracks and cracks terminating at the interface without asymptotic enrichment}

\author[a]{Sundararajan Natarajan\fnref{label1}\corref{cor1}}
\author[a]{Chongmin Song}
\author[b]{Salim Belouettar}

\address [a]{School of Civil and Environmental Engineering, The University of New South Wales, Sydney, NSW 2052, Australia.}
\address[b]{Centre de Recherche Public Henri Tudor, 29, Avenue John F Kennedy, L-1855 Luxembourg-Kirchberg, Luxembourg.}

 \fntext[1]{School of Civil and Environmental Engineering, The University of New South Wales, Sydney, NSW 2052, Australia. Tel: +61(2)9385 5030. Email: s.natarajan@unsw.edu.au; sundararajan.natarajan@gmail.com}

\begin{abstract}
In this paper, we extend the recently proposed extended scaled boundary finite element method (xSBFEM)~\cite{natarajansong2013} to study fracture parameters of interfacial cracks and cracks terminating at the interface. The approach is also applied to crack growth along the interface and crack deflecting into the material within the context of linear elastic fracture mechanics. Apart from the stress intensity factors, the T-stress can be computed directly from the definitions, without any requirement of path independent integrals. The method aims at improving the capability of the extended finite element method in treating crack tip singularities of cracks at interfaces. An optimum size of the scaled boundary region is presented for multimaterial junctions. The proposed method: (1) does not require special numerical integration technique; (2) does not require a priori knowledge of the asymptotic fields and (3) the stiffness of the region containing the crack tip is computed directly. The robustness of the proposed approach is demonstrated with a few examples in the context of linear elastic fracture mechanics. A discussion on the crack growth along the interface and crack deflecting into the material is also presented.
\end{abstract}

\begin{keyword}
scaled boundary finite element method \sep extended finite element method \sep interfacial crack \sep T-stress \sep stress intensity factor \sep enrichment functions \sep numerical integration.
\end{keyword}

\end{frontmatter}

\section{Introduction}

The finite element method (FEM) reformulates the continuous boundary and initial value problems into equivalent variational forms. The FEM requires the domain to be subdivided into non-overlapping regions called `elements' and local polynomial representation is used for the fields within the element. The FEM with piecewise polynomials are inefficient to deal with internal discontinuities such as material interfaces or singularities. The material interfaces can be captured by a conforming mesh and the singularities can be captured by using quarter-point elements~\cite{barsoum1977} or by locally enriching the finite element approximation space with singular functions~\cite{strangfix1973}. In case of the quarter point elements, the singularity is captured by moving the element's mid-side node to the position one quarter of the way from the crack tip. This was considered to be a major milestone in applying the FEM for linear elastic fracture mechanics (LEFM). Tracey~\cite{tracey1971} and Atluri \textit{et al.,}~\cite{atlurikobayashi1975} proposed a new element that has an inverse square root singularity near the crack. The main advantage of this element is that the stress intensity factors (SIFs) can be computed more accurately. Strang and Fix in their well known book~\cite{strangfix1973}, proposed to add singular functions to approximate the displacement field near the singularity. Benzley~\cite{benzley1974} developed an arbitrary quadrilateral element with singular corner node by `enriching' a bilinear quadrilateral element with singular terms. All the above mentioned approaches still relies on a conforming finite element mesh. In an effort to overcome the limitations of the FEM, meshfree methods were introduced. Treatment of evolving discontinuities in meshfree methods is more straightforward because it does not require conforming mesh or mesh adaptation as the discontinuities evolve~\cite{rabczukgracie2010,rabczukbelytschko2006,rabczukzi2008a,zirabczuk2007,bacnguyen-xuan2013,chau-dinhzi2012}. Three dimensional crack propagation have been presented in~\cite{rabczukbelytschko2007,zhuangaugarde2012,rabczukzi2010,bordasgzi2008,rabczukzi2008a,rabczukbordas2007}.On another albeit related front, in the past decade reducing the meshing burden has spurred interest in the research community, whilst still retaining the finite element framework~\cite{melenkbabuvska1996,belytschkoblack1999,garciafancello2000,rabczukbelytschko2004,rabczukzi2008,nguyen-xuanliu2013}. Of particular note, in this study, we focus our attention on enrichment methods, especially, on the methods which are built on the finite element framework, for example, the extended FEM (XFEM). However, the idea presented here can be coupled with other methods, such as the meshless methods~\cite{belytschkoorgan1995} and the boundary element method~\cite{bird2012} to name a few.

Based on the seminal work of Babu\v{s}ka~\textit{et al.,}~\cite{melenkbabuvska1996}, Belytschko's group in 1999~\cite{belytschkoblack1999} introduced the XFEM that can model crack propagation and strong discontinuities with minimal remeshing. Since its inception, the XFEM has emerged as a versatile tool to handle internal discontinuous surfaces or moving boundaries independent of the underlying finite element discretization. Built on a finite element framework, the method relies on a modification (augmentation, enrichment) of the finite element spaces to capture the internal discontinuities. The success of the XFEM when applied to moving boundary problems, esp, crack growth, relies on the a priori knowledge of the set of functions that span the asymptotic fields ahead of the crack tip. The requirement for a priori knowledge of the asymptotic fields hinders the application of the XFEM directly to heterogenous materials for which the asymptotic fields do not exist in closed form or are very complex. When the required functions are not known a priori, these functions can be computed numerically~\cite{duartekim2008,waismanbelytschko2008,menkbordas2010,mousavigrinspun2011,mousavigrinspun2011a}. Nevertheless, all the approaches adopted in the literature depend on the following two steps:
\begin{itemize}
\item If the enrichment functions are unknown a priori, the enrichment functions are computed numerically, for example, the work of Menk and Bordas~\cite{menkbordas2010}, Mousavi \textit{et al.,}~\cite{mousavigrinspun2011} and Zhu~\cite{zhu2012}.
\item The numerically computed functions are then used as enrichment functions within the XFEM/GFEM framework.
\end{itemize}
In addition to the requirement of the knowledge of asymptotic fields, augmenting these functions to the FE approximation basis leads to other difficulties such as, blending problem, numerical integration of enrichment functions, ill-conditioning and additional unknowns. Recent research has been directed towards improving the method by alleviating some of the aforementioned difficulties (for example, hybrid crack element~\cite{xiaokarihaloo2007}, spider XFEM~\cite{chahinelaborde2008}, hybrid analytical XFEM~\cite{rethoreroux2010}, stable GFEM~\cite{babuvskabanerjee2012}, etc.,) whilst the other focus has been in developing new enrichment functions when XFEM is applied to heterogeneous materials (for example, enrichment functions for orthotropic materials~\cite{asadpouremohammadi2007,asharimohammadi2011,hattorirojas-diaz2012}, for functionally graded materials~\cite{dolbowgosz2002}, for crack at bimaterial interfaces~\cite{sukumarhuang2004,liuxiao2004} and cracks terminating at the interface~\cite{bouhalashao2013}, to name a few). However, in spite of these advancements, to accurately evaluate the fracture parameters (such as the stress intensity factors and the T-stress), special techniques, such as path-independent integrals, are necessary~\cite{kimpaulino2003,courtingardin2005,bergerkarageorghis2007,morais2007,passieuxgravouil2011}.

\subsection{Approach} In this paper, we study the singular stress states that exist in solids when discontinuities are present in the geometry and/or the mechanical properties of the material using the scaled boundary formulation. The scaled boundary formulation is then combined with the XFEM to study the fracture parameters, for example, the stress intensity factors (SIFs) and the T-stress. The salient feature is that the combined method does not require a priori knowledge of the asymptotic expansions of the displacements. Moreover, the stiffness of the region containing the crack tip can be computed directly, unlike, the method described in~\cite{menkbordas2010,mousavigrinspun2011,duartekim2008,zhu2012}, which requires a two step procedure.

\subsection{Outline} The paper commences with a brief introduction to the extended finite element method (XFEM) in Section \ref{xfemoverview} and the scaled boundary finite element method (SBFEM) in Section \ref{sbfem}. The coupling of the XFEM and the SBFEM is discussed in Section \ref{couplexfemsbfem}. Section \ref{gsiftress} discusses the computation of the stress intensity factors (SIFs) and the T-stress in an isotropic and in a heterogeneous medium. Section \ref{numresults} presents the displacement and stress modes and the order of singularity for various configuration using the scaled boundary formulation. The effectiveness of the coupled scaled boundary and XFEM is demonstrated with a few problems from linear elastic fracture mechanics. Also, the formulation is used to model crack growth along the interface and into the material. The concluding remarks are given in the last section.

\section{Overview of the extended finite element method}

\input{XFEMintro}

\section{Overview of the scaled boundary formulation}
\label{sbfem}
The scaled boundary finite element method (SBFEM) is a semi-analytical computational technique developed by Wolf and Song~\cite{wolfsong2001}. It reduces the governing partial differential equations to a set of ordinary differential equations and is suitable for solving linear elliptic, parabolic and hyperbolic partial differential equations. In the SBFEM, the coordinates are defined by a radial-circumferential like coordinate system. Numerical solutions are sought around the circumferential direction using the conventional FEM, whilst in the radial direction, the solution is defined by smooth analytical functions. A scaling centre $O$ is selected at a point from which the whole boundary of the domain is visible (see \fref{fig:sbfemdescrip}). The radial coordinate $\xi$ is defined from the scaling centre $(\xi=0)$ to the boundary $(\xi=1)$. Each edge on the boundary is discretised using 1D finite elements with local coordinate $\eta$, defined in the range $-1\leq\eta\leq1$. The coordinate of a point on the element $\mathbf{x}_{\eta}(\eta)=[\begin{array}{cc} x(\eta) & y(\eta)\end{array}]^{\mathrm{T}}$ are expressed as
\begin{equation}
\mathbf{x}_{\eta}(\eta)=\mathbf{N}(\eta)\mathbf{x}_{\mathrm{b}}\label{eq:lineEleGeom}
\end{equation}
where $\xx_b$ are the coordinates of the nodes on the boundary and $\mathbf{N}(\eta)$ is the shape function matrix. In this study, standard 1D Gauss-Lobatto-Lagrange shape functions are used. By choosing the origin of the Cartesian coordinate system as the scaling center, the geometry of the domain is described as:
\begin{equation}
\mathbf{x}(\xi,\eta)=\xi\mathbf{x}_{\eta}(\eta)\label{eqn:sbccartesian}
\end{equation}
\begin{figure}[htpb]
\centering

\subfigure[]{\scalebox{0.6}{\input{./Figures/sbfemnocrk.pstex_t}}}
\subfigure[]{\scalebox{0.6}{\input{./Figures/sbfemcrk.pstex_t}}}
\caption{Coordinate transformation in the scaled boundary finite element method: (a) domain without a crack. The scaling centre $O$ can be placed anywhere as long as it satisfies the visibility criteria and (b) domain with a crack and a material interface. The scaling center $O$ is placed at the crack tip. The material interface within the domain need not be discretized, while on the boundary a conforming mesh is built with a node at the material interface.}
\label{fig:sbfemdescrip}
\end{figure}
Using \Eref{eqn:sbccartesian}, the polar coordinates $r$ and $\theta$ are expressed in the scaled boundary coordinates $\xi$ and $\eta$ as
\begin{align}
r(\xi,\eta) & =\xi r_{\eta}(\eta)=\xi\sqrt{x_{\eta}^{2}(\eta)+y_{\eta}^{2}(\eta)}\label{eq-polar-r}\\
\theta(\eta) & =\arctan\dfrac{y_{\eta}(\eta)}{x_{\eta}(\eta)}\label{eq-polar-theta}
\end{align}
where $r_{\eta}(\eta)$ is the radial coordinate on the boundary. The angle $\theta(\eta)$ is a single-valued function in its principal value ($-\pi<\theta\leq\pi$). The element number and the local coordinate $\eta$ can be regarded as a discrete representation of the angle $\theta$.

\paragraph*{Displacement approximation} The displacements of a point (for example, point $A$ see \fref{fig:sbfemdescrip}) is approximated by:
\begin{equation}
\mathbf{u}(\xi,\eta)=\mathbf{N}(\eta)\mathbf{u}(\xi)\label{eqn:dispapprox}
\end{equation}
 where $\mathbf{u}(\xi)$ are the radial displacement functions. Substituting \Eref{eqn:dispapprox} in the definition of strain-displacement relations, the strains $\boldsymbol{\varepsilon}(\xi,\eta)$ are expressed as:
\begin{equation}
\boldsymbol{\varepsilon}(\xi,\eta)=\mathbf{L}\mathbf{u}(\xi,\eta)\label{eqn:sbfemstrain}
\end{equation}
 where $\mathbf{L}$ is a linear operator matrix formulated in the
scaled boundary coordinates as
\begin{equation}
\mathbf{L}=\mathbf{b}_{1}(\eta)\frac{\partial}{\partial\xi}+\xi^{-1}\mathbf{b}_{2}(\eta)\label{eqn:Loperator}
\end{equation}
 with
\begin{align}
\mathbf{b}_{1}(\eta) & =\frac{1}{|\mathbf{J}(\eta)|}\left[\begin{array}{cc}
y_{\eta}(\eta)_{,\eta} & 0\\
0 & -x_{\eta}(\eta)_{,\eta}\\
-x_{\eta}(\eta)_{,\eta} & y_{\eta}(\eta)_{,\eta}
\end{array}\right] \nonumber \\
\mathbf{b}_{2}(\eta) & =\frac{1}{|\mathbf{J}(\eta)|}\left[\begin{array}{cc}
-y_{\eta}(\eta) & 0\\
0 & x_{\eta}(\eta)\\
x_{\eta}(\eta) & y_{\eta}(\eta)
\end{array}\right]\label{eq:b2}
\end{align}
 In \Eref{eq:b2}, $\mathbf{J}(\eta)$ is the Jacobian on the boundary \cite{wolfsong2001,deekswolf2002}, given by:
\begin{align*}
\mathbf{J}(\eta)= & \left[\begin{array}{cc}
x_{\eta}(\eta) & y_{\eta}(\eta)\\
x_{\eta}(\eta)_{,\eta} & y_{\eta}(\eta)_{,\eta}
\end{array}\right]
\end{align*}

Using ~\Eref{eqn:sbfemstrain} and the Hooke's law $\boldsymbol{\sigma}=\mathbf{D}\boldsymbol{\varepsilon}$, the stresses $\boldsymbol{\sigma}(\xi,\eta)$ is expressed as
\begin{equation}
\boldsymbol{\sigma}(\xi,\eta)=\mathbf{D}\left(\mathbf{B}_{1}(\eta)\mathbf{u}(\xi)_{,\xi}+\xi^{-1}\mathbf{B}_{2}(\eta)\mathbf{u}(\xi)\right)\label{eqn:sbfemstress}
\end{equation}
 where $\mathbf{D}$ is the material constitutive matrix and
\begin{align}
\mathbf{B}_{1}(\eta) & =\mathbf{b}_{1}(\eta)\mathbf{N}(\eta) \nonumber \\
\mathbf{B}_{2}(\eta) & =\mathbf{b}_{2}(\eta)\mathbf{N}(\eta)_{,\eta}
\end{align}
Substituting~\Erefs{eqn:dispapprox}, (\ref{eqn:sbfemstrain}) and (\ref{eqn:sbfemstress}) for the displacement, strain and stress fields, respectively, in the virtual work statement results in \cite{wolfsong2001,deekswolf2002}
\begin{equation}
\begin{split}\delta\mathbf{u}{}_{\mathrm{b}}^{{\rm T}}\left((\mathbf{E}_{0}\xi\mathbf{u}(\xi)_{,\xi}+\mathbf{E}_{1}^{{\rm T}}\mathbf{u}(\xi))|_{\xi=1}-\mathbf{f}\right)\\
-\int\limits _{0}^{1}\delta\mathbf{u}(\xi)^{{\rm T}}\left(\mathbf{E}_{0}\xi^{2}\mathbf{u}(\xi)_{,\xi\xi}+(\mathbf{E}_{0}+\mathbf{E}_{1}^{\mathrm{T}}-\mathbf{E}_{1})\xi\mathbf{u}(\xi)_{,\xi}-\mathbf{E}_{2}\mathbf{u}(\xi)\right)d\xi=0
\end{split}
\label{eqn:sbfemvirtual}
\end{equation}
where $\mathbf{f}$ is the equivalent boundary nodal forces and $\mathbf{u}_{\mathrm{b}}$ is the nodal displacement vector. By considering the arbitrariness of $\delta\mathbf{u}(\xi)$, the following ODE is obtained:
\begin{equation}
\mathbf{E}_{0}\xi^{2}\mathbf{u}(\xi)_{,\xi\xi}+(\mathbf{E}_{0}+\mathbf{E}_{1}^{\mathrm{T}}-\mathbf{E}_{1})\xi\mathbf{u}(\xi)_{,\xi}-\mathbf{E}_{2}\mathbf{u}(\xi)=0\label{eqn:governODEsbfem}
\end{equation}
where $\mathbf{E}_{0},\mathbf{E}_{1}$ and $\mathbf{E}_{2}$ are coefficient matrices given by \cite{wolfsong2001,deekswolf2002}
\begin{align}
\mathbf{E}_{0} & =\int_{\eta}\mathbf{B}_{1}(\eta)^{{\rm T}}\mathbf{D}\mathbf{B}_{1}(\eta)|\mathbf{J}(\eta)|d\eta,\nonumber \\
\mathbf{E}_{1} & =\int_{\eta}\mathbf{B}_{2}(\eta)^{{\rm T}}\mathbf{D}\mathbf{B}_{1}(\eta)|\mathbf{J}(\eta)|d\eta,\nonumber \\
\mathbf{E}_{2} & =\int_{\eta}\mathbf{B}_{2}(\eta)^{{\rm T}}\mathbf{D}\mathbf{B}_{2}(\eta)|\mathbf{J}(\eta)|d\eta.\label{eqn:coeffmat}
\end{align}
They are evaluated element-by-element on the element boundary and assembled similar to the standard FE procedure of assemblage. For a 2 noded line element, the coefficient matrices given by~\Eref{eqn:coeffmat} can be written explicitly. A simple routine is given in~\cite{natarajansong2013}. The code can easily be extended to treat higher order elements.

\paragraph*{Computation of the stiffness matrix}
\Eref{eqn:governODEsbfem} is a homogeneous second-order ordinary differential equation. Its solution is obtained by introducing the variable $\boldsymbol{\chi}(\xi)$
\begin{align}
\boldsymbol{\chi}(\xi)= & \left\{ \begin{array}{c}
\mathbf{u}(\xi)\\
\mathbf{q}(\xi)
\end{array}\right\} \label{eq:chi}
\end{align}
 where $\mathbf{q}(\xi)$ is the internal load vector
\begin{align}
\mathbf{q}(\xi)= & \mathbf{E}_{0}\xi\mathbf{u}(\xi)_{,\xi}+\mathbf{E}_{1}^{{\rm T}}\mathbf{u}(\xi)
\end{align}
The boundary nodal forces are related to the displacement functions by:
\begin{equation}
\mathbf{f}=\mathbf{q}(\xi=1)=(\mathbf{E}_{0}\xi\mathbf{u}(\xi)_{,\xi}+\mathbf{E}_{1}^{{\rm T}}\mathbf{u}(\xi))|_{\xi=1}\label{eqn:nodalforce}
\end{equation}
This allows~\Eref{eqn:governODEsbfem} to be transformed into a first order ordinary differential equation with twice the number of unknowns in an element as:
\begin{equation}
\xi\boldsymbol{\chi}(\xi)_{,\xi}=-\mathbf{Z}\boldsymbol{\chi}(\xi)\label{eq:first order}
\end{equation}
 where $\mathbf{Z}$ is a Hamiltonian matrix
\begin{equation}
\mathbf{Z}=\left[\begin{array}{cc}
\mathbf{E}_{0}^{-1}\mathbf{E}_{1}^{\mathrm{T}} & -\mathbf{E}_{0}^{-1}\\
\mathbf{E}_{1}\mathbf{E}_{0}^{-1}\mathbf{E}_{1}^{\mathrm{T}}-\mathbf{E}_{2} & -\mathbf{E}_{1}\mathbf{E}_{0}^{-1}
\end{array}\right]\label{eq:Hamiltonian matrix}
\end{equation}
An eigenvalue decomposition of $\mathbf{Z}$ is performed. The blocks of eigenvalues and transformation matrices necessary are:
\begin{align}
\mathbf{Z}\left[\begin{array}{c}
\boldsymbol{\Phi}_{\mathrm{u}}\\
\boldsymbol{\Phi}_{\mathrm{q}}
\end{array}\right] & =\left[\begin{array}{c}
\boldsymbol{\Phi}_{\mathrm{u}}\\
\boldsymbol{\Phi}_{\mathrm{q}}
\end{array}\right]\boldsymbol{\Lambda}_{\mathrm{n}}\label{eq:eigen decomp}
\end{align}
In \Eref{eq:eigen decomp}, $\boldsymbol{\Lambda}_{\mathrm{n}}=\mathrm{diag}\left(\lambda_{1},\,\lambda_{2},\,...,\lambda_{n}\right)$ contains the eigenvalues with negative real part. $\boldsymbol{\Phi}_{\mathrm{u}}$ and $\boldsymbol{\Phi}_{\mathrm{q}}$ are the corresponding transformation matrices of $\boldsymbol{\Lambda}_{\mathrm{n}}$. They represent the modal displacements and forces, respectively. The general solution of \Eref{eq:first order} is given by:

\begin{align}
\mathbf{u}(\xi)= & \boldsymbol{\Phi}_{\mathrm{u}}\xi^{-\boldsymbol{\Lambda}_{\mathrm{n}}}\mathbf{c}\\
\mathbf{q}(\xi)= & \boldsymbol{\Phi}_{\mathrm{q}}\xi^{-\boldsymbol{\Lambda}_{\mathrm{n}}}\mathbf{c}
\label{eqn:raddispfunc}
\end{align}
where $\mathbf{c}$ are integration constants that are obtained from the nodal displacements $\mathbf{u}_{\mathrm{b}}=\mathbf{u}(\xi=1)$ as:
\begin{align}
\mathbf{c}= & \boldsymbol{\Phi}_{\mathrm{u}}^{-1}\mathbf{u}_{\mathrm{b}}\label{eq:int constants}
\end{align}
The complete displacement field of a point defined by the sector covered by a line element on the element is obtained by substituting~\Eref{eqn:raddispfunc} into \Eref{eqn:dispapprox} resulting in:
\begin{align}
\mathbf{u}(\xi,\eta)= & \mathbf{N}(\eta)\boldsymbol{\Phi}_{\mathrm{u}}\xi^{-\boldsymbol{\Lambda}_{\mathrm{n}}}\mathbf{c}\label{eq:dispapprox final}
\end{align}
Taking the derivative of $\mathbf{u}(\xi)$ with respect to $\xi$ and substituting into~\Eref{eqn:sbfemstress} the stress field $\boldsymbol{\sigma}(\xi,\eta)$ can be expressed as:
\begin{align}
\boldsymbol{\sigma}(\xi,\eta)= & \boldsymbol{\Psi}_{\sigma}(\eta)\xi^{-\boldsymbol{\Lambda}_{\mathrm{n}}-\mathbf{I}}\mathbf{c}\label{eq:stress field complete}
\end{align}
where the stress mode $\boldsymbol{\Psi}_{\sigma}(\eta)$ is defined as:
\begin{align}
\boldsymbol{\Psi}_{\sigma}(\eta)= & \mathbf{D}\left(-\mathbf{B}_{1}(\eta)\boldsymbol{\Phi}_{\mathrm{u}}\boldsymbol{\Lambda}_{\mathrm{n}}+\mathbf{B}_{2}(\eta)\boldsymbol{\Phi}_{\mathrm{u}}\right)\label{eq:stress mode}
\end{align}
The stiffness matrix of an element is obtained by first substituting \Eref{eq:int constants} into~\Eref{eqn:raddispfunc} at $\xi=1$. This results in:
\begin{align}
\mathbf{f}= & \boldsymbol{\Phi}_{\mathrm{q}}\boldsymbol{\Phi}_{\mathrm{u}}^{-1}\mathbf{\mathbf{u}_{\mathrm{b}}}\label{eq:equil poly}
\end{align}
From \Eref{eq:equil poly}, the stiffness matrix $\mathbf{K}$ can be identified to be given by the expression
\begin{equation}
\mathbf{K}=\boldsymbol{\Phi}_{\mathrm{q}}\boldsymbol{\Phi}_{\mathrm{u}}^{-1}\label{eqn:sbfemkmat-b}
\end{equation}

\begin{remark}
No a priori knowledge of the asymptotic solution and special integration technique is required to compute the stiffness matrix. The SBFEM solution includes the stress singularity as well. The details are given in Section \ref{gsiftress}.
\end{remark}

\begin{remark}
The SBFEM can be formulated over arbitrary regions, for example over polygonal elements~\cite{ooisong2012,natarajanooi2013}.
\end{remark}

\section{Coupling of the XFEM \& the SBFEM} \label{couplexfemsbfem}
Recently, the authors proposed the extended SBFEM (xSBFEM)~\cite{natarajansong2013}. The main idea of the xSBFEM is to approximate the non-smooth behaviour in the close vicinity by the scaled boundary formulation instead of augmenting the approximation basis with asymptotic expansion. The nodes whose nodal support is completely cut by the crack interior or the material interface are treated within the XFEM framework and enriched with a Heaviside function or 'abs' enrichment function. The scaling center is chosen as the crack tip (see \fref{fig:couplexfemsbfem}). This local modification and enrichment strategy introduces four types of elements:

\begin{itemize}
\item \textit{Split elements} are elements completely cut by the crack interior or by the material interface. Their nodes are enriched with the Heaviside function to represent the jump across the crack face or by the 'abs' function to represent the weak discontinuity.
\item \textit{Split-blending elements} are elements neighbouring split elements. They are such that some of their nodes are enriched with the weakly discontinuous function and others are not enriched at all.
\item \textit{Scaled boundary elements} are elements in the close proximity of the crack tip treated by the semi-analytical formulation.
\item \textit{Standard elements} are elements that are in neither of the above categories. None of their nodes are enriched.
\end{itemize}

\begin{remark}
It is noted that no special treatment is required for the elements that share nodes between the $\Omega^{\rm fem}$ and $\Omega^{^{\rm sbfem}}$.
\end{remark}

\begin{figure}[htpb]
\centering
\includegraphics[scale=0.5]{./Figures/couplingxfemsbfem}
\caption{Coupling of XFEM and SBFEM}
\label{fig:couplexfemsbfem}
\end{figure}

The displacement approximation can be decomposed into the standard part $\uu^h_{\rm std}$ and an enriched part $\uu^h_{\rm xfem}$ as given by \Eref{eqn:uS2} without asymptotic enrichments. In the proposed approach, the standard part is further decomposed into the standard FEM $\uu^h_{\rm fem}$ and into a semi-analytical part $\uu^h_{\rm sbfem}$. The displacement approximation is then given by:
\begin{equation}
\renewcommand{\arraystretch}{2}
\uu^h(\xx) = \left\{ \begin{array}{lr} \sum\limits_{I \in \mathcal{N}^{\rm fem}} N_I(\xx) \qq_I + \sum\limits_{J \in \mathcal{N}^c} N_J(\xx) H(\xx) \aaa_J & \xx \in \Omega^{^{\rm fem}} \cup \Omega^{^{\rm xfem}} \\
 \sum\limits_{K \in \mathcal{N}^{\rm sbfem}} N_K(\eta)  \uu_{bK}(\xi) & \xx = f(\xi,\eta) \in \Omega^{^{\rm sbfem}} \end{array} \right.
 \label{eqn:dispApprox}
\end{equation}
where $\mathcal{N}^{\rm fem}$ is the set of all nodes in the FE mesh that belongs to $\Omega^{^{\rm fem}}$, $\mathcal{N}^{\rm sbfem}$ is the set of nodes in the FE mesh that belongs to $\Omega^{^{\rm sbfem}}$ and $\mathcal{N}^{\rm xfem}$ is the set of nodes that are enriched with the Heaviside function $H$. $N_I$ and $N_J$ are the standard finite element functions, $\qq_I$ and $\aaa_J$ are the standard and the enriched nodal variables associated with node $I$ and node $J$, respectively. $\uu_{bK}$ are the nodal variables associated with the nodes on the boundary of $\Omega^{^{\rm sbfem}}$ and $N_K$ are the standard FE shape functions defined on the scaled boundary coordinates in $\Omega^{^{\rm sbfem}}$. For the nodes that on the boundary of $\Omega^{^{\rm fem}}$ and $\Omega^{^{\rm sbfem}}$, no special coupling technique is required. As the nodal unknown coefficients are required to be continuous across the boundary, the unknown coefficients from the SBFEM and the FEM are the same and are assembled in the usual way. A similar procedure is followed when assembling the $\Omega^{^{\rm fem}}$ and $\Omega^{^{\rm xfem}}$ to the global stiffness and to the global force vector. For more details, interested readers are referred to~\cite{natarajansong2013}.

\section{Computation of the stress intensity factors and the T-Stress} \label{gsiftress}
\input{gsif}

\section{Results and discussion}
\label{numresults}
In this section, we first employ the SBFEM to estimate the order of the singularity and the stress and the displacement modes for known benchmark problems in the context of linear elastic fracture mechanics. Later, we will adopt the combined XFEM and the SBFEM to compute the stress intensity factors and the T-stress for cracks aligned to the material interfaces. The stress distribution ahead of the crack and the T-stress for a crack terminating at an interface is also studied. The results of which are compared with a finite element solution. In this study, unless otherwise mentioned, bilinear quadrilateral elements are used. An extension to other arbitrary polytopes is possible as illustrated in~\cite{natarajanooi2013}. The results from the present approach are compared with available numerical results in the literature, for example, with the XFEM~\cite{yuwu2012} and with the BEM~\cite{sladeksladek1997,kimvlassak2006}.

\subsection{Computation of displacement \& stress modes and the order of singularity}
To test and to demonstrate the convergence and the accuracy of the SBFEM in estimating the strength of the singularity and the corresponding displacement and stress modes, we consider three examples.

\paragraph{Crack in an isotropic medium} \fref{fig:isoprobdescription} shows a crack in an isotropic medium with Young's modulus $E$ and Poisson's ratio $\nu$. The displacement fields are given by~\cite{williams1957}:
\begin{align}
u_x(r,\theta) &= \frac{K_I}{2\mu} \sqrt{ \frac{r}{2\pi}} \cos \frac{\theta}{2} \left[ \kappa -1 + 2 \sin^2\frac{\theta}{2} \right] + \frac{K_{II}}{2\mu} \sqrt{ \frac{r}{2\pi}} \sin \frac{\theta}{2} \left[ \kappa + 1 + 2\cos^2 \frac{\theta}{2} \right] \nonumber \\
u_y(r,\theta) &= \frac{K_I}{2\mu} \sqrt{ \frac{r}{2\pi}} \sin \frac{\theta}{2} \left[ \kappa +1 - 2 \cos^2\frac{\theta}{2} \right] + \frac{K_{II}}{2\mu} \sqrt{ \frac{r}{2\pi}} \sin \frac{\theta}{2} \left[ \kappa - 1 - 2\sin^2 \frac{\theta}{2} \right]
\end{align}
where $\kappa$ is the Kolosov constant and $\mu$ is the shear modulus. The reference solution is obtained along a circle of unit radius $r=$ 1 and with stress intensity factors $K_I$ is set to unity and $K_{II}$ to zero. The displacement modes, the stress modes and the order of singularity are computed by the scaled boundary formulation. A typical scaled boundary finite element mesh is shown in \fref{fig:isoprobdescription}. It is noted that only the boundary of the domain needs to be discretized and the solution within the domain is represented analytically.

\begin{remark}
In our study, we have employed spectral element to discretize the boundary. We employ Lagrange interpolants within each element, where the nodes of these shape functions are placed at the zeros of Legendre polynomials (Gauss Lobatto points) mapped from the reference domain $[-1,1]\times[-1,1]$ to each element.
\end{remark}

In Table \ref{table:isocrksingularity}, the convergence of the order of singularity with increasing number of elements is shown. With decreasing element size, it is seen that the numerically estimated order of singularity converges to the theoretical value 0.5. Also, the influence of the order of the shape functions on the convergence of the order of singularity is shown. As expected, the theoretical value of 0.5 can be achieved by either increasing the order of the shape functions or by decreasing the element size. \fref{fig:isocrkres} shows the relative error in the displacement norm for each component (i.e., $u_x$ and $u_y$) with the increasing number of elements along the boundary. It can be seen that increasing the number of elements leads to decreasing error.

\Eref{eq-stress-bd-slu} is the parametric equation for the stress field in the polar coordinates $r$ and $\theta$. The terms $\left(r_{\eta}^{\lambda_{i}+1}(\eta)\boldsymbol{\psi}_{\sigma i}(\eta)\right)$ in \Eref{eq-stress-bd-slu} together with $\theta(\eta)$ in \Eref{eq-polar-theta} are the stress modes describing the angular stress distribution at a constant radial coordinate $r$. \fref{fig:isomatsingstressD} shows the stress distribution ahead of the crack tip for mode I fracture problem. Each of the stress mode is normalized with its value of $\sigma_{yy}$ at $\theta=$ 0$^\circ$. A total of 8 third order elements were used to discretize the boundary. The results from the scaled boundary formulation are compared with the analytical solutions (represented by the superscript $E$). A very good agreement is observed.

\begin{figure}[htpb]
\centering
\subfigure[]{\scalebox{0.5}{\input{./Figures/isocrk.pstex_t}}}
\subfigure[]{\includegraphics[scale=0.5]{./Figures/isomeshcrk}}
\caption{Crack in an isotropic medium: (a) geometry and (b) SBFEM Mesh detail}
\label{fig:isoprobdescription}
\end{figure}

\begin{table}[htpb]
\centering

\caption{Convergence of the order of singularity for an isotropic plate with a crack.}
\begin{tabular}{lrrrr}
\hline
Number of & \multicolumn{4}{c}{Order of shape functions}\\
\cline{2-5}
Elements & $p=$2 & $p=$ 3 & $p=$ 4 & $p=$ 5 \\
\hline
8 & 0.50826515 & 0.50007949 & 0.49995729 & 0.50000291\\
40 & 0.50035216	 & 0.49999979	 & 0.50000000	& 0.50000000\\
80 & 0.50008862	 & 0.49999999	 & - & - \\
160 & 0.50002219 & 0.50000000 & - & - \\
\hline
\end{tabular}
\label{table:isocrksingularity}
\end{table}

\begin{figure}
\centering
\newlength\figureheight
\newlength\figurewidth
\setlength\figureheight{8cm}
\setlength\figurewidth{10cm}
\subfigure[Relative error in displacement, $u_x$]{\input{./Figures/isocrkux.tikz}}
\subfigure[Relative error in displacement, $u_y$]{\input{./Figures/isocrkuy.tikz}}
\caption{Crack in an isotropic medium: Convergence in the relative error in the displacement.}
\label{fig:isocrkres}
\end{figure}

\begin{figure}
\centering

\setlength\figureheight{8cm}
\setlength\figurewidth{10cm}
\input{./Figures/strDIsoCrkAngle.tikz}
\caption{Angular distribution of singular stress ahead of a crack tip for a crack in an isotropic medium under mode I loading conditions, i.e., $K_I=$ 1 and $K_{II} = $ 0.}
\label{fig:isomatsingstressD}
\end{figure}

\paragraph{Crack aligned to a bi-material interface} The second test case is a crack aligned to the bimaterial interface. \fref{fig:bimatprobdescription} shows the geometry and the corresponding finite element mesh. Again, only the boundary of the domain is discretized. The Lagrange shape functions of order $p=$ 9 and a total of eight elements are used to discretize the domain in the circumferential direction. The order of the singularity computed from the semi-analytical approach is compared with the analytical expression  given by, $0.5 \pm i \epsilon$, where $\epsilon$ is the oscillatory index given by:
\begin{equation}
\epsilon = \frac{1}{2\pi} \mathrm{log} \left( \frac{1 - \beta}{1+\beta} \right)
\label{eqn:osciindex}
\end{equation}
where $\beta$ is the second Dundur's parameter
\begin{equation}
\beta = \frac{\mu_1(\kappa_2-1) + \mu_2(\kappa_1-1)}{\mu_1(\kappa_2+1) + \mu_2(\kappa_1+1)}
\end{equation}
and,
\begin{equation}
\kappa_i = \left\{ \begin{array}{cr} (3-\nu_i)/(1+\nu_i) & \textup{plane stress} \\ 3 -4\nu_i & \textup{plane strain} \end{array} \right.
\end{equation}
where $\mu_i, \nu_i$ and $\kappa_i, (i=1,2)$ are the shear modulus, Poisson's ratio and the Kolosov coefficient, respectively. \fref{fig:bimatsing} shows the influence of the ratio of Young's modulus on the order of the singularity. The angular stress distribution ahead of the crack aligned to a bimaterial interface is shown in \fref{fig:bimatsingstressD}.  The ratio of Young's modulus is $E_1/E_2=$ 10. The singular stress fields are obtained for a crack length $a$ = 0.5 and $K_I/p\sqrt{\pi a}=$ 1.123 and $K_{II}/p\sqrt{\pi a}=$ -0.123. A very good agreement is observed.

\begin{figure}[htpb]
\centering
\subfigure[]{\scalebox{0.5}{\input{./Figures/bimatcrk.pstex_t}}}
\subfigure[]{\includegraphics[scale=0.5]{./Figures/bimeshcrk}}
\caption{Crack aligned to a bi-material interface: (a) geometry and (b) SBFEM Mesh}
\label{fig:bimatprobdescription}
\end{figure}

\begin{figure}
\centering

\setlength\figureheight{8cm}
\setlength\figurewidth{10cm}
\input{./Figures/bimatsinguorder.tikz}
\caption{Bimaterial interface crack: order of singularity $(\lambda)$ as a function of the Young's modulus mismatch.}
\label{fig:bimatsing}
\end{figure}

\begin{figure}
\centering

\setlength\figureheight{8cm}
\setlength\figurewidth{10cm}
\input{./Figures/StrDBimatInterCrk.tikz}
\caption{Angular distribution of singular stress ahead of a crack tip for a crack aligned to the bimaterial interface. The ratio of Young's modulus is $E_1/E_2=$ 10. The singular stress fields are obtained for a crack length $a$ = 0.5 and $K_I/p\sqrt{\pi a}=$ 1.123 and $K_{II}/p\sqrt{\pi a}=$ -0.123.}
\label{fig:bimatsingstressD}
\end{figure}

\paragraph{Triple junction} As a last example, we consider the two specific cases of a triple junction: (a) a fully bonded triple junction and (b) a triple junction with a crack aligned to one of the interfaces. \fref{fig:crkatint} shows the problem description and the corresponding SBFEM mesh. Note that only the boundary of the domain needs to be discretized with 8 elements (on the boundary) with Lagrange shape functions of order $p=$ 9 are used. All materials are assumed to be homogeneous and isotropic. The computed order of singularity for various combination of material parameters are shown in \fref{fig:tjunsingular}. The results shown in \fref{fig:tjunsingular} indicate the existence of significant singularities over the full range of $E_2/E_1$. The singularities are more severe for higher values of $E_2/E_1$ than for lower values. \fref{fig:tjunsingular} also shows the case when the crack is present at the interface of material 3 and material 2. Again, a full range of the order of singularity is shown for various combinations of material parameters. For certain values of $E_3/E_2$, multiple real roots are obtained. The present formulation is efficient in capturing the different orders of singularity.

\begin{figure}[htpb]
\centering
\subfigure[]{\scalebox{0.5}{\input{./Figures/trimatnocrk.pstex_t}}}
\subfigure[]{\includegraphics[scale=0.5]{./Figures/triplejnmeshdetail}}
\subfigure[]{\scalebox{0.5}{\input{./Figures/trimatcrk.pstex_t}}}
\subfigure[]{\includegraphics[scale=0.5]{./Figures/triplejnmeshdetailcrk}}
\caption{(a-b) a fully bonded triple junction and SBFEM mesh and (c-d) a triple junction with a crack aligned to one of the interfaces.}
\label{fig:crkatint}
\end{figure}

\begin{figure}
\centering

\setlength\figureheight{8cm}
\setlength\figurewidth{10cm}
\subfigure[$E_3/E_1=$ 10, $E_1=$ 1]{\input{./Figures/TripleNoCrk.tikz}}
\subfigure[$E_1/E_2=$ 10, $E_2=$ 1]{\input{./Figures/TripleDisbond.tikz}}
\caption{Order of singularities for a triple junction: (a) a fully bonded junction and (b) a crack along the interface between the material 2 and material 3.}
\label{fig:tjunsingular}
\end{figure}

\input{bimaterialcrack}

\section{Conclusion}
In this paper, we studied the strength of the singularity of cracks aligned to bimaterial interfaces and cracks terminating at the interface by employing a semi-analytical approach. Later, the stress intensity factors and the T-stress was computed by using a combined XFEM and a semi-analytical approach. The salient features of this technique are: (i) no a priori knowledge of the asymptotic fields is required and (ii) the stiffness matrix of the region enriched for a particular crack is obtained directly without a special numerical integration technique. This eliminates the need to compute the enrichment functions as in~\cite{duartekim2008,menkbordas2010,mousavigrinspun2011,zhu2012}, whilst employing the Heaviside functions allows the discontinuous surfaces to be presented independent of the underlying mesh.The crack growth along the interface and the deflection of the crack into a material is also discussed. The proposed technique can readily be coupled with the existing XFEM/GFEM code and can be applied to study fracture behaviour in heterogeneous materials.

\section*{Acknowledgements}
S Natarajan would like to acknowledge the financial support of the School of Civil and Environmental Engineering, The University of New South Wales for his research fellowship for the period September 2012 onwards.

\section*{References}
\bibliographystyle{elsarticle-num}
\bibliography{SinguStress}

\end{document}


TITLE: Find the number of points in $(-\infty,\infty)$ for which $x^2-x\sin x-\cos x=0.$
QUESTION [0 upvotes]: Find the number of points in $(-\infty,\infty)$ for which $x^2-x\sin x-\cos x=0.$

How should i attempt this question?I have no idea.Please help.Thanks.

REPLY [2 votes]: Here we have to calculate number of real values of $x$ in $x^2=x\sin x+\cos x$
So we will draw graph of $f(x) = x^2$ and $f(x) = x\sin x+\cos x$
First one is upward parabola passing through origin and for drawing graph of
$f(x)=x\sin x+\cos x\;,$ we will put $\displaystyle x=0\;,\frac{\pi}{2}\;,\pi,\frac{3\pi}{2},2\pi$ and $f(-x) = f(x)$
Means $f(x)$ is an even function .

So here these two graph intersect each other at $2$ distinct point
So we get number of real solution of the equation is $\bf{Two}$


TITLE: Evaluate integration with parameter
QUESTION [1 upvotes]: I want to evaluate the integration $$I(\epsilon)=\int_0^\infty \frac{1}{\sqrt {x(x+\epsilon)(1+x^5)}}\,dx \quad(\epsilon\to0)$$
Although $I(0)$ is divergent, one may still derive an approximate expression  for $I(\epsilon)$ at small $\epsilon$ or a numerical expression with $\epsilon$ being a parameter.
Using Mathematica's LogLinearPlot to draw a logarithmatic graph,I found that $I(\epsilon)$ is linear to $\epsilon$'s logarithm. And by linear fitting I made a guess: When $\epsilon$ goes to $0$, there will be the relation $$I(\epsilon)= -\ln\epsilon+\textrm{const.}$$
But how can I figure out the constant above? By assymptotic expansion or something?
Or the observation and assumption I made are wrong as a whole? Is there any other way to evaluate the approximate or numerical expression of $I(\epsilon)?$

REPLY [1 votes]: Integrate by parts, writing $\displaystyle\frac{1}{\sqrt{x(x+\epsilon)}}=\frac{d}{dx}\ln\color{gray}{\frac{1}{2}}\left(x+\frac{\epsilon}{2}+\sqrt{x(x+\epsilon)}\right)$. As a result, you get $$\lim_{\epsilon\to+0}\left(\ln\epsilon+\int_0^\infty\frac{dx}{\sqrt{x(x+\epsilon)(1+x^a)}}\right)=2\ln 2+\frac{a}{2}\int_0^\infty\frac{x^{a-1}\ln x}{(1+x^a)^{3/2}}\,dx=2\ln 2+\frac{A}{2a}$$ (for real $a>0$) where $A=\displaystyle\int_0^\infty\frac{\ln t\ dt}{(1+t)^{3/2}}=4\ln 2$ (integrated elementarily).
Thus the limit is equal to $2(1+1/a)\ln 2$. Put $a=5$ to get your "$\mathrm{const.}$".