**1. Introduction**

The fracture phenomenon is a fundamental research topic in the field of solid mechanics, and a misunderstanding of the fracture mechanism may lead to a poor evaluation of a structure's integrity. Several experimental and theoretical studies have been conducted in this field covering this vast physical phenomenon. The cracking issue is among the problems that arise in this regard, so the study of the various kinds of discontinuities and defects that appear on the surface of structures is of major importance in the field of design and modeling for ensuring the reliability of engineering structures. These defects may appear in the form of notches, cracks, inclusions, holes, corrosion, and other types of material degradation. The existence of a defect on an engineering structure can cause material and human damage; this is caused by the propagation of cracks either in terms of direction or propagation rates. Therefore, it is necessary to predict the crack propagation rate to be able to estimate the critical load, and the acceptable length of the crack to ensure the stability of structures. The initiation and propagation of the crack requires a specific study of the defect, i.e., it is necessary to adopt a criterion of rupture in the context of fracture mechanics. Since the beginning of this research field, many methods have been developed to give a general vision of the fracture process in order to find adequate solutions and improve the strength of structures. With this technological progress, the simulation of the physical phenomena tends to become purely numerical; of course, experiments play a very important role, but due to their difficulty, and to save time and take advantage of the data-processing tools that exist today, we had to adapt to and develop in the world of digitalization. Numerical simulation and discontinuity problems still do not agree, since researchers are always looking to improve numerical solutions to be able to reflect reality.

**Citation:** Montassir, S.; Moustabchir, H.; Elkhalfi, A.; Scutaru, M.L.; Vlase, S. Fracture Modelling of a Cracked Pressurized Cylindrical Structure by Using Extended Iso-Geometric Analysis (X-IGA). *Mathematics* **2021**, *9*, 2990. https://doi.org/10.3390/ math9232990

Academic Editor: David Carfì

Received: 13 October 2021 Accepted: 16 November 2021 Published: 23 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In this sense, various numerical methods are developed that cite the classical finite element method (CFEM) [1,2], element-free Galerkin methods (EFG) [3], the extended finite element method (X-FEM) [4–6], the phase field numerical manifold method (PFNMM) [7], the boundary element method [8], and the thermo-mechanical peridynamic model (TMPM) [9]. Despite the huge use of the CFEM, however, it suffers from a certain lack of treatment of the cracking problems; it imposes that the mesh conforms with the crack or to a surface defect, and adopting a specific mesh creates difficulties at the simulation stage. The appearance of the enrichment approach [10–12] overcomes the shortcomings of the conventional method and is able to simulate all types of notches, cracks, holes, and other defects, regardless of the defect shape and mesh type. The Lagrange polynomials were used to make the interpolation with these new techniques, i.e., the geometry and the solution of the problem are approximated by these polynomials. Since the basis functions used are not the same as those used in the design, discretization errors appeared in the analysis stage [13]. The real issue occurred when we wanted to move from modeling to simulation by the process of the finite element method. Despite the integration of simulations in modeling software, such as CATIA, most software is still not good in this area, compared to software that is already simulation-integrated, such as ABAQUS. That is why there are still gaps between these two parts. This subject has become an area of interest for researchers in approximating conical shapes, as it is known that CAD software constructs this type of shape by B-spline curves. The advantage of these curves is that they are able to reproduce all types of conical shapes and free surfaces in an exact way. In this field, Hughes et al. [14] started to develop a new approach to eliminate the shortcomings of the classical method, and due to their research, iso-geometric analysis (IGA) was developed. On the other hand, industrial computational software has not yet integrated this functionality without some works that have made the implementation of these functions in specific fields of application, for example in LS-DYNA [15–17]. They were among the first in this field; precisely, they integrated IGA to study shells. In Altair RADIOSS [18], IGA was implemented and proved on industrial benchmarks. In ABAQUS [19–21], an application of IGA was implemented for linear elasticity problems. IGA has been employed in a variety of engineering fields, including the mechanics of vibration [22,23], fluid mechanics [24], electromagnetic problems [25], the medical field [26], and digital image correlation [27]. The issues addressed are diverse, including nonlinear mechanics [28], shell analysis [29–34], contact problems [35,36], fluid– structure interactions [24], the optimization of structural design [37], buckling failure [38], and crack problem analysis [39]. From this research, IGA has indicated the effectiveness of its results and it can sweep the world of numerical computation. Therefore, it can become an alternative method to the classical finite element method in the future.

Over the last few years, the mathematical formulation of the IGA method has been revised. In addition, it was expanded to fix issues concerning discontinuities, e.g., the crack, under the partition of the unity finite element method [40]. This new concept is known as the extended iso-geometric analysis (X-IGA). It has been utilized by a range of authors to address crack-growth issues in the field of fracture mechanics [41–47]. Two main functions have been added: the Heaviside jump function is used to enrich the displacement fields around the crack surfaces, while crack-tip enrichment functions are used to model the singularity at the crack tip. Most of the existing research works that use this innovative approach are limited to 2-D plate cracking problems [48,49] and cracking problems in the cantilever beam [50]. Therefore, the strategy has only been implemented with simple domain geometries. As for cracked shells, research is still in progress because of the difficulty of modeling the cracks and approximating the geometry with the NURBS functions. Furthermore, it requires a mathematical background. That is what made this kind of model interesting. In this field, X-IGA has been used for a cracked 2-D pipe. However, the method has been applied through FORTRAN language [51]. For the present investigation, we will study the efficiency and accuracy of X-IGA for cracks with 2-D pipe geometry given by CAD curves, with a special focus on ensuring that all stages of the calculation are done in MATLAB code. It is noted that mesh generation is used in MATLAB independently

of another meshing software. The strategy presented in this implementation follows the philosophy used in the traditional FE codes, and also benefits from having a MATLAB routine that allows for the discretization of partial differential equations based on NURBS and B-spline. Hence, the specific aims of this paper are to formulate the X-FEM concept in the framework of NURBS and to enrich the solution in MATLAB code to adapt it for a cracked pipe under pressure. To validate this strategy, a P264GH steel gas pipe has been used with uniform inner pressure. The mechanical characteristics data of this model are taken from an experimental study [52]. Based on the above literature review, stress intensity factors (SIF) are widely used to characterize fracture mechanics. Therefore, we evaluated the SIF with this implementation and the obtained results were compared with the CFEM available in ABAQUS software and the X-FEM approach, which was implemented in the ABAQUS user-defined element (UEL) [53]. Finally, to verify the results obtained by this strategy of X-IGA application, we made a comparison with [51].

The present study begins with a brief discussion of the IGA concept, including an assessment of the B-spline and NURBS basis functions. The concept of discontinuity inside a continuum formulation is outlined to provide a background for the following discussion on the X-IGA. Later, the methodology for implementing this improved approach is introduced. The paper concludes by providing case studies of a pipe under inner uniform pressure. The obtained results of the present study were evaluated against the CFEM and X-FEM solutions.

#### **2. Outline on B-splines, NURBS, and IGA Concepts**

Piecewise polynomial functions with a specified degree of continuity are known as B-splines. In a knot vector, vector Ξ = *ξ*1, *ξ*2,..., *ξn*+*p*+<sup>1</sup> , with *ξi*+<sup>1</sup> ≥ *ξi*; *i* is the index of knots, *n* is the number of control points, and *p* is the polynomial degree. A collection of coordinates in parametric space is used to create univariate B-Spline shape functions. The Cox-de Boor recursion model on the appropriate knot vectors determines the *i*th B-splines basis functions *Ni*,*p*(*ξ*) [14].

For a polynomial order *p* = 0,

$$N\_{i,0}(\xi) = \left\{ \begin{array}{c} 1 \text{ if } \xi\_i \le \xi \le \xi\_{i+1} \\ 0 \text{ otherwise} \end{array} \right\} \tag{1}$$

For *p* ≥ 1,

$$N\_{i,p}(\xi) = \frac{\mathfrak{z} - \mathfrak{z}\_i}{\mathfrak{z}\_{i+p} - \mathfrak{z}\_i} N\_{i,p-1}(\xi) + \frac{\mathfrak{z}\_{i+p+1} - \mathfrak{z}\_i}{\mathfrak{z}\_{i+p+1} - \mathfrak{z}\_{i+1}} N\_{i+1,p-1}(\xi) \tag{2}$$

The B-splines shape function's first derivative with respect to the parameter is:

$$\frac{d}{d\xi} = \left(N\_{i,p}\left(\xi\right)\right) = \frac{p}{\mathfrak{z}\_{i+p} - \mathfrak{z}\_i} N\_{i,p-1}\left(\xi\right) - \frac{p}{\mathfrak{z}\_{i+p+1} - \mathfrak{z}\_{i+1}} N\_{i+1,p-1}\left(\xi\right) \tag{3}$$

An example of a cubic B-spline basis function is illustrated in Figure 1. An open knot vector was used in this example due to the multiplicity of the first and the last knots, which are equal to *p +* 1. With IGA, the construction of the geometrical model is done with B-spline functions of open vectors. This type of vector allows for an interpolation of the basis function at the end; it is more than suitable for enforcing the boundary conditions.

**Figure 1.** An example of cubic B-splines basis functions with knot vector Ξ = (0, 0, 0, 0.5, 1, 1, 1).

#### *2.1. Curve and Surface Building with B-spline*

A linear combination of the shape functions and coefficients denoted by the control points constructs a B-spline curve:

$$\mathbb{C}(\xi) = \sum\_{i=1}^{n} N\_{i,p}(\xi) B\_i \tag{4}$$

*Ni*,*<sup>p</sup>* and *Bi* are the *i*th B-spline function and control points, respectively.

A B-spline surface is defined by a tensor product and parameterized by two-knot vectors, knot vector = *ξ*1, *ξ*2,..., *ξn*+*p*+<sup>1</sup> , *η*1, *η*2,..., *ηn*+*p*+<sup>1</sup> as:

$$S(\mathcal{J}, \eta) = \sum\_{i=1}^{n} \sum\_{j=1}^{m} N\_{i,p}(\mathcal{J}) M\_{j,q}(\eta) B\_{i\prime} \tag{5}$$

where *Ni*,*p*(*ξ*) and *Mj*,*q*(*η*) are univariate shape functions and B is the 3-vector of control point coordinates.

#### *2.2. Non-Uniform Rational B-splines (NURBS)*

NURBS is the generalization of the B-spline functions:

$$R\_{i,p}(\xi) = \frac{N\_{i,p}(\xi) \, w\_i}{\sum\_{l=1}^{n} N\_{l,p}(\xi) w\_l} \tag{6}$$

where *wi* > 0 is the weighting parameter and *Ni*,*p*(*ξ*) is the B-spline basis function. For *wi* = 1, NURBS functions are transformed into B-splines functions. The Rhino Python Editor performs the generation of the weights.

The NURBS shape function's first derivative with respect to the parameter is:

$$\frac{d}{d\xi}R\_{i,p}(\xi) = w\_i \frac{\mathcal{W}(\xi)N\_{i,p}'(\xi) - \mathcal{W}'(\xi)N\_{i,p}(\xi)}{\left(\mathcal{W}(\xi)\right)^2} \tag{7}$$

where *W*(*ξ*) = ∑*<sup>n</sup>* <sup>î</sup>=<sup>1</sup> *N*î,*p*(*ξ*) *wi* and *W* (*ξ*) = ∑*<sup>n</sup>* <sup>î</sup>=<sup>1</sup> *N* î,*p*(*ξ*)*wi*.

#### *2.3. Curve and Surface Building with NURBS*

To construct the control points for the NURBS geometry, Equation (8) is utilized:

$$(B\_i)\_j = \frac{(B\_i^w)\_j}{w\_i}, \ j = 1, \dots, k \tag{8}$$

where (*Bi*)*<sup>j</sup>* is the *j th* element of the vector *Bi* and *wi* the *i th* weight. By using Equation (6) in combination with Equation (8), the NURBS curve was constructed as:

$$C(\mathcal{J}) = \sum\_{i=1}^{n} R\_{i,p}(\mathcal{J}) B\_i \tag{9}$$

By analogy with the B-spline, NURBS surfaces are constructed from a tensor product through two-knot vector arrays. It is introduced as:

$$S(\boldsymbol{\xi}, \boldsymbol{\eta}) = \sum\_{i=1}^{n} \sum\_{j=1}^{m} \frac{N\_i^p(\boldsymbol{\xi}) \mathcal{M}\_i^q(\boldsymbol{\eta}) w\_{i,j}}{\sum\_{i=1}^{n} \sum\_{j=1}^{m} N\_i^p(\boldsymbol{\xi}) \mathcal{M}\_i^q(\boldsymbol{\eta}) w\_{i,j}} \tag{10}$$

The reformulation (10) can be done by setting up in the following form:

$$S(\mathcal{J}, \boldsymbol{\eta}) = \sum\_{i=1}^{n} \sum\_{j=1}^{m} R\_{i,j}^{p,q} B\_{i,j} \tag{11}$$

where *N<sup>p</sup> <sup>i</sup>* (*ξ*) and *<sup>M</sup><sup>p</sup> <sup>i</sup>* (*ξ*) are shape functions, respectively. Then, *p* and *q* are the order of the basis function in the two directions and m and n are the numbers of control points in the two directions.

#### *2.4. Governing Equation*

Let us briefly review the concept of discontinuity inside a continuum formulation. The improved displacement field term is introduced, taking into account the displacement field computation at the discontinuity. The concept of the partition of unity has been used for Lagrangian basis functions; the properties of this technique are also suitable for the B-spline and NURBS functions, which form the basis of the iso-geometric approach. The full displacement field can be expressed as the amount of two subfields by using this property:

$$\mathcal{U}^{\hbar}(\mathbf{x}) = \sum\_{I=1}^{N} \phi\_I(\mathbf{x}) \left( \alpha\_I + \sum\_{J=1}^{M} \beta\_{IJ} \mathcal{K}\_J(\mathbf{x}) \right) \tag{12}$$

In Equation (12), *φ<sup>I</sup>* represents the standard basis function, *KJ* is the improved basis function with m expressions, and the standard and the improved nodal degree of freedom are *α<sup>I</sup>* and *βI J*, respectively. To develop the approximation, *M* will assume the value 2, and Equation (12) can be represented as:

$$ML^\hbar(\mathbf{x}) = \sum\_{I=1}^N \Phi\_I(\mathbf{x}) \left[ a\_I + H(\mathbf{x}) a\_I + \sum\_{\beta=1}^4 F\_\beta(\mathbf{x}) b^\beta{}\_I \right] \tag{13}$$

In Equation (13), *KJ* has been decomposed into two different enrichment terms, the Heaviside function *H*(*x*) and the crack-tip function *Fβ*(*x*), in order to capture the discontinuity and the singular fields. *aI* and *b<sup>β</sup> <sup>I</sup>* introduce the improved nodal degrees of freedom.

To follow the crack, the level set methodology has been applied. This allows for representing the discontinuity as a moving interface. The signed distance, which defines the location of an arbitrary point with respect to the interface, is still the most popular:

$$\left|\varphi(\mathbf{x}) = \right| \left|\mathbf{x} - \mathbf{\hat{x}}\right| \left|\operatorname{sign}(n\_{\Gamma\_d}.(\mathbf{x} - \mathbf{\hat{x}})) \right| \tag{14}$$

where *x* is a point in a mesh element, *x*ˆ is the closet point to x on the discontinuity, and *n*Γ*<sup>d</sup>* is the interface normal.

The position of the crack in a domain is defined by the values of the level set function that is represented as follows:

$$\rho(\mathbf{x}) = \begin{array}{l|l} < 0 \; \text{if } \mathbf{x} \in \Gamma\_d^- \\ = 0 \; \text{if } \mathbf{x} \in \Gamma\_d^+ \\ > 0 \; \text{if } \mathbf{x} \in \Gamma\_d^+ \end{array} \tag{15}$$

When the body's forces are not present, the elastostatics equation, in its strongest form, is:

$$
\stackrel{\rightarrow}{div}\stackrel{\rightarrow}{\sigma} = \stackrel{\rightarrow}{0} \text{in } \Sigma \tag{16}
$$

With the appropriate set of boundary conditions: Essential boundary conditions:

$$
\overrightarrow{\boldsymbol{\mu}} = \overrightarrow{\boldsymbol{\mu}\_D} \text{ on } \Gamma\_{\boldsymbol{\mu}} \tag{17}
$$

Natural boundary conditions:

$$
\stackrel{\rightarrow}{\sigma}\stackrel{\rightarrow}{n} = \stackrel{\rightarrow}{T\_D} \text{ on } \Gamma\_t \tag{18}
$$

Crack surface: <sup>=</sup>

$$
\stackrel{\rightarrow}{\sigma} \stackrel{\rightarrow}{\dot{n}}\_{\Gamma\_d} = \stackrel{\rightarrow}{0} \text{ on } \Gamma\_d \tag{19}
$$

where *n* is the normal vector, *n*Γ*<sup>d</sup>* is the normal vector regarding the interface (see Figure 2), *σ* is the stress tensor, and *u* is the displacement vector.

**Figure 2.** Domain Σ with discontinuity Γ*d*.

#### *2.5. Variational Formulation*

As illustrated in Figure 2, a 2-D domain has been considered in this work with conventional boundary conditions: the Dirichlet Γ*<sup>u</sup>* and the Neumann Γ*<sup>t</sup>* boundary. The crack face presents further traction-free boundaries Γ*c*.

By using the virtual work theory, the weak form of the problem can be established from the equilibrium equation, and is represented as:

$$\mathcal{L}\_{c\eta} = \int\_{\Sigma} (\sigma(u) : \varepsilon(q)) d\Sigma - \int\_{\Gamma} (T\_D.q \mid) d\Gamma \tag{20}$$

The stress, the strain tensor, and the traction vector are described by *σ*, *ε*, and *TD*, respectively. *u* and *q* denote, respectively, the displacement vector and the virtual displacement.

#### **3. Extended Iso-Geometric Analysis (X-IGA)**

*3.1. X-IGA Formulation for Cracks*

X-IGA aims to evaluate fractures in an engineering component without having to re-mesh it. In the framework of the X-FEM approach, the enrichment of the shape functions (B-spline, NURBS) is ensured by a Heaviside and crack-tip function since they constitute a partition unit (PU). The former function was introduced to insert a discontinuity and the latter to treat singularities at the crack tip. The approximation of the solution is given as [54]:

$$\mathrm{L}\boldsymbol{I}^{\mathrm{H}}(\boldsymbol{\xi}) = \sum\_{\mathrm{I}\in\mathrm{N}\_{\mathrm{sland}}} \mathrm{R}\boldsymbol{\eta}(\boldsymbol{\xi})\boldsymbol{u}\boldsymbol{I} + \sum\_{\mathrm{J}\in\mathrm{N}\_{\mathrm{CrSpli}}} \mathrm{R}\boldsymbol{\eta}(\boldsymbol{\xi})\boldsymbol{H}(\boldsymbol{\xi})\boldsymbol{a}\boldsymbol{\eta} + \sum\_{\mathrm{K}\in\mathrm{N}\_{\mathrm{CrUp}}} \mathrm{R}\boldsymbol{\chi}(\boldsymbol{\xi})(\sum\_{\mathrm{a}=1}^{4} \mathrm{F}\_{\boldsymbol{\mathfrak{a}}}(\boldsymbol{\xi})\boldsymbol{b}\_{\mathrm{K}}^{\mathrm{R}}) \tag{21}$$

where *RI* represents shape functions *uI*, and *aJ* and *b<sup>α</sup> <sup>K</sup>* define the standard and the further DOFs, respectively. *Nstand* includes the standard nodes of the mesh, *NCrSplit* includes the nodes of elements that have been split by the crack faces, and *NCrTip* includes the nodes of the crack-tip elements. The parameter coordinates are represented by *ξ*. The Heaviside function, wherein the quantities on both parts of the split element are different, is denoted by *H*(*ξ*). The purpose of the crack-tip function, among others, is to increase the accuracy of the results, and it is denoted by *Fα*. These enrichment functions are characterized by the following equations [55]:

$$H(\xi) = \text{sign}\left(\boldsymbol{\varrho}(\xi)\right) = \begin{cases} -1 \text{ if } \boldsymbol{\varrho}(\xi) < 0 \\\ +1 \text{ if } \boldsymbol{\varrho}(\xi) > 0 \end{cases} \tag{22}$$

$$F\_a(\frac{r}{2}) = F\_a(r, \theta) = \left\{ r^{\frac{1}{2}} \left[ \sin \frac{\theta}{2}, \cos \frac{\theta}{2}, \sin \theta \sin \frac{\theta}{2}, \cos \frac{\theta}{2} \sin \theta \right] \right\} \tag{23}$$

The polar coordinates of the crack tip are represented by *r*, *θ*.

After substituting Equation (19) in the equilibrium equation, the final phase of the processing stage is to solve the linear algebra system:

$$K^{cur} \mathcal{U}^{cur} = F^{cur} \,, \tag{24}$$

where *Kenr* denotes an enriched stiffness matrix, *Fenr* denotes a force vector, and *Uenr* denotes an enriched displacement vector as:

$$\mathcal{U}^{err} = \begin{Bmatrix} \mathcal{U} \mathcal{U} \ b\_1 \ b\_2 \ b\_3 \ b\_4 \end{Bmatrix}^T,\tag{25}$$

where *U* is the DOF vector for the IGA normal, d is the DOF for Heaviside enrichment, and *b*1, *b*2, *b*<sup>3</sup> , and *b*<sup>4</sup> are the DOF vectors for the crack-tip enrichment functions. *Kenr* and *Fenr* are constructed as follows from the element stiffness matrix:

$$K^{\text{err}} = \begin{bmatrix} K^{\text{uu}} & K^{\text{ua}} & K^{\text{ub}\_1} & K^{\text{ub}\_2} & K^{\text{ub}\_3} & K^{\text{ub}\_4} \\ K^{\text{au}} & K^{\text{aa}} & K^{\text{ab}\_1} & K^{\text{ab}\_2} & K^{\text{ab}\_3} & K^{\text{ab}\_4} \\ K^{\text{b}\_1u} & K^{\text{b}\_1a} & K^{\text{b}\_1b\_1} & K^{\text{b}\_1b\_2} & K^{\text{b}\_1b\_3} & K^{\text{b}\_1b\_4} \\ K^{\text{b}\_2u} & K^{\text{b}\_2u} & K^{\text{b}\_2b\_1} & K^{\text{b}\_2b\_2} & K^{\text{b}\_2b\_3} & K^{\text{b}\_2b\_4} \\ K^{\text{b}\_3u} & K^{\text{b}\_3u} & K^{\text{b}\_3b\_1} & K^{\text{b}\_3b\_2} & K^{\text{b}\_3b\_3} & K^{\text{b}\_3b\_4} \\ K^{\text{b}\_4u} & K^{\text{b}\_4u} & K^{\text{b}\_4b\_1} & K^{\text{b}\_4b\_2} & K^{\text{b}\_4b\_3} & K^{\text{b}\_4b\_4} \end{bmatrix} \tag{26}$$

$$F^{\text{err}} = \left\{ F^u \; F^u \; F^{b\_1} \; F^{b\_2} \; F^{b\_3} \; F^{b\_4} \right\}^T. \tag{27}$$

The number of control points *Nstand*, *NCrSplit*, *NCrTip*, and the number of enrichment basis functions *Ne f* determine the size of each *Kenr* and *Fenr*. Each component can be written as follows:

$$\mathcal{K}\_{ij}^{rs} = \int\_{\Sigma\_r} (\mathcal{B}\_i^r)^T \mathcal{C} \mathcal{B}\_j^s \, d\Sigma \, (r, s = u, \, d, \, b\_1, b\_2, \, b\_3, \, b\_4) \tag{28}$$

$$F\_i^u = \int\_{\Gamma\_u} R\_i^T T\_D d\Gamma\_\prime \tag{29}$$

$$F\_i^a = \int\_{\Gamma\_t} R\_i^T H T\_D d\Gamma\_\prime \tag{30}$$

$$\, \_iF\_i^{b\_\mathbf{a}} = \int\_{\Gamma\_t} \mathcal{R}\_i^T f\_{b\_\mathbf{a}} \, T\_D d\Gamma \, (\mathbf{a} = 1, \ 2, \ 3, \ 4) \tag{31}$$

$$B^u\_i = \begin{bmatrix} \frac{\partial R\_i}{\partial x} & 0\\ 0 & \frac{\partial R\_i}{\partial y} \\ \frac{\partial R\_i}{\partial y} & \frac{\partial R\_i}{\partial x} \end{bmatrix} \tag{32}$$

$$B\_i^a = \begin{bmatrix} \frac{\partial R\_i}{\partial x} H & 0\\ 0 & \frac{\partial R\_i}{\partial y} H \\ \frac{\partial R\_i}{\partial y} H & \frac{\partial R\_i}{\partial x} H \end{bmatrix} \tag{33}$$

$$B\_i^{b\_R} = \begin{bmatrix} \frac{\partial R\_i f\_{b\_R}}{\partial x} & 0\\ 0 & \frac{\partial R\_i f\_{b\_R}}{\partial y} \\ \frac{\partial R\_i f\_{b\_R}}{\partial y} & \frac{\partial R\_i f\_{b\_R}}{\partial x} \end{bmatrix} (\alpha = 1, 2, 3, 4) \tag{34}$$

The Heaviside function (Equation (22)) is represented by H, while the crack-tip enrichment function (Equation(23)) is represented by *fb<sup>α</sup>* .

#### *3.2. Construction of 2-D Pipe with X-IGA Concept*

A patch with an internal interface is used to model a 2-D pipe problem (Figure 3); this interface is the result of coinciding control points in the circumferential direction at the beginning and end of the patch. To solve the problem correctly, we must confine the control values of these overlapping control points so that each pair of control values for the corresponding coincident control points is the same. The master–slave approach is a simple solution to this problem. This method is implemented in the present study by the global renumbering of DOFs in the physical space, where coincident DOFs are numbered by the same indices. The global numbering is saved in an array and given as an additional input argument.

**Figure 3.** 2-D pipe with the iso-geometric analysis approach.

#### *3.3. Enrichment Topology for Control Points*

In IGA, every basis function is linked to its corresponding control point in a specific way. The intersection of each supported domain with the crack face leads to enrichment by the Heaviside function. The domain support, which contains the crack tip, will be enriched by the singular function. According to [56], there are two ways to enrich the crack tip: with either topological or geometrical enrichment. In this work, topological enrichment has been employed. Figure 4 shows a schematic representation of this concept.

**Figure 4.** The enrichment concept used in NURBS modelling [57].

According to the Figure 4, *Cs* denotes the crack-face control points, while *CT* denotes the crack-tip-enriched control points, and *Ci* denotes the standard control point. For the purpose of selecting enriched control points, the level set method has been used. We applied the procedure that has been used by [58]. Initially, the level set values of the crack at the mesh's vertices are computed according to these level sets, and the formulation determines the elements intersected by the crack and the crack-tip element.

#### **4. Numerical Integration in the Elastic Field**

For numerical integration, the standard Gaussian quadrature method cannot be applied explicitly to the XIGA since it contains various discontinuous elements. To assess the integration rule of the crack's split and tip elements in this study, the triangular subdomain methodology is used (Figure 5). This technique has been successfully implemented by X-FEM [42], as seen in Figure 6. Each sub-triangle element has a defined number of integration points.

**Figure 5.** Partitioning the crack face (**a**) and crack tip (**b**) using the sub-triangulation technique [42].

**Figure 6.** Sub-triangles technique around the crack and the distribution of the Gauss point.

#### **5. Fracture Parameter Evaluation**

The stress intensity factor (SIF) is a crucial criterion in crack formation and growth research. As a result, when a numerical simulation is used to solve fracture mechanics issues, one of the objectives is to quantify the SIF. There are a range of ways to calculate this parameter, which is used in so much of the interaction integral method [55]:

$$\mathcal{M} = \int\_{\Gamma} \left[ \sigma\_{ij}^{act} \frac{\partial u\_i^{aux}}{\partial \mathbf{x}\_1} + \sigma\_{ij}^{aux} \frac{\partial u\_i^{act}}{\partial \mathbf{x}\_1} - \mathcal{W}^M \delta\_{1j} \right] \frac{\partial q}{\partial \mathbf{x}\_j} d\Gamma\_{\prime} \tag{35}$$

where *M* is the interaction integral with the aux and act index defining the auxiliary and the actual state, respectively. The stress and the displacement are represented by *σij* and *ui*, respectively. *W<sup>M</sup>* represents the interaction work; it can be described in the following form:

$$\mathcal{W}^M = \frac{1}{2} \left( \sigma\_{ij}^{act} \varepsilon\_{ij}^{aux} + \sigma\_{ij}^{aux} \varepsilon\_{ij}^{act} \right), \tag{36}$$

where *εij* and *q* define the strain and an arbitrary function, whose values are defined as:

$$q = \begin{cases} 1, \text{|at the crack tip} \\ 0, \text{|along the contour} \end{cases} \tag{37}$$

The stress intensity factor with these two modes (*I*, *II*) and the *J*-integral are related by:

$$J = \frac{1}{E'} \left( K\_I^2 + K\_{II}^2 \right). \tag{38}$$

Based on this correlation, the following equation has been derived:

$$M = \frac{2}{E'} \left( \mathcal{K}\_I^{\text{act}} \mathcal{K}\_I^{\text{aux}} + \mathcal{K}\_{II}^{\text{act}} \mathcal{K}\_{II}^{\text{aux}} \right) \tag{39}$$

where

$$E' = \begin{cases} \begin{array}{c} \text{E for plane stress} \\ \frac{\text{E}}{1 - \nu^2} \text{ for plane strain} \end{array} \tag{40}$$

*E* represents the Young's Modulus, *v* is the Poisson's ratio, and *M* represents the interaction integral.

The crack opening is the field of interest of this study. Therefore, the SIF in mode *I* can be obtained by choosing *Kaux <sup>I</sup>* =1 and *<sup>K</sup>aux I I* = 0. Finally, the SIF is defined as:

$$K\_I = \frac{E'}{2}M\tag{41}$$

#### **6. Process of Implementing a X-FEM Code in ABAQUS**

The idea of implementing the X-FEM technique is based on the fact that the ABAQUS software does not include the stress intensity factor computation for a 2-D crack. The two enriched functions that correspond to the Heaviside and the crack-tip functions that appeared in Equation (13) were implemented through a user-defined element (UEL). The implementation process requires three phases: pre-processing, processing, and postprocessing. In this study, to make the implementation easier, the input file (XFEM.inp) has been generated by ABAQUS, and then the interaction between the crack and the mesh was constructed to apply Equation (14). Therefore, the enriched elements and nodes have been determined. The user-defined element (UEL) in ABAQUS is used to program the processing stage (UEL.for) to compute the stresses and strains. The last stage consists of calculating the Mode-I SIF, KI, by a post-processing code; indeed, the interaction integral method explained in section five (Fracture Parameter Evaluation) has been programmed in FORTRAN. The details for the file descriptions are based on [53], who introduce the X-FEM implementation in ABAQUS. Thus, we have adopted the principle of this implementation for our problem. The ABAQUS command is used to run the main file (XFEM.inp) and the user subroutine file (UEL.for): Abaqus job = X-FEM user = UEL.for.

The resulting information of this simulation is stored in (XFEM.fil), and the ABAQUS output conventions are included in this file. Then, we used an external subroutine for computing the stress intensity factors. This subroutine is compiled through the ABAQUS command "ABAQUS make job = interaction\_integral", and then run with the command "ABAQUS interaction\_integral". The implementation process is described in Figure 7.

**Figure 7.** X-FEM implementation process in a Finite Element Code.

### **7. Process of Implementing a X-IGA Code in MATLAB**

In this section, the main steps of the XIGA implementation code to numerically model a structure with a pre-existing discontinuity (a crack) have been summarized in Figure 8. To understand the implementation, it is necessary to know two major things: the three steps of a finite element code (pre-processing, processing, and post-processing), and the identification of a crack using the level set method (LSM). The flowchart starts with the introduction of the input parameters, which contain the geometrical data and the material properties. Then, the elasticity matrix is integrated. In addition, the data required by the NURBS functions, such as the polynomial order, control points, and node vectors, are introduced to build the NURBS model. The next step is to introduce the crack data, length, and coordinates of the crack points with the level set method to determine the position of the crack and to select the enrichment points. Then, the Heaviside and crack-tip functions are imposed on the nodes according to the technique that precedes this step. Then, the boundary conditions, nodal force vector, and stiffness matrix are calculated. The stress, strain, and displacement values are the output of this process. These data can be put to use in the interaction integral to compute the stress intensity factor.

**Figure 8.** Flowchart of the implementation.

#### **8. Numerical Results and Discussions**

The iso-geometric analysis is extended and used in this section to analyze the cracked model under mechanical loading; a specimen from the industrial field is chosen. Pressure pipes are commonly used and the performance of these systems is still under investigation. In a 2-D linear static analysis case, an isotropic and homogeneous pipe has been studied. The plane stress is thought to be the stress distribution state. For this analysis, The P264GH material is used and its mechanical and chemical characteristics are described in Tables 1 and 2, respectively.


**Table 1.** Mechanical characteristics for the P264GH [59].

**Table 2.** Chemical characteristic of P264GH—% by weight [60].


To examine the reliability, accuracy, and efficiency of this approach, an external axial crack was studied. The results of the analysis were compared with the results of standard ABAQUS software using the classical finite element method (CFEM), which is based on the integral contour and the X-FEM method using a subroutine UEL code.

In each of these parametric directions *ξ* and *η*, the degree of the NURBS polynomial is two (*p* = *q* = 2). In the first direction *ξ*, the knot vector is open and with interior duplication, and in the second direction *η*, the knot vector is open and without internal duplication.

The integration is done along each Gauss point direction (*p* + 1) × (*q* + 1), and as the sub-triangles approach was used in this case, 13 Gauss quadrature points were imposed for each sub-triangle. For each numerical method, different mesh sizes were examined. It is important to be aware that the crack was represented as a straight segment.

In order to figure out the fracture parameter, i.e., to estimate the stress state near a crack tip, the stress intensity factor (SIF), KI, and the J-Integral are extracted. For the three techniques used in this work, the interaction integral method was implemented. The X-IGA technique was implemented in MATLAB code. The ABAQUS software was used for the CFEM and X-FEM to extract the KI, but for the X-FEM, it is important to mention that the software does not support the computation of this parameter in the 2-D domain, which led to the use of a subroutine UEL.

### *8.1. Two-Dimensional Pipe with an Axial Crack under Uniform Pressure*

In the present study, an isotropic and homogenous 2-D pipe including a 2 mm straight edge crack with uniform pressure distribution and an outer and inner radius of Ro=20 mm and Ri=10 mm, respectively, is examined. Three models were used in this study, which included 320, 480, and 640 element numbers for the X-IGA method, and 1297, 3844, and 85,942 element numbers for the CFEM, with a step size of 1, 0.5, and 0.1, respectively, and 470, 1265, and 2747 element numbers for X-FEM with a step size of 1, 0.45, and 0.25 around the crack, respectively. The NURBS geometry is represented by using several patches (subdomains) with an internal interface. In the circumferential direction, the control points coincide with each other on the patch, and this implies the creation of interface terms. The NURBS geometry is illustrated in Figure 9a.

**Figure 9.** Geometry construction: (**a**) the NURBS model, (**b**) the CFEM model, (**c**) X-FEM model.

The FEM model has been represented by a finite element mesh on ABAQUS software. The most important thing is that the modeling of the crack requires a specific treatment, i.e., it requires building a particular mesh around the crack. A ring of a triangle, as represented in Figure 9b, is formed at the crack tip, along with concentric layers of structured quads [61].

For the X-FEM method (Figure 9c), the meshing was done in a simple way since the crack was modeled independently of the meshing. It should be noted that, for this technique, a half tube was treated because the implementation is heavy when using UEL subroutines, so to overcome this problem, a half tube was analyzed with the symmetry boundary conditions.

To perform the numerical simulation, we present Table 3, which summarizes the crack length, crack position, and pressures applied.

**Table 3.** Data details for the numerical simulation.


Figure 10 illustrates the distribution of the Von Mises stress by the three techniques for *<sup>a</sup> <sup>t</sup>* = 0.3, and a pressure of 2.5 MPa [52] has been applied. The accuracy of the stress and strain distribution in a geometry takes a very important place in fracture mechanics, especially when it concerns cracking problems.

**Figure 10.** The distribution of Von Mises stress: (**a**) with FEM analysis, (**b**) with X-FEM analysis, (**c**) with the improved technique (X-IGA).

The zone of interest is the crack tip where the degree of damage of the defect has been known. Obtaining a regularity of stresses in this region is a priority for numerical methods since the calculation of fracture parameters, such as the stress intensity factor, as well as the angle of deviation of the crack propagation, is based on the value of stresses at the crack tip.

By comparison of the three results of Figure 10, there is a similarity in the results of the numerical simulation. However, some errors can appear when using these numerical techniques. The errors' origins are detailed as follows: with the CFEM analysis, the field of strains and stresses become singular at the crack tip, even with the integration of the singularity in the model to improve the accuracy of the results. With the implementation of the X-FEM method via user subroutine UEL in ABAQUS, a mesh sensitivity affects the simulation results. For the present study, there is no singularity at the crack tip and there is a minor sensitivity of mesh. The next section of this work shows a calculation of different meshes for all the different methods. The X-IGA method implemented in MATLAB can be an alternative to these numerical methods. The results presented in Figure 10 can support this conclusion because 4112 elements were used for the full tube with the CFEM method, 2747 elements for the half tube with the X-FEM method, and only 384 elements for the XIGA method. For this reason, the evaluation of the stresses around an existing crack on a pipe can be made by the present study with large elements and with a weak error, which appears in the solution discretization.

#### *8.2. Evaluation of the Fracture Parameter*

In order to present the accuracy of the X-IGA technique, as well as the regularity of the stress distribution around the crack tip, the SIFs were extracted and the value of the mode I (*KI*) was calculated in this study, since the degree of damage that corresponds to the opening of the crack is more severe than with the other modes [62]. To ensure that this study is inscribed in linear elasticity, a pressure of 2.5 MPa [52] was applied. Three models are evaluated for different mesh sizes. Therefore, the result of CFEM, X-FEM, and X-IGA are compared. Moreover, the interaction integral approach, which is known as M-integral, was used for calculating the SIFs.

For this comparison, the depth of the crack was varied from the thickness of the model a/t; the thickness is *t* = 10 mm and *a* = 2, 3, ... , 8 mm. Figures 11 and 12 illustrate the results of the computation. The comparison between the results obtained by the X-IGA and CFEM methods shows a good similarity. It is observed that fine meshes must be used to discretize the areas around the crack tips for the CFEM, whereas the meshes are coarser for the XIGA. It is clear to observe that only 640 elements are required by the present study to obtain the result of 1297 elements with the CFEM method. Results obtained by implementing X-IGA in MATLAB are in good agreement with the CFEM method implemented in ABAQUS. The mentioned technique did not require much effort in dealing with the mesh in the crack tip, whereas using ABAQUS required a specific meshing technique, due to the singularity at the crack tip, and this can affect the numerical results.

As for the comparison between the X-IGA and X-FEM, the results show a good similarity. It is also observed that, with X-FEM, a fine mesh must be used without a particular treatment around the crack, since the crack is modelled independently of the mesh. It is clear to observe that the present study requires just 640 elements for obtaining the result of 2747 elements using X-FEM. The same enrichment functions for overcoming the problem of singularity are used for both methods and only the shape functions are modified. Therefore, there is a similarity of results between the two numerical methods. When the depth of the crack approaches the inside of the pipe, the SIFs values become maximized. It is clear to see that whatever technique is employed to simulate the crack, the SIFs gradually increase with the crack length.

**Figure 11.** The comparison between the CFEM and X-IGA method for a 2-D cracked pipe.

**Figure 12.** The comparison between the X-FEM and X-IGA method for a 2-D cracked pipe.

By these results, implementing X-IGA in MATLAB shows the possibility to calculate the fracture mechanics parameters for a cracked pipe under uniform pressure with a large mesh size, compared to the other approaches (Figure 12). Therefore, it can be an alternative way to evaluate the damage of a cracked pipe.

In addition to the SIF, the J-integral has specific importance when it comes to the numerical stress analysis of cracks. To give more physical meaning to the analysis, and to validate the strategy used in the application of X-IGA to address the cracking pipe problem, we evaluated the J-integral with various crack lengths, as shown in Figure 13. Here, we used 640 elements for the present study, and 2747 and 85,642 elements for X-FEM and CFEM, respectively. The J-integral value increased gradually with the crack length and the results obtained are similar to other numerical results.

**Figure 13.** J-integral value obtained by the present study, X-FEM, and CFEM.

This problem is also analyzed by [51]; they used a UEL subroutine in the ABAQUS software to evaluate the stress intensity factor. With the same element number and for *a* = 5 mm and *p* = 2.5 MPa, we evaluated the performance of the present study and the results of both implementations are compared with Folias solutions [63], as illustrated in Table 4.


**Table 4.** Comparative study of the present investigation and implementation using Fortran for *a* = 5 mm.

Table 3 proves the significance of the present study implemented in MATLAB. It is interesting to note that, in the present study, the error for a model with 470 elements is 3.645%, while for the study implemented in Fortran the error is 8.306%. Additionally, with 767 elements, the error for the present study is 6.01%, while for the study implemented in Fortran the error is 9.454%. It is observed that the X-IGA implementation strategy that was followed in the present study had a higher accuracy than the strategy followed by [51].

In addition, the accuracy of the present study can be realized by computing the error of SIFs (%) with various crack lengths, as illustrated in Figure 14. It is noted that the maximum error of the present study is 12.25%, while the study implemented in Fortran [51] is 27.83%.

**Figure 14.** Error of SIFs with various crack length for the present study and the reference [51].

#### *8.3. Effect of Pressure on the Fracture Parameter Calculation*

Finally, in order to check the efficiency of the present study for the calculation of the fracture parameter in the 2-D pipe domain, the inner pressure was varied with the initial crack length *a* = 4 mm, so that the pressure did not exceed 4.5 MPa [64], keeping the analysis in the elastic domain. The results of the SIFs and the J-Integral of CFEM, X-FEM, and X-IGA are presented in Figure 15. It is observed that for the three-analyses technique, the stress intensity factor increases with the increase of inner pressure. It is also observed that the results for both the X-FEM and X-IGA methods became similar each time the pressure was increased; this is due to the use of enrichment functions at the crack tip. The singularity at the crack tip and the mesh dependency for the CFEM method makes their results inaccurate.

**Figure 15.** (**a**) SIFs and (**b**) J-Integral values for different pressures by the CFEM, X-FEM, and EX-IGA methods.

From these results, the X-IGA method can be used in the fracture parameters computation on cracked cylindrical structures under uniform pressure, the X-IGA method has been well validated for the calculation of cracked plates [48], and this study is an extension of the research work that has already been done in this field. The comparison between the most-known methods in the field of numerical computation will allow us to justify the role of this new technique: it can replace the existing methods in the current calculation codes in the future, and the validation of the efficiency of this technique by the different research works will convince the users in the industrial field. This is due to several reasons; the most important reason, for the industrial sector, is the cost of calculation. Instead of designing the model and approximating it by several meshing processes, it is enough to use the geometry directly by the same shape functions, named B-spline or NURBS; after that, the model will be reproduced precisely. It has been observed that X-IGA just needs a small element for a crack in a pipe, compared to other techniques, and that with these elements, the same results have been obtained, so if there is a complex geometry, the calculation procedure will be very fast compared to the current method.

#### **9. Conclusions**

The aim of this investigation was to implement the X-IGA technique into a MATLAB code for modeling cracked pipelines. The main theoretical approach was the NURBS principle, which provided a higher-order continuity in numerical modeling.

An external crack in a two-dimensional pipe subjected to a uniform pressure has been studied. The accuracy of this technique has been examined by deriving the stress intensity factors and the J-integral. For several mesh sizes and for different inner pressures, SIFs and the J-integral were extracted by X-IGA analysis using MATLAB code and its accuracy was validated with the enrichment technique (X-FEM) using FORTRAN language and the conventional finite element method (X-FEM) using ABAQUS software. It has been shown that, when using the X-IGA analysis:


Other problems can be addressed by the provided technique, such as crack-growth problems and dynamic fracture analyses in pipelines, which are planned for future research. **Author Contributions:** All the authors conceived the framework and structured the whole manuscript, checked the results, and completed the revision of the paper. The authors have equally contributed to the elaboration of this manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors thank Catalin Iulian Pruncu from the University of Strathclyde, UK, for proofreading assistance and overall guidance, which made drafting this work possible.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

