An Application of Isogeometric Analysis and Boundary Integral Element Method for Solving Nonlinear Contact Problems

A novel method to solve nonlinear contact between two bodies with plane strain behavior is presented in this paper. This method is based on the Boundary Integral Equation (BIE) and the Isogeometric Analysis (IGA). Unlike works that divide the boundary into elements, this method evaluates it as a single element, reducing the degrees of freedom in the solution. Moreover, the Particle Swarm Optimization Algorithm (PSO) was used in order to estimate the deformation when two bodies are in contact and there is penetration between them. The obtained results were compared with the Finite Element Method (FEM) and the Hertz contact equation.


Introduction
This paper presents an alternative solution to overcome the drawbacks of the boundary element method (BEM). The conventional BEM solves the continuum mechanics problems by using the reciprocal work theorem also known as Betti's theorem, the Navier-Cauchy equations, and the divergence theorem. The divergence theorem simplifies the continuum mechanics equations solution by modeling only at the boundaries [1]. Such modeling consists of dividing the boundary of the body under analysis into a discrete set of functions, denominated "elements" with nodes, and describing the displacements and tractions by polynomial functions. Since BEM is a consequence of the Betti's theorem, there are two sets of displacements and tractions. One set is known in advance through entities called traction and displacement kernels, which are valid for any equilibrium geometry. The other set belongs to the body under analysis which is the problem to be solved. Naturally, to obtain a unique solution, boundary conditions such as prescribed tractions and displacements must be applied. The calculated traction and displacement kernels depend on the normal of the surface, the geometric variable, and the interpolation points that describe the boundary. However, the discretization process produces a lack of continuity on the boundary.
An alternative to solve this problem is to use the Non-Uniform Rational B-Splines functions (NURBS) from the CAD model instead of using interpolation functions. This concept was called Isogeometric Analysis (IGA), and it was proposed by Hughes et al. [2] in 2005. Initially, the IGA was used with the FEM framework, but eventually, it was combined with BEM and with the finite cell method (FD).
Simpson et al. [3,4] proposed the Isogeometric Boundary Element Method (IGA-BEM) in 2012 and established the guidelines for combining IGA with BEM, e.g., how the elements are generated control points of the NURBS are modified repeatedly, until a correct area is found. With the proposed method, the solution domain is closed (not discrete), by defining the boundary with a single NURBS, and with fewer analysis points in contrast to traditional FEM. This paper is arranged as follows: Section 2 presents the basis to understand the isogeometric boundary element method, such as the NURBS functions and conventional BEM. Section 3 sets the guidelines to implement the modified IGA-BEM, and later, it details the methodology to be implemented in conjunction with the optimization algorithm. The following section presents the results and properties of the analyzed bodies, as well as a comparison with FEM and BEM. Finally, Section 5 concludes the work and defines future work with the theory presented.

Modeling Geometric with NURBS
A NURBS is a parametric curve that can generate lines, conics, and circles accurately. Due to this versatility, it has been widely adopted in computer-aided design programs (CAD). Mathematically, a NURBS curve is defined by Equation (1), where n is the control points number, ξ symbolizes the function parameter, B represents the control points, w the weights, N denotes the B-spline basis functions, Equations (2) and (3), with order k (degree k − 1). Figure 1 shows the basis function as a function of ξ.
The parameter ξ gets its values from the knot vector Ξ = [ξ min . . . ξ max ], which goes from a minimum ξ min to a maximum value ξ max , and its length is n + 1 + k. If C is a closed curve, the values at the ends of the knot vector are repeated k times. The weights, w, allow local control in the shape of the curve; they move the curve closer or farther to the corresponding control point.
Appl. Sci. 2020, 10, x 3 of 20 This paper is arranged as follows: Section 2 presents the basis to understand the isogeometric boundary element method, such as the NURBS functions and conventional BEM. Section 3 sets the guidelines to implement the modified IGA-BEM, and later, it details the methodology to be implemented in conjunction with the optimization algorithm. The following section presents the results and properties of the analyzed bodies, as well as a comparison with FEM and BEM. Finally, Section 5 concludes the work and defines future work with the theory presented.

Modeling Geometric with NURBS
A NURBS is a parametric curve that can generate lines, conics, and circles accurately. Due to this versatility, it has been widely adopted in computer-aided design programs (CAD). Mathematically, a NURBS curve is defined by Equation (1), where is the control points number, symbolizes the function parameter, represents the control points, the weights, denotes the B-spline basis functions, Equations (2) and (3), with order (degree − 1). Figure 1 shows the basis function as a function of .
The parameter gets its values from the knot vector Ξ = [ … ], which goes from a minimum to a maximum value , and its length is + 1 + . If is a closed curve, the values at the ends of the knot vector are repeated times. The weights, , allow local control in the shape of the curve; they move the curve closer or farther to the corresponding control point. The curve continuity can be controlled with the repeated values in the knot vector . On the other hand, a NURB curve is times derivable, and its first derivate is shown in Equation (4).
where: The curve continuity can be controlled with the repeated values in the knot vector ξ. On the other hand, a NURB curve is k times derivable, and its first derivate is shown in Equation (4).
T ij and U ij represent the traction and displacement kernels from the Kelvin solution, which gives us tractions and displacements on any point of the surface Q when a load is applied to an interior point, p, and is applicable to any geometry. The variables u i , t i , and f i are the displacements, tractions, and forces to be found. Equations (9) and (10) represent the displacement and traction kernels for the two-dimensional elastostatic problem. These kernels depend on the values of the shear modulus (µ), the Poisson coefficient (v), and the distance between the load point (p) and field point, Q, r(p, Q), described by Equation (11). For each load point and field point, displacement and traction kernels are calculated.
As Equations (8)- (11) have the load point, p, in the interior of body and field point, Q, in the boundary body, it is necessary to move the load point, p, to the boundary. Thus, Equation (7) is transformed into Equation (12), with the absence of body forces: Appl. Sci. 2020, 10, 2345 where C(P) is called "Jump Term". This element moves the internal load point (p) to the boundary, and now the load point is denoted as P. The numerical way to solve Equation (12) is to divide the body contour under analysis into "elements". Each element is composed of "nodes" (Figure 2). Shape functions describe the geometry and displacement and traction variables. To generate the kernels, a node is taken as the load point, P, and the field point, Q, is generated by each element using shape functions. The diagonal of these kernels corresponds to the point P being equal to the point Q, or P equal to any node of the element. When this happens, the values of the kernels are singular or very close to being singular. A correct approximation of these values is detailed in [25]. Then, a method as the Gaussian quadrature is applied to integrate the kernels on each element. Given these conditions, Equation (12) can be rewritten as follows: where is the local coordinate to describe the geometry, denotes the displacement, the traction at the element, is the total number of elements, represents the quadratic shape functions, and denotes the Jacobian element transformation to convert the variables from the boundary Γ to the local coordinate. The integral range [−1,1] is part of the Gaussian quadrature.

The Proposed Isogeometric Boundary Element Method
The guidelines to use IGA-BEM were established by [3]. In such a way, Equation (13) remains the same, except that the parameter replaces the term and the shape functions are the B-spline basis functions with order : However, unlike the method proposed by [3], our model generates the body boundary through a single NURBS, that is, with one "macro-element". In such a way, it is not necessary to integrate the To generate the kernels, a node is taken as the load point, P, and the field point, Q, is generated by each element using shape functions. The diagonal of these kernels corresponds to the point P being equal to the point Q, or P equal to any node of the element. When this happens, the values of the kernels are singular or very close to being singular. A correct approximation of these values is detailed in [25]. Then, a method as the Gaussian quadrature is applied to integrate the kernels on each element. Given these conditions, Equation (12) can be rewritten as follows:  (13) where ζ is the local coordinate to describe the geometry, u e i denotes the displacement, t e i the traction at the e element, n e is the total number of elements, S represents the quadratic shape functions, and J e denotes the Jacobian element transformation to convert the variables from the boundary Γ to the local coordinate. The integral range [−1, 1] is part of the Gaussian quadrature.

The Proposed Isogeometric Boundary Element Method
The guidelines to use IGA-BEM were established by [3]. In such a way, Equation (13) remains the same, except that the parameter ξ replaces the term ζ and the shape functions are the B-spline basis functions N with order k: C(P)u i (P) + However, unlike the method proposed by [3], our model generates the body boundary through a single NURBS, that is, with one "macro-element". In such a way, it is not necessary to integrate Appl. Sci. 2020, 10, 2345 6 of 19 the kernels. Nevertheless, it is required to integrate the Jacobian because the element represents the whole-body contour. The result is the body perimeter.
In addition to the boundary conditions, three sets define the problem to be analyzed: • Parameter values ξ to get the nodes P and Q on the boundary and generate the values of the kernel.

•
The control points (B) that define the shape of the body.

•
Parameter values ξ to calculate the Jacobian of the macro-element.
The proposed integral boundary equation is given by Equation (15), where η represents the ξ values that are used to calculate the Jacobian, which are different from the ξ values to find the P and Q nodes.
In conventional BEM, the jump term C was calculated indirectly, using rigid body considerations [25]. However, as Simpson et al. [4] described, these considerations are no longer valid, because of the nature of NURBS so the jump term has to be calculated explicitly [29]. However, because the geometry of the body-under-analysis is represented by a single element, the jump term equations are reduced to Equation (16). This matrix only affects the diagonal of the traction kernel when the load point, P coincides with the field point, Q.
On the other hand, because the boundary variable, Γ, is a function of ξ, the term of Jacobian transformation, J(η), is described by the following: The general flowchart of the isogeometric boundary method proposed is shown in Figure 3. The first step is to read the properties of the material-under-analysis and the NURBS geometry data. That is, the Young modulus (E), the Poisson coefficient (ν), the control points (B), the knot vector (Ξ), and the curve order (k). Then, a quantity N of load/field points is chosen. The coordinates of the load points on boundary are generated by replacing a value of ξ in the NURBS function.
The next step is to calculate the kernels for each combination of load points, P, and field points, Q. Every Kernel is saved in a global matrix of displacement (U) and traction (T) kernels with size (N × 2) × (N × 2). When the load point matches with the field point, the calculation of that traction kernel is skipped and replaced by the jump term of Equation (14). In the case of the displacement kernel, a tolerance is added or subtracted from the field point, to prevent the term (1/r) from going to infinity. The obtained global kernels are multiplied by the integral of the body boundary Equation (17); if a more exact value of the integral is desired, more ξ values may be added and may be different from the ξ values for the load and field points.
Following the flowchart in Figure 3, the next step is to apply the boundary conditions to generate a unique solution. Thus, it is necessary to accommodate the unknown variables, either displacements or tractions, on the left side, and those known on the right side, to arrive at the system of equations of the form Ax = b. The kernels would be the matrix "A," and the displacement and the tractions prescribed would be "b". Finally, by solving the equation system, the unknown values of displacement and tractions are obtained.
The next equations represent some additional elements to compute the kernels; the dy(ξ) dξ and dx(ξ) dξ terms are gotten by Equation (5).

Application: Contact between Two Bodies
The contact is analyzed by following the unilateral equations described in the Sognirini problem (contact between a rigid and an elastic body) and the Winkler-Westergaard problem (contact between Appl. Sci. 2020, 10, 2345 8 of 19 elastic bodies). One of these equations, Equation (23), establishes that the normal displacements of one body, with respects to the other, are less than or equal to zero; therefore, there is no penetration between both bodies.
u(x)· n| Γ ≤ 0 (23) If the contact is analyzed with the proposed modified IGA-BEM, there would be two bodies that share a common area that, after applying a load, reach equilibrium, as observed in Equations (24) and (25): Γ 1 and Γ 2 represent the boundaries of body 1 and body 2, respectively, that are not in contact. Γ AB represents the region of the common boundary between the two bodies, that is, the contact zone, (Figure 4). Since this region is unknown, an optimization algorithm will be applied to find it. The flowchart of the contact methodology can be seen in Figure 5.
(contact between a rigid and an elastic body) and the Winkler-Westergaard problem (contact between elastic bodies). One of these equations, Equation (23), establishes that the normal displacements of one body, with respects to the other, are less than or equal to zero; therefore, there is no penetration between both bodies.
If the contact is analyzed with the proposed modified IGA-BEM, there would be two bodies that share a common area that, after applying a load, reach equilibrium, as observed in Equations (24) and (25): Γ and Γ represent the boundaries of body 1 and body 2, respectively, that are not in contact. Γ represents the region of the common boundary between the two bodies, that is, the contact zone, ( Figure 4). Since this region is unknown, an optimization algorithm will be applied to find it. The flowchart of the contact methodology can be seen in Figure 5. The first stage of the flowchart consists of executing the complete diagram of Figure 3. The second stage is to find the intersection between the two bodies (points and ). If these elements had the same Young's modulus and the applied force is large enough, the nature of the IGA-BEM would create penetration between both bodies. This situation is not real; what physically happens is that these elements are deformed without penetration. The deformation value depends on the applied force and the materials properties. The first stage of the flowchart consists of executing the complete diagram of Figure 3. The second stage is to find the intersection between the two bodies (points A and B). If these elements had the same Young's modulus and the applied force is large enough, the nature of the IGA-BEM would create penetration between both bodies. This situation is not real; what physically happens is that these elements are deformed without penetration. The deformation value depends on the applied force and the materials properties.
The third stage is to modify the control points, using the PSO algorithm to create new body geometries with the same area as the originals but fixing the contact area. To ensure that the contact area is correct, the same load is reapplied. If the contact zone changes, then the cycle is repeated; if it does not change, or the difference between one iteration and another is minimal, then the contact algorithm is stopped. The details of the PSO algorithm are described in the next section.
The third stage is to modify the control points, using the PSO algorithm to create new body geometries with the same area as the originals but fixing the contact area. To ensure that the contact area is correct, the same load is reapplied. If the contact zone changes, then the cycle is repeated; if it does not change, or the difference between one iteration and another is minimal, then the contact algorithm is stopped. The details of the PSO algorithm are described in the next section.

The Particle Swarm Optimization
An optimization algorithm is a method to find the maximum or minimum of a given function. For the contact problem of Figure 4, the objective function is shown in Equation (26), where corresponds to the perimeter of the deformed body, and corresponds to the body perimeter before applying the load. The algorithm varies the NURBS control points that define the bodies, and the contact zone is assumed to be a straight line ( Figure 6).

The Particle Swarm Optimization
An optimization algorithm is a method to find the maximum or minimum of a given function. For the contact problem of Figure 4, the objective function is shown in Equation (26), where I corresponds to the perimeter of the deformed body, and I or corresponds to the body perimeter before applying the load. The algorithm varies the NURBS control points that define the bodies, and the contact zone is assumed to be a straight line ( Figure 6).
Appl. Sci. 2020, 10, x 10 of 20 There is a wide variety of optimization algorithms, but those algorithms based on nature stand out. For example, the genetic algorithms, the particle swarm optimization, the ant colony optimization, and the differential evolution, among others. In this paper, the particle swarm algorithm was chosen because it is easy to implement and has fast convergence [30].
The PSO is inspired by how animal herds behave. The flowchart of this method is shown in Figure 7. The first step is to set the population size or "swarm" ( ). Each "individual" of the community has a set of optimization variables or design variables that represent their position ( ). There is a wide variety of optimization algorithms, but those algorithms based on nature stand out. For example, the genetic algorithms, the particle swarm optimization, the ant colony optimization, and the differential evolution, among others. In this paper, the particle swarm algorithm was chosen because it is easy to implement and has fast convergence [30].
The PSO is inspired by how animal herds behave. The flowchart of this method is shown in Figure 7. The first step is to set the population size or "swarm" (N swarm ). Each "individual" of the community has a set of optimization variables or design variables that represent their position (X). These elements move in the solution space by adjusting their speeds (V). However, the individual speed depends on their previous position and speed values. Moreover, it depends on the best individual location (L) and the best position of all individuals (G). The algorithm runs until the maximum number of iterations (N i ) is reached. The last best global position values correspond to the design optimized variables.
The function to be optimized is the initial perimeter of the bodies, and the design variables are the control points (B) that are outside of the contact zone. The next position X i and velocity V i values of each particle are calculated by using Equations (27) and (28), respectively, where the best particle position is L i , and the best position of all particles is G N .
To update the best local and global values, the function to be optimized must be analyzed. This means the perimeter of the body for each X i value must be calculated. The elements of X i are the new control points. If the perimeter is higher or less than the original, then a penalization is applied by adding or subtracting an amount from the obtained perimeter value, Equations (29) and (30), where k represents a proportional constant to the difference between the value of I and I or . To update the best local position, L i , the penalized I is compared with the best saved local perimeter; if the new I value is better than the saved, the best local position L i is updated. The same procedure is done to update the best global location.

Results
The methodology described in Section 3.1 was programmed to solve a plane strain problem, and its results were compared with FEM and the conventional BEM. Figure 8 represent the problem: A load was applied on the top of the cylinder, and the opposite side was fixed. Being a plane problem, it is not necessary to create the cylinder model; instead, only the cross-section is required. The

Results
The methodology described in Section 3.1 was programmed to solve a plane strain problem, and its results were compared with FEM and the conventional BEM. Figure 8 represent the problem: A load was applied on the top of the cylinder, and the opposite side was fixed. Being a plane problem, it is not necessary to create the cylinder model; instead, only the cross-section is required. The cylinder dimensions and the material properties can be found in Table 1, likewise the load value. The results of this problem are presented in Sections 4.1-4.3.

Results
The methodology described in Section 3.1 was programmed to solve a plane strain problem, and its results were compared with FEM and the conventional BEM. Figure 8 represent the problem: A load was applied on the top of the cylinder, and the opposite side was fixed. Being a plane problem, it is not necessary to create the cylinder model; instead, only the cross-section is required. The cylinder dimensions and the material properties can be found in Table 1, likewise the load value. The results of this problem are presented in Sections 4.1-4.3.  The proposed modified IGA-BEM was used to analyze the contact between two cylinders. The contact problem was also analyzed as a plane strain problem. Both cylinders have the material The proposed modified IGA-BEM was used to analyze the contact between two cylinders. The contact problem was also analyzed as a plane strain problem. Both cylinders have the material properties defined in Table 1. The load was applied as shown in Figure 9, and a displacement restriction was placed at the bottom of the cylinders. A symmetric behavior in the deformation is assumed as the contact zone is the same for both circles. The results of this problem are presented in Sections 4.4 and 4.5.
Appl. Sci. 2020, 10, x 12 of 20 properties defined in Table 1. The load was applied as shown in Figure 9, and a displacement restriction was placed at the bottom of the cylinders. A symmetric behavior in the deformation is assumed as the contact zone is the same for both circles. The results of this problem are presented in Sections 4.4 and 4.5.

Plane Strain Problem-Solve with FEM
A commercial FEM software was used to solve the problem presented in Figure 8. First, a load was set on the upper side, and the load was distributed between two neighboring nodes. Next, displacement restrictions on the nodes were added on the opposite side. A fine mesh was made,

Plane Strain Problem-Solve with FEM
A commercial FEM software was used to solve the problem presented in Figure 8. First, a load was set on the upper side, and the load was distributed between two neighboring nodes. Next, displacement restrictions on the nodes were added on the opposite side. A fine mesh was made, generating 2566 nodes and 827 elements (Figure 10a). The result of this problem is shown in Figure 10b.

Plane Strain Problem-Solve with FEM
A commercial FEM software was used to solve the problem presented in Figure 8. First, a load was set on the upper side, and the load was distributed between two neighboring nodes. Next, displacement restrictions on the nodes were added on the opposite side. A fine mesh was made, generating 2566 nodes and 827 elements (Figure 10a). The result of this problem is shown in Figure  10b.

Plane Strain Problem-Solve with Conventional BEM
The strain problem of Figure 8 was solved with a program based on BEM and the algorithm proposed by [25]. To solve the problem, the following data was used: 71 nodes, 36 elements distributed on the boundary (Figure 11a), 6 integration points for Gaussian Quadrature, quadratic shape functions to describe tractions and displacements, and the data from Table 1. The result can be found in Figure 11b. To generate the internal displacements, 13 internal points were chosen using the guidelines described by [25].

Plane Strain Problem-Solve with Conventional BEM
The strain problem of Figure 8 was solved with a program based on BEM and the algorithm proposed by [25]. To solve the problem, the following data was used: 71 nodes, 36 elements distributed on the boundary (Figure 11a), 6 integration points for Gaussian Quadrature, quadratic shape functions to describe tractions and displacements, and the data from Table 1. The result can be found in Figure 11b. To generate the internal displacements, 13 internal points were chosen using the guidelines described by [25].

Plane Strain Problem Solve with the Isogeometric-Boundary Integral Element Method Proposed
To create the circle with NURBS, the curve order, the control points, and the knot vector were first calculated. Those values are shown in Table 2. As seen in Section 3.2, the load points have to be generated on the boundary. Therefore, 24 values of were chosen, to generate the coordinates of the load/field points, using Equation (1). Table 3 specifies the additional elements to use the proposed method. The used load/field points and control points are shown in Figure 12a. The deformed cylinder is shown in Figure 12b. Like the conventional BEM, the internal displacements are calculated by adding internal points.

Plane Strain Problem Solve with the Isogeometric-Boundary Integral Element Method Proposed
To create the circle with NURBS, the curve order, the control points, and the knot vector were first calculated. Those values are shown in Table 2. As seen in Section 3.2, the load points have to be generated on the boundary. Therefore, 24 values of ξ were chosen, to generate the coordinates of the load/field points, using Equation (1). Table 3 specifies the additional elements to use the proposed method. The used load/field points and control points are shown in Figure 12a. The deformed cylinder is shown in Figure 12b. Like the conventional BEM, the internal displacements are calculated by adding internal points.

Contact Problem Solve with Hertz Theory
According to [31], the width of the rectangular contact area between two cylinders can be found by Equation (28) when they have the same elasticity modulus ( = = ) and the same Poisson In Equation (31), represents the load per unit length ( / ), is the relative curvature ( ⋅ + ) and, is the width of the rectangular of contact area. More detailed information about Hertz theory can be found at [13]. The length of the cylinder is 100 mm.
If the values described in Table 1 and in the previous paragraph are substituted in Equation (31), the width of the contact area is = .056 . Following Hertz's theory, the maximum Hertzian contact pressure, the maximum shear stress, and its depth can be found. These values are shown in Table 4.

Contact Problem Solve with Hertz Theory
According to [31], the width of the rectangular contact area between two cylinders can be found by Equation (28) when they have the same elasticity modulus (E = E 1 = E 2 ) and the same Poisson coefficient (v = 0.3). In Equation (31), p represents the load per unit length (P/L), K D is the relative curvature ( D 1 ·D 2 D 1 +D 2 ) and, b is the width of the rectangular of contact area. More detailed information about Hertz theory can be found at [13]. The length of the cylinder is 100 mm.
If the values described in Table 1 and in the previous paragraph are substituted in Equation (31), the width of the contact area is b = .056 mm. Following Hertz's theory, the maximum Hertzian contact pressure, the maximum shear stress, and its depth can be found. These values are shown in Table 4.

Contact Problem Solve with FEM
The contact problem was also solved with FEM software. The parameters used in the simulation are described in Table 5. A fine mesh was chosen to obtain greater results' precision ( Figure 13a). To avoid penetration between the bodies, the augmented Lagrange method was chosen. The behavior of this algorithm was briefly described in the introduction section. Figure 13b shows the deformed bodies, and the color map represents the displacements in Y. The maximum deformation is located in the load application area. The contact zone has a length of approximately 0.65 mm, as seen in Figure 14.

Contact Problem Solved with Proposed IGA-BEM.
To use the proposed contact algorithm from Section 3.2, some control points were added near to the common point (0, 0). The NURBS data to create the two cylinders in contact are in Table 4. As seen, weights were also changed to maintain the circular shape of the cylinders. Table 5 specifies the data to solve the problem with the proposed IGA-BEM. The load points were generated by selecting 24 values of distributed between the minimum and maximum value of this variable. Unlike the

Contact Problem Solved with Proposed IGA-BEM.
To use the proposed contact algorithm from Section 3.2, some control points were added near to the common point (0, 0). The NURBS data to create the two cylinders in contact are in Table 4. As seen, weights were also changed to maintain the circular shape of the cylinders. Table 5 specifies the data to solve the problem with the proposed IGA-BEM. The load points were generated by selecting 24 values of distributed between the minimum and maximum value of this variable. Unlike the

Contact Problem Solved with Proposed IGA-BEM.
To use the proposed contact algorithm from Section 3.2, some control points were added near to the common point (0, 0). The NURBS data to create the two cylinders in contact are in Table 6. As seen, weights were also changed to maintain the circular shape of the cylinders. Table 7 specifies the data to solve the problem with the proposed IGA-BEM. The load points were generated by selecting 24 values of ξ distributed between the minimum and maximum value of this variable. Unlike the problem in Section 4.3, the control points increased 9-13 ( Figure 6). Finally, the number of points to calculate the Jacobian was set to 50, to obtain a more accurate value; however, this value is only calculated once. Figure 15 shows the two cylinders deformed after applying the load. Figure 16 is a close-up to the contact area and shows the interference between the bodies. Table 8 details the input parameters of the optimization algorithm. The values of C1 and C2 were set according to the recommendation of Clerc [31].   problem in Section 4.3, the control points increased 9-13 ( Figure 6). Finally, the number of points to calculate the Jacobian was set to 50, to obtain a more accurate value; however, this value is only calculated once. Figure 15 shows the two cylinders deformed after applying the load. Figure 16 is a close-up to the contact area and shows the interference between the bodies. Table 6 details the input parameters of the optimization algorithm. The values of C1 and C2 were set according to the recommendation of Clerc [31].  Curve order = 3 Table 5. Data to analyze the contact problem with the proposed method.
Load/field points number 24

Elements number 1
Control points number (for NURBS geometry) 13 Jacobian points number 50     Figure 17 shows the bodies in contact without penetration, after having implemented the algorithm proposed in Section 3.2. As can be seen, in Figure 18, the contact area is a straight line. This zone has a length of approximately 62 mm.     Figure 17 shows the bodies in contact without penetration, after having implemented the algorithm proposed in Section 3.2. As can be seen, in Figure 18, the contact area is a straight line. This zone has a length of approximately 62 mm.  Figure 17 shows the bodies in contact without penetration, after having implemented the algorithm proposed in Section 3.2. As can be seen, in Figure 18, the contact area is a straight line. This zone has a length of approximately 62 mm.

Discussion
As seen in the previous section, the three methods generated similar deformed bodies when solving the plane strain problem. With a fine mesh, the maximum displacement in the Y-axis with FEM was −0.094 mm, while with the conventional BEM, it was −0.023 mm, and with the proposed IGA-BEM method, it was −0.073 mm. The maximum displacement was seen where the P load was applied. To solve the problem, the FEM software used 2566 nodes or analysis points, while BEM used 71 analysis points, and the proposed IGA-BEM used 24 analysis points. Internal points were added in BEM and proposed IGA-BEM to obtain the internal displacements. In contrast, given the nature of FEM, the analysis did not require adding additional points.
Regarding the contact problem, Table 9 shows a comparison of the contact area widths obtained with the Hertz model, the FEM software, and the proposed IGA-BEM of contact. Taking the analytical solution as a reference, the error obtained with FEM was 16%, while the proposed method was 10%. In both cases, interference between the bodies was eliminated after applying the load. With FEM, the problem was solved by using a total of 2184 nodes, equivalent to 4368 Degrees of Freedom (DOF). With the proposed IGA-BEM, the problem was solved using 24 analysis points (48 DOF), but the control points that define the bodies' shape increased 9-13. On the other hand, the deformed figure was obtained after the PSO algorithm converged in a value gap.
Although the proposed method solves the problems with fewer degrees of freedom compared to FEM, it requires more robustness to avoid sensitivity problems in the calculation of n x and n y ; The same issue was also observed with the conventional BEM.

Conclusions
The proposed IGA-BEM got consistent results, while solving the problem of Figure 8, compared to those obtained from the proven methods FEM and BEM. Although the displacement values were different, the deformation generated was similar in all three methods. However, fewer degrees of freedom were required in proposed IGA-BEM, because a single NURBS was used to model the completed body. Moreover, the load/field points used in the conventional BEM were separated from the points, to represent the geometry.
Regarding the nonlinear contact problem, the contact-zone widths obtained with FEM and with the proposed contact method varied 16% and 10%, respectively, compared to the one obtained analytically, with the Hertz method. Nevertheless, the DOF used by FEM to solve the contact problem were 4368, while the proposed method only used 48 DOF.
Decreasing degrees of freedom is very significant, for example, when solving multiphysics problems, where it is required to invert matrices. Furthermore, as seen in [32], the propagation error is tied to the size of the matrix: The more elements the matrix has, the higher the error.
Although the results are promising, future work is required to prove that the methodology is valid in other situations where displacements are not symmetrical or additional factors, such as friction, are considered.