Anisotropic elastic applications in composite materials using the isogeometric boundary element method

Abstract This paper describes an isogeometric analysis of the boundary element method, called IGABEM, applied to anisotropic 2D plane elastic problems. The Lekhnitskii's anisotropic fundamental solution is used and the singular integral terms is regularized. For the weak singularity kernel, the Telles transform is used. On the other hand, in the strong singularity term, the Singularity Subtraction Technique - SST is used. The shape functions used are the Non-Uniform Rational B-Spline - NURBS. Thus, the same mathematical representation of Computer Aided Design - CAD is used in the implemented computational code. This avoids the generation of meshes and provides an exact representation of most complex geometries. As a reflection of this isoparametric concept, errors from geometric interpolation are excluded, improving numerical results. Furthermore, to increase the numerical efficiency of the code, the NURBS are decomposed into Bézier curves. To evaluate the accuracy of the formulation, complex problems using high-order isogeometric boundary elements are analyzed. Their results are compared to analytical solutions showing good agreement. Graphical Abstract


INTRODUCTION
As it is well known, Computer-Aided Design (CAD) programs use non-uniform rational B-Splines (NURBS) to model boundaries of solids.NURBS can exactly represent circles, ellipses, spheres, cylinders, and many other engineering geometries, which is not possible to obtain with the use of other functions such as, for example, polynomials, even those of higher order.Among the difficulties in working with high-order polynomials, we can highlight the globality: when changing any point of the curve, the entire curve will be affected.Therefore, when working with polynomials, it is usually necessary to divide a curve into elements.This behavior is not observed when using NURBS, because each control point only affects a portion of the whole curve (Cottrell et al., 2009).The NURBS geometric representation is typically smooth, whereas the representation for polynomials is typically continuous but not smooth.The ability to efficiently use a smooth basis in analysis has computational advantages over polynomials.
In this paper, NURBS are used not only to model geometry but also to approximate displacements and tractions throughout the boundary.NURBS curves are given through the Bézier extraction operator, which decomposes a set of NURBS basis functions for the Bernstein polynomial (Borden et al., 2011).This will allow the generation of continuous Bézier elements of class 0 C , according to Beer and Duenser (2019), giving local representation of the basis functions, in addition to establishing a well-defined element from where collocation points are defined.The proposed Bézier elements provide a structure of elements for isogeometric analysis (Hughes et al., 2005), which is incorporated into the BEM, giving the idea of the isogeometric analysis boundary element method (IGABEM).Thus, the isoparametric concept is preserved, once that geometry, displacements, and tractions are interpolated by the same functions.Applications are presented for twodimensional linear anisotropic elasticity problems, whose fundamental solution is written using through the formalism of Lekhnitskii (1968), in which the constitutive equation in terms of the flexibility tensor for the Airy stress is used.These fundamental solutions are present in several works and in different areas dealing with anisotropic elasticity, some of which are fracture mechanics (Rajapakse and Xu, 2001), (Albuquerque et al., 2004); bending problems (Pastorino et al., 2021); stress concentration factor (Özaslan et al., 2019).However, when using the fundamental solutions of Lekhnitskii, depending on the position of the collocation points of the discretized geometric region, singularities arise.These singular integrals are represented by the displacement (weak singularity) and traction (strong singularity) kernels.For weak singularity, the Telles transformation (Telles, 1987) was used, which develops a new transformation for the traditional Gaussian quadrature.And for strong singularity, the Singularity Subtraction Technique (SST) is used, as presented by Guiggiani (1998).The SST allows integration over high-order curved boundary elements, and this idea fits perfectly with the NURBS curves.In addition, according to Cordeiro and Leonel (2020), the SST is a methodology that can be applied to cores in a generic way since it does not impose any formal restrictions on the type of core to be integrated.
Therefore, this article presents an Isogeometric Analysis of the Boundary Element Method, with the anisotropic fundamental solution integrated by the Singularity Subtraction Technique -SST.As far as the authors know, the IGABEM for anisotropic materials has not yet been presented in the literature.Furthermore, the NURBS are decomposed into Bézier curves without the loss of continuity properties, which helps to reduce the cost of the developed computational algorithm.As a result, this formulation has superior accuracy when compared to codes that use the conventional MEC.Different numerical examples of simple and complex structures are presented to demonstrate the efficiency of the proposed methodology.

FUNDAMENTALS OF ANISOTROPIC ELASTICITY
If considered an infinitesimal element within the domain Ω the equilibrium of forces is expressed in the form The traction at any point on the domain boundary Ω is given by: where j n is the normal vector of the boundary.According to Hwu (2010), it is necessary to have the compatibility conditions for the deformation fields to guarantee the existence of a continuous displacement field of single value for a continuous body.Thus, the compatibility equation for infinitesimal deformation, given the symmetry of the deformations, for the two-dimensional case is satisfied if its derivatives obey the following relationship: Constitutive relationships for linear elastic material in an infinitesimal element, are related by generalized Hooke's law, as follows: The terms ij σ and ij  are the second order stress and strain tensors, respectively.Each tensor is represented by nine index components , 1, 2,3 i j = . In Equation (3), the linearity coefficient is the fourth order tensor with 81 elements, called a tensor of elastic constants, which satisfies Green's symmetry conditions ijkl For plane stress or plane strain problems, Equation ( 4) can be written as: where,   are given in terms of material compliance coefficients ij a (Lekhnitskii, 1968) as:

Formulation of the boundary integral equation
The boundary integral equation for anisotopic elastic problems can be derived from Betti's reciprocal theorem and Somigliana Identity, as can be seen in Gaul et al. (2013).Therefore, for a bidimensional anisotropic elastic problem of boundary Γ and considering a null body forces, a boundary integral equation can be written in the form: The right hand side integral represents an integral in the sense of Cauchy main value, and ij c is the free coefficient term.In the case of this work, as the collocation points are always in a smooth part of the boundary, we have were ij δ is the Kronecker delta.
Fundamental solutions for plane anisotropic elasticity used in Equation ( 8) are based on the Airy stress function as seen in Lekhnitskii (1968).For 2D problems in plane strain state, the fundamental solution of displacement ij U and the fundamental solution of traction ij T , which represent the displacement and the traction at the field point z , in the direction j , when the concentrated force is acting on the direction i at the source point 0 z , are given through a mapping in the complex plane, as: Based on Equation ( 9), fundamental solutions of displacement and traction are written as follows: Latin American Journal of Solids and Structures, 2023, 20(1), e471 4/23 where

[ ]
Re stands for the real part of the terms.Elements   , ji g , and jk A can be seen in Sollero & Aliabadi (1993).

REPRESENTATION OF NURBS THROUGH BÉZIER EXTRACTION
Some geometric shapes are not possible to draw exactly using conventional polynomials.To solve this problem, rational functions, known as NURBS, are used as basis functions.With this basis is possible to achieve all forms of geometries, from the simplest to the most complex.In NURBS basis functions, in addition to the vector of nodes and the degree of the curve, it is necessary to introduce a weight vector.Thus, NURBS basis are defined through a weighted average between the B-Spline basis functions and may have a different weight for each control point.
From the algebraic point of view, a NURBS curve of grade p is defined as: In Equation ( 12), the term i P are the control points, and the set , ( ) R t is the rational basis written as, ( , i p t φ is the basis B-Spline of degree p for a vector of non-uniform knot written as { } 0, 0..., 0, 1, ...,1 u = , ( ) W t , is the weight function.Equation (13) can also be written in the matrix form as: Important characteristics for Equation ( 14) can be consulted in Hughes et al. (2005).

NURBS in the form of Bézier curves
The Bézier extraction operator for NURBS breaks down a set of NURBS base functions into Bernstein's polynomial.This allows the generation of continuous 0 C Bézier elements according to Beer et al. (2020), giving local representation of basis functions in addition to establishing a well-defined element where collocation points are located.The work of Borden et al. (2011) presented details of this methodology.
To build NURBS by Bézier curves, two steps are required: 1) Bézier decomposition, where all inner nodes are repeated by insertion of nodes up to multiplicity equal to p and 2) Bézier extraction operator.

Insertion of knots
The refinement process is performed inserting new nodes . New nodes are added at the midpoint between the old ones.In order to avoid Latin American Journal of Solids and Structures, 2023, 20(1), e471 5/23 change in the continuity of the curve, 2 m + new control points are inserted.Furthermore, there is no change in geometric or parametric properties of the curve.Therefore, for each inserted node in the existing node vector, a new control point is added.This is done through:

(
) According to Borden et al. (2011), the continuity of basis functions is reduced by one with each node insertion.However, by defining the new control point, the continuity of the curve is maintained.

Bézier decomposition
To illustrate the Bézier decomposition, consider Figure 1.In Figure 1(a) and 1(b) the NURBS curves are shown through the gray dashed line, constructed from the control points and their basis functions.The NURBS of Figure 1(a) was built using the basis functions obtained via conventional NURBS, as shown in Figure 1(c), while the NURBS of Figure 1(b) was constructed using the basis functions obtained via Bézier decomposition by the Bernstein polynomial, as shown in Figure 1(d).Figure 1 shows a cubic curve, whose associated node vector is { } 0,0,0,0,1, 2,3, 4, 4, 4, 4 = U .The process of decomposing the Bézier curve is carried out by inserting repeated nodes

{ }
1,1, 2, 2,3,3 = U in all internal nodes, through Equation ( 15), until there is a multiplicity equal to the degree of the curve.According to Borden et al. (2011), it is important to highlight that those internal knots should have multiplicity to form Bézier elements.However, p multiplicity is sufficient to represent Bernstein's polynomials, which in this context are also called Bézier basis functions.With each node that is inserted, the continuity of the basis functions is reduced by one degree at the knot location.Note that knots internal to element, 1, 2 and 3 ξ = , Figure 1(c), appear only once in the knot vector and, therefore, the maximum possible continuity level is represented.

Bézier extraction operator
The Bézier extraction operator C is constructed considering the n control points and the m knots necessary to perform the decomposition of the vector and j α according to Equation (15).Equation ( 15) can be rewritten in matrix form to represent the sequence of control points created by inserting knots, as: ( ) By repeating the form of Equation ( 16) m times, the final set of control points Latin American Journal of Solids and Structures, 2023, 20(1), e471  6/23 In this way, Bézier elements are obtained through the final form of the decomposition, and the relationship between the new Bézier control points and the initial NURBS control points is written as: Since, inserting a knot does not cause geometric or parametric changes in the curve, then, given the set of Bernstein basis functions defined by the knot vector ( ) After matrix multiplication, Equation ( 19) is written as, where C is called the Bézier extraction operator and depends only on the knot vector.Note that B is represented as class 0 C , and C maps the Bernstein polynomial basis piecewise into NURBS, performing the linear combination of the Bézier basis functions.Thus, when writing CB , this term reverts to the continuity of the original NURBS.
Inserting Equation ( 20) into Equation ( 14), the NURBS basis functions using the Bézier extraction operator, becomes: In this format, Equation ( 21) is ready to be used as a shape function in the discretization of IGABEM.

ISOGEOMETRIC BOUNDARY ELEMENTS
In IGABEM, the definition of curves is based on patch which are the knots in the 1 1 interval are the knots span, ξ being the Gaussian element,  (Simpson et al. 2013).From this connectivity between elements, it is possible to establish isogeometric approximations for the geometric element, and the geometry of the problem is described using the basis functions as follows: Displacements and tractions are written as: u and e c t represent the displacements and traction at collocation points that may not belong to the domain (Simpson et al. 2013).
Using Equations ( 23) and ( 24), in Equation ( 8), we have: where le j u and le j t represent the components of vectors e c u and e c t , and J is the Jacobian of the transformation.In conventional BEM, the Kronecker delta property of the basis functions guarantees that the basis functions are interpolating.In the IGABEM, on the other hand, it cannot guarantee that this is true, and the displacement must be given as , where ξ ′ is the local coordinate of the collocation point in the element e .
The Bézier extraction operator promotes the construction of elements very similar to the one used in the conventional BEM.Each patch of the parametric space corresponds to an element independent of the adjacent elements, which, in turn, corresponds to the parametric subintervals adjacent to the initial one.Thus, in this paper, Equation ( 25) is applied to each control point.Finally, a non-sparse and non-symmetric matrix algebraic system is obtained as follows: where H and G are influence matrices.Note that the control points, except for rectilinear geometries, do not belong to the NURBS curve.And so, the boundary conditions cannot be applied directly to Equation (26).To solve this problem, an E transformation matrix for NURBS is generated according to Cabral et al., (1990).This matrix uses basis functions to relate the values at the control points to the values at the collocation points, in the form: where u and t are the vectors containing the displacements and tractions at the collocation points, c u and c t are vectors containing values at the control points.Now, Equations ( 27) and ( 28) can be applied to Equation ( 26): Reorganizing the known and unknown terms of Equation ( 29) we arrive at a linear equation system of the type Where x is the solution vector of unknown tractions and displacements.

Strategy for the generation of collocation points
As commented, control points may not be necessarily on the curve, but on the control polygon.This is not the best choice for collocation points in the BEM.In Li & Qian (2011), different forms of distribution for the collocation points were tested.Based on their results, the Greville coordinates is the best choice, with stable and accurate results for the isogeometric method.Thus, the Greville abscissa used here is defined as the mean value of knots: When the curve is smooth, Greville coordinates are perfect to be use as collocation points.Otherwise, when there are corners, in order to avoid collocation points at corners, it is proposed to change the positions of the first and last collocation point: Latin American Journal of Solids and Structures, 2023, 20(1), e471 8/23 ( ) where s is a coefficient that defines how much the collocation point moves.As shown by Wang & Benson (2015), an optimal value is 0.5 s = .

Singular integrals
In the conventional BEM, only one basis function associated with an element is non-zero at a given collocation point.This is not true for IGABEM, since more than one basis functions can be different from zero at the collocation point.Therefore, there is only one singular integral in each element.While in NURBS, the basis functions do not have this characteristic, as several of them can be non-null in the collocation points.For each collocation point there are k singular integrals.This characteristic related to the analysis of singularities in IGABEM is a crucial point to solve the BIE.
The fundamental displacement solution containing the weak singularity was treated in this work with the Telles transform.Details about this method will not be given here, once it can be found in Telles (1987).On the other hand, for traction equation, the Subtraction Singularity Technique (SST) seen in Guiggiani (1998) is used.It is based on the expansion of the kernel in series for the treatment of singularities.In general, this technique is formulated based on the substitution of the singular integral, creating two new terms: a new regular kernel and a new singular integral, in such a way that the latter can be integrated analytically.
The first step of this method is to consider the expansion of k z , which is the complex mapping shown in Equation ( 9), in Taylor series around the singular term 0 ξ , which becomes: Expanding the derivative in Equation ( 27), it becomes: Now, substituting Equation (34) into Equation (33), results in: It is worth remembering that the unit vector n which is normal to the boundary Γ , and its components in 1 x and in 2 x directions are written as 1 x n and 2 x n in the form: From Equation (36) one arrives at: So when replacing Equation (37) in Equation ( 35), it is possible to find the following equation: Latin American Journal of Solids and Structures, 2023, 20(1), e471 9/23 where [ ] Through Equation ( 38), we arrive at the final expression for the expansion of k z around 0 ξ , writing as: ) The expansion written in this form as presented by Cordeiro & Leonel (2020), when replaced in the strongly singular integral, can isolate singularities of kernels of integrals.Therefore, for the regularization process of the strongly singular term, expanded auxiliary kernels for the boundary will be added and subtracted.Then, the strongly singular term is written as: Note that, in the integral on the left hand side of Equation ( 41 where When analyzing Equation (41) after replacing Equations ( 42) and ( 43), the strong singularity is noted.To regularize the integral between square brackets of Equation ( 43), Marczak & Creus (2002) used the first-order term of the expansion of the basis function.On the other hand, Cordeiro & Leonel (2020) considered only the constant term of the As it has been shown that only one term is necessary for an adequate representation, then only the constant expansion term is used in this work, and with that, the integral between square brackets of Equation ( 41) is solved numerically by the Gauss Legendre quadrature.On the other hand, the other integral on the right side is still singular and must be solved analytically in the sense of the Cauchy Principal Value -CPV through a limit analysis.So * 0

( , )
T K ξ ξ is written as: From Equation ( 40), it is possible to find the relationship: Thus, it is possible to remove the term that is replaced in Equation ( 44), which is now written as: Therefore, with Equation ( 46) it is possible to calculate the singular integral in a regular way using the limit [ ] CPV 1 / ( ) ln( 1) ln( 1) It is necessary to draw attention to an important fact in the work of Cordeiro & Leonel (2020).Their final equation to regularize the strongly singular term, similar to Equation ( 46), has a Jacobian 0 ( ) J ξ on the right side of the equation.However, when using Equation (45) in Equation ( 46), this Jacobian is canceled out and at this point this equation differs from theirs.

Treatment of body forces
In order to efficiently apply the boundary element method without volume integration for problems involving body forces, such as gravity and centrifugal force, the present authors used the modified boundary conditions method in the context of anisotropic plane elasticity.This method was used by Caicedo & Portela (2017) and by Oliveira and Portela (2019) to regularize singular stress fields at the tip of cracks in applications of the finite element method and meshless methods, respectively.In the case of this work, the method is the same, but the objective is different from that presented in the literature.The particular solutions used for domain integrals are those presented by Deb and Banerjee (1990).In the work of Deb and Banerjee (1990), an extra integral was generated in the boundary integral equation, which led to a new vector in the matrix equation.However, in this work, no extra integrals are generated, only the boundary conditions are modified.
By decomposing the elastic field a homogeneous part H and a particular part P , stresses and displacement can be written as: ( ) where are respectively the homogeneous solutions for the stress and displacement fields that is, disregarding the body forces of the original problem, P ij σ and P i u are respectively the stress and displacement particular solutions of the original problem, representing the elastic field with the body forces.
Considering that the body forces is due to the rotation of the body (centrifugal load), the motion equation is given by: Particular solutions for this equation can be found in Deb and Banerjee (1990).For displacements, we have: And for stress, we have: ( ) Components 1 c and 2 c are represented in the form of Equations ( 55) and ( 56), for the plane stress state and for plane strain state, according the Equations ( 6) and ( 7), respectively.ρω where ρ is the density of the material and ω is the angular velocity.

NUMERICAL RESULTS
In order to assess the proposed method, four applications were analyzed, with different mesh refinements and curve orders.The first application is a problem of elastic analysis in an orthotropic rotating disk, whose analytical solution was presented by Chang (1974).The second application is the determination of the stress concentration factor in an infinity plate, whose analytical solution was presented by Lekhnitskii (1968).The third application is the calculation of the stress distribution in a rectangular plate.In this one, the solution was the same used in Nuismer & Whitney (1975).Finally, the fourth application is the calculation of the tangential stress distributions along the circumference of a hole, whose analytical solution was presented by Lekhnitskii (1968).

Elastic analysis on a rotating disk under constant angular velocity.
Consider an orthotropic rotating disk, of radius R , as shown in Figure 2. The disk is in a plane stress state, at a constant angular velocity ω .The direction of the longitudinal and transversal elastic module 1 E and 2 E respectively, are also shown in Figure 2. Equations representing stress tensor components were given by Chang (1974).When evaluated in polar coordinates, stress tensor components are written as: where r σ is the radial stress, t σ is the tangential stress (circumferential), and rt σ is the shear stress.Finally, the equations for displacement vector components are written as: ( ) where the radial displacement is written through Equations ( 58) and ( 59) in the form   Table 1 presents elastic constants and the geometric information used in the simulation of the problem.
Table 1.Mechanical properties of the material and geometric information: problem 01.

Description of variables Values of variables Description of variables Values of variables
Longitudinal elastic modulus The quarter of the disk was modeled with different amount of Source Points -SP and Integration Points -IP, as can be seen in Figure 4.These quantities are defined by each NURBS, such as 04SP and 08IP.This means the total geometry has 12 source points.Numerical results for displacements are presented for the 1 N edge in Figure 4(a) and for the 3 N edge in Figure 4(b).The circumferential stress is shown in Figure 4(c) and the radial stress in Figure 4(d), for the 1 N edge.Numerical results show that the mesh refinement with 4 source points already presents good agreement with analytical solutions.On the other hand, the number of integration points has strong influence on the numerical results, as shown in Figure 4.In order to identify the influence of mesh refinement on the numerical error, the study of convergence for displacements and stresses was carried out using a relative error measurement: Latin American Journal of Solids and Structures, 2023, 20(1), e471 14/23 where ε is the Euclidean norm of the difference between analytical and numerical solutions and n is the total number of NURBS source points.Results for the errors are shown in Figure 5, with displacements at edge 1 and at edge 3 N being shown in Figure 5(a), and circumferential stress and radial stress error in Figure 5(b), for the edge 1 N .These errors represent global errors computed 0,at edges 1 N , 2 N , and 3 N .

Stress concentration factor in a circular hole in a plate of infinite dimensions.
Consider a plate of infinite dimensions with a circular hole as shown in Figure 6.The aim of this example is the computation of the Stress Concentration Factor -SCF on a rectangular plate of an anisotropic material, with constant thickness and under an infinite with load σ in the y direction.This problem is a good example to illustrate the efficiency of higher-order IGABEM, mainly for modeling the curved hole boundary around which stresses vary rapidly.According to Lekhnitskii (1968) the stress concentration factor for an orthotropic material on a plate of infinite width is given by:  , 2 where ij A are laminate stiffness matrix terms.For the anisotropic SCF, Konish & Whitney (1975) was used, as: Latin American Journal of Solids and Structures, 2023, 20(1), e471 15/23 where i µ are roots of the characteristic equation for an anisotropic material.In addition to the maximum SCF given by Equation (63), it is also possible to write the equation for the minimum SCF: Consider a symmetrical laminate of graphite/epoxy with elastic properties given in Table 2. Plate dimensions are much larger than the hole diameter, so a uniform uniaxial stress-tension region completely surrounds the hole.Due to the symmetry of the problem, the modeling can be performed considering the geometry shown in Figure 7 with the specified dimensions, which is the representation of a quarter of the model shown in Figure 6.The geometry was built with NURBS of degree 2 p = with unidirectional fiber reinforced lamina with fiber angles varying at (0 , 45 , 90 ) θ °°°= .In Figure 7 is shown point A where the maximum stress occurs and point B where the minimum stress occurs.Finally, to obtain the numerical solution, the boundary conditions given by equation ( 65) are considered.
Table 3 shows the numerical results for the maximum SCF (point A of Figure 7) and the minimum SCF (point B of Figure 7).In the analysis of three laminate fiber directions, results are in good agreement when compared to analytical solutions, as can be seen through the percent error shown in the table, and also with quadratic discontinuous boundary element method (shape functions are Lagrangian polynomials).

Description of variables Values of variables
Longitudinal elastic modulus     Figure 8 shows displacements for the infinite plates for the three analyzed cases.The black line is the undeformed plate, and the blue line is the deformed plate.

Stress distributions along the main direction.
The analytical solution for the distribution of stresses along the principal direction of a circular hole in an orthotropic plate of infinite dimensions was presented by Lekhnitskii (1968).For the same problem, in Nuismer and Whitney (1975) we see a regression polynomial obtained by adding second, fourth, sixth and eighth order terms, which very well simulates the Lekhnitskii analytical solution.Terms of this polynomial were determined from the use of the exact value of the SCF given by Equation ( 62) for the orthotopic case or by Equation ( 63) for the anisotropic case.The polynomial equation is written as: ( ) For applications in this section, geometric configurations similar to those applied in the problem presented in section 5.2, see Figure 9, are considered, but with radius 1 r = , width of the plate   for the given materials.However, in Konish & Whitney (1975), authors point out that this solution is restricted to geometries where the plate dimensions are much larger with respect to the hole size.
To measure the magnitude of the errors between the analytical and the numerical solution, errors are plotted for each type of material, as shown in Figures 10 (b), (d), and (f).This error analysis was performed using Equation (67).As it can be seen in Figure 10, the behavior of errors is as expected, i.e., they are larger near the hole, decreasing with the increase of the distance from the hole.This application is based on a geometry similar to the problem in Figure 6, being the plate of infinite dimensions with a circular hole.However, now the stress σ is applied in the x direction, as shown in Figure 11.Note that in this application an angle ϕ is considered, which describes the direction of the fiber in the laminate.Notice that the maximum SCF is now located at point B and the minimum SCF at point A. Thus, the analytical equation for the tangential stress a the hole θ σ in Figure 11 is given as: x y the coordinates of the fibers and ( , ) x y of the laminate.
Constants k and n are given by Lekhnitskii (1968).Note that θ is the polar angle measured from the x axis and E θ is the Young's modulus in the tangent direction to the opening boundary written according to Equation ( 69), which is related to the elastic constants in the principal direction by: For computational simulation of this problem, consider the geometry shown in Figure 12 given by five NURBS curves.Information on material type, source points, control points and NURBS degree are detailed in Tables 5, 6, and 7. Boundary conditions necessary for the numerical solution of the problem are given in Equation (70).In these examples, boron/aluminium, carbon/epoxy and carbon/phenolic laminates are used, whose elastic constants were taken from Daniel et al., (2006).Figures 13(a       All numerical results for the different laminates given in Figures 13 -15 show good agreement with the analytical solution.However, it is worth mentioning two situations that we must be aware when simulating this type of problem.The first is the degree of anisotropy 1 2 / E E , which when it is large, it can increase the difficulty of the numerical results to agree to the analytical ones, as the anisotropic degree larger in Figure 13 with respect to that used in Figure 15, for example.And the second situation is that in the regions close to the hole, the numerical results always present a little more difficulty in following the analytical solution as well, as shown in the error analysis of Figures 15, for example.This can be justified, because close to the boundaries there are sudden variations of the stress.

CONCLUSION
In this paper, the IGABEM with Bézier decomposition was successfully applied to solve 2D anisotropic plane elastic problems.The Lekhnitskii's anisotropic fundamental solution was used and the singular integral terms, which arise in the boundary integral equation of displacement and traction, were regularized.For the weak singularity kernel, the Telles transformation was used.For the strong singularity term, the singularity subtraction technique was used.The convergence and numerical precision of the proposed formulation was demonstrated through the analysis of different problems.Stress distribution in a rotating disk was computed by a different approach, using modified boundary conditions.Numerical solutions computed using IGABEM have shown excellent agreement with analytical solutions available in the literature, even with a small amount of degrees of freedom.
tensor, i b is the vector of body forces, and  is the density of the material.

Figure 1 .
Figure 1.In (a) and (b) NURBS cubic curve and its control points.In (c) and (d) the basis functions of each curve.
, it is necessary to construct a term of connectivity between elements, in which a set of local basis functions are related to global functions in the form of the local basis function ( ) c , the element number ( ) e and the number of the global basis function are related by a connectivity function, see

t
are vectors of the geometric coordinates, displacement coefficients and traction coefficients respectively, associated with the control point corresponding to the basis function c for the Bézier curve e .It is important to note that the term coefficient is used, since NURBS basis functions do not necessarily obey the Kronecker Latin American Journal of Solids and Structures, 2023, 20(1), e471 7/23 delta property.Therefore e c and ( ) J ξ is the parameter that transforms the dimensionless coordinate to real coordinates.While the terms of the right-hand integrals are the original kernels of the integral, that is: is written in the form of the expansion of the Taylor series around 0 ξ ξ = , obtained through the linear term of the expansion of Equation (40).

Figure 2 .
Figure 2. Representation of the orthotropic rotating disk.

3 N
the problem can then be simulated by modeling a quarter of the disk, as shown in Figure3.The model geometry was built with three NURBS curves ( 1 N , 2 N , and) as shown in Figure3.Boundary conditions (bc) are given by Equation (60).

Figure 3 .
Figure 3.A quarter of the orthotropic rotating disk.

Figure 6 .
Figure 6.Plate with hole in infinite plane with uniform stress in y axis.

Figure 7 .
Figure 7. Illustration of a quarter of the plate with the hole pulled in the y direction.

Figure 8 .
Figure 8. Original plate (black) and deformed (blue) with fiber orientations in (a) 0 θ 30 x = , and the height of the plate 30 y = .The boundary conditions are the same as shown in Equation (65).The laminate considered in this application is a crossply laminate with stacking sequence [ ] 0 / 90 S .Numerical results were obtained based on three types of symmetrical laminates: graphite/epoxy, boron/epoxy, and aramid/epoxy as shown in Table4.

Figure 9 .
Figure 9. Illustration of the collocation points for a quarter of the plate with a hole pulled in the y direction.

Table 4 .Figure 10
Figure 10 shows the numerical results for ( ,0) y x σ for the three types of materials.A total of 105 collocation points with 8 control points were used.For plotting the graphs, only the values in (1, 0) y σ

Figure 11 .
Figure 11.Hole in infinite plane and uniform stress in x , being ( ´, ´)

Figure 12 .
Figure 12.Illustration of a quarter of the plate with a hole pulled in the x direction.
), 14(a), and 15(a) show different variations of tangential stresses around the circumference of the circular hole for the laminates.Whereas, the errors shown in Figures13(b), 14(b) and 15(b) were obtained through Equation (67).Furthermore, in Figures13(c), 14(c), and 15(c) the deformations of each laminate are shown.

Table 3 .
Numerical results for SCF maximum (point A) and SCF minimum (point B) shown in Figure7.

Table 7
NURBS data and elastic properties for carbon / phenolic (8H Fabric) laminate.