Modelling a Coupled Thermoelectromechanical Behaviour of Contact Elements via Fractal Surfaces

A three-dimensional coupled thermoelectromechanical model for electrical connectors is here proposed to evaluate local stress and temperature distributions around the contact area of electric connectors under different applied loads. A micromechanical numerical model has been developed bymerging together the contact theory approach, whichmakes use of the so-called roughness parameters obtained from experimental measurements on real contact surfaces, with the topology description of the rough surface via the theory of fractal geometry. Particularly, the variation of asperities has been evaluated via the Weierstrass-Mandelbrot function. In this way the micromechanical model allowed for an upgraded contact algorithm in terms of effective contact area and thermal and electrical contact conductivities. Such an algorithm is subsequently implemented to construct a global model for performing transient thermoelectromechanical analyses without the need of simulating roughness asperities of contact surfaces, so reducing the computational cost. A comparison between numerical and analytical results shows that the adopted procedure is suitable to simulate the transient thermoelectromechanical response of electric connectors.


Introduction
Several engineering applications involve connections where the electrical contact relies on a relatively weak pressure, able to ensure a partial adhesion between two elements.The real contact surfaces are not flat but include many asperities [1]; contacts occur in a number   of small surfaces named aspots, so that the effective contact area   is much smaller than the apparent, macroscopic, contact one   (Figure 1).The effective contact area can also be reduced due, for example, to natural oxidation and the presence of other superficial contaminant films, that can be present in the contact zone.A direct consequence of the fact that the electrical contact is established only through these a-spots is that contact resistivity increases, producing the constriction resistance.The total contact resistance can be generally estimated by means of a statistical approach linking the microscopic contact physics to the meso-scale modelling [2][3][4]; with this approach the total constriction resistance is defined as a function of parameters such as surface roughness, material hardness, normal pressure   , and temperature .
Numerical simulations of electrical contacts can be conducted in agreement with the contact theory [5] which, in general, is able to evaluate only the apparent contact area   , possibly overestimating the electrical and the thermal connection; for example, in [6,7] the electromechanical contact has been evaluated by considering a micro-macro approach where   has been statistically defined.
To improve the knowledge of the effective interconnection between two bodies in contact, a micromechanical numerical model is here developed to simulate the irregular contact surface.Micromechanical models require a homogenized material response which takes into account microstructural heterogeneity [8].Analytical models of composite structures at nanoscale can be found in [9][10][11][12][13][14].
More specifically, in the present work the micromechanical model is characterized by means of three-dimensional contact surfaces obtained via the fractal theory [15,16] which makes use of roughness parameters resulting from experimental measurements on real contact surfaces to simulate the surface asperities (Figure 2).Asperities heights, peaks distribution, and so forth are dependent on the contact surface preparation and the load history; in fact residual strains under elastoplastic cycles can modify the contact surfaces.The fractal theory has been followed to represent the initial asperities for a virgin material and to evaluate the effective electrical contact area; the material is modelled to behave elastoplastically and able to undergo large strains.
Previous three-dimensional models have been developed [17] to evaluate the global effects in terms of stress and temperature distributions around the contact area; the analyses were supported by experimental tests to retrieve the thermoelectromechanical characteristics of interest.
Here an analytical-numerical comparison is additionally carried out to validate this procedure, considering a classic case of a cylinder and a hemisphere in contact under different compression loads.
Such problems may undergo instability [18][19][20], which would require a local investigation of the stress field in large strains [21][22][23][24]; nevertheless this topic is not addressed in the present micro-macro approach, under the assumption that the compressive load does not generate buckling.

The Analytical Model
On the macroscopic scale two contact elements pressed together apparently touch each other on area   that, due to surface roughness, is much larger than area   of effective mechanical contact.In fact, on such a scale contact occurs between the asperities of the two surfaces in a number   of small surfaces per macro-unit area, named a-spots, which both increase with the applied contact force  [25].Current density lines are forced to pass through a-spots, causing a local increase in resistance   .For a single circular a-spot of radius   and sufficiently thick elements,   is analytically known.A more complex question is determining the total constriction resistance corresponding to   and given by many a-spots whose number and radius depend on the compressive load , the surface roughness, and the material properties.
A statistical approach that allows linking the macroscopic model needed for numerical simulation with the physics of microscopic contacts is extensively presented in [2][3][4] and it has been taken as reference for the present work.The statistical characterization of the surface roughness is considered as representative quantities: the mean plane between peaks and valleys of the asperities, the mean absolute asperity slope , and the root mean square (RMS) of the surface roughness   , whereas the contact between the two surfaces is expressed by the mean plane distance .As long as relatively weak contact forces are considered,  and   can be assumed to be constant, so that only  varies with the applied load .
When the load increases, the asperities in contact are crushed and the mean distance  between two bodies decreases (Figure 3).We define  max as the maximum mean distance between the two bodies ( =  max when the two bodies are in contact but the external load is equal to 0).If the external load is applied, the mean distance  is evaluated by where  spot is the spot stiffness in the contact surface.The spot stiffness is evaluated in [4,17] via the relationship where  1 ,  2 are two material parameters.
Starting from these conditions, the following expression is derived for the a-spot mean radius, their number, and the total contact surface In order to relate  to , in agreement with most of the available literature, a plastic deformation of the asperities has been assumed, because a-spots are very small so that weak forces can produce very high local pressures over them, highly exceeding the yield limit [3,4].However it must also be noted that hypothesis of elastic deformations is even reported in literature, considering that asperities may behave plastically at first but when the a-spots enlarge enough, loads are elastically supported.As far as the linear assumption is accepted, the plastic condition allows expressing the contact force as  =     , with   the contact micro-hardness, whereas   is related to the macroscopic apparent pressure   = /  , so that   Considering the Vickers hardness test,   can be found starting by a Vickers micro-hardness measure where  V is the mean indentation diagonal; that is, in a similar way Considering the plurality of a-spots and the proximity among them, each current tube can expand far from its contact point to a radius  (Figure 1).The total contact resistance   can be defined as the inverse of the contact conductance   = 1/ℎ  that can be evaluated as in which   is the electrical conductivity.Moreover, as the contact surfaces generally present some level of oxidation and often also stray deposits of insulating substances, a film resistance has been added to the constriction resistance from (7).As it strongly depends on the level of cleanness of the surfaces and on the used metal, values taken from literature have been considered.

The Thermoelectromechanical Model
Under an electrical flux, bodies' configurations are subjected to modifications caused by electrical energy dissipation.Correspondingly, a coupled thermoelectromechanical model has been accounted for.
The electric field is governed by Maxwell's equation where Γ is the surface area, Ω is the volume, J is the electrical current density vector (per unit area), and   is the internal volumetric current source (per unit volume).J can be defined via Ohm's law by considering the current density as a function of the electrical conductivity matrix   =   (), with  being temperature, and the gradient of the electrical potential where x is the position vector in the current configuration and  is the electrical potential.The thermal energy generated by the electrical current can be defined by Joule's law that evaluates the dissipated electrical energy   as The energy released in the form of internal heat is where  is the heat energy generated during dissipation and  is the energy conversion factor.Thermal variations within a body involve a new mechanical configuration.The stiffness matrix for a coupled thermoelectromechanical model results as where k  , k  , k  are the mechanical, thermal, and electrical parts, respectively, whereas k  , k  , k  , k  are the coupling terms.At present the piezoelectric behaviour is not considered and the constitutive relation can be defined as where the stress tensor T is dependent on the logarithmic elastic strain E  , obtained by subtracting the thermal strain tensor E  to total logarithmic strain E, and the constitutive tensor C. Considering an elastoplastic material, the constitutive tensor C is dependent on both the position vector in the reference configuration X and temperature .

The Contact Algorithm
In order to simulate the thermoelectric interconnection when two bodies are touching, the contact theory has been applied [5].
Considering two bodies, Ω 1 and Ω 2 (Figure 4), for example, representative of connection components, two surfaces can be identified, Γ 1 (with Γ 1 ∈ Ω 1 , named slave) and Γ 2 (with Γ 2 ∈ Ω 2 named master), where contact is possible.The closed contact condition is achieved and the two bodies are in contact if the contact surface Γ  = Γ 1 ∩ Γ 2 ̸ = 0. Contact is defined when the fundamental conditions are set: (i) Nonpenetration conditions: where u  (with  equal to , master, and , slave) are the displacement vectors, X  are the position vectors in the reference configuration,  is the gap function (the distance between two points in contact), and n is the normal vector (Figure 4).(ii) Action-reaction conditions: where t  are the stress vectors.(iii) Kuhn-Tucker conditions: where   is the normal pressure:   = t ⋅ n.Generally, friction effects between two bodies in contact are dependent on the roughness of the contact surface; in our case, if one considers the adopted material asperity model, friction is a direct consequence of the model itself and it comes from the transversal asperity connections, so that a normal condition only has been considered (Figure 5).
The normal contact is characterized by a normal pressure, defined based on the penalty method; that is, where   is the gap function along the normal direction and   is the penalty coefficient.
In the developed three-dimensional numerical model here presented, characterized by tetrahedral elements, master and slave surfaces have been defined on the contact faces of the elements (Figure 6).The closed contact condition has been considered in the contact pair, defined through a slave node and a master point.On these surfaces only the gap function  is evaluated.
The master point x  in a contact pair is chosen as the point with the minimum distance from the slave node x  (minimal distance rule); generally it does not coincide with a master node but with a generic point belonging to surface Γ 2 .
The minimum distance between slave nodes and master points can be defined by considering that the master point is obtained by the orthogonal projection of the slave node onto the master surface (see Figure 7), where the normal direction n is defined by and a  are the tangent vectors to the master surface in x  (Figure 4) [26].
A "master element" is characterized by the elements faces, and the tangent vector at a generic point x  can be consequently defined taking into account the shape functions derivatives in the master element.A master point can generally be represented by referring to two different reference systems: a global reference system (g.r.s.), x  = x  ( 1 ,  2 ,  3 ), and a local curvilinear reference system (l.r.s.) belonging to the master surface, x  = x  () = x  ( 1 ,  2 ).The shape Master Slave  functions of a master element related to point x  can be defined as in which  is the number of element nodes,   is the shape function for node , and x  = (  1 ,   2 ,   3 ) are the coordinates of the master element nodes referred to as the g.r.s.

Thermoelectromechanical Conditions
The contact surfaces shown before must additionally transfer thermal and electrical potential fluxes.
The heat flux per unit area  defined between a contact pair has been assumed as where the thermal conductivity ℎ is dependent on the temperature   (with  equal to , master, and , slave) and the normal contact pressure   .The electric flux density  has been defined as with Δ being the electric potential between master and slave surfaces and ℎ  the gap of electrical conductance or electrical conductivity in the contact zone, which is dependent on the average temperature in the contact zone  avg =   −   and on the normal pressure   .By assuming that the effective contact area   can be represented as a circular area, the contact resistance   can be obtained following Yovanovich's resistance relation [3,4] and considering that, in general, the contact size is much larger than the mean free path of electrons; correspondingly the contact resistance can be obtained as where   is the resistivity of the contact surface  and  is the contact diameter.Experimental tests have been conducted on samples in aluminium alloy with assumed resistivity equal to 3.4⋅10 −8 Ωm [27].A thermoelectrical conductance due to the interstitial gas layers [16] has been here neglected.
Modelling and Simulation in Engineering

Roughness Tests
To evaluate the electrical resistance between two contact bodies under different compression load levels, experimental tests have been conducted with an aluminium cylinder and a hemisphere (Figure 8(a)).In this simplified case the apparent contact surface   remains circular and can be also analytically estimated.The effective contact area, that is, the zone where the thermoelectric connections are defined, is dependent on contact surface asperities; the surface roughness in the contact area has consequently been characterized.
For each contact surface of cylinders and hemispheres, 20 roughness profiles have been evaluated, in 20 different positions (the measurements have been taken in the neighbourhood of the contact zone between the cylinder and the hemisphere).
Profile lengths are equal to 4.0 mm considering a cut-off length of 0.8 mm.The sampling measures have been taken every 2 m.
Typical examples of roughness profiles evaluated in different positions of the sample are shown in Figure 9.
The RMS surface roughness   in the contact zone has been defined as where  , is the RMS value in the contact surface .The RMS for the two contact surfaces has been obtained based on the roughness tests, reaching an average value of   = 0.105 m.
The mean absolute asperity slope  has been evaluated equal to 13%.

The Fractal Surface
The roughness parameters obtained via the roughness tests are dependent on the resolution of available measuring instruments because of the multiscale character of the roughness and its nonstationary features.The fractal theory can be adopted [27] to evaluate the invariant roughness parameters for all scale levels.At this purpose the numerical representations of the effective contact area have been conducted in agreement with the fractal approach.A micromechanical model has been developed to evaluate where the thermoelectrical connection occurs, taking into account the material asperity in the contact surface.
So, starting from the theory of the fractal geometry, the topology description of the rough surface has been possible.The asperity height (, ) in a plane (, ) has been calculated with the modified Weierstrass-Mandelbrot (WM) function [16]  (, ) where the mechanical characteristics are , sample length, , fractal roughness, , fractal dimension; considering a twodimensional surface, the limits of parameter  are 2 <  < 3.  represents the extent of space occupied by the rough surface;  is the scaling parameter (equal to 1.5);  is the number of superimposed ridges used in constructing the surface profile;  is a frequency index; and  , is a random phase with interval [0-2].
The only unknown variables in (24) are  and , which can be experimentally estimated.

Fractal Surface Generation.
To define the fractal surface two different codes have been developed: the first written in C++ language is able to define the WM coefficients after reading different 1D roughness profiles experimentally obtained, with roughness surface tests at different positions in the same area of the analysed contact surface.This code evaluates the surface conditions with different parameters such as maximum and minimum height peak present in the profiles, average height of asperities, profile slope, and asperity density distribution.
These surface characteristics are necessary to obtain an equivalent fractal surface developed via a second code, written in Fortran 90, where the three-dimensional WM surface (Figure 10) is defined.
After adopting such a method for generating a fractal surface, a solid geometry has been then defined and afterwards a numerical model has been created (Figures 10(b) and 10(c)).The irregular surface has been meshed via a triangular discretization.
As shown in (24), the fractal surface can be defined provided that a "random phase" parameter, representative of an independent random variable, is known.Hence, a specific Fortran function has been called to generate randomly a phase parameter as shown in Box 1.

Mechanical Characteristics
An aluminium alloy has been tested.The mechanical characteristics in terms of elastic modulus  and yield stress have been obtained at room temperature via tensile tests.A thermal variation of such parameters has been accounted for in agreement with Eurocode 9 prescriptions [28].When temperature increases, the elastic modulus and yield stress decrease (Figure 11(a)).
For the sake of simplicity, a bilinear temperature-dependent elastoplastic relation has been assumed for the numerical analyses (Figure 11(b)).
The hardness parameters have been experimentally defined by considering the Vickers measure; namely, a value  V = 315 MPa at room temperature was obtained for the tested alloy.

Numerical Analyses
Two different models have been carried out in the following: a micromechanical model, to evaluate the effective contact area between two contact surfaces (and consequently upgrading the contact algorithm), and a global model where the thermoelectromechanical characteristics have been modified in agreement with the micromechanical results.

The Micromechanical Model.
To evaluate the effective contact zone, where the real electrothermomechanical connections are established, fractal interfaces associated with the contact algorithm have been considered.In the real contact surfaces the irregular asperities generate stress concentrations that can be estimated through the fractal surfaces themselves.
The WM equation used to generate the roughness surface considers a random phase parameter that allows the creation of different fractal surfaces with the same mechanical characteristics (derived by experimental tests).The contact surfaces used in the numerical models have been built by taking an average value of 10 different surfaces obtained via the WM equation to simulate the cylinder contact surface and 10 for the hemispheric one.The maximum in-plane extension of the fractal surface has been assumed 1.0 × 1.0 mm 2 and the average asperity height is equal to about 0.2 m.
The independent variables  and  have been calculated by taking into account that maximum, minimum, and average asperity heights must be the same, when comparing fractal surfaces and experimental measurements.
The micro-model of the cylinder in contact with the hemisphere consists of 123610 nodes and 513833 tetrahedral elements; circular macro asperities have been created to simulate the effective contact surface, by means of surface preparation techniques.
During compression of the two electrical joint surfaces, the different asperity heights come into contact at different times and the effective contact area depends on the external load as shown in Figure 12.
Figures 12(b) and 12(c) show that contact occurs first in the circular macro asperities and, subsequently, contact involves the internal part of the surface.As shown in Figure 13(a), the joint stiffness increases when surface asperities enter into contact but stiffness decreases at higher temperatures.The relation in terms of contact stiffness, temperature, and gap function has been used to evaluate the mechanical penalty coefficient in the global model.
When two asperities are touching, the transferred mechanical stress rapidly increases above the yield stress (Figure 13(b)).
As already stated, when compression occurs, at the contact surface the effective contact area   changes with temperature, but as reported also in [27], the relation between the applied load and the effective contact area results to be linear, with decreasing slope at increasing temperatures (see Figure 14).At thermal melting point (about 550 ∘ C) the curve becomes almost horizontal.
If the effective contact area is equal to zero the two surfaces are not in contact and the contact resistance   is infinite.In agreement with (22),   decreases as   increases (Figure 15(a)).
As evidenced in Figure 14, by varying normal pressure and temperature, the effective contact area   changes, which allows obtaining   law in terms of pressure and temperature variations (Figure 15(b)).
Hence, the contact algorithm adopted in the macro scale model has been modified in light of the thermoelectromechanical relations coming from the micromechanical model just illustrated.In the global one, roughness contact surfaces have not been represented to reduce the computational cost, and the contact area itself gives the apparent contact area   related to the pair, once the effective contact area   is evaluated following the relationship shown in Figure 16.
Consequently, the contact resistance   has been defined for each contact pair.

Global Model and Comparison of the Results.
The global model has been carried out as illustrated in Guarnieri et al. [29], where a hemisphere and a cylinder are put in contact at different compressive load levels and different electrical potentials.These types of shapes ensure a closed solution for the analytical approach, because the resulting contact surface is always circular.The contact characteristics have been compared considering the analytical formulation and the previously described fractal characterization.
The schematic view of contact between the two surfaces is shown in Figure 8, where the maximum diameter for the  samples  is assumed to be equal to 90 mm and the material is conceived to behave thermoelastoplastically.
A current density  = 1.0 A/mm 2 is applied to the hemisphere after the load application, at an initial temperature of 20 ∘ C. The three-dimensional numerical model consists of 730000 tetrahedral elements (Figure 17(a)) and linear shape functions.Due to symmetry, only 1/4 of the geometry has been considered.
After the application of the compressive load, the sphere touches the cylinder in the contact zone and the contact pressure increases following the elastoplastic relationship.A permanent indentation occurs if the elastic field is exceeded as shown in Figure 17(b) [30].
The geometric configuration of the sample allowed for obtaining a circular apparent contact surface at different load levels (see Figure 18), which is proved in the following to correspond to that analytically evaluated [1].
In the contact algorithm the electrical conductance ℎ  at the contact pair  has been obtained via the apparent contact area  , and the effective contact area  , , based on Yovanovich's formulation [2,4] where   is the average electrical conductance of the two bodies in contact ) . ( has been calculated, via the micromechanical models, as a function of contact pressure and temperature.In the analytical and numerical models   has been assumed to be equal to 2.9 × 10 4 S/mm.Similarly, the thermal conductance has been evaluated and the thermal conductivity ℎ has been assumed to be equal to 0.237 W/(K⋅mm).This assumption allowed evaluating the variation of the electrical and the thermal conductance at different contact positions, under different contact pressures.Analytical and numerical results have been compared by considering the contact surface at 20 ∘ C.
A comparison in terms of real contact area   and ℎ V is reported in Table 1 and in Figure 19.
As shown, the numerical and the analytic results for   are pretty similar at different load levels, which confirms that fractal surfaces are good to evaluate the real contact surface in these types of problems.
A slight discrepancy is encountered in the first load steps: such difference can be explained by the fact that the analytical method is based on a statistical approach and the correspondent results are averaged, whereas the fractal approach takes into account the effective contact area found through the micromechanical model during the overall compressive load histories; that is, the analytical method is not able to locally catch the real contact area at low load values.Indeed, the deviation between the results gradually decreases when the load increases, due to the fact that the contact spots defined by the fractal model are comparable to those obtained via the statistical approach.In fact, the resulting contact areas must finally coincide, getting closer to the real one.
According to the micromechanical model, the real contact area   is consequently a function dependent on load and temperature,   =   (, ).This allows evaluating the      The global model is also able to evaluate the displacement variations and the stress fields (Figures 20(a) and 20(b)) when an electrical current is applied at different compressive loads.
The contact discontinuity generates a sudden deviation in terms of electrical potential and temperature (Figure 21); when the contact area increases, at higher loads, the peak of temperature and electrical potential decreases, being locally higher than the thermal and the electrical contact conductivities.

Conclusions
A three-dimensional coupled thermoelectromechanical model for electrical connectors has been here proposed to evaluate local stress and temperature distributions around the contact area of the connectors under different applied loads.The electrical resistance arising in the contact area varies with the effective contact area   , where thermal and electrical contacts are established.It can be assumed that   depends on the surface roughness and the contact pressure.
A micromechanical numerical model has been developed by merging together the contact theory approach [5], which makes use of roughness parameters obtained from experimental measurements on real contact surfaces, with the topology description of the rough surface via the theory of fractal geometry.Particularly, the asperities variation has been evaluated via the Weierstrass-Mandelbrot function [16].
In this way the micromechanical model has allowed for obtaining an upgraded contact algorithm in terms of effective contact area   , as well as thermal and electrical contact conductivities, qualifying a small region of the total contact zone.
Such an algorithm has been subsequently implemented in a global model for performing transient thermoelectromechanical analyses without the need of simulating roughness asperities of contact surfaces, so reducing the computational costs.
A comparison between numerical and analytical results proved that the adopted procedure is suitable to simulate the transient thermoelectromechanical response of electric connectors.Stress tensor t:

Nomenclature
Contact stress vector u: Displacement vector x: Position vector in actual configuration

Figure 1 :
Figure 1: Schematic view of a-spot area.

Figure 3 :
Figure 3: Connection between two contact surfaces.

Figure 5 :
Figure 5: Normal contact definition in the fractal surface.

Figure 6 :
Figure 6: Master and slave surfaces in finite elements.

Figure 7 :
Figure 7: Definition of the normal direction in a point and location of a master point.

Figure 10 :
Figure 10: Example of a generated fractal surface (a); solid geometry (b); mesh of the numerical model (c).

Figure 12 :Figure 13 :Figure 14 :
Figure 12: Effective contact surface at 79% of the applied load (a); at 89% of the applied load (b); at 100% of the applied load (c).

Figure 15 :
Figure 15: Contact resistance   versus a-dimensional contact area (a);   versus normal pressure (b).

Figure 16 :
Figure 16: Relation between temperature, contact area, and normal pressure.

Figure 19 :
Figure 19: Apparent and real contact area at different load levels (20 ∘ C) (a); comparison between analytical and numerical electrical conductivities (b).

Figure 21 :
Figure 21: Comparison between electrical potential versus position at different loads levels (a); comparison between temperature versus position at different loads levels (b).

a
: Tangent vector along  direction   : Apparent contact area   : Effective/real contact area   : a-spot radius : Contact point radius C: Constitutive tensor   : Material parameters ( = 1, 2) : Sample diameter : Fractal dimension  max : Maximum mean distance between two bodies  V : Mean indentation diagonal E: Total logarithmic strain E  : Electrical potential gradient E  : Elastic logarithmic strain E  : Thermal logarithmic strain : Gapfunction : Fractalroughness   : Gap function along normal direction ℎ: Thermalconductivity ℎ  : Electrical contact conductance   : Contact microhardness  V : Vickers microhardness : Current density J: Electrical current density k  : Stiffness matrix parameters   : Penalty coefficient : Sample length : Meanasperityslope : Number of superimposed ridges used in constructing the surface profile n: Normalvector   : Shape function   : Numbersofa-spots : Compressiveload   : Contact pressure   : Dissipated energy : Heatfluxperunitarea : Heat energy release   : Internal volumetric current source   : Contact resistance : Temperature T:

Table 1 :
Comparison between analytical and numerical results.