Extended isogeometric analysis using analysis-suitable T-splines for plane crack problems

The brutal fracture problem has a great importance in some industrial fields, such as in aeronautics, aerospace and nuclear. This phenomenon has very serious consequences, lead to the need for analysing more and better understanding the mechanical behaviour of structures, with taking into account the effects of the fort and the weak discontinuities, especially in critical conditions. It is a scientific challenge that represents an important issue, analytically as well as numerically. The complexity of the analytical solutions even for simple cases imposes to setting up effective numerical methods to model the mechanical behaviour. Various studies have been developed using different methods, such as extended finite element method (XFEM) [1], boundary element method (BEM) [2] , element free Galerkin method (EFGM) [3], and other methods [4, 5]. Recently, a large field was opened by Hughes et al. [6] offers the possibility of exploiting computer aided design (CAD) tools in the analysis methods using isoparametric concept. This novel alternative method called isogeometric analysis (IGA). The basic idea of this method is use the computational geometry technologies as shape bases to describe exactly the geometry and also to approximate the unknown fields. Following this discovery, several researches in various fields have been investigated by this method, including: structural dynamics [7], composite materials [8], fluid– structure interaction [9], electromagnetic problems [10], contact problems [11], turbulent flow [12], aero-dynamics [13] and thermomechanical problems [14]. However, in fracture mechanics problems, Benson et al. [15] and De Luycker et al. [16] have proposed extended isogeometric analysis (XIGA) for modelling cracks. In this method the general principle of the XFEM is used into IGA by including the asymptotic and signed distance enrichment functions. Therefore this method has the advantages of both XFEM and IGA, which are summarized by the ability to represent complex geometries independently of any discontinuities and without explicit meshing for obtain solution with higher orders. There are many CAD basis functions can be used in isogeometric analysis, the most widely used of them are Non Uniform Rational B-splines (NURBS) due to their properties like continuity, smoothness, variation diminishing, convex hull and possibility of using Knot insertion and degree elevation refinements. They have the ability to describe exactly all conic sections but they are not well for all complex geometries, even for multiple patches NURBS generate a complicated mesh of control points and this leads to produce superfluous control points. In order to handle these disadvantages Sederberg et al. [17] proposed T-splines as a generalized tool of NURBS, in which the index space (T-mesh) cans have T-junctions and locally refined [18]. Therefore the major advantages of this technique are local refinement and ability to represent complex geometries with minimal number of control points compared with those used in NURBS. According to their ability in engineering design, T-splines have been used by the analysts to serves as basis functions for isogeometric analysis in many advanced searches. However T-splines bases are not always valid to use in analysis for different geometric configurations, because there are not rules ensure these bases to be linear independent and form a partition of unity. Li et al. [19] introduced analysis-suitable T-splines, wherein for any choice of knot vectors the blending functions are linearly independent. Like NURBS bases, analysis-suitable Tspline bases have the important properties of the analysis basis functions. Moreover, they provide an efficient algorithm which allows making highly localized refinement [20]. In this paper, plane crack problem is analysed numerically using XIGA, which is enhanced by adopting analysis-suitable T-splines in order to approximate the solution, construct the geometry and make the local refinement around the discontinuities. The asymptotic functions are used to model the crack tip, and the M-integral is used to evaluate the stress intensity factors.


Introduction
The brutal fracture problem has a great importance in some industrial fields, such as in aeronautics, aerospace and nuclear.This phenomenon has very serious consequences, lead to the need for analysing more and better understanding the mechanical behaviour of structures, with taking into account the effects of the fort and the weak discontinuities, especially in critical conditions.It is a scientific challenge that represents an important issue, analytically as well as numerically.
The complexity of the analytical solutions even for simple cases imposes to setting up effective numerical methods to model the mechanical behaviour.Various studies have been developed using different methods, such as extended finite element method (XFEM) [1], boundary element method (BEM) [2] , element free Galerkin method (EFGM) [3], and other methods [4,5].Recently, a large field was opened by Hughes et al. [6] offers the possibility of exploiting computer aided design (CAD) tools in the analysis methods using isoparametric concept.This novel alternative method called isogeometric analysis (IGA).The basic idea of this method is use the computational geometry technologies as shape bases to describe exactly the geometry and also to approximate the unknown fields.Following this discovery, several researches in various fields have been investigated by this method, including: structural dynamics [7], composite materials [8], fluidstructure interaction [9], electromagnetic problems [10], contact problems [11], turbulent flow [12], aero-dynamics [13] and thermomechanical problems [14].However, in fracture mechanics problems, Benson et al. [15] and De Luycker et al. [16] have proposed extended isogeometric analysis (XIGA) for modelling cracks.In this method the general principle of the XFEM is used into IGA by including the asymptotic and signed distance enrichment functions.Therefore this method has the advantages of both XFEM and IGA, which are summarized by the ability to represent complex geometries independently of any discontinuities and without explicit meshing for obtain solution with higher orders.
There are many CAD basis functions can be used in isogeometric analysis, the most widely used of them are Non Uniform Rational B-splines (NURBS) due to their properties like continuity, smoothness, variation diminishing, convex hull and possibility of using Knot insertion and degree elevation refinements.They have the ability to describe exactly all conic sections but they are not well for all complex geometries, even for multiple patches NURBS generate a complicated mesh of control points and this leads to produce superfluous control points.In order to handle these disadvantages Sederberg et al. [17] proposed T-splines as a generalized tool of NURBS, in which the index space (T-mesh) cans have T-junctions and locally refined [18].Therefore the major advantages of this technique are local refinement and ability to represent complex geometries with minimal number of control points compared with those used in NURBS.
According to their ability in engineering design, T-splines have been used by the analysts to serves as basis functions for isogeometric analysis in many advanced searches.However T-splines bases are not always valid to use in analysis for different geometric configurations, because there are not rules ensure these bases to be linear independent and form a partition of unity.Li et al. [19] introduced analysis-suitable T-splines, wherein for any choice of knot vectors the blending functions are linearly independent.Like NURBS bases, analysis-suitable Tspline bases have the important properties of the analysis basis functions.Moreover, they provide an efficient algorithm which allows making highly localized refinement [20].
In this paper, plane crack problem is analysed numerically using XIGA, which is enhanced by adopting analysis-suitable T-splines in order to approximate the solution, construct the geometry and make the local refinement around the discontinuities.The asymptotic functions are used to model the crack tip, and the M-integral is used to evaluate the stress intensity factors.

Analysis-suitable T-splines
Analysis-suitable T-splines are a subset of Tsplines which have a restricted T-mesh topology; it is essentially based on the positions of the T-junctions in the Tmesh.And that by representation extensions in all Tjunctions, if there are intersections between T-junction extensions, then the extended T-mesh is not analysis suitable T-splines.A T-junction extension is composed by extension of two line segments, the first is from the Tjunction until one edge or vertex in the opposite direction of a missing edge, and the second is from the T-junction until two edges or vertices in the other direction, an example is shown in Fig. 1.

Building an analysis-suitable T-spline
For a given analysis-suitable T-mesh and net of control points (Pα), an analysis-suitable T-spline surface is given by: where Ξα and Hα which correspond to the anchor α (e.g., Fig. 2 shows local knot vectors of the anchor β), her formula is written in terms of B-spline bases Ni,p and Mj,q as: k is the number of the anchors in T-mesh, wα is a set of control point weights, p and q are orders of the basis functions corresponding to the parametric directions ξ and η, respectively.The B-spline basis functions in one parametric direction are given by and

Local refinement of analysis-suitable T-splines
In order to make a local refinement in T-splines, knots must be inserted to selected local knot vectors, so some basis functions are subdivided [18].For example, in the case of cubic orders (as in the present study) and one inserted knot z to a local vector, two basis functions are produced.We can combine them linearly to form the original basis function as: where For more than one inserted knot, the application of these equations is repeated until derive all coefficients.
The original basis functions B 1 and the generated basis functions B 2 of a T-spline space can be combined using the refinement operator M in linear system as: where N1 and N2 are the column vectors of basis functions of the original and the refined T-spline space, respectively.A special algorithm of refinement was introduced by Scott et al. [20]  For more detail refer you to [20].

Governing equation
Consider a cracked body Ω with outer boundary Γ subjected to a uniform body forces f b , traction forces f t applied at Γt and displacement conditions applied at Γu, in the state of equilibrium.The crack boundary Γc is considered to be traction free Fig. 3.
The strong form of the equilibrium equation and the boundary conditions are defined as [1]: = on The variational formulation of the boundary value problem is defined as: where σ is the stress tensor, ε is the strain tensor and n is the unit outward normal.
Fig. 3 An arbitrary cracked body with boundary conditions

Discretization
The Extended isogeometric analysis (XIGA) [15,16] uses the same methodology of the extended finite element method (XFEM) for modelling the discontinuities but with basis functions derived from the geometry like in isogeometric analysis [6].For the crack problems, XIGA provides the possibility of modelling the crack independently of the mesh and within a geometry presented exactly.Uncommonly, in this study T-splines are adopted where R is the analysis-suitable T-spline basis function Eq. ( 2), H is the Heaviside function used to represent crack discontinuity, it equals 1 and -1 above and below the crack, respectively.Fl are the asymptotic enrichment functions used to reproduce the singular field near the crack tips.u, a and b are vectors of degrees of freedom corresponding to classical (ns), crack face (ncf) and crack tip (nct) control points, respectively.
For isotropic domains, the crack-tip enrichment function is defined by four components in each point as: XIGA discretization uses the same XFEM procedure to discrete Eq. ( 13) where K is the global stiffness matrix, f is applied forces vector and u is the displacement vector.For one finite element, they are defined as: C is the matrix of elastic constants and B is the matrix of shape function derivatives which is given by:

Stress intensity factor computation
For approximate this parameter we used the interaction Integral, is based on the definition of an auxiliary state to extract mixed mode stress intensity factors (SIF) as following: where J act and J aux are the actual and auxiliary states J integrals, and J is the J integral value for the superposition state, her general form is written as [21]: and M is the interaction integral defined as where κ is the strain energy density, ℓ is an arbitrary contour surrounding the crack tip, nj is the jth component of the outward unit normal to ℓ, δ1j is the Kronecker delta, g is a smoothly function varies from 1 on the crack-tip to 0 on the edge of the contour ℓ, and A is the area inside ℓ.
where state 1 and state 2 defined as  

Numerical examples
In this section, two numerical examples are simulated by programming all the above techniques and methods into one code.The geometric models of the below examples are constructed using cubic order in both directions.First, a rectangular plate with double edge cracks is chosen in order to study the convergence and the domain independence in the computations of SIF, and then a disk with center crack is described exactly and analyzed for different crack lengths and angles in order to verify the accuracy of the proposed approach.Four types of finite elements are distinguished in these examples according to their positions with the crack, the standard element contains 3×3 Gauss points, the element that having tip enriched control points contains 5 × 5 Gauss points and the sub-triangle technique is used for the tip-element by 7 Gauss points in each triangle, however the split element contains 6 × 6 Gauss points for the double edge crack problem and the sub-triangle technique is used by 7 Gauss points in each triangle for the cracked disk problem.

Double edge cracked specimen under uniaxial tension
Consider a tensile plate of width 2w = 10 and height 2h = 20 with double edge cracks of length a = 3, as shown in Fig. 4. The analytical normalized SIF solution of this problem is given by: () where TI (a / w) is the analytical formula corresponding to the mode I, which can be computed as [22] 5) are used to study the convergence of the proposed approach with normalized M-integral radius equal to 1. Table 1 shows the obtained values with their errors, and Table 2 compares the results for different domain radius to study the domain independence in the case of 503 elements (Fig. 5, c).
According to Table 1, analysis suitable T-splines give us precise results for different numbers of control points and that attributed to the local refinement property.We can observe from Table 2 that the SIF values are almost not sensitive at all to the radius of the M-integral.A disk with an inclined central crack subjected to compression load is considered for both fracture modes, as depicted in Fig. 6.Different inclined angles are tested, namely, φ = 0°, 15°, 30°, 45°, 60°, 75° and 90° for different crack lengths.The analytical SIFs of this problem can be obtained by the following equations: where σ * is the characteristic stress (σ * = P0 / (πDt)), γ is crack length ratio (γ = a / D), t is disk thickness, TI and TII are the geometric functions corresponding to the mode I and II, respectively (all the geometric functions are taken from [23]).
In order to evaluate the mixed mode stress intensity factors in the case of a = 30 mm, we used a mesh consists of 1337 control points and 1216 elements, as shown in Fig. 7, a.The crack tip and the crack face enriched control points are defined like in Fig. 8.The obtained results are compared with those of the reference [23] in Fig. 9.For different crack length ratios we used a uniform mesh (see Fig. 7, b) to show the variations of KI and KII in Fig. 10.According to Fig. 9, the results obtained from the present approach are in good agreement with the analytical results and that probably due to the ability of the T-splines to describe such geometries exactly.The mode I stress intensity factor reduces steadily with an increase in crack angle (Fig. 10, a), whereas the mode II stress intensity factor increases and reaches its maximum value at (φ = 45° for γ = 0.2, 0.3 and φ = 30° for γ = 0.4, 0.5, 0.6, 0.7 and φ = 15° for γ = 0.8), and then decreases (Fig. 10, b.Therefore the crack angle of the maximum value of the mode II stress intensity factor depends on the crack length ratio.

Conclusion
In the present study, extended isogeometric analysis has been enhanced by analysis-suitable T-spline for modeling and analyzing cracks in plane problems.A compatible refinement algorithm has been contributed to make the local refinement and avoid the emergence of superfluous control points.The obtained stress intensity factors have been compared with analytical solutions, a good agreement proved the accuracy and efficiency of the proposed method.

Fig. 1 aFig. 2
Fig. 1 a -T-mesh; b -extended T-mesh (black points are Tjunctions); c -Analysis-suitable T-mesh (formed by adding the red edges) and in order to make the local refinement of analysis-suitable T-spline spaces.It consists of the follow-Ξβ Hβ β ing steps: -Create the refined T-mesh T2 from the original analysis-suitable T-mesh Ts1.-Form the extended T-mesh of T2.-If the extended T-mesh of T2 has intersecting Tjunction extensions, one edge must inserted into T2 in such a way that reduce the number of the intersections.-Repeat step 3 until the extended T-mesh has no Tjunction extensions intersect.-Compute the refinement matrix M.
in XIGA using analysis-suitable T-splines to approximate the displacement in any point ζ = (x1, x2) as follows ,θ) are local crack tip polar coordinates.
Young's modulus and ν is Poisson's ratio).The mode I and II stress intensity factors can be obtained as follows:

Fig. 4
Fig. 4 Finite rectangular plate with double edge cracks Different mesh configurations consist of 119, 251, 503, 655 and 1021 elements (all shown in Fig.5) are used to study the convergence of the proposed approach with normalized M-integral radius equal to 1. Table1shows the obtained values with their errors, and Table2compares the results for different domain radius to study the domain independence in the case of 503 elements (Fig.5, c).According to Table1, analysis suitable T-splines give us precise results for different numbers of control points and that attributed to the local refinement property.We can observe from Table2that the SIF values are almost not sensitive at all to the radius of the M-integral.

Table 1
Convergence of the SIF for various control nets