BULETINUL INSTITUTULUI POLITEHNIC DIN IAŞIPublicat deUniversitatea Tehnică …Gheorghe Asachi” din IaşiTomul LIV (LVIII), Fasc. 3, 2011SecţiaCONSTRUCŢII. ARHITECTURĂPENALTY BASED ALGORITHMS FOR FRICTIONALCONTACT PROBLEMSBYANDREI-IONUŢŞTEFANCU*, SILVIU-CRISTIAN MELENCIUCand MIHAI BUDESCU“Gheorghe Asachi” Technical University of Iaşi,Faculty of Civil Engineering and Building ServicesReceived: June 21, 2011Accepted for publication: August 29, 2011Abstract. The finite element method is a numerical method that can be successfully used to generate solutions for problems belonging to a vast array of engineering fields: stationary, transitory, linear or nonlinear problems. For the linear case, computing the solution to the given problem is a straightforward process, the displacements are obtained in a single step and all the other quantities are evaluated afterwards. When faced with a nonlinear problems, in this case with a contact nonlinearity, one needs to account for the fact that the stiffness matrix of the systems varies with the loading, the force vs. stiffness relation being unknown prior to the beginning of the analysis. Modern software using the finite element method to solve contact problems usually approaches such problems via two basic theories that, although different in their approaches, lead to the desired solutions. One of the theories is known as the penalty function method, and the other as the Lagrange multipliers method. The hereby paper briefly presents the two methods emphasizing the penalty based ones. The paper also underscores the influence of input parameters for the case of the two methods on the results when using the software ANSYS 12.Key words: finite element method; pure penalty methods; Lagrange multipliers method; H-adaptive meshing.*Corresponding author: e-mail: stefancu@ce.tuiasi.ro120Andrei-IonuţŞtefancu, Silviu-Cristian Melenciuc and Mihai Budescu1. IntroductionThe finite element method (FEM) is a numerical method that can be applied to obtain solutions to problems belonging to a variety of engineering disciplines: stationary problems, transitory problems, linear or nonlinear stress/strain problems. Heat transfer, fluid flow, or electromagnetism problems can also be solved using FEM.The first ones to publish articles in this field are Alexander Hrennikoff (1941) and Richard Courant (1943). Although their approaches are different, they have a common feature – discretization of continua into a series of discrete sub-domains called elements. Olgierd Zienkiewicz summarizes the work of his predecessors in what was to be known as FEM (1947). In 1960 Ray Clough is the first one to use the term finite element. Zienkiewicz and Cheung published in 1967 the first book entirely devoted to the finite element method.In this context, in 1971, the first version of the ANSYS software is released. Currently, the ANSYS software package is able to solve static, dynamic, heat transfer, fluid flow, electromagnetics problems, etc. ANSYS is a market leader for more than 20 years (Moavenim, 1999). In the present study version 12.0.1 of the software was used (ANSYS Workbench 2.0, 2009).2. Finite Element Formulations of Contact ProblemsThere are two basic theories that, although different in their approaches, offer the desired solutions to body contact problems: the penalty function method and the Lagrange multipliers method.The main difference between them is the way they include in their formulation the potential energy of contacting surfaces.The penalty function method, due to its economy, has received a wider acceptance. The method is very useful when solving frictional contact problems, while the Lagrange method, based on multipliers, is known for its accuracy.The main drawback of the Lagrange method is that it may lead to ill-converging solutions while the penalty formulation may lead to inaccurate ones.In the following the pure penalty, the augmented Lagrange methods will be presented.2.1. The Penalty MethodThe penalty method involves adding a penalty term to enhance the solving process. In contact problems the penalty term includes the stiffness matrix of the contact surface. The matrix results from the concept that one body imaginary penetrates the another (Wriggers et al., 1990).Bul. Inst. Polit. Ia şi, t. LVII (LXI), f. 3, 2011 121The stiffness matrix of the contact surface is added to the stiffness matrix of the contacting body, so that the incremental equation of the Finite Element becomes[K b + K c ]u = F , (1)where: K b is the stiffness matrix of contacting bodies; K c – stiffness matrix of contact surface; u – displacement; F – force.The magnitude of the contact surface is unknown (Stein & Ramm, 2003), therefore its stiffness matrix, K c , is a nonlinear term. The total load and displacement values are∑Δ=F F tot , (2)∑Δ=u u tot, (3)where: F tot is the force vector; u tot – displacement vector.To derive the stiffness matrix, the contact zone (encompassing the contact surface) is divided into a series of contact elements. The element represents the interaction between the surface node of one body with the respective element face of the other body. Fig. 1 shows a contact element in a two dimensional application. It is composed of a slave node (point S ) and a master line, connecting nodes 1 and 2. S 0 marks the slave node before the application of the load increment, and S marks the node after loading.Fig. 1 – Contact element – penalty method formulation.Given the nature of the numerical simulations presented afterwards only the sliding mode of friction will be presented. In this case, the tangential force acting at the contact surface equals the magnitude of the friction force, hence the first variation of the potential energy of a contact element ist n n d t t n n t t n n c g g k g g g k g f g f δμδδδδ)sgn( +=+=Π, (4)122 Andrei-Ionu ţ Ştefancu, Silviu-Cristian Melenciuc and Mihai Budescuwhere: k n represents penalty terms used to express the relationship between the contact force and the penetrations along the normal direction; k t – penalty terms used to express the relationship between the contact force and the penetrations along the tangential direction; g n – penetration along the normal direction; g t – penetration along the tangential direction;n n n f k g =, (5))()sgn(n n d t t g k g f μ−=. (6)2.2. The Augmented Lagrange Multiplier MethodIn the case of classical Lagrange Multiplier Method the contact forces are expressed by Lagrange multipliers. The augmented Lagrange method involves the regularization of classical Lagrange method by adding a penalty function from the penalty method (Simo & Laursen, 1992). This method, unlike the classical one, can be applied to sticking friction, sliding friction, and to africtionless contactFig. 2 – Contact element – Lagrangian methods.The contact problem involves the minimization of potential П by equating to zero the following expression:kg g g u u T T b 21)(),(+Λ+Π=ΛΠ, (7) where ⎥⎥⎦⎤⎢⎢⎣⎡⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧=Λk t k n t n t n Tλλλλλλ,...,,2211, (8) with:n λ – Lagrange multiplier for the normal direction;t λ– Lagrange multiplier for the tangential direction;Bul. Inst. Polit. Ia şi, t. LVII (LXI), f. 3, 2011 12312k n n n 12k t t t g g g g ,,...,g g g ⎡⎤⎧⎫⎧⎫⎧⎫⎪⎪⎪⎪⎪⎪=⎢⎥⎨⎬⎨⎬⎨⎬⎢⎥⎪⎪⎪⎪⎪⎪⎩⎭⎩⎭⎩⎭⎣⎦. (9)3. Parametric Analysis of Frictional ContactIn order to illustrate the way that the contact algorithms may influence the results a parametric analysis is performed. The purpose of this analysis is to exemplify how various input parameters can alter the results.3.1. Finite Elements Formulations of Contact Problems in ANSYSThe Finite Elements (FE) software ANSYS, for the penalty method, assumes that contact force along the normal direction is written as follows:cont.penetr.cont.F x K Δ=Δ, (10)where: K cont. is the contact stiffness, defined by real constant FKN for the 17x contact elements (in the current analysis the 174 contact element is used); x penetr. – distance between two existing nodes on separate contact bodies; F cont. – contact force. ANSYS automatically chooses the real constant FKN as a scale factor of the stiffness of the underlying elements. This value can be modified by the user (via FKN – a scale factor).Given the fact that the augmented Lagrange method is actually a penalty method with penetration control, the contact force is computed according to eq. (10), the only difference being the contact stiffness formulationpenetr.cont.1x K i i +=+λλ, (11)where λi is a Lagrange multiplier.Although the Lagrange multipliers are condensed out at the element level, one can think regarding this method as the same as a regular penalty one except that the contact stiffness is “updated” per contact element (Imaoka, 2001).Similar to the normal direction, a real constant – FKT models the tangential stiffness of the contact.3.2. Adaptive Solutions in ANSYSIn order to overcome the influence of the meshing upon the final results of the analysis and to improve the accuracy of the solution an adaptive solution will be used.124 Andrei-Ionu ţ Ştefancu, Silviu-Cristian Melenciuc and Mihai BudescuIn ANSYS the desired accuracy of a solution can be achieved by means of adaptive and iterative analysis, whereby h -adaptive methodology is employed.The h -adaptive method begins with an initial FE model that is refined over various iterations by replacing coarse elements with finer ones in selected regions of the model. This is effectively a selective remeshing procedure.The criterion for which elements are selected for adaptive refinement depends on geometry and, for the current analysis, on a 10% allowable difference between the maximum values of the frictional (obtained in two consecutive runs with different meshes).The user-specified accuracy is achieved when convergence is satisfied as follows:)in ...,3,2,1( ,1001R n i E i i i −=<⎟⎠⎞⎜⎝⎛−+φφφ, (11) where:φis the result quantity; E – expected accuracy (10% for this case); R – the region on the geometry that is being subjected to adaptive analysis (entire geometry in this case); i – the iteration number.The results are compared from iteration i to iteration i + 1. Iteration in this context includes a full analysis in which h -adaptive meshing and solving are performed.For this case of adaptive procedures, the ANSYS product identifies the largest elements, which are deleted and replaced with a finer FE representation (ANSYS, 2009).The overall results show a good behavior of the model. Only two iteration are performed in order to satisfy reach the expected accuracy of the solution.Table 1 h-Adaptive Methodology Convergence HistoryIteration Frictional Stress, [MPa]Change, [%] Nodes Elements 1 0.130042,518 352 2 0.12739–2.0571 14,158 8,2343.3. The Model Used in the Parametric AnalysisThe model used, represented in Fig. 3, comprises two solids made up of nonlinear structural steel materials. The larger solid has its lower surface fixed while at the upper end interacts with the smaller solid via a frictional contact (coefficient of friction 0.2). A normal pressure of 0.5 MPa is applied on top of the smaller solid, and displacement is applied on the left hand side face.Bul. Inst. Polit. Iaşi, t. LVII (LXI), f. 3, 2011 125Fig. 3 – The model used in the parametric analysis. A – frictionalcontact; B – applied pressure; C – fixed support; D – applieddisplacement.Input parameters:a) FE formulation (P1); this parameter can take two values: 0 for augmented Lagrange method and 1 for pure penalty method;b) the normal contact stiffness factor FKN (P2) that varies between 0.01 and 1;c) the tangent contact stiffness factor FKT (P3) that varies between 0.01 and 1.Output parameters: a number of output parameters have been monitored, such as: maximum (P5) and minimum (P8) normal elastic strain, maximum (P9) and minimum (P10) shear elastic strain, maximum (P12) and minimum (P13) normal stress, maximum (P14) and minimum (P15) shear stress, maximum (P11) frictional stress, maximum (P6) penetration, analysis run time (P7), maximum stiffness energy (P16).3.4. Results of the Parametric AnalysisThe parametric analysis provides a wide range of information regarding the dependence of the output parameters on the input ones. Based on the relevance of the results only a limited amount of them will be presented.The local sensitivity chart allows one to appreciate the impact of the input parameters on the output ones. This means that the output is computed based on the change of each input independently of the current value of each input parameter. The larger the change of the output, the more significant is the input parameter that was varied (ANSYS, 2009). Since the local sensibilities can only be computed for continuous parameters (P1 is a discrete one) the sensibility chart will be presented for the pure penalty and augmented Lagrange method individually.126Andrei-IonuţŞtefancu, Silviu-Cristian Melenciuc and Mihai BudescuFig. 4 – Local sensibility chart – Lagrangian method.Fig. 5 – Local sensibility chart – pure penalty method.Bul. Inst. Polit. Iaşi, t. LVII (LXI), f. 3, 2011 127It can be seen, from Fig. 5, that the pure penalty method is less sensitiveto contact normal and tangent stiffness than the augmented Lagrange method.The only output parameter influenced by the contact stiffness is the analysisrun-time.In Figs. 4 and 5 only sensitivities of the three parameters (P5, P6, P7)have been presented because the sensitivities of the other are zero or almostzero. Based on this the variation of P5, P6 and P7, with P2 and P3, arepresented in what follows.a bFig. 6 – Variation of P5 with P2 and P3: a – augmented Lagrangemethod; b – pure penalty method.a bFig. 7 – Variation of P6 with P2 and P3: a – augmented Lagrangemethod; b – pure penalty method.As it can be seen from Figs. 6 and 7 there isn’t an exact pattern of thevariation of the maximum normal elastic strain (P5), maximum penetration128Andrei-IonuţŞtefancu, Silviu-Cristian Melenciuc and Mihai Budescu(P6) or analysis run time (P7) with the normal (P2) and tangent (P3) contactstiffness factor.Give the kinematic nature of the problem, and the contact type(frictional) it can be observed from Fig. 8 that the analysis run time (the timeneeded to compute a solution for the given problem) is tangent stiffnessdependent.a bFig. 8 – Variation of P7 with P2 and P3: a – augmented Lagrangemethod; b – pure penalty method.4. ConclusionsGiven the high nonlinear characteristic of the frictional contacts an extraattention is necessary to be paid to contact algorithms and their inputparameters. In such case an h-adaptive solution is recommended to be usedbecause such approach can “fade out” the influence of contact parameters onmost of the output parameters.If working circumstances require fulfilling certain limitations, accuracyconditions may be enforced, thus improving the confidence level of the finalsolution. One must keep in mind though, that an increased number of accuracyconvergence conditions leads to prohibitive analysis run time.REFERENCESClough R.W., The Finite Element Method in Plane Stress Analysis. Proc. 2nd ASCEConf. on Electron. Comp., Pittsburg,Pennsylvania, 1960.Courant R., Variational Methods for the Solution of Problems in Equilibrum and Vibra-tions. Bul. of the Amer. Mathem. Soc., 49,1-23 (1943).Hrennikoff A., Solution of Problems of Elasticity by the Frame-Work Method. ASME,J. of Appl. Mech., 8, 169–175 1(941).Bul. Inst. Polit. Ia şi, t. LVII (LXI), f. 3, 2011 129 Imaoka S., Sheldon’s ANSYS Tips and Tricks: Understanding Lagrange Multipliers .available on-line at /ansys/tips_sheldon/STI07_Lagrange_Mul-tipliers.pdf , 2001.Moavenim A., Finite Element Analysis – Theory and Applicaion with ANSYS . PrenticeHall, Upper Saddle River, New Jersey, 1999.Simo J.C., Laursen T.A., An Augmented Lagrangian Treatment of Contact ProblemsInvolving Friction . Comp. a. Struct., 42, 1, 97-116 (1992).Stein E., Ramm E., Error-Controlled Adaptive Finite Elements in Solid Mechanics . J.Wiley a. Sons, NY, 2003.Wriggers P., Vu Van T., Stein E., Finite Element Formulation of Large DeformationImpact-Contact Problems with Friction . Comp. a. Struct., 37, 3, 319-331 (1990).Zienkiewicz O.C., Cheung Y.K., The Finite Element Method in Continuum andStructural Mechanics . McGraw Hill, NY, 1967.Zienkiewicz O.C., The Stress Distribution in Gravity Dams . J. Inst. Civ. Engng., 27,244-271 (1947). * * * ANSYS Workbench 2.0 Framework Version: 12.0.1. ANSYS Inc., 2009. * * * Design Exploration . ANSYS Inc., 2009. * * * Theory Reference, Analysis Tools, ANSYS Workbench Product Adaptive Solutions . ANSYS Inc., 2009.ALGORITMI BAZA ŢI PE METODA COREC ŢIILOR UTILIZA ŢI ÎNREZOLVAREA PROBLEMELOR DE CONTACT CU FRECARE(Rezumat)Metoda elementului finit este o metod ă numeric ă ce poate fi aplicat ă cu succes pentru a ob ţine solu ţiile problemelor dintr-o multitudine de discipline inginere şti: probleme sta ţionare, probleme tranzitorii, liniare sau neliniare. În cazul liniar g ăsirea solu ţiei unei probleme date este un proces simplu. Deplas ările sunt ob ţinute într-un singur pas de analiz ă, tensiunile şi deforma ţiile fiind evaluate ulterior. În cazul problemelor neliniare – în acest caz neliniaritate de contact – trebuie s ă se ţin ă cont de faptul c ă matricea de rigiditate a sistemului variaz ă func ţie de înc ărcare, rela ţia for ţă vs . rigiditate nefiind cunoscut ă a priori . Programele moderne, ce folosesc metoda elementului finit pentru a rezolva probleme de contact, abordeaz ă de obicei astfel de probleme prin intermediul a dou ă teorii care, de şi diferite în abord ările lor, conduc la solu ţia dorit ă. Una dintre teorii este cunoscut ă sub numele de metoda corec ţiilor, iar cealalt ă ca metoda multiplicatorilor Lagrange. În lucrare se prezint ă pe scurt cele dou ă metode, accentul punându-se pe metodele bazate pe corec ţii. Lucrarea eviden ţiaz ă, de asemenea, influen ţa parametrilor de intrare caracteristici algoritmilor de rezolvare a problemelor de contact asupra rezultatelor atunci când se utilizeaz ă pachetul software ANSYS 12.。
