Numerical modeling of contact problems with the finite element method utilizing a B-Spline surface for contact surface smoothing

The basic aim of this work is to present a method to treat structural mechanics problems dealing with tridimensional contact and large elastic deformations. The formulation presented in this article offers a B-Spline based discretization for the contact surface, which solves the continuity problems presented by the classic Lagrangian contact element formulation. A B-Spline surface is utilized to discretize the contact surface, and the BSpline basis functions are used to distribute the contact forces between the contact surfaces nodes. The finite element utilized is an 8-node hexahedron, and the formulation is continuum mechanics based, written in the current configuration, utilizing a neo-Hookean material model. The contact restrictions are enforced by the augmented Lagrangian method. Key-words Contact mechanics, B-Spline, Surface smoothing, Augmented Lagrangian method, Friction.


INTRODUCTION
Finite element contact simulations exhibiting large deformations has been growing as a subject of study in the last few decades.Many works have been produced about it, and one can cite WRIGGERS in 1989and 1996, BENSON (1990) and LAURSEN AND SIMO in 1993 as authors giving important contributions to the field, as well as CURNIER in 1999, andBANDEIRA in 2001, who offered further research on the subject.
The classical way to discretize a contact surface on a tridimensional problem is through the Lagrangian contact element.The Lagrangian contact element for a hexahedron finite element mesh is a 4-node, linear element, which results in a flat surface.When discretizing a curved surface, the Lagrangian contact element, being a linear element, leads to a facetization of the contact surface, i. e., the curved surface is represented by several flat surfaces, side by side, with an angle between them.This facetization leads to a geometric discontinuity of the contact surface, and therefore, generates discontinuity of the normal vector.This in turn, leads to convergence problems, especially when large sliding occurs.
To tackle this issue, many approaches have been proposed.One may consider, between many options, the Mortar method, described in (PUSO, LAURSEN and SOLBERG, 2008), the Hermitian element (WRIGGERS, 2006) and (ALAIN BATAILLY, 2013), the smoothing of the contact surface using B-Splines or NURBS surfaces (CALLUM and CORBETT, 2014), a combination of NURBS and Mortar presented by (TEMIZER, 2011), or the isogeometric analysis as proposed by (M.E.MATZEN, 2012) and (DE LORENZIS, WRIGGERS, ZAVARISE, 2012), which does away with the classical Lagrangian element, and utilizes a new contact element based on B-Spline or NURBS basis function.
This paper describes an isogeometric approach, where a B-Spline patch is used to discretize a contact surface over the master body.The difference between this paper and other works which utilize the isogeometric approach is that no smoothing of the slave surface is done.The smoothing is done only on the master surfaces.Each node of the slave surface is tested individually, and the slave surface conforms to the master surface.The authors believe that when contact happens, the slave surface will assume the master surface's shape, and thus will also be smoothed.
In the first section of this article, a brief description of the contact problem, the virtual work applied to the balance of linear momentum equation and the Neo-Hookean material constitutive equation are presented.Then, the equation for contact forces is shown.Afterwards, the B-Spline concept is briefly described, and the method for discretizing the contact surface using the B-Spline is detailed.The relevant equations for the contact pressures, their linearization and finally the relevant vectors for the discretization of the contact surfaces are presented.
In the last section numerical examples are shown, utilizing the formulation presented in this paper.The results obtained are compared to the commercial software ANSYS and afterwards the results are briefly discussed.

PROBLEM DESCRIPTION
It is shown, in Figure 1, a tridimensional contact problem between two deformable solids as described in (SIMO, LAURSEN, 1992).The group of material points, that belong to the bodies in the reference position, are called B and B .Additional configurations can be obtained by the transformations φ : B × I → R and φ : B × I → R in the interval t = [0, T].The bodies B and B , subject to contact, are called slave and master.In each of these bodies, there is a boundary where contact can be expected.These boundaries are called Г and Г .Г is the slave boundary, belonging to the body B , and Г the master boundary, belonging to the body B .These boundaries are in the reference configuration, and are subject to the same transformations φ and φ that the bodies B and B are subject.A point belonging to body B is denoted by X in the reference configuration and x on the current configuration, and the same applies for body B , exchanging the index.To solve the contact problem, it is necessary to calculate the gap between the bodies B 1 and B 2 at the points where contact might occur.This is done by searching in Г 2 for a point x 2 , closest to a point x 1 in Г 1 , by the means of the function The minimum for this function can be calculated by finding the roots of the first derivate.

VIRTUAL WORK PRINCIPLE
The equilibrium equation, derived from the balance of linear momentum is defined by, where  is the first Piola Kirchoff tensor, ρ is the body's density and  ̅ the body forces for unit of volume.
After applying virtual displacements, and bringing the tensors from the reference configuration to the current configuration, one obtains In this equation,  is Cauchy's tensor and ρ t is density for current unit of volume.This equation represents equilibrium between the internal work and the external work defined at the current configuration.
It's necessary to include the material behavior in the finite element simulation, thus an equation for a material model is necessary.The Neo-Hookean material model is adopted on this article.This material model is adequate for simulating large tridimensional displacements.Therefore, Cauchy's tensor on the equation (3.2) is obtained through the deformation energy function for a Neo-Hookean material, given by  = (J − 1) + ( − ).Moreover, since the problem being simulated is a contact problem, the equations need to depict the behavior of both bodies, , as exhibited on (LAURSEN, 2002).Equation (3.2) is then expanded to include both bodies, To satisfy the impenetrability condition and to represent the contact forces, a contact energy term is added to equation (3.6).
The definition of П depends upon the restriction method utilized.In this paper, the augmented Lagrangian method was utilized, as described in (BANDEIRA, et al., 2010), (BERTSEKAS, 1995) and (LUENBERGER, 1984).The contact equation is defined, as shown in (BANDEIRA, 2001), (WRIGGERS, 2008) and (LAURSEN and SIMO, 1993a), by П = ∫ t δg dA + ∫  δ dA .(3.8) The terms t and  represent the normal and tangencial contact pressures, respectively.g represents the gap function, which notates the penetration between the bodies, and  is the tangencial displacement function, which represents how much one of the bodies slides over the other.
The first integral on equation (2.8) accounts for the work done by normal forces, while the second integral accounts for the work done by the tangential forces.If the problem doesn't include friction, the second integral is null.
Equation (3.7) then can be discretized and solved by the finite element method.In this paper, the isoparametric hexahedral element was used.More details on the finite element discretization can be found on (BATHE, 1996), (ZIENKIEWICZ, TAYLOR, ZHU, 2005), and(COOK, 2002).

NODE TO SURFACE FORMULATION
In this section, the node-to-surface contact formulation is expressed, for more details, refer to (WRIGGERS, 2008) and (BANDEIRA, 2001).The contact formulation is independent of the finite element utilized for its discretization, hence it can be used for a B-Spline based finite element.
The equation for the contact forces as found in BANDEIRA is defined in (3.8) and is repeated here for convenience, The first term accounts for the normal forces and the second one to the tangential forces.When the augmented Lagrangian method is used to enforce the contact restriction, the term t is defined by where λ is the Lagrangian multiplier, ξ is the penalty factor and g is the gap function.In turn, g which is defined by where  is the point of the slave surface currently analyzed, or rather, a node, and  is the projection of the slave point on the master surface.
The term δg is given by δg The second term of equation ( 4.1) accounts for the friction forces, parallel to the contact surface.On that equation,  depends on the friction model utilized.In this paper, the Coulomb friction model was deemed sufficient.In the Coulomb friction model, the term  depends on whether the friction is on the stick case, where the two bodies don't move relative to each other, or the slip case, where one body slides over the other.Depending on which is the current case,  is calculated differently.The following functions are used to calculate it.
The value of Φ defines whether the current node is on a stick or slip case, and then t is defined accordingly.t is defined by and Φ is given by ϕ For more details regarding the node-to-surface formulation utilizing the augmented Lagrangian method to enforce contact restrictions, one can refer to (BANDEIRA, 2001).

CONTACT ELEMENT DISCRETIZATION
The solids bodies involved in the contact problem must be discretized into a finite element mesh.In this paper, the bodies are discretized by 8-node hexahedron elements, and the contact surface is discretized using a B-Spline based finite element.
Using a B-Spline based finite element for the discretization of the contact surfaces gives a few advantages over the Lagrangian 4-node contact element, such as a smooth surface and a higher degree of continuity between elements, avoiding the numerical problems associated with the node-surface formulation.The definition for B-Spline surfaces is briefly shown in the next session, and it's followed by the procedure to discretize the contact surface using a B-Spline based finite element.

B-Spline surface
A B-Spline surface is constructed based on the interpolation of a number of spatial points.These points are called control points, and their spatial coordinates are interpolated using a number of functions called base functions to construct a tridimensional surface.
The B-Spline parametric surface is constructed using the base functions.In turn, the base functions are built based on a parameterization of space given by the knot vectors, i. e. where  is a control point, N and M are the base functions, one increasing in a different direction to the other.k and l are the order of the base functions.The B-Spline base function is given by . (5.3) The base functions are built recursively, with higher orders needing more calculations and thus processing power on a numerical analysis.An order of 1 reduces the base function to a Lagrangian base function.
B-Spline surfaces enjoy  continuity, where k is the order of the base function, and m is the multiplicity of an inner knot i.In this paper, other than an open knot vector, no knot vector manipulation was used, so the multiplicity of any inner knot is 1.
More detailed information about B-Spline and NURBS curves can be found on (LES PIEGL, 1997) and (ROGERS, 2001).

B-Spline based finite elements
When the classic Lagrangian 4-node contact element is used, multiple contact elements are utilized to form a contact surface.When there's contact between a slave node and the master surface, the resulting contact forces are distributed between the slave node, and the four nodes belonging to the single contact element where contact was detected.
On this paper, B-Spline based isogeometric analysis is utilized.A single B-Spline surface, or patch, is utilized to represent the whole contact surface on the master body.As such, a B-Spline surface is constructed using the nodes of the master's surface contact area.
More details about isogeometric analysis can be found on (AUSTIN COTTRELL, 2009) and a bi-dimensional isogeometric analysis can be seen on (DE LORENZIS, 2011).
In case contact between a slave node and the master surface is found, the resulting contact forces are distributed among the slave node and contact surface nodes, using the B-Spline base functions as criteria for the distribution.As a general rule, the further away from the contact point the node is, smaller is the contact force acting on that node.
For the numerical examples executed, the methodology found in (BANDEIRA, 2001) was followed, but utilizing a B-Spline based finite element surface in place of a Lagrangian element contact surface mesh.As such, only the master surface is discretized with the B-Spline finite elements, the slave surface is not discretized, but each of its nodes is checked for contact against the master surface.
The knot vectors used to construct the B-Spline surface patch are constructed dynamically, based on the number of nodes taken as control points in the appropriate direction.

DISCRETIZATION OF CONTACT FORCES
Discretization of the contact work described in equation (4.1) can be written in vector form as W (, ) = ∑ A  ⋅  , (6.1) where A represents the contact finite element area, ns is the total number of slave nodes,  represents the derivative of the spatial position of the element nodes and  represents the contact stresses.Expressing this in a vector representation  = t  + t  = t  + t  + t  (6.2) with the vectors being defined by where Q / is defined as Q / = N , (t)M , (u), (6.3) with "a" and "b" being the index number of the control points that construct the surface patch in each direction, and "k" and "l" are the selected base function's order, omitted for brevity.The linearization of equation ( 6.1) is necessary to solve the problem for the displacements, so it's given by  ∆ = ∑ A ( )  ∆ (6.4) with  representing the stiffness matrix of the element.To create the stiffness matrix contribution due to contact forces, a discretization of the equation is necessary.This discretization, in vector form, is given by  = ϵ  + t   + a   −  n ⋅ a , (6.5) and the additional vectors necessary to write the equation are defined by

NUMERICAL EXAMPLES
Three numerical examples were considered to test the B-Spline based finite element contact algorithm.The first one was the Parisch model (PARISCH, 1989), which was executed to compare directly with the results obtained by the Lagrangian finite element under large displacements.The second one was an authors created example, two beams subject to a significant bending moment in which both show large displacements.Finally, the third example shows two parallel beams, and it was designed to test the B-Spline finite element when the contact surface is generated after the body is significantly deformed.Lastly, the final example was designed to exhibit an amount of sliding over the contact surface.

Parisch
In Figure 2, one can see the example created by (PARISCH, 1989) and referenced in (BANDEIRA, 2001).It consists of two blocks, shaped like square based prisms, one supporting the other.The inferior block has all its base nodes restricted in all directions, while none of the superior block's nodes has any restrictions.The superior block is supported only by contact forces.There are 50 nodes and 13 elements on this example, and more geometry information can be seen on Figure 2.
This example was executed without any tangential forces (no friction).The normal penalty parameter was ε = 10 , the blocks elasticity modules were 10 Pa for the superior block and 10 0 Pa for the inferior block, and it needed only a single load increment to reach convergence.The second example can be observed on Figure 5. Created by the authors, the example consists of two beams crossed and superposed, with a 0.2 mm gap between them.On the free end of the superior beam, at the edge, a load totaling 9000N is applied.This example has 288 nodes and 120 elements.More information about the example geometry can be observed on Figure 5.This example was simulated with friction, the contact penalty parameters were ε = 10 for normal penalty and ε = 10 for the tangential penalty.The superior beam has an elasticity modulus of 4.10 Pa, while the inferior beam has 10 Pa.The friction coefficient used was μ = 0.2, and three load increments were applied.
The resulting displacements can be observed in Figure 6 and Figure 7.In Figure 8, one can observe the results obtained on ANSYS for this example, with the same loads, number of nodes and elements.
Numerical modeling of contact problems with the finite element method utilizing a B-Spline surface for contact surface smoothing    Illustrated in Figure 13, this example consists in two parallel beams, restricted on both ends, with forces applied perpendicularly.There are 384 and 180 elements in this example.The loads were applied in a way to arch the beams before contact happens, so the contact surfaces are curved.
The superior beam has double the width of the inferior one; otherwise they're the same size.Information about restricted nodes, load positioning and details about the dimensions can be observed on Figure 8.The elasticity modulus of both beams is 10 Pa.Friction was used and the friction coefficient is μ = 0.2.30 load increments were utilized, and the normal and tangential penalty values are ε = 10 and ε = 10 , respectively.The deformations can be observed on Figure 14 and Figure 15.The last example shown, observed in Figure 21, consists of a beam fixed on one end, suspended over another beam which is fixed on both ends.There's a 0.5 mm gap between them, and there are twelve 150N loads applied on the nodes of the free end of the superior beam.Details about the loads positions and the example's geometry can be observed on Figure 12.This example consists of 441 nodes and 200 elements.The elasticity modulus for both beams is 10 Pa, the normal penalty utilized was ε = 5. 10 and 30 load increments were necessary for convergence.
The displacement results can be observed on Figure 22 and Figure 23.The ANSYS results from the same mesh can be observed on Figure 24.In this paper, problems exhibiting large displacements, curved contact surfaces, as well as large sliding were presented.In all cases presented, convergence was attained, and the displacements obtained were consistent with the results generated in ANSYS, as one can observe on table 1.The Parisch numerical example isn't constrained enough to run on ANSYS, and was left out of the comparisons.
Numerical modeling of contact problems with the finite element method utilizing a B-Spline surface for contact surface smoothing Latin American Journal of Solids and Structures, 2018, 15(8), e77 22/23 The stress values found on this paper diverged from the ones found on ANSYS by a fairly significant amount.The reason for this discrepance is a difference in the Neo-Hookean energy strain equation of the material.The equation used in this paper was developed by Ciarlet-Simo, as described on (WRIGGERS, 2008) On equation 8.1, J is the jacobian, and μ , λ are the Lamé parameters.The Neo-Hookean energy strain equation used by ANSYS, as found on its help file, is given by W (I , J) = (J − 1) + μ ( − ), (8.2) where b is the left Cauchy Green tensor, and K is given by K = λ + μ .(8.3) The ANSYS stress results can still be very valuable for a qualitative comparison.One can observe that the stress distributions displayed on the results found on this paper and the ones derived from ANSYS are consistently similar.
It was shown that B-Spline based finite element and the discretization technique utilized on this paper are valid forms of smoothing the contact surface and solving the facetization problem.It provides consistent stresses across the contact surfaces and forbids penetration, especially with better discretized meshes.
The implementation of the B-Spline based finite element is relatively simple, and further studies can be done to compare performance and convergence between it and other existing smoothing techniques.

Figure 1 :
Figure 1:Notation for the contact problema with friction (3.3) Equation (3.2) can be divided into two parts, internal and external forces.The first term represents the internal forces, and the second and third terms represent the external forces, Ξ = [ξ , ξ , ξ … ξ ]. (5.1) represents the knots.The knots define the limits of the base functions domain, by dividing the parametric space into discrete elements.The knot vector has  elements,  being the sum of the number of control points plus the order k of the base function.The order of the base function represents the number of intervals between knots that composes the function's domain.An open knot vector has its first and last entries repeated k times.The open knot vector makes the curve or surface to behave the same when close to its borders as it behaves in the middle.A point in the B-Spline surface is given by (t, u) = ∑ ∑  , N , (ξ)M , (η) , (5.2)

Figure 5 :
Figure 5: Crossed Beams example Numerical modeling of contact problems with the finite element method utilizing a B-Spline surface for contact surface smoothing Latin American Journal of Solids and Structures, 2018, 15(8), e7711/23

Figure 9 :
Figure 9: Crossed Beams maximum main stress

Figure 11 :
Figure 11: Crossed Beams ANSYS maximum main stress

Figure 17 :
Figure 17: Parallel Beams maximum main stress

Figure 19 :
Figure 19: Parallel Beams ANSYS maximum main stress

Table 1 :
Result comparison