THERMAL DESIGN OF SLIDING SYSTEMS BY

                  Paolo Decuzzi1, Tommaso Soranno1 and Przemyslaw Zagrodzki2
  Center of Excellence in Computational Mechanics - Politecnico di Bari, Via Re David 200, 70125,
BARI; ITALY; e-mail:;
  Raytech Composites Inc., Crawfordsville, IN 47933; USA;

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.

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.

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
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.

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

                                              q fi = K i             i=1,2                                         (3)

The initial conditions are

                                 Ti ( x, y; t = 0 ) = Ti (x, y )                                                   (4)

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

                                        ∆t 2          ∆t 2                                                          (10)

              qtfi+∆t − qtfi
being q =
        fi                     the heat flux at time tm = t + ∆t/2. Introducing auxialiary matrices
% 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)

and Ti n ,t +∆t the solutions of the non - homogeneous part of (11)

                                                         ATi n ,t +∆t = Fi q tm
                                                          i                  fi

Equation (14) can be easily inverted giving

                                                  Ti n ,t +∆t = Ai−1Fi q tm

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 
                                      
                                                           Bi ,11     Bi ,12   qtm 
                                      n                =                    
                                     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

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 )
                                                           h            h

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 
                                                                   
                            fi    =   , q(2) =   ,………., q(fin ) =  
                                    .          .                  . 
                                    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 
                                         
                                                              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

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)

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

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
            -4                                                                                                2 inc
            -6                                                                                                5 inc
            -8                                                                                                10 inc
                                                                                                              100 inc
           -10                                                                                                500 inc
           -12                                                                                                1000 inc

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



                                 -1,6                                           SUB1

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


                          -1,5                                                                                SUB1


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

                             1   6      11      16     21   26   31     36   41         46   51
                       -10              SUPERFICIE

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)
                       10              gri100




                             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)


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.

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.

[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.

To top