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.

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 [1921].

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.For 2.For The derivative of the B-Spline function is written as follows:

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.

and 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 and are knots vectors (nodes).

2.1.3. Surface NURBS

The surface NURBS functions are represented by the following equation:

is the control points and 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 displacement for a typical point can be written in generalized form

is NURBS basic function, is the nodal shift in the classical method, is the Heaviside function, and 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. and 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 byThe functions of enrichments in the crack tip are defined by

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.

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, COORDS(MCRD,NNODE),...)
CRead real and integer properties set at the ABAQUS input file
E= PROPS!!! Module Young
Nu= PROPS!!! coeifficentpoisson
Typ_Prb= JPROPS!!! Plane strain/stress
C==============================================================================
CDefining Gauss points coordinates and Gauss Gauss weights
C==============================================================================
CALL Gauss ( NbPtInt,DimenProb,GaussPtsCoord,0)
C==============================================================================
CStiffness 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==============================================================================
CCalculation 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)
1DetJac,NDOFEL,DimenProb,JELEM,NumPatch)
CCrack tip derivates Functions Psi calculation at Gauss and nodal points
CALL Psi_CrackTip(NbrCrackPs,CoorCrack,dRdx,Xcoord,Ycoord,DerPsi,
1DerPsinodal)!!! Eq (9)
CDerPsi: Function Derivatives at the end of crack in gauss points
CDerPsinodal: Function Derivatives at the end of crack in nodal points
CCoorCrack: Coordinates of crack points
CHeaviside Function H calculation at Gauss and nodal points
CALLHeaviside_Func(NnoeudX,Distance,NNODE,ix,dRdx,Heav,Heavnodal)
!!! Eq (8)
CHeav: Function Heaviside in point x,y
CHeavnodal: Fonction Heaviside in elements nodes
C==============================================================================
CAssembly of matrix [B] Derivatives the functions NURBS in coordinates points x, y
C==============================================================================
B=0
Iter=1
CLoop over the nodes
DO i= 1,NNODE
CThe 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)
cThe 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)
cThe 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,2j+2+Iter)=DerPsi(i,1,j)-DerPsinodal(i,1,j)
B(2,2j+3+Iter)=DerPsi(i,2,j)-DerPsinodal(i,2,j)
B(3,2j+2+Iter)=DerPsi(i,2,j)-DerPsinodal(i,2,j)
B(3,2j+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==============================================================================
CAssembking 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 UEL

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:

is weight function, in the crack tip, and in the contour around the crack

for plane strain and for plane stress.

For mode I, we have and .

For mode II, we have and .

Finally, we get

The theoretical Stress Intensity Factor (KI), for the case of a two-dimensional plate with a crack, is written as follows:where

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 Figure 2. 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].

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.

Figure 4 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 5 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.

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. Figure 6 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.

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.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.