VIEWS: 5 PAGES: 10 CATEGORY: Guides POSTED ON: 5/26/2010 Public Domain
THERMAL DESIGN OF SLIDING SYSTEMS BY DOMAIN DECOMPOSITION TECHNIQUES Paolo Decuzzi1, Tommaso Soranno1 and Przemyslaw Zagrodzki2 1 Center of Excellence in Computational Mechanics - Politecnico di Bari, Via Re David 200, 70125, BARI; ITALY; e-mail: p.decuzzi@poliba.it; tsoranno@hotmail.com 2 Raytech Composites Inc., Crawfordsville, IN 47933; USA; e-mail: pzagrodzki@raybestosproducts.com SOMMARIO La progettazione di freni e frizioni multidisco richiede l’analisi teorico, numerica e sperimentale di un sistema caratterizzato da una notevole complessità geometrica, dalla interazione fra deformazioni termoelastiche e scambio termico, dalla rapida evoluzione delle variabili termomeccaniche durante le manovre di innesto e/o frenatura. Tutto questo comporta elevati tempi di calcolo. In questo lavoro, si presenta un nuovo approccio per la soluzione numerica di problemi termomeccanici per sistemi di slittamento basato su: (i) l’impiego di sistemi di riferimento locali (Formulazione Lagrangiana), e (ii) la decomposizione del problema in sotto-problemi (Decomposizione del Dominio di Calcolo). La simulazione numerica è realizzata adoperando il software ABAQUS. Si sono considerati due casi: uno strato sottile in movimento rispetto ad una sorgente di calore sinusoidale distribuita sulla sua superficie; due strati sottili in moto relativo con una distribuzione sinusoidale di flussi termici all’interfaccia. L’evoluzione nel tempo del campo di temperatura è stato determinato e confrontato con risultati di codici standard realizzati in ABAQUS. Si è mostrato che il metodo proposto porta a differenze contenute (2-4%) dopo brevi periodi di simulazione. Questo lavoro rappresenta il primo passo verso una parallelizzazione del calcolo termoelastico per una progettazione ottimale di freni e frizioni. ABSTRACT The complexity of the analysis of sliding systems as clutches and brakes is related to their 3D geometry, thermo-elastic coupling and transient operating conditions. This all leads to extreme computational costs for the analysis of practical problems. In this work, a new approach is presented for the numerical simulation of thermal problems in frictional systems resting on two main ideas: (i) the use of local reference systems fixed to each domain (Lagrangian Formulation) and (ii) the decomposition of the problem into a collection of smaller problems (Domain Decomposition Technique). The numerical simulation is performed using a commercially available softwer ABAQUS. Two test cases have been considered: a moving single layer with a stationary sinusoidal heat flux applied to its surface; two layers in relative motion with a sinousoidal heat flux at the contact interface, representing non-uniform contact conditions. The time evolution of the temperature field has been predicted and compared with a standard convective-diffusive solutions provided by ABAQUS. It is shown that the two approaches are in excellent agreement giving percentage difference of just 2-4% after relatively short simulation periods (about 0.01 sec). This work is the first step towards a parallelization of the thermo-elastic coupled problem for the optimal design of clutches and brakes for automotive applications. 1. INTRODUCTION In clutches and brakes, the frictional heat generated at the interface, which is proportional to the local contact pressure and relative sliding speed, induces thermoelastic deformations affecting the surface profile and consequently the interface mechanical pressure. Such a system with thermomechanical coupling between surface displacements, frictional heat and pressure can be unstable, depending on the sliding velocity and material properties, leading to a progressive concentration of pressure and temperature over small areas (hot spots) [1, 2]. The analysis of this well documented phenomenon (Barber 1969 [3], Burton 1972 [4], Zagrodzki 1990 [5], Lee&Barber 1993 [6], Decuzzi et al. 2001 [7]), named Frictionally Excited Thermoelastic Instability, relays on a fully coupled analysis where the evolution of the temperature field is firstly defined by solving the Fourier’s equation and subsequently the elastic field is derived by solving the thermoelastic governing equations. Geometry of clutches and brakes can sometimes be considered two- dimensional (plane or axisymmetric), as most of the numerical and analytical analyses presented in the literature assume. In many real systems, however, the problem is fully three dimensional. Also clutches and brakes usually operate with sliding speeds rapidly changing during engagement in a fairly general manner. These features lead to extreme computational costs for the analysis of most of practical problems. Usually the increase in computational time has been tackled by employing more sophisticated processors and optimized codes. More recently, parallel computing or ‘grid computing’ are becoming more and more common (Foster, 2003 [8]). These systems relay on the sharing of the computational work among several ‘slave’ machines connected to one or more ‘master’ machines. Such a very efficient approach can be followed only after ‘parallelizing’ the problem and computational codes. Mutlidisk clutches and brakes are inherently multi-domain systems in that they are made up of multiple friction and metal disks each of them constituting one sole computational sub-domain. This feature has inspired and guided this work. A new approach for the numerical solution of thermal problems in frictional sliding systems is presented, as the first step towards a fully coupled thermomechanical analysis. Two main ideas are described and their feasibility is checked: (i) the use of local reference systems fixed over each domain (disk) (Lagrangian Formulation) and (ii) the decomposition of the problem into a collection of smaller problems (Domain Decomposition Technique). Point (i) relays on the fact that by writing the governing equations of heat conduction with respect to a fixed frame of reference a diffusive-convective (d-c problem) problem is transformed into a purely diffusive problem (d- problem), for which standard and more effective algorithm are available in the literature. Point (ii) relayes on the natural decomposition of multidisk systems in sub-problems where each single disk is treated as an independent element. The solution in each domain (disk) is related to the adjacent domain by satisfying appropriate boundary conditions which implies temperature continuity and energy conservation at the contacting interfaces. In this work, the strategies described at point (i) and (ii) have been employed for estimating the temperature distribution for a multidisk clutch (or brake), being this the first step towards the parallelization of a coupled thermoelastic problem. 2. FORMULATION Two different geometries have been considered (Figure1). In CaseI, a single frictional layer is moving with speed V with respect to a fixed sinusoidal heat flux applied along x (the sole upper layer in Figure1). In CaseII, two layers are included moving with a relative speed V, one made up of frictional material and the other one of metal. At their interface a sinusoidal distribution of heat flux, fixed with respect to the metal layer, is applied. In both cases, the frictional layer has a thickness of 0.67 mm, whilst the metal layer is 2.75 mm thick. Their length along x is of 30 mm. From Figure1a, the frictional layer has a very refined mesh just beneath the sliding interface where the temperature gradient is larger. The bias ratio has been choosen as suggested by Zagrodzki et al. (2001) [9], where the existance of a thermal boundary layer has been pointed out. A uniform mesh is instead considered for the metal layer. Frictional Layer Sliding Interface Metal Layer Figure1: A frictional layer (upper layer) in sliding contact with a metal layer (lower layer). At the frictional interface the mesh is more dense. 2.1 The governing equations Refering to the 2D two-layer system made up of a friction layer (1) and a metal layer (2), two different frames of reference are considered, each fixed to the respective layer, for which the transient Fourier equation takes the form 1 ∂Ti ( x, y; t ) ∇ 2Ti ( x, y; t ) = i=1,2 (1) ki ∂t with the following boundary conditions at the sliding interface ∂T1 ∂T T1 = T2 = To ; K1 + K 2 2 = q f 1 + q f 2 = q f = fVp on Γ c (2) ∂y ∂y where qf, the frictional heat generated at the sliding interface, is partitioned between the two layers being qfi the heat flux entering the layer i at the sliding interface ∂Ti q fi = K i i=1,2 (3) ∂y The initial conditions are Ti ( x, y; t = 0 ) = Ti (x, y ) (4) 0 i=1,2 where Ti = Ti ( x, t ) is the temperature distribution of the layer i, ki is the coefficient of thermal diffusivity and Ki is the coefficient of thermal conductivity. Notice that the Fourier equation for the friction layer written with respect to the stationary frame fixed to the metal layer would have been ∂T2 ( x, y; t ) ∂T ( x, y; t ) k 2 ∇ 2T2 ( x, y; t ) = +V 2 (5) ∂t ∂x where ∂T2 / ∂t is the diffusive term and V ⋅ ∂T2 (x, y; t ) / ∂x is the convective term, which disappears once the equation is written with respect to (x2,y2). 2.2 The Domain Decomposition Algorithm The crucial point in the decomposition is the definition of the boundary conditions linking adjacent layers (sub-domains). The temperature field Ti ( x, y; t ) in each domain can be considered as the superposition of a homogeneous solution Θ i ( x, y; t ) and a non-homogeneous solution Θ i ( x, y; t ) . h n The homogeneous solution is temperature field determined for homogeneous boundary conditions and non - zero initial conditions (no heat penetrating the layer through the sliding interface) that is 1 ∂Θih ∂Θih Θih ( x,0) = Ti o ( x ) ∇2 Θih − = 0; = 0 on Γc (6) k ∂t ∂y The non-homogeneous solution Θin ( x, t ) is the temperature field determined for non-homogeneous boundary conditions and zero initial conditions 1 ∂Θin ∂Θin Θin ( x,0) = 0 (7) ∇2Θin − = 0; = q fi on Γc k ∂t ∂y A finite element discretization of the Fourier’s equation (1) leads to the following matrix form ⋅ (8) Ci Ti + LiTi + Fi q fi = 0 where Ti is the temperature vector which elements are temperatures at nodes; Ci thermal capacity matrix; Li the thermal conductivity matrix; q fi the heat flux vector and Fi the matrix transforming the vector qfi into nodal concentrated heat sources. Applying a central difference algorithm to the system of differential equation (8), it results Ti t +∆t − Ti t T t +∆t − Ti t qt +∆t + qtfi (9) Ci + Li i = − Fi fi ∆t 2 2 where Ti t is the solution at time t and Ti t +∆t is the solution at time t+∆t. The above equation can be rewritten as Ci Li t +∆t Ci Li + ⋅ Ti = + ⋅ Ti − Fi q fi tm ∆t 2 ∆t 2 (10) qtfi+∆t − qtfi being q = tm fi the heat flux at time tm = t + ∆t/2. Introducing auxialiary matrices 2 % C L C L Ai = i − i and Ai = i + i , it finally follows ∆t 2 ∆t 2 % ATi t +∆t = ATi t − Fi q tm (11) i i fi Recalling the definition of homogeneous and non-homogeneous problems given above, the general solution of the problem can be written as Ti t +∆t = Ti h ,t +∆t + Ti n ,t +∆t (12) being Ti h ,t +∆t the solutions of the homogeneous part of (11) % ATi h ,t +∆t = ATi t (13) i and Ti n ,t +∆t the solutions of the non - homogeneous part of (11) ATi n ,t +∆t = Fi q tm i fi (14) Equation (14) can be easily inverted giving Ti n ,t +∆t = Ai−1Fi q tm fi (15) being Ai−1Fi = Bi . Considering the nodes at the interface separately from the subsurface nodes, the non – homogenous temperature field can be rewritten as t +∆t Ti ,1 n Bi ,11 Bi ,12 qtm n = fi (16) Ti ,2 Bi ,21 Bi ,22 0 noticing that the heat fluxes generated at the subsurface nodes are zero. Thus picking up the sole interface temperature field (row 1) Ti ,1,t +∆t = Bi ,11q tm n fi (17) and invoking at the sliding interface, firstly, temperature continuity T1t +∆t = T2t +∆t (18) and, secondly, energy balance it follows: q tm = ( B1,11 − B2,11 ) −1 ( −T1,1,t +∆t + T2,1,t +∆t + B2,11q tm ) f1 h h f (19) The above equation allows for the estimation of the nodal heat fluxes entering body 1; whilst the heat fluxes entering body two are determined by (2). Notice that no iteration is required. The solution of the homogeneous problem Ti h ,t +∆t can be easily derived by using for instance an ABAQUS code (or any other Finite Element, Finite Volume or Finite Difference method) solving a simple purely d- problem. More efficiently, the solution of the non-homogenous problem Ti n ,t +∆t , which is needed to estimate the matrix Bi ,11 and Bi ,21 , can be derived by a Green’s Function approach. The Fourier’s equation for the non-homogenous problem is solved for n independent heat flux vectors, being n the number of surface nodes, 1 0 0 0 1 0 q(1) fi = , q(2) = ,………., q(fin ) = fi (20) . . . 0 0 1 and the resulting temperature field does coincide with the matrix Bi ,11 and Bi ,21 . The general solution of the non-homogenous problem is then derived as t +∆t Ti ,1 n B n = i ,11 qtm fi (21) Ti ,2 Bi ,21 It should be noticed that the matrix Bi ,11 are expected to be diagonally dominated, since the contribution of the concentrated heat flux decreases rapidly as the distance from the point of application increases. 2.3 The Moving Frame of Reference Figure 2 shows two layers, stationary and moving with speed V, in sliding contact. In Figure 2a locations of nodes at the interface at time t+∆t at which the solution has to be determined, are shown. Temperature at t, which is known, constitutes the initial conditions for the sought solution. Figures 2b shows retrospecively node locations at t. Because of the relative motion, a different set of nodes in layer 2 matches the nodes in layer 1 at t than at t+∆t . However, node numbering is only a matter of convention and it can be changed even from step to step. We will move the mesh in part 2 at the beginning of each step so that always the same node numbers in part 2 match the nodes in part 1 through the simulation. In standard finite element codes for solid mechanics, solutions at t and at t+∆t use by default the same node numbering system. In order to realize our idea in a standard code, temperature T2 in part 2 at time t is being shifted in each time step by distance traveled by the part during that time step. The shift is being performed just before staring a new step In the central difference algorithm, Eq.(12), used for discretization of the problem in time, the vector of heat fluxes q tm fi is defined at the mid-time tm = t + ∆t/2. At that time the nodes of layer 2 are halfway between the nodes of layer 1 as illustrated in Figure 2c. Consequently, the heat flux entering layer 2, defined by Eq.(7), is determined at mid-points rather than at the nodes of that layer. Figure2a Node locations at time Figure2b Node locations at time Figure2c Node locations at time t+∆t t t+∆t/2 3. RESULTS Two different cases have been considered: CaseI where a single layer is moving with respect to a sinusoidal heat flux and CaseII where two layers are sliding with relative speed V with a sinusoidal distribution of heat flux at the interface. In CaseI the idea of a moving frame of reference has been employed transforming a d-c problem into a d- problem; whilst in CaseII the idea of a moving frame is used in conjuction with domain decomposition. The solution derived using the presented approach are compared with solutions obtained by modelling a d-c problem in ABAQUS. The Petrov-Galerkin (P-G) algorithm has been employed which is stable for sufficiently small Courant numbers V ⋅ ∆t Cu ≤ 1 where Cu = (23) ∆x Consequently, discretization along x (∆x) and time discretization (∆t) are not independent: as the sliding speed V increases smaller and smaller time steps ∆t are rquired. The material properties employed in the present analysis are listed in Table1. Metal Layer (Steel) Frictional Layer Young’s modulus, E (GPa) 200 0.27 Poisson’s ratio, ν 0.3 0.12 coeff. thermal expansion, α (1/K) 12x10-6 14x10-6 thermal conductivity, K (W/mK) 42 0.241 specific heat, cp (J/Kg K) 452 1610 density, ρ (Kg/m3) 7800 1125 Table 1: Material properties used in the present analysis 3.1 CaseI: single layer model The single layer solution derived using the P-G algorithm obtained in ABAQUS (Td-c) has been compared with the solution obtained employing the idea of a moving frame of reference (Td). In particular the percentage difference DIF% = (Td − f − Td ) / Td − f × 100 ) has been estimated as a function of depth y, material properties and sliding speed V. In Figure3, the parameter DIF% has been plotted along the x-direction at different depths: on the surface (y=0, blue line); subsurface at y=0.1 mm (violet line) and at y=0.3 mm (yellow line). It is clearly shown that the percentage difference is larger at the surface and progressively reduces as the depth increases. Subsequently, it has been estimated the variation of DIF% with the simulation time. Results are shown in Figure4, from which it is clearly observed a reduction of the percentage difference as the number of iterations increase (∆t = 10-5 sec): a 0.4% difference is estimated after 1000 steps (10-2 sec). The difference is even smaller for 1500 time steps for which a 0.35% difference has been measured. It has been also observed that the effect of the material properties is generally negligible, whilst as the sliding speed V increases the percetage difference increases too and a more refined mesh is required. It can be concluded that the idea of a moving frame of reference can be effectively used. 1 6 11 16 21 26 31 36 41 46 0,3 0,2 0,1 0 -0,1 Sup -0,2 sub1 -0,3 sub2 -0,4 Figure3 The percentage difference DIF% is plotted at different depths along x (y=0 azure; y=0.1 blue; y=0.3 mm violet line) 1501 1506 1511 1516 1521 1526 1531 1536 1541 1546 1551 16 14 12 10 8 6 4 2 0 -2 -4 2 inc -6 5 inc -8 10 inc 100 inc -10 500 inc -12 1000 inc -14 -16 Figure4: The percentage difference DIF% is plotted as the simulation time grows (1 increment = ∆t = 10-5 sec) 3.2 CaseII: double layer model The double layer solution derived using a standard d-c model generated in ABAQUS is compared with the solution derived by applying both the moving frame of reference and domain decomposition ideas. The percentage difference DIF% is again plotted along the x direction at different depths for the frictional layer (Figure5) and the metal layer (Figure6). The percentage difference is shown to be smaller than 2% in absolute value just after 0.01 sec of simulation (103 steps) at the surface. Clearly, for the metal layer (Figure6) where there is no migration speed, the percentage difference is in phase regardless of the depth, whilst a phase shift is observed for the frictional layer (Figure5). The percentage difference is negative meaning that the present approach overestimates the temperature given by the d-c solution. 1 6 11 16 21 26 31 36 41 46 51 0 -0,4 -0,8 -1,2 SUPERFICIE -1,6 SUB1 SUB2 -2 Figure5: The percentage difference DIF% for the friction layer after 0.01 sec: (y=0, violet; y=0.1 mm blue; y=0.3 mm azure line) 1 6 11 16 21 26 31 36 41 46 51 0 -0,5 -1 SUPERFICIE -1,5 SUB1 SUB2 -2 -2,5 Figure6 The percentage difference DIF% for the metal layer after 0.01 sec: (y=0, violet; y=0.1 mm blue; y=0.3 mm azure line) As the number of simulations increases, it is estimated a progressive increase of the temperature difference, and in particular at the surface of the frictional layer a percentage difference of 4% is measured after 0.03 seconds of simulation (Figure7). Since the difference between the two approaches is strictly related to the time and space discretization, in order to reduce such a percentage difference several runs for different mesh refinements have been done. In particular, along x, a mesh with of 50, 90 and 100 elements have been considered, at which different time steps have been associated satisfying the Courant number condition. The percentage difference variation for different level of discretizations have been plotted in Figure8, where the yellow line is for 50 elements (the discretization used so far), the blue line is for 90 elements and the azure line is for 100 elements. As the number of elements along x increases, the percentage difference decreases going from 1.63% to 0.5%. It is not convenient going any further in discretization in that an increase of the number of elements that is a decrease of ∆x would lead to smaller and smaller time steps ∆t.. The approach proposed gives again results in good agreement with those estimated by P-G algorithm in ABAQUS. 1 6 11 16 21 26 31 36 41 46 51 4 2 0 -2 -4 -6 -8 -10 SUPERFICIE SUB1 -12 SUB2 -14 Figure7 The percentage difference DIF% for the metal layer after 0.03 sec: (y=0, violet; y=0.1 mm blue; y=0.3 mm azure line) 14 gri91 12 gri51 10 gri100 8 6 4 2 0 0 0,02 0,04 0,06 0,08 0,1 0,12 tempo [sec] Figure8: The maximum percentage difference DIF% on the surface as the number of elements along x increases (Friction Layer) 4. CONCLUSIONS A new procedure has been presented for the thermoelastic design of brakes and clutches based on domain decomposition and moving frame of referances. The procedure has been verified in ABAQUS code, but it could be implemented using any Fnite Element, Finite Volume and Finite Difference method and software. Two cases have been considered: (i) a single layer sliding with respect to a fixed sinuosidal heat flux has been used to investigate the effectivness of transforming a diffusive-convective problem into a purely diffusive problem by means of a reference motion; (ii) a two-layer model which has been used to test the idea of domain decmpositioning and code parallelization. For case (i), it has been demonstrated that the difference between the proposed procedure and the classical P-G model used in ABAQUS is smaller than 0.5% after 1000 iteration steps. For case (ii), by comparing the proposed approach with the standard P-G algorithm, it has been observed a difference between the two even of 12% after 0.01 sec of simulation and this difference tends to increase with time differently from the single layer case. By increasing the discretization along x, however, it has been observed a steady reduction of the percentage difference with the refiniment of the mesh: the precentage difference has reduced to 4 % for 100 elements against the 12% meaured for 50 elements. The described approach can be used for running efficient parallel computation of thermal problems and constitues the starting point for the definition of a more general thermomechanical parallel code for clutches and brakes optimal design. ACKNOWLEDGMENTS Dr. Paolo Decuzzi acknowledges financial support from the Center of Excellence for Computational Mechanics and Raybestos Products Company, and wish to express his sincere gratitude to Dr. J. Macy of Raytech Corporation. Ing. T. Soranno acknowledges financial support from the Center of Excellence for Computational Mechanics. REFERENCES [1] Anderson, A.E. and Knapp, R.A., (1990), Hot spotting in automotive friction systems, Wear, Vol. 138, pp. 319-337; [2] Takezaki K. and Kubota M., (1992), Thermal and Mechanical Damage of Paper Wet Friction Material Induced by Non-Uniform Contact, SAE 922095; [3] Barber, J.R., (1969), Thermoelastic instability in the sliding of conforming solids, Proceedings of the Royal Society London, A Vol. 312, pp. 381-394; [4] Burton, R.A., Nerlikar V. and Kiliparti, S.R., 1972, Thermoelastic instability in seal-like configuration, Wear, Vol. 24, pp. 177-188; [5] Zagrodzki P., (1990), Analysis of thermomechanical phenomena in multidisk clutches and brakes, WEAR, 140 (2): 291-308; [6] Lee K. E and Barber J.R., (1993), The effect of shear traction on frictionally-exited thermoelastic instability, Wear, Vol.160 , pp.237-242. [7] Decuzzi P., Ciavarella M. and Monno G., (2001), Frictionally Excited Thermoelastic Instability in Multi-Disk Clutches and Brakes, Journal of Tribology, OCTOBER 2001, Vol. 123, 865 – 871; [8] Foster I., The Grid: Computing without bounds, Scientific American, April 2003; [9] Zagrodzki P., Lam K.B., Al Bahakali E. , Barber J.R., (2001) Nonlinear transient behavior of sliding system whit frictionally exited thermoelastic instability, ASME Journal of Tribology,Vol. 123, pp.699-708.