# Application of the Extended Isogeometric Analysis (X-IGA) to Evaluate a Pipeline Structure Containing an External Crack

Authors: S. El Fakkoussi; H. Moustabchir; A. Elkhalfi; C. I. Pruncu
Journal: Journal of Engineering
 (2018)

---

## Abstract

This work proposes a novel strategy for a two-dimensional problem that includes the approach of extended isogeometric analysis (X-IGA) in order to detect the behavior of a crack in pipeline structures. The nonrational B-Spline uniform function (NURBS) was used for the approximation of the solution fields (displacements) taking into account its geometry constrains. The modeling of the X-IGA was implemented under Abaqus/Standard software via subroutine (UEL) where the Stress Intensity Factor (KI) was extracted. The results permit detecting with accuracy the fracture toughness of a pipeline structure containing an external crack that can be submitted to critical pressures. To validate the performances of the novel strategy a careful comparison with existing literature and analytical and numerical computation methods was performed.


---

## Body

## 1. Introduction

The mechanical evaluation of a metal structure that contains cracks still represents an important research topic in the aeronautical, nuclear, and naval fields. Further preventive detection/prediction of a crack remains a major concern for researchers and engineers. An early detection makes it possible to minimize the cost of the corrective maintenance and to avoid important material and human damage. Because of this challenge, many researchers have devoted their efforts to study the phenomenon of crack propagation in metal structures. They proposed modern and reliable tools to quantify the damage rate near the crack. The literature highlights the existence of some methods based on the Stress Intensity Factor (KI) concept that permits evaluating with accuracy the fracture of the structures. This factor can be determined analytically, semianalytically, and by numerical computation methods. The numerical methods are based on the classical finite element method (FEM) and the extended finite element method (X-FEM) [1]. The X-FEM method is considered a powerful tool in terms of its precision to calculate the KI factor; it makes it possible to model the crack without taking into consideration the mesh near the crack tip. To enhance its robustness, in the X-FEM an enrichment function is added with respect to the conventional finite element formulation. However, this method has some limitations when an elliptical crack is simulated and is unable to model the semielliptic or circular cracks. The KI factor generates instabilities on the mentioned above cases [2, 3]. To avoid this problem that generates errors in the modeling of displacement and its geometry, Hughes et al. [4] proposed a Computer Aided Design (CAD) approach based on the finite elements procedure. In this routine the Lagrange form functions with the nonrational B-Spline uniform function (NURBS) is replaced. This approach is known as isogeometric analysis (IGA). These functions permit checking the interpolation conditions; however, the NURBS functions are always positive.The new IGA approach has been well implemented in several fields such as the heat transfer [5, 6], vibration analysis [7, 8], fluid flow analysis [9, 10], and wave propagation [11]. Besides, the IGA method showed promising results in the composite structures [12, 13], in thin-shell structures [14], and for crack propagation simulations [15]. In the fracture mechanics, the IGA approach was validated by comparison to the extended finite elements method (X-FEM) and the classical finite elements method (FEM). Sanches et al. [16] proposed a new function named i-Spline while imposing some modifications to the B-Spline function. The proposed method was implemented in two-dimensional problems and was integrated in the finite elements calculation using Matlab software.Bornemann and Cirak [17] identified a new technique that can be implemented in the B-Spline function by using the finite element method. The new techniques are applied to obtain convergence on the linear and geometrically nonlinear, two- and three-dimensional problems, using Matlab and Abaqus/Standard software. Choi and Cho [18] calculated the KI factor for a semielliptic crack using the NURBS function in the finite element method (MEF) formulation. Several researchers started to implement the NURBS function in the Abaqus/Standard software via subroutine (UEL) in the problems case without crack [19–21].In the modeling of some discontinuities problem like cracks, the IGA method has been coupled with the X-FEM method to generate a novel approach known as X-IGA. The new approach is proposed to understand the crack propagation phenomenon and fatigue of metal structures. Ghorashi et al. [22] investigated stationary and propagation of cracks using the combination of IGA approach and X-FEM method. They validated the results for a two-dimensional plane crack on a plate. Later, they compared the Stress Intensity Factor (KI) obtained with the analytical values and with the values of X-FEM method. Jung and Taciroglu [23] have used the natural cubic spline function in the X-FEM formulation to evaluate a semielliptic crack propagation on a two-dimensional plate. Bhardwaj et al. [24] studied the benefits of using the NURBS functions by applying the X-FEM method with the main objective of validating the X-IGA method. Based on the X-IGA approach, Singh et al. [25] used Bézier extraction configured on the T-Spline function to treat cracks in the two-dimensional specimen. On the other hand, Tinh Quoc Bui [26] used an extended isogeometric analysis (X-IGA) for the simulation of two-dimensional fracture mechanics problems in electromagnetic piezoelectric materials. A robust research on understanding of the thermal buckling analysis of functionally graded plates with internal defects was considered later by [27]. Bui et al. [28] have developed an innovative dynamic extended isogeometric analysis (X-IGA) for transient fracture of magnetoelectroelastic crack (MEE) solids under coupled electro-magneto-mechanical loading. The X-IGA approach permits obtaining a greater accuracy with respect to X-FEM methods [29] and conventional FEM approach [30], in terms of errors convergence and modeling strategy. Therefore, the novel X-IGA approach can be a robust tool to evaluate the fracture of metal structures.It is noticed that the IGA and X-IGA approaches are implemented only in the Matlab software. So far, from our knowledge, this was partially implemented in Abaqus/Standard software and this is done only for structures without cracks. In that reason, we propose to extend the novel X-IGA approach 0068 in the evaluation of the structure containing a crack and presenting a particular case to calculate the Stress Intensity Factor (KI) using the X-IGA approach for a pipeline containing an external crack. To validate our model a comparison with the Stress Intensity Factor (KI) results obtained via an analytical method and numerical methods (FEM and X-FEM) was implemented. The proposed method implemented under the Abaqus/Standard software offers a modern approach to investigate the fracture mechanisms for industrial site where the complex problems are encountered. The coupling between the two powerful methods IGA and X-FEM provides a new strategy capable of modeling cracked structures without remeshing and to model complex geometries. This approach will allow reducing considerable the computing time with direct benefit of improving the cost.


## 2. Theoretical Background

In this section, IGA and X-IGA as well as the analytical solution for SIF finite plate and an annular tube were briefly reviewed. More details of IGA and X-IGA and the LEFM are referred to the literature [4, 24, 33].


### 2.1. Isogeometric Analysis (IGA)

#### 2.1.1. B-Spline Function

The B-Spline functions are written in the form of Cox-de Boor recursion formula [34]:1.
Forp=0(1)Ni,0ξ=1ξi≤ξ≤ξi+10otherwise2.
Forp≥1(2)Ni,pξ=ξ-ξiξi+p-ξiNi,p-1ξ+ξi+p+1-ξξi+p+1-ξi+1Ni+1,p-1ξThe derivative of the B-Spline function is written as follows:(3)dNi,pξdξ=pξi+p-ξiNi,p-1ξ-pξi+p+1-ξi+1Ni+1,p-1ξ


#### 2.1.2. Nonuniform Function Rational B-Spline (NURBS)

NURBS are a generalization of B-Splines based not on polynomials but on rational functions [34]. The advantage of using the NURBS function and not the B-Spline function is that the former is able to represent complex shapes such as conic sections circles and elliptical.(4)Ri,pξ=ωA·Ni,pξ∑B=1nωA·NB,pξ(5)RAp,qξ,η=ωA·Ni,pξ·Mj,qη∑i=1n∑j=1mωA·Ni,pξ·Mj,qηN i , p ( ξ ) and Mj,q(η) are the basic functions, B-Splines with p and q are the order of the interpolations function, n and m are the number of control points and are the weights corresponding to the control points, and T=ζ1,ζ2,…,ζn+p+1 and H=η1,η2,…,ηn+q+1 are knots vectors (nodes).


#### 2.1.3. Surface NURBS

The surface NURBS functions are represented by the following equation:(6)Sζ,η=∑A=1n·mPA·RAp,qξ,ηP A is the control points and RAp,q(ξ,η) is the bivariate NURBS basic functions.


### 2.2. Extended Isogeometric Analysis (X-IGA)

The advantage of using the X-FEM method [1] together with IGA approach [4] is that the first method is able to model the crack in the structures using the enrichment functions and the second one is able to model geometries. The coupling between two methods is made by using NURBS functions.The approximation of displacementu(ξ) for a typical point ξ=(ζ,η) can be written in generalized form(7)uξ=∑inRiξui+∑jncfRjξHξ-Hξjaj+∑knctRkξ∑α=14ψαξ-ψαξkbαkR i is NURBS basic function, ui is the nodal shift in the classical method, H(ξ) is the Heaviside function, and ψα(r,θ) are the enrichment functions near the crack tip and are additional degree of freedom for Heaviside function and enrichment function in of the crack tip. ncf and nct are, respectively, the number of the points of the enriched controls for the Heaviside functions and the functions of crack tip.The Heaviside function is defined by(8)Hξ=+1φξ≻00φξ=0-1φξ≺0The functions of enrichments in the crack tip are defined by(9)ψjr,θ=rsinθ2,rcosθ2,…rsinθ2sinθ,rcosθ2sinθ


### 2.3. IGA/X-FEM Implementation in Abaqus/Standard Software

Duval et al. [35] implemented the IGA method, in the Abaqus/standard software. The elastic domain was considered for the crack case. They used a subroutine UELMAT written in Fortran language as a main code programming. There the NURBS functions were added as secondary file in the form (∗ NB). This file contains all the necessary input information of the NURBS functions such as the vector knot, the weight functions, and the order of the function. The present work will apply this strategy to implement the IGA approach coupled with the X-FEM method in the Abaqus/Standard software. Subsequently, the functions NURBS are implemented in Abaqus/Standard software via the UEL subroutine (UEL) that allows modifying the formulation of the finite elements. Some other parameters of the method (X-FEM) were added in this code, such as the functions of enrichment at the crack tip and the Heaviside function. The main code (UEL) is detailed in Algorithm 1 at the end of this document. To facilitate the addition of extra information to the main code, we have used the subroutine uexternaldb. This later on allows adding the necessary information of the crack and information of the NURBS function. These parameters are used later in the main program (UEL). The details of program implantation are summarized in Figure 1.Algorithm 1
SUBROUTINE UEL(RHS,AMATRX , SVARS , ENERGY,NDOFEL,NRHS,NSVARS, 1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE...)  INCLUDE 'ABA_PARAM.INC'  DIMENSION RHS(MLVARX, ∗ ), AMATRX(NDOFEL,NDOFEL), PROPS( ∗ ), 1 SVARS(NSVARS), ENERGY ( 8 ) , COORDS(MCRD,NNODE),...)  C  Read real and integer properties set at the ABAQUS input file  E  = PROPS ( 1 )  !!! Module Young Nu  = PROPS ( 2 )  !!! coeifficent poisson Typ_Prb  =  JPROPS ( 1 )  !!! Plane strain/stress  C ============================================================================== C  Defining Gauss points coordinates and Gauss Gauss weights C ==============================================================================  CALL Gauss ( NbPtInt,DimenProb,GaussPtsCoord,0)  C ============================================================================== C  Stiffness matrix and force vector are initialized to zero C ==============================================================================  DO i = 1, NDOFEL RHS(i,1) = zero DO j = 1, NDOFEL AMATRX( j,i ) = zero ENDDO ENDDO  C ============================================================================== DO Gauss integration point loop C ==============================================================================  C  Calculation of the matrix (elasticity)  CALL CALC_D( Typ_Prb,D,E,Nu )  CALL Shap_Funct_NURBS ( COORDS,dRdx,NNODE,R,GaussPtsCoord (2:,n), !!! Eq ( 5)  1  DetJac,NDOFEL,DimenProb,JELEM,NumPatch )  C  Crack tip derivates Functions Psi calculation at Gauss and nodal points  CALL Psi_CrackTip ( NbrCrackPs,CoorCrack,dRdx,Xcoord,Ycoord,DerPsi , 1  DerPsinodal )  !!! Eq ( 9 )    C  DerPsi  : Function Derivatives at the end of crack in gauss points C  DerPsinodal  : Function Derivatives at the end of crack in nodal points C  CoorCrack  : Coordinates of crack points C  Heaviside Function H calculation at Gauss and nodal points  CALL Heaviside_Func ( NnoeudX,Distance,NNODE,ix,dRdx,Heav,Heavnodal ) !!! Eq ( 8 )  C  Heav  : Function Heaviside in point x,y C  Heavnodal : Fonction Heaviside in elements nodes  C ============================================================================== C  Assembly of matrix [B] Derivatives the functions NURBS in coordinates points x, y C ==============================================================================  B=0 Iter =1  C  Loop over the nodes  DO i= 1,NNODE  C  The matrix [B] in the standard case part without crack  B(1,Iter)  = dRdx (1,i,1) B(2,Iter+1)  = dRdx (2,i,1) B(3,Iter)  = dRdx (2,i,1) B(3,Iter+1)  = dRdx (1,i,1)  c  The matrix B in the case of crack compute by heaviside function  IF  ( NodesInfo (i) = 1 )  THEN !!! nodes enriched by heaviside functions  B(1,2+Iter)  = ( Heav-Heavnodal (i)) ∗ dRdx(1,i,1) B(2,3+Iter)  = ( Heav-Heavnodal (i)) ∗ dRdx(2,i,1) B(3,2+Iter)  = ( Heav-Heavnodal (i)) ∗ dRdx(2,i,1) B(3,3+Iter)  = ( Heav-Heavnodal (i)) ∗ dRdx(1,i,1) c  The matrix B in the case of crack compute by the functions (Psi) at the points of the crack  ELSEIF  ( NodesInfo (i) = 2 )  THEN !!! nodes enriched by Crack tip functions DO j= 1,4 B(1,2 ∗ j+2+Iter)= DerPsi (i,1,j)- DerPsinodal (i,1,j) B(2,2 ∗ j+3+Iter) = DerPsi (i,2,j)- DerPsinodal (i,2,j) B(3,2 ∗ j+2+Iter) = DerPsi (i,2,j)- DerPsinodal (i,2,j) B(3,2 ∗ j+3+Iter)= DerPsi (i,1,j)- DerPsinodal (i,1,j ) ENDDO END IF  Iter =Iter+12  ENDDO !!! i End loop over the element nodes  C ============================================================================== C  Assembking RHS and AMATRIX C ==============================================================================  AMATRX = Ke RHS = − Fint + Fext  C ============================================================================== ENDDO !!! end intergartion gauss C ==============================================================================    c Calculate and store in SVARS the values of strains, stresses, strain energy, dRdX ,  ...   CALL SVARS_Sub (JTYPE,JELEM,SVARS,NSVARS,U,NDOFEL,...)  RETURN END SUBROUTINE UELFigure 1
Implementation methodology of X-IGA approach in Abaqus/Standard software.The projection of the results such as the stresses and strains determined by the X-IGA approach is based on the technique of Bézier composition, technique that is detailed in [36]. The results extraction and the main fracture mechanic variables are determined using Fortran algorithm. It will help to release (10) and to calculate later the Stress Intensity Factor (KI) as stated in (13).


### 2.4. Calculation of the Stress Intensity Factor (KI)

The computation of the Stress Intensity Factor (KI) is based on the integral interaction method [37]. The mathematical formulation of this method for a two-dimensional problem is determined in the following way:(10)M1,2=∫Aσ1ijδu2iδx1+σ2ijδu1iδx1-W1,2δijδqδxjdAq is weight function, q=1 in the crack tip, and q=0 in the contour around the crack(11)W1,2=12σ1ijε2ij+σ2ijε1ij=σ1ijε2ij=σ2ijε1ij(12)M1,2=2E∗K1IK2I+K1IIK2IIE ∗ = E for plane strain and E∗=E/1-μ2 for plane stress.For mode I, we haveK(2)I=1 and K(2)II=0.For mode II, we haveK(2)I=0 and K(2)II=1.Finally, we get(13)K1Ι=M1,Ι·E∗2The theoretical Stress Intensity Factor (KI), for the case of a two-dimensional plate with a crack, is written as follows:(14)KI=Fσπcwhere(15)F=1,12-0.23cL+10,55cL2-21.72cL3+30.39cL4


## 2.1. Isogeometric Analysis (IGA)

### 2.1.1. B-Spline Function

The B-Spline functions are written in the form of Cox-de Boor recursion formula [34]:1.
Forp=0(1)Ni,0ξ=1ξi≤ξ≤ξi+10otherwise2.
Forp≥1(2)Ni,pξ=ξ-ξiξi+p-ξiNi,p-1ξ+ξi+p+1-ξξi+p+1-ξi+1Ni+1,p-1ξThe derivative of the B-Spline function is written as follows:(3)dNi,pξdξ=pξi+p-ξiNi,p-1ξ-pξi+p+1-ξi+1Ni+1,p-1ξ


### 2.1.2. Nonuniform Function Rational B-Spline (NURBS)

NURBS are a generalization of B-Splines based not on polynomials but on rational functions [34]. The advantage of using the NURBS function and not the B-Spline function is that the former is able to represent complex shapes such as conic sections circles and elliptical.(4)Ri,pξ=ωA·Ni,pξ∑B=1nωA·NB,pξ(5)RAp,qξ,η=ωA·Ni,pξ·Mj,qη∑i=1n∑j=1mωA·Ni,pξ·Mj,qηN i , p ( ξ ) and Mj,q(η) are the basic functions, B-Splines with p and q are the order of the interpolations function, n and m are the number of control points and are the weights corresponding to the control points, and T=ζ1,ζ2,…,ζn+p+1 and H=η1,η2,…,ηn+q+1 are knots vectors (nodes).


### 2.1.3. Surface NURBS

The surface NURBS functions are represented by the following equation:(6)Sζ,η=∑A=1n·mPA·RAp,qξ,ηP A is the control points and RAp,q(ξ,η) is the bivariate NURBS basic functions.


## 2.1.1. B-Spline Function

The B-Spline functions are written in the form of Cox-de Boor recursion formula [34]:1.
Forp=0(1)Ni,0ξ=1ξi≤ξ≤ξi+10otherwise2.
Forp≥1(2)Ni,pξ=ξ-ξiξi+p-ξiNi,p-1ξ+ξi+p+1-ξξi+p+1-ξi+1Ni+1,p-1ξThe derivative of the B-Spline function is written as follows:(3)dNi,pξdξ=pξi+p-ξiNi,p-1ξ-pξi+p+1-ξi+1Ni+1,p-1ξ


## 2.1.2. Nonuniform Function Rational B-Spline (NURBS)

NURBS are a generalization of B-Splines based not on polynomials but on rational functions [34]. The advantage of using the NURBS function and not the B-Spline function is that the former is able to represent complex shapes such as conic sections circles and elliptical.(4)Ri,pξ=ωA·Ni,pξ∑B=1nωA·NB,pξ(5)RAp,qξ,η=ωA·Ni,pξ·Mj,qη∑i=1n∑j=1mωA·Ni,pξ·Mj,qηN i , p ( ξ ) and Mj,q(η) are the basic functions, B-Splines with p and q are the order of the interpolations function, n and m are the number of control points and are the weights corresponding to the control points, and T=ζ1,ζ2,…,ζn+p+1 and H=η1,η2,…,ηn+q+1 are knots vectors (nodes).


## 2.1.3. Surface NURBS

The surface NURBS functions are represented by the following equation:(6)Sζ,η=∑A=1n·mPA·RAp,qξ,ηP A is the control points and RAp,q(ξ,η) is the bivariate NURBS basic functions.


## 2.2. Extended Isogeometric Analysis (X-IGA)

The advantage of using the X-FEM method [1] together with IGA approach [4] is that the first method is able to model the crack in the structures using the enrichment functions and the second one is able to model geometries. The coupling between two methods is made by using NURBS functions.The approximation of displacementu(ξ) for a typical point ξ=(ζ,η) can be written in generalized form(7)uξ=∑inRiξui+∑jncfRjξHξ-Hξjaj+∑knctRkξ∑α=14ψαξ-ψαξkbαkR i is NURBS basic function, ui is the nodal shift in the classical method, H(ξ) is the Heaviside function, and ψα(r,θ) are the enrichment functions near the crack tip and are additional degree of freedom for Heaviside function and enrichment function in of the crack tip. ncf and nct are, respectively, the number of the points of the enriched controls for the Heaviside functions and the functions of crack tip.The Heaviside function is defined by(8)Hξ=+1φξ≻00φξ=0-1φξ≺0The functions of enrichments in the crack tip are defined by(9)ψjr,θ=rsinθ2,rcosθ2,…rsinθ2sinθ,rcosθ2sinθ


## 2.3. IGA/X-FEM Implementation in Abaqus/Standard Software

Duval et al. [35] implemented the IGA method, in the Abaqus/standard software. The elastic domain was considered for the crack case. They used a subroutine UELMAT written in Fortran language as a main code programming. There the NURBS functions were added as secondary file in the form (∗ NB). This file contains all the necessary input information of the NURBS functions such as the vector knot, the weight functions, and the order of the function. The present work will apply this strategy to implement the IGA approach coupled with the X-FEM method in the Abaqus/Standard software. Subsequently, the functions NURBS are implemented in Abaqus/Standard software via the UEL subroutine (UEL) that allows modifying the formulation of the finite elements. Some other parameters of the method (X-FEM) were added in this code, such as the functions of enrichment at the crack tip and the Heaviside function. The main code (UEL) is detailed in Algorithm 1 at the end of this document. To facilitate the addition of extra information to the main code, we have used the subroutine uexternaldb. This later on allows adding the necessary information of the crack and information of the NURBS function. These parameters are used later in the main program (UEL). The details of program implantation are summarized in Figure 1.Algorithm 1
SUBROUTINE UEL(RHS,AMATRX , SVARS , ENERGY,NDOFEL,NRHS,NSVARS, 1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE...)  INCLUDE 'ABA_PARAM.INC'  DIMENSION RHS(MLVARX, ∗ ), AMATRX(NDOFEL,NDOFEL), PROPS( ∗ ), 1 SVARS(NSVARS), ENERGY ( 8 ) , COORDS(MCRD,NNODE),...)  C  Read real and integer properties set at the ABAQUS input file  E  = PROPS ( 1 )  !!! Module Young Nu  = PROPS ( 2 )  !!! coeifficent poisson Typ_Prb  =  JPROPS ( 1 )  !!! Plane strain/stress  C ============================================================================== C  Defining Gauss points coordinates and Gauss Gauss weights C ==============================================================================  CALL Gauss ( NbPtInt,DimenProb,GaussPtsCoord,0)  C ============================================================================== C  Stiffness matrix and force vector are initialized to zero C ==============================================================================  DO i = 1, NDOFEL RHS(i,1) = zero DO j = 1, NDOFEL AMATRX( j,i ) = zero ENDDO ENDDO  C ============================================================================== DO Gauss integration point loop C ==============================================================================  C  Calculation of the matrix (elasticity)  CALL CALC_D( Typ_Prb,D,E,Nu )  CALL Shap_Funct_NURBS ( COORDS,dRdx,NNODE,R,GaussPtsCoord (2:,n), !!! Eq ( 5)  1  DetJac,NDOFEL,DimenProb,JELEM,NumPatch )  C  Crack tip derivates Functions Psi calculation at Gauss and nodal points  CALL Psi_CrackTip ( NbrCrackPs,CoorCrack,dRdx,Xcoord,Ycoord,DerPsi , 1  DerPsinodal )  !!! Eq ( 9 )    C  DerPsi  : Function Derivatives at the end of crack in gauss points C  DerPsinodal  : Function Derivatives at the end of crack in nodal points C  CoorCrack  : Coordinates of crack points C  Heaviside Function H calculation at Gauss and nodal points  CALL Heaviside_Func ( NnoeudX,Distance,NNODE,ix,dRdx,Heav,Heavnodal ) !!! Eq ( 8 )  C  Heav  : Function Heaviside in point x,y C  Heavnodal : Fonction Heaviside in elements nodes  C ============================================================================== C  Assembly of matrix [B] Derivatives the functions NURBS in coordinates points x, y C ==============================================================================  B=0 Iter =1  C  Loop over the nodes  DO i= 1,NNODE  C  The matrix [B] in the standard case part without crack  B(1,Iter)  = dRdx (1,i,1) B(2,Iter+1)  = dRdx (2,i,1) B(3,Iter)  = dRdx (2,i,1) B(3,Iter+1)  = dRdx (1,i,1)  c  The matrix B in the case of crack compute by heaviside function  IF  ( NodesInfo (i) = 1 )  THEN !!! nodes enriched by heaviside functions  B(1,2+Iter)  = ( Heav-Heavnodal (i)) ∗ dRdx(1,i,1) B(2,3+Iter)  = ( Heav-Heavnodal (i)) ∗ dRdx(2,i,1) B(3,2+Iter)  = ( Heav-Heavnodal (i)) ∗ dRdx(2,i,1) B(3,3+Iter)  = ( Heav-Heavnodal (i)) ∗ dRdx(1,i,1) c  The matrix B in the case of crack compute by the functions (Psi) at the points of the crack  ELSEIF  ( NodesInfo (i) = 2 )  THEN !!! nodes enriched by Crack tip functions DO j= 1,4 B(1,2 ∗ j+2+Iter)= DerPsi (i,1,j)- DerPsinodal (i,1,j) B(2,2 ∗ j+3+Iter) = DerPsi (i,2,j)- DerPsinodal (i,2,j) B(3,2 ∗ j+2+Iter) = DerPsi (i,2,j)- DerPsinodal (i,2,j) B(3,2 ∗ j+3+Iter)= DerPsi (i,1,j)- DerPsinodal (i,1,j ) ENDDO END IF  Iter =Iter+12  ENDDO !!! i End loop over the element nodes  C ============================================================================== C  Assembking RHS and AMATRIX C ==============================================================================  AMATRX = Ke RHS = − Fint + Fext  C ============================================================================== ENDDO !!! end intergartion gauss C ==============================================================================    c Calculate and store in SVARS the values of strains, stresses, strain energy, dRdX ,  ...   CALL SVARS_Sub (JTYPE,JELEM,SVARS,NSVARS,U,NDOFEL,...)  RETURN END SUBROUTINE UELFigure 1
Implementation methodology of X-IGA approach in Abaqus/Standard software.The projection of the results such as the stresses and strains determined by the X-IGA approach is based on the technique of Bézier composition, technique that is detailed in [36]. The results extraction and the main fracture mechanic variables are determined using Fortran algorithm. It will help to release (10) and to calculate later the Stress Intensity Factor (KI) as stated in (13).


## 2.4. Calculation of the Stress Intensity Factor (KI)

The computation of the Stress Intensity Factor (KI) is based on the integral interaction method [37]. The mathematical formulation of this method for a two-dimensional problem is determined in the following way:(10)M1,2=∫Aσ1ijδu2iδx1+σ2ijδu1iδx1-W1,2δijδqδxjdAq is weight function, q=1 in the crack tip, and q=0 in the contour around the crack(11)W1,2=12σ1ijε2ij+σ2ijε1ij=σ1ijε2ij=σ2ijε1ij(12)M1,2=2E∗K1IK2I+K1IIK2IIE ∗ = E for plane strain and E∗=E/1-μ2 for plane stress.For mode I, we haveK(2)I=1 and K(2)II=0.For mode II, we haveK(2)I=0 and K(2)II=1.Finally, we get(13)K1Ι=M1,Ι·E∗2The theoretical Stress Intensity Factor (KI), for the case of a two-dimensional plate with a crack, is written as follows:(14)KI=Fσπcwhere(15)F=1,12-0.23cL+10,55cL2-21.72cL3+30.39cL4


## 3. Numerical Simulation

To validate the novel proposed strategy in the Abaqus/Standard software, we have used a two-dimensional plate with an axial crack. The dimension of this plate is L = 5 mm and the crack length a = 3 mm. The configuration of the investigated plate was sketched in Figure2. Then, a half-tube characterized by an inner radius and outer radius Ri = 10 mm, Re = 20 mm, and a crack length a = 5 mm was introduced to validate the novel approach. The major details of a tube containing an axial crack and its boundary conditions were introduced in Figure 3. The model studied is a static linear one, where the objective was to extract the values of the Stress Intensity Factor (KI). The NURBS functions used in this study are of the cubic order (p = q = 3). The number of elements used for the two-dimensional plate case is, respectively, 121, 441, and 676. In addition, the number of elements used for the case of a half-tube are, respectively, 470 and 767. The mechanical properties of the material investigated in this study are presented in the Table 1 [31].Table 1
Mechanical properties of the material [31].
Module d’Young E= 207000 [MPa]  Poisson coefficient ν =0.3  Yield strength Re= 340 [MPa]  Ultimate tensile strength Rm= 440 [MPa]  Elongation at break A=35 %Figure 2
A two-dimensional plate with an axial crack.Figure 3
Annular tube with an axial crack and its boundary conditions.


## 4. Results and Discussions

The execution of X-IGA approach was performed under the Abaqus/Standard software via UEL subroutine and the postprocessing results were made by Fortran, allowing in this way computing the integral interaction equation (10) and the KI Factor (13). The obtained results are introduced in Table 2. The studied cases contain different meshes: a mesh of 121, 441, and 691 elements. Then, a comparison with the analytical, finite element method (FEM), and X-FEM methods was conducted. The obtained results are presented in Table 2. These results show a very good agreement with the other numerical results FEM and X-FEM that do not exceed an error of 0.270%. For the case of a half-tube problem with an external crack under pressure, we identified that the results of the new applied strategy show good agreement with literature results without exceeding an error of 1.03%. Details of the results are shown in the Table 3.Table 2
Comparison of the KI factors between the proposed method of X-IGA and the analytical, FEM, and X-FEM methods for the case of a crack with a = 3 mm.
KI [MPa.mm 0.5 ] Analytical FEM Method X-FEM Method X-IGA 121 elements X-IGA 441 elements X-IGA 676 elements 123.62 126.40 122.03 123.95 123.89 123.75 Error ( % ) 2.19 1.29 0.270 0.218 0.105Table 3
Comparison of the KI factors between the proposed method of X-IGA and the analytical, FEM, and X-FEM methods for the case of a crack with a = 5 mm.
KI [MPa.mm 0.5 ] Analytical FEM Method X-FEM Method X-IGA 470 elements X-IGA 767 elements 12.88[32] 13.79 13.69 13.18 13.015 Error ( % ) 6.59 5.91 2.27 1.03Figure4 shows the comparison of von Mises stresses for a half-tube with an external crack under pressure. The results obtained after simulation show a good agreement between the X-FEM method and the new X-IGA approach. The difference is due to the problem of the singularity near the crack tip, whereas the FEM method is unable to capture the constraints near the crack tip. On the other hand, we found that the constraints in the X-IGA method are close to the X-FEM method because of the enrichment functions that was already added to avoid the singularity problem and the constraints near the crack tip.Figure 4
The von Mises stress for the case of X-IGA, X-FEM, and FEM.
(a) (b) (c)Figure5 shows the variation of the Stress Intensity Factor values as a function of the pressure applied on the half-tube. The results of the X-IGA, X-FEM and analytical methods are close to each other and they vary linearly with the pressure. Instead, the FEM method results are far inaccurate compared to other methods, due to the weakness of the used method.Figure 5
Variation of the Stress Intensity Factors (KI) with applied pressure for a half-tube case.The variation of the Stress Intensity Factor values as a function of the crack length with respect to the thickness of the half-tube was performed as well. Figure6 shows that the depth of the crack increases while the KI factor also increases up to when it reaches a max value of 22.5 for a ration 0.8. The X-IGA results agree very well with to the analytical and the X-FEM method values.Figure 6
The variation of the Stress Intensity Factors with crack length for different crack sizes.The X-IGA approach implemented under the Abaqus/Standard software was successfully validated by corroborating its results against analytical and numerical computations models. The confidence generated by the proposed method can help to solve major fracture mechanics challenges and provides a better accuracy for an economical cost.


## 5. Conclusions

The present study introduces the novel X-IGA approach that was implemented under the Abaqus/Standard software via UEL subroutine. The postprocessing results, in terms of the Stress Intensity Factors (KI), were performed in the Fortran language. The proposed approach was validated by two different problems considering a linear elasticity case. The modeling was performed considering a two-dimensional plate containing an axial crack and a two-dimensional half-tube with an external crack under pressure. The obtained results were compared with the analytical and numerical methods, such as the FEM and X-FEM. The comparison between these results, in terms of the Stress Intensity Factors (KI), and the new proposed ones shows a very good agreement for the two cases. The comparison between their errors demonstrated that the new proposed X-IGA approach generates lower errors compared to other methods. The numerical errors values obtained by X-IGA do not exceed a value of 0.270% when using two-dimensional plate and 1.03% when using the half-tube.The novel proposed strategy permits using the X-IGA approach in computational software providing greater approximation for complex geometry and displacement solution. Besides, the use of the novel strategy permits solving the most challenged computation problems compared to the other methods, X-FEM and FEM, by using smaller elements with a high order.For future research, we are planning to implement the novel strategy (X-IGA) in order to solve a three-dimensional semielliptic crack problem and later to extend the existing strategy in the problem solving of fatigue crack with a more complex configuration.



---

*Source: 4125765-2018-10-17.xml*