Numerical Simulation of cracked orthotropic materials using extended isogeometric analysis

In the present study, extended isogeometric analysis (XIGA) is used to analyse cracks in orthotropic media. NURBS and T-splines geometric technologies are used to define the geometry and the solution. Knot insertion and order elevation are used in NURBS models, while a new local refinement algorithm is applied to T-spline models. In XIGA, the basic idea of the extended finite element method (X-FEM) is used along with isogeometric analysis for modelling discontinuities by including enrichment functions. Special orthotropic crack tip enrichments are used to reproduce the singular fields near a crack tip, and fracture properties of the models are defined by the mixed mode stress intensity factors (SIFs), which are obtained by means of the interaction integral (M-integral). Results of the proposed method are compared with other available results.


Introduction
Composite materials are not immune to manufacture defects such as voids, inclusions and cracks. The crack problems greatly influence in the macroscopic response of composites materials, and this lead to the need for better understanding and analyzing the behavior of these materials under applied loads especially in critical conditions.
The most widely geometric function used in CAD is the NURBS, due to its simplicity and its efficiency in geometrical representations, especially all conical sections, which can be represented exactly, such as, circles, cylinders, spheres and ellipsoids. NURBS basis functions have very useful mathematical properties, like partition of unity, compact support, and high degree of continuity. Also, NURBS are characterized by the multitude of refinement methods (h-refinement, p-refinement, krefinement). In addition, there are many effective algorithms that can be programmed numerically to generate NURBS objects. As in CAD, NURBS was also used for IGA in most researches, because their mathematical properties are always compatible with analysis conditions regardless of the used geometric shape. Despite all the advantages mentioned there are some disadvantages observed during the analysis, which cannot be avoided even by using multiple patches, where NURBS generates a complicated mesh, which lead to produce superfluous control points.
Sederberg et al. [62] proposed T-splines as more generalized tools than NURBS in order to handle their disadvantages. T-splines are characterized by their robustness in the generation of adjacent patches, and also by the local refinement property of the mesh even for a single patch. In addition, T-splines provide several algorithms that can be used to refine geometries locally. Recently, T-splines were introduced as additional modules in some CAD software, like Maya and Rhino. As part of the IGA, Tsplines are considered as applicable analytical tools that have proven to be efficient in several studies. From the analysis point of view, T-splines allow to use a minimum number of degrees of freedom compared to the number that can be used in the NURBS, and this leads to simplify the computation. In some domains, T-splines are considered as a better tool for IGA, like in the fracture mechanics, contact mechanics and mechanics of bi-materials, generally in any domain where the property of local refinement is necessary and the geometry is complex. However T-spline bases are not always valid to be used in analysis for different geometric configurations, because the linear independent and the partition unity properties are not always ensured. Li et al. [63] introduced analysis-suitable T-spline, where for any choice of knot vectors the blending functions are linearly independent. Like NURBS bases, analysis-suitable T-spline bases have the properties of the analysis basis functions. Moreover, they provide an efficient algorithm, which allows making highly localized refinement [64].
In this paper, cracks in orthotropic media are analysed using XIGA, whereas NURBS and T-splines are used for modelling the crack and construct the geometry. Orthotropic crack tip enrichment functions [65][66][67] are implemented to accurately calculate the displacement and stress fields, and mixed mode stress intensity factors (SIFs) are numerically evaluated using the interaction integral (M-integral) [68] to determine fracture properties of the domain. An edge and central crack are investigated by considering the effect of fibre orientation.

Fracture mechanics for 2D orthotropic materials
Consider a cracked elastic homogeneous orthotropic body Ω subjected to traction forces ft applied at Γt and displacement conditions applied at Γu in the absence of body forces, as shown in Figure 1. The partial differential equation of the stress functions of this problem can be obtained using equilibrium and compatibility conditions [ where (x,y) is local Cartesian coordinates, φ is the stress function and cij is the components of the compliance matrix, which is computed in terms of Young moduli (E1, E2), Poisson ratios (ν12,ν13 and ν23) and shear modulus (μ12, μ13 and μ23) as follows: The general characteristic equation of the partial differential equation is: The four roots of this characteristic equation are complex and conjugated two-two (s1, ̅ 1 and s2, ̅ 2). These roots were used by [1] to derive the displacement and the stress fields in the vicinity of the crack tip for an infinite domain. The stress components for pure mode I are: and for pure mode-II: The displacements for pure mode-I are where Re represents the real part of a complex number, KI and KII are the stress intensity factors for modes I and mode II, respectively, and the constant values di and ei are computed as: 2 11 12 16 1, 2 The relation between energy release rate G and stress intensity factors was expressed by Sih et al. [1] for homogeneous composites as follows: 12

Non uniform rational B-splines (NURBS)
where n is the number of basis functions, w is a set of positive weights, ξ is a parametric coordinate of the first direction and Ni are the B-spline basis functions corresponding to the knot vector Ξ={ξ1,ξ2,......,ξn+p+l}, they are given by: and  (14) For the two parametric directions ξ and η, bivariate NURBS basis functions are given by:  (15) where Mj are B-spline basis functions of order q corresponding to the second knot vector H={η1,η2,......,ηm+q+1} of which m is the number of basis functions.

T-splines
In order to construct a T-spline surface, the T-mesh, which is a mesh of rectangular elements defined by the lines corresponding to knot values, as shown in Figure 2, must be defined according to the basis function orders p and q to extract the local knot vectors Ξα and Hα, and compute the blending functions Rα for each anchor. Then, in a similar way as for NURBS, we can define the T-spline surface by: where the blending functions given by: (17) and P is the control points, α is the anchor index, k is the number of the anchors in T-mesh and B is the bivariate local basis function, which is given by: For more than one inserted knot, the application of these equations is repeated until all coefficients are derived. In refined T-spline space, the relation between the original basis functions B 1 and the generated basis functions B 2 can be expressed in linear system form using the refinement operator M 1 2  N MN (22) where N1 and N2 are the column vectors of blending functions of the original and the refined T-spline spaces, respectively.

Extended isogeometric analysis (XIGA)
In XIGA, the displacement approximation for cracks at a particular point ζ=(ξ,η), can be written in generalized form by extending the IGA approximation: where Q is the T-spline basis function (R), Eq. 17 or NURBS basis function (R) Eq. 15. H is the Heaviside function used for modelling the crack face, it takes the value 1 above the crack and -1 below the crack, F are the crack-tip enrichment functions derived from the analytical solution of the displacement field around the crack tip, ui, aj and bk are the displacement vectors correspond to ns, ncf and nct control points, respectively. The crack-tip enrichment functions are obtained from the analytical solution of the displacement field in the vicinity of the crack-tip as [66]  

Numerical simulations
In this section, a circular orthotropic plate with a central crack and a rectangular orthotropic plate with an edge crack are simulated in the plane stress state. The corresponding boundary conditions are shown in Figure 3. First, in the circular plate (E1=0.1 GPa, E2=1 GPa, μ12=0.5 GPa, v12=0.03) different inclined cracks are examined, where the orthotropic axes (1,2) coincide with the Cartesian axes (x,y). The problem is analysed using 793 control points and 688 elements (Figure 4). Variation of mode I and II SIF is presented in Figure 6.    In Figure 6 and Table 1 the obtained SIFs are in good agreement with the other methods. T-splines provide a minimal number of nodes due to their ability in describing exactly such geometries. Also, the results of the edge crack problem (Figure 7)  of freedom compared to XFEM and NURBS based XIGA. The mode I SIF changes their trend in 45° and the mode II SIF changes their trend in 30°, so in orthotropic materials the stress intensity factors are not only depend on the crack angle, but they also depend on the direction of elasticity axes. For more details about the results that obtained by other numerical methods for this problem, the reader may refer to [23,69].

Conclusion
Static fracture behaviour for a crack in orthotropic materials is analysed numerically in this study using the Extended Isogeometric Analysis (XIGA). NURBS and T-splines bases are used to define the exact geometry and approximate the solution throughout the analysis. Orthotropic crack-tip enrichment functions, which can be applied to all types of orthotropic materials, are implemented. SIFs for both mode I and II are evaluated using the M-integral to determine fracture properties of the domain. The quality of the obtained results regarding their good agreement with the results of literature demonstrates the accuracy of the present approach and the robustness of the developed code. The inclination angle of the orthotropic axes generates and affects the evolution of mixed mode SIF even if the crack is only subjected to mode I loading, therefore the mode II must also be evaluated.