Abstract
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 BSpline based discretization for the contact surface, which solves the continuity problems presented by the classic Lagrangian contact element formulation. A BSpline 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 8node hexahedron, and the formulation is continuum mechanics based, written in the current configuration, utilizing a neoHookean material model. The contact restrictions are enforced by the augmented Lagrangian method.
Keywords
Contact mechanics; BSpline; Surface smoothing; Augmented Lagrangian method; Friction
1 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 1989 and 1996, BENSON (1990)D.J. BENSON, J. O. H. A single sufrace contact algorithm for the postbuckling analysis of shell structures. Comput. Methods Appl. Mech. Eng 78 (1990) 141163, 1990. and LAURSEN AND SIMO in 1993 as authors giving important contributions to the field, as well as CURNIER in 1999, and BANDEIRA in 2001BANDEIRA, A. A. Análise de problemas de contato com atrito em 3D. Tese (Doutorado)  Departamento de Engenharia de Estruturas e Fundações, Escola Politécnica, Universidade de São Paulo, São Paulo, 2001. 275., 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 4node, 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, 2008PUSO, M. A.; LAURSEN, T. A.; SOLBERG, J. A segmenttosegment mortar contact method for quadratic elements and large deformations. Comput Methods Appl Mech Eng, 197, 2008. 555566.), the Hermitian element (WRIGGERS, 2006WRIGGERS, P. Computational Contact Mechanics. Second Edition. ed. Hannover, Germany: SpringerVerlag Berlin Heidelberg, 2006.) and (ALAIN BATAILLY, 2013ALAIN BATAILLY, B. M. N. C. A comparative study between two smoothing strategies for the simulation of contact with large sliding. Comput Mech, p. 581601, 2013.), the smoothing of the contact surface using BSplines or NURBS surfaces (CALLUM and CORBETT, 2014CALLUM J. CORBETT, R. A. S. NURBSenriched contact finite elements. Computer methods in applied mechanics and engineering, 2014.), a combination of NURBS and Mortar presented by (TEMIZER, 2011I. TEMIZER, P. W. T. J. R. H. Contact treatment in isogeometric analysis with NURBS. Comput. Methods Appl. Mech. Engrg., p. 11001112, 2011.), or the isogeometric analysis as proposed by (M.E. MATZEN, 2012M.E. MATZEN, T. C. M. B. A point to segment contact formulation for isogeometric, NURBS based finite elements. Elselvier, 2012.) and (DE LORENZIS, WRIGGERS, ZAVARISE, 2012DE LORENZIS, L.; WRIGGERS, P.; ZAVARISE, G. A mortar formulation for 3D large deformation contact using NURBSbased isogeometric analysis and the augmented Lagrangian method. Comput Mech, 49, 2012. 120.), which does away with the classical Lagrangian element, and utilizes a new contact element based on BSpline or NURBS basis function.
This paper describes an isogeometric approach, where a BSpline 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 NeoHookean material constitutive equation are presented. Then, the equation for contact forces is shown. Afterwards, the BSpline concept is briefly described, and the method for discretizing the contact surface using the BSpline 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.
2 PROBLEM DESCRIPTION
It is shown, in Figure 1, a tridimensional contact problem between two deformable solids as described in (SIMO, LAURSEN, 1992SIMO, J. C.; LAURSEN, T. A. An augmented Lagrangian Treatment of Contact Problems Involving Friction. Computers & Structures, 42, 1992. 97116.). The group of material points, that belong to the bodies in the reference position, are called
The displacement u is defined as,
with
To solve the contact problem, it is necessary to calculate the gap between the bodies
The minimum for this function can be calculated by finding the roots of the first derivate.
3 VIRTUAL WORK PRINCIPLE
The equilibrium equation, derived from the balance of linear momentum is defined by,
where
After applying virtual displacements, and bringing the tensors from the reference configuration to the current configuration, one obtains
In this equation,
It’s necessary to include the material behavior in the finite element simulation, thus an equation for a material model is necessary. The NeoHookean 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 NeoHookean material, given by
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,
and
Moreover, since the problem being simulated is a contact problem, the equations need to depict the behavior of both bodies, , as exhibited on (LAURSEN, 2002LAURSEN, T. A. Computational Contact and Impact Mechanics. [S.l.]: Springer, 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
The terms
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, 1996BATHE, K.J. Finite Element Procedures. First. ed. New Jersey: PrenticeHall, Englewood Cliffs, 1996.), (ZIENKIEWICZ, TAYLOR, ZHU, 2005ZIENKIEWICZ, O. C.; TAYLOR, R. L.; ZHU, J. Z. The Finite Element Method for Solid and Structural Mechanics. 6. ed. Oxford: Elsevier ButterworthHeinemann, 2005.), and (COOK, 2002COOK, M. P. W. Concepts and Applications of Finite Element Analysis. 4th. ed. [S.l.]: John Wiley & Sons, INC., 2002.).
4 NODE TO SURFACE FORMULATION
In this section, the nodetosurface contact formulation is expressed, for more details, refer to (WRIGGERS, 2008WRIGGERS, P. Nonlinear Finite Element Methods. Hannover, Germany: SpringerVerlag Berlin Heidelberg, 2008.) and (BANDEIRA, 2001BANDEIRA, A. A. Análise de problemas de contato com atrito em 3D. Tese (Doutorado)  Departamento de Engenharia de Estruturas e Fundações, Escola Politécnica, Universidade de São Paulo, São Paulo, 2001. 275.). The contact formulation is independent of the finite element utilized for its discretization, hence it can be used for a BSpline 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
where
where
The term
The second term of equation (4.1) accounts for the friction forces, parallel to the contact surface. On that equation,
and
The value of
and
For more details regarding the nodetosurface formulation utilizing the augmented Lagrangian method to enforce contact restrictions, one can refer to (BANDEIRA, 2001BANDEIRA, A. A. Análise de problemas de contato com atrito em 3D. Tese (Doutorado)  Departamento de Engenharia de Estruturas e Fundações, Escola Politécnica, Universidade de São Paulo, São Paulo, 2001. 275.).
Going back to (4.1),
with
and the following terms defined as follows,
and
4.1 Linearization
The solution for the equation (4.1) is discovered using NewtonRaphson’s method, and as such, the derivatives of equation (4.1) are necessary. The derivative is given by
For the first integral, the term
with
and
where
In the second term of equation (4.14),
for the stick case, with
For the slip case,
with
Following that,
and
5 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 8node hexahedron elements, and the contact surface is discretized using a BSpline based finite element.
Using a BSpline based finite element for the discretization of the contact surfaces gives a few advantages over the Lagrangian 4node contact element, such as a smooth surface and a higher degree of continuity between elements, avoiding the numerical problems associated with the nodesurface formulation. The definition for BSpline surfaces is briefly shown in the next session, and it’s followed by the procedure to discretize the contact surface using a BSpline based finite element.
5.1 BSpline surface
A BSpline 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 BSpline 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.
A point in the BSpline surface is given by
where
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.
BSpline surfaces enjoy
More detailed information about BSpline and NURBS curves can be found on (LES PIEGL, 1997LES PIEGL, W. T. The NURBS Book. 2nd. ed. [S.l.]: Springer, 1997.) and (ROGERS, 2001ROGERS, D. F. An Introduction to NURBS. [S.l.]: Morgan Kaufman Publishers, 2001.).
5.2 BSpline based finite elements
When the classic Lagrangian 4node 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, BSpline based isogeometric analysis is utilized. A single BSpline surface, or patch, is utilized to represent the whole contact surface on the master body. As such, a BSpline surface is constructed using the nodes of the master’s surface contact area.
More details about isogeometric analysis can be found on (AUSTIN COTTRELL, 2009J. AUSTIN COTTRELL, T. J. R. H. Y. B. Isogeometric Analysis. 1st. ed. [S.l.]: Wiley, 2009.) and a bidimensional isogeometric analysis can be seen on (DE LORENZIS, 2011L. DE LORENZIS, I. T. P. W. G. Z. A large deformation frictional contact formulation using NURBSbased isogeometric analysis. INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, 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 BSpline 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, 2001BANDEIRA, A. A. Análise de problemas de contato com atrito em 3D. Tese (Doutorado)  Departamento de Engenharia de Estruturas e Fundações, Escola Politécnica, Universidade de São Paulo, São Paulo, 2001. 275.) was followed, but utilizing a BSpline based finite element surface in place of a Lagrangian element contact surface mesh. As such, only the master surface is discretized with the BSpline 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 BSpline surface patch are constructed dynamically, based on the number of nodes taken as control points in the appropriate direction.
6 DISCRETIZATION OF CONTACT FORCES
Discretization of the contact work described in equation (4.1) can be written in vector form as
where
Expressing this in a vector representation
with the vectors being defined by
where
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
with
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
and the additional vectors necessary to write the equation are defined by
7 NUMERICAL EXAMPLES
Three numerical examples were considered to test the BSpline based finite element contact algorithm. The first one was the Parisch model (PARISCH, 1989PARISCH, H. A Consistent Tangent Stiffness Matrix for ThreeDimensional NonLinear Contact Analysis. International Journal for Numerical Methods in Engineering, 28, 1989. 180312.), 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 BSpline 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.
7.1 Parisch
In Figure 2, one can see the example created by (PARISCH, 1989PARISCH, H. A Consistent Tangent Stiffness Matrix for ThreeDimensional NonLinear Contact Analysis. International Journal for Numerical Methods in Engineering, 28, 1989. 180312.) and referenced in (BANDEIRA, 2001BANDEIRA, A. A. Análise de problemas de contato com atrito em 3D. Tese (Doutorado)  Departamento de Engenharia de Estruturas e Fundações, Escola Politécnica, Universidade de São Paulo, São Paulo, 2001. 275.). 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
The resulting deformations can be observed on Figure 3. The results obtained were very close to the values obtained by Parisch. The vertical coordinate of the node where the load was applied is 21.95 on the Parisch article, and the value found on this paper was 21.94, a negligible difference.
Figure 4 shows the maximum and minimum main stresses acting on the Parisch example. ANSYS doesn’t run this example, because it considers it not constrained enough.
7.2 Crossed Beams
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
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.
The obtained results for maximum and minimum main stresses for the crossed beams example are shown in figures 9 and 10.
Figures 11 and 12 show the maximum and minimum stresses obtained by ANSYS.
7.3 Parallel Beams
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
The deformations can be observed on Figure 14 and Figure 15.
In Figure 16, one can observe the ANSYS results, again with the same loads and discretization.
Figures 17 and 18 show the maximum and minimum main stresses obtained.
Figures 19 and 20 show the maximum and minumum main stresses obtained on ANSYS.
7.4 Supported Beam
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
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 figures 25 and 26 one can see the maximum and minimum main stresses for the Supported Beam example.
In figures 27 and 28, the ANSYS maximum and minimum main stresses are shown.
8 CONCLUSION
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.
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 NeoHookean energy strain equation of the material. The equation used in this paper was developed by CiarletSimo, as described on (WRIGGERS, 2008WRIGGERS, P. Nonlinear Finite Element Methods. Hannover, Germany: SpringerVerlag Berlin Heidelberg, 2008.),
On equation 8.1,
where b is the left Cauchy Green tensor, and K is given by
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 BSpline 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 BSpline based finite element is relatively simple, and further studies can be done to compare performance and convergence between it and other existing smoothing techniques.
REFERENCES
 ALAIN BATAILLY, B. M. N. C. A comparative study between two smoothing strategies for the simulation of contact with large sliding. Comput Mech, p. 581601, 2013.
 BANDEIRA, A. A. Análise de problemas de contato com atrito em 3D. Tese (Doutorado)  Departamento de Engenharia de Estruturas e Fundações, Escola Politécnica, Universidade de São Paulo, São Paulo, 2001. 275.
 BANDEIRA, A. A. et al. Algoritmos de otimização aplicados à solução de sistemas estruturais nãolineares com restrições: uma abordagem utilizando os métodos da Penalidade e do Lagrangiano Aumentado. Exacta (São Paulo. Impresso), 8, 2010. 345361.
 BATHE, K.J. Finite Element Procedures. First. ed. New Jersey: PrenticeHall, Englewood Cliffs, 1996.
 BERTSEKAS, D. P. Nonlinear programming. Belmont: Athena Scientific, 1995.
 CALLUM J. CORBETT, R. A. S. NURBSenriched contact finite elements. Computer methods in applied mechanics and engineering, 2014.
 COOK, M. P. W. Concepts and Applications of Finite Element Analysis. 4th. ed. [S.l.]: John Wiley & Sons, INC., 2002.
 D.J. BENSON, J. O. H. A single sufrace contact algorithm for the postbuckling analysis of shell structures. Comput. Methods Appl. Mech. Eng 78 (1990) 141163, 1990.
 DE LORENZIS, L.; WRIGGERS, P.; ZAVARISE, G. A mortar formulation for 3D large deformation contact using NURBSbased isogeometric analysis and the augmented Lagrangian method. Comput Mech, 49, 2012. 120.
 I. TEMIZER, P. W. T. J. R. H. Contact treatment in isogeometric analysis with NURBS. Comput. Methods Appl. Mech. Engrg., p. 11001112, 2011.
 J. AUSTIN COTTRELL, T. J. R. H. Y. B. Isogeometric Analysis. 1st. ed. [S.l.]: Wiley, 2009.
 L. DE LORENZIS, I. T. P. W. G. Z. A large deformation frictional contact formulation using NURBSbased isogeometric analysis. INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, 2011.
 LAURSEN, T. A. Computational Contact and Impact Mechanics. [S.l.]: Springer, 2002.
 LAURSEN, T. A.; SIMO, J. C. A ContinuumBased Finite Element Formulation for the Implicit Solution of Multibody, Large Deformation Frictional Contact Problems. Int. J. Num. Meth. Engng., 36, 1993a. 34513485.
 LES PIEGL, W. T. The NURBS Book. 2nd. ed. [S.l.]: Springer, 1997.
 LUENBERGER, D. G. Linear and Nonlinear Programming. 2. ed. Reading. ed. Massachusetts: AddisonWesley Publishing Company, 1984.
 M.E. MATZEN, T. C. M. B. A point to segment contact formulation for isogeometric, NURBS based finite elements. Elselvier, 2012.
 PARISCH, H. A Consistent Tangent Stiffness Matrix for ThreeDimensional NonLinear Contact Analysis. International Journal for Numerical Methods in Engineering, 28, 1989. 180312.
 PUSO, M. A.; LAURSEN, T. A.; SOLBERG, J. A segmenttosegment mortar contact method for quadratic elements and large deformations. Comput Methods Appl Mech Eng, 197, 2008. 555566.
 ROGERS, D. F. An Introduction to NURBS. [S.l.]: Morgan Kaufman Publishers, 2001.
 SIMO, J. C.; LAURSEN, T. A. An augmented Lagrangian Treatment of Contact Problems Involving Friction. Computers & Structures, 42, 1992. 97116.
 WRIGGERS, P. Computational Contact Mechanics. Second Edition. ed. Hannover, Germany: SpringerVerlag Berlin Heidelberg, 2006.
 WRIGGERS, P. Nonlinear Finite Element Methods. Hannover, Germany: SpringerVerlag Berlin Heidelberg, 2008.
 ZIENKIEWICZ, O. C.; TAYLOR, R. L.; ZHU, J. Z. The Finite Element Method for Solid and Structural Mechanics. 6. ed. Oxford: Elsevier ButterworthHeinemann, 2005.
Publication Dates

Publication in this collection
2018
History

Received
02 Feb 2018 
Reviewed
16 Apr 2018 
Accepted
20 Apr 2018