Evaluation of contact force distribution along a curve, based on measured electric potentials

The paper presents the first results of an exploratory research work, aiming the experimental evaluation of the mechanical contact between conforming surfaces of metallic bodies. The proposed procedure—that is based on the measurement of electric potentials—is able to determine the actual contact pattern and estimate the force distribution on the opposing surfaces. In the present paper, the contact properties along a curve are evaluated, but the basic principles supporting the solution of the more general surface contact problem are also laid down. The goal of this research work is to provide a reliable basis for the investigation of contact-related effects on the modal properties of structures and machines. The distribution of electric contact resistivity can also be assessed using this method, opening the way to the refined characterization of electric switches and connectors. During the application of the method, measured electric voltage values are processed by a neural network trained on a data set generated by finite element software. The resulting electric contact conductivity distribution is used to determine the contact force distribution, taking into account the equilibrium of forces and moments. The results are also checked by available alternative experimental methods.

the efficient solution of the equations of motion poses challenges even if rigid multi-body systems are considered [4]. Moreover, one cannot ignore the fact that the contact phenomenon is inherently nonlinear. Thus, the contact stiffness and damping must be characterized by parameters that physically differ from the parameters used in linear models. Based on pre-determined equivalent requirements, modal stiffness and damping can be defined, but these parameters cannot entirely reflect the possible complex behaviour related to the contacts.
While the stiffness can be modelled fairly well based on geometrical and material properties, the prediction of the damping of materials and structures is difficult, since it is related to several physical phenomena [5]. In engineering structures or machines, the material damping of metals is usually negligible compared to the socalled structural damping [6]. Structural damping is strongly related to the relative motion between structural elements; thus, the corresponding models are mostly based on the well-known Coulomb friction model [7]. Although the Coulomb friction force is primarily related to the normal force between the bodies, the distribution of the force may have a large influence on the measurable natural frequencies and damping properties of the mechanical system [8,9]. In view of this problem, we propose a new experimental method for the evaluation of the mechanical contact between metallic machine or structure parts, by the estimation of the contact force distribution.
Several models were introduced in the scientific literature in order to describe the contact of solid bodies and the unavoidable friction between them. Elastic contact problems were rigorously solved by Hertz [10], and his results still provide a reliable basis for the estimation of contact stress if smooth surfaces touch each other. However, macroscopically smooth surfaces are rough on finer scales and-as it was pointed out by Bowden and Tabor [11]-the real contact area is only a small fraction of the apparent area. As a consequence, the normal stress does not characterize the contact of uneven surfaces properly, because the stress is calculated by dividing the force by the nominal area of the surface where the force is acting. Thus, the normal force is one of the most important variables in several contact models.
Although the surface roughness influences the contact properties, it is rather difficult to characterize the surface quality properly. Single numbers, like the average roughness R a or peak to valley roughness R T , are often used for the quality control of manufacturing processes, but they are less relevant if one wants to describe the small-scale interactions between the surfaces [12]. Thus, several researchers used multi-scale or probabilistic models. The multi-scale or fractal models of contact can be traced back to Archard [13], while Greenwood and Williamson [14] developed a very successful model that is based on the height distribution of asperities. There are several extensions of the model that revise or relax certain assumptions, e.g. [15,16]. Nevertheless, Greenwood and Williamson pointed out that the real local contact areas are sometimes evenly distributed on the apparent contact area, while clustering can be observed in other cases. Thus, the evaluation of the contact pattern or contact patch has significance in the characterization of the contact.
It is worth mentioning that the interaction of contacting bodies can be described on an atomic or molecular level, too (see, e.g. [17,18]). However, for the application of the contact theories in the engineering practice, one must choose a proper macroscopic scale for the characterization of the contact. On a large scale, even the application of machinist's blue (soft liquid or hard paint) can reveal the contacting areas or highlight the points of excessive wear-but only after the separation of surfaces.
As the computational capacity of computers reached a sufficiently high level, researchers started to use numerical simulation techniques for the calculation of the dependence of contact properties on the normal force, based on measured or generated surface topography and the theories of elasticity and plasticity [19][20][21].
However, two surfaces cannot be aligned with the accuracy of the surface topology measurements in the engineering practice. Even small differences in position lead to inaccuracies in the theoretical predictions about the contact force distribution. It means that the information about the contact pattern loses its reliability once the bodies are no longer in contact. This is especially true if the damping properties are to be characterized, since the friction force-that is related to the real contact area-was found to increase in time [22]. Thus, in case of well-machined surfaces-e.g. the mating surfaces of the tool-holder and the main spindle in milling machines [23]-a more refined method needs to be developed that can provide information about the actual contact, i.e. when the parts are not separated.
The authors are aware of only a few experimental methods for the check of contact calculations while the bodies are still in contact, but each of these methods has serious deficiencies. A large group of methods, like the photoelastic coating technique [24] or the application of strain gauges [25], is based on the measurement of the strain field on a free surface not far from the contact area. If the stiffness of the contacting bodies is large and the magnitude of the applied force is given-which is typical in the case of machine parts-the achievable strains are rather small, leading to inaccurate results (see Sect. 5.1). Moreover, the application of these methods is time consuming, and a relatively large area must be covered with strain gauges or special coatings, which may be impossible in practice. Another measurement method is the application of the X-ray technique [26]. Unfortunately, this method is capable only to estimate the area of real contact, and it does not provide information about the contact forces. Moreover, the application of this method may require extremely high power X-ray tubes for the analysis of thick metal bodies like machine parts. To overcome these difficulties, we propose a novel method that is based on the measurement of electric contact resistance.
In Sect. 2, the main steps of the method are described. Since the proper determination of the contact resistances is of high importance from the point of view of the method, Sect. 3 deals with the theoretical analysis of this electrical problem. Practical considerations about the application of the method are presented in Sect. 4, where the developed measurement procedure is described, focusing on the 1D contact evaluation between two tubes. In Sect. 5, the method is tested in special cases to demonstrate its reliability. Then, the force distribution along a curve is determined for various values of the resultant compressing force. Finally, the results are summarized in Sect. 6.

Description of the method
Our method exploits the monotonic dependence of electrical contact resistance on the normal force. This relationship was examined by several authors (see, e.g. [27,28]). The connection between these two physical quantities is often utilized for the evaluation of the contact resistance, based on mechanical properties. Two typical examples are the design of electric switches and arc-welding processes [29,30]. However, the inverse problem-evaluation of mechanical contact based on measured electrical resistance-is not examined in the scientific literature to the best knowledge of the authors. A possible reason is the large statistical deviation of resistance measurements. The electrical contact resistance varies in time due to physical and chemical reactions between the contacting bodies, making the repeatability of the resistance measurements poor. Due to the-usually insulating-oxide layer on the surface of metals, electric conductivity is assured only in domains where the oxide layer breaks. Thus, the conducting area is only a small part of the load-bearing area, presumably at small spots dispersed in the area of mechanical contact [28].
On fine scales (typically in the range of micro-and nanometers in case of machine parts), the electric resistance depends on the stochastic pattern of cracks on the oxide layer. On larger scales (in the range of millimeters), however, one can take into account that the major cause of the contact resistance is that the current lines are constricted to the areas of electrical contact [29], as it is illustrated in Fig. 1. Since the constriction of current lines over the mechanical contact domains has a major contribution to the contact resistance, we assume that the macroscopic pattern of contact resistance coincides with the pattern of the mechanical contact. This is especially true if the normal force is relatively large [28].
The proposed measurement procedure-about which a utility model application has been accepted [31]-is based on the following steps:  • First, the nominal contact area is partitioned into a discrete number of smaller regions, as illustrated in Fig. 2.
The physical characteristics of the contact (distribution of electrical resistance and contact force) can be defined according to the resolution of the partition. • To determine the distribution of contact resistance, a constant current is generated between pairs of points on opposite sides of the contact surface, using electrodes touching the contacting bodies. Meanwhile, the electric voltage is detected between neighbouring points using another pair of electrodes. This kind of voltage measurement can be repeated at several point-pairs and several different allocations of the current source's electrodes. In the most favourable case, voltage can be measured at each of the pre-defined regions, as illustrated in Fig. 2. This situation corresponds to the contact of two thin sheets of metal. It is more common that only a part of the regions can be reached, e.g. measurements can be conducted only at the perimeter of the contact curve-see, for example, points A, B, C and D in Fig. 2.
Since the generated electric current is kept constant and the resistance of the voltmeter is large, the contact resistance of the electrodes does not influence the results. Thus, the measured voltage values provide information about the properties of the mechanical contact. The detailed explanation of voltage measurement is presented in Sect. 4.2, while the developed experimental device is described in Sect. 4.3. • Performing the aforementioned measurement for several pairs of points, a voltage distribution-and the elements of the impedance matrix, the so-called effective resistance distribution (see Sect. 4.2)-can be obtained. Since the effective resistance is simply the ratio of the measured voltage and current values, it depends on the path of the current lines in the material and the specific resistivity of the material. Thus, it is not equal with the contact resistance. However, as it is derived in Sect. 3 using the theory of linear electric circuits, the contact resistance can be determined if a sufficient number of voltage measurements is accomplished.
To determine the contact resistance distribution, the measured voltage and current values are filtered (see Sect. 4.4) and can be processed by two alternative methods. The calculation of the so-called second discrete derivatives of the resistances (Sect. 4.5) may provide a quick way to evaluate the contact pattern along a curve. The other approach is based on the finite element analysis of the contacting bodies (Sect. 4.6) and the training of a neural network (Sect. 4.7). Supplementing these numerical tools with the results of preliminary experiments, and so providing an estimate about the order of magnitude of typical contact resistances, the reliability of the method can be checked before building a full test rig. • As is described in Sect. 4.8, there are several alternative approaches for the estimation of the contact force distribution. Even if the force distribution is not determined, the contact pattern can be characterized by the distribution of the contact resistance. The monitoring of this pattern may reveal the imperfections in the contact that can be useful during the assembly of machine parts. The application of this approach is illustrated in Sect. 5.1. The relation between the contact force and contact resistance is rather difficult to establish, due to the sensitivity of the contact resistance on various-practically uncontrollable-effects, like the oxidation of the surfaces and the local plastic deformation of materials. Thus, the same surface quality and cleanliness must be assured along the contact surface, but this requirement is usually fulfilled in case of industrial-grade machine parts. According to our assumption, the repeatability of the force distribution measurement can be achieved by compensating the temporal variation of the resistance. The elaborated compensation and calibration procedure exploits the equilibrium equations of forces and moments and a special Fourier series expansion. According to the proposed method, a constant known current is led through two metallic bodies that are in contact with each other. The current flows into one of the bodies, while the same current is led out from the other body. As a consequence, a potential field develops that satisfies Laplace's equation and the boundary conditions. The goal is the characterization of the contact resistance on a discrete partition of the contact surface between the bodies, based on the measured voltages.
To analyse the electrical system comprising the contacting bodies and the attached current generator, we set up a lumped parameter model such that the electrical contact is represented by the finite number n 0 of contact resistors with resistances R ci , i 1, 2, . . . , n 0 . Our goal is the determination of these contact resistances.
The 2n 0 terminal points of the contact resistors define the nodes of an electrical network. As an example, see Fig. 3, where the branches corresponding to contact resistors are illustrated by dashed lines for n 0 5. These branches will be referred to as contact branches.
Besides the contact resistances, internal (material) resistances can also be defined between each (noncontact) pairs of nodes. These resistances (see the branches drawn by solid lines in Fig. 3) are considered to be known, since their values can be accurately determined by the finite element analysis of the monolithic metallic bodies.
We restrict ourselves to the cases when both the current generator terminals and the voltage measurement electrodes are attached to so-called contact pairs-these are the nodes that are directly connected by a contact resistor. One of the possible scenarios is illustrated in Fig. 3.
Based on Kirchhoff's circuit laws [32], it is possible to set up the system of governing equations for the obtained network. There are altogether n 2n 0 nodes, so the number of independent nodal equations is e n n − 1. On the left hand sides of these equations the sums of unknown currents appear that enter (− sign) or leave (+ sign) the nodes. The right hand sides are typically zeros, except for the nodal equations of those nodes where the current generator is connected to the network.
The number of branches within one of the bodies is b 0 n 0 (n 0 − 1)/2, that is, there are n 0 (n 0 − 1) branches within the two bodies, each representing the material resistance between two nodes. With n 0 contact branches, the total number of branches is b n 2 0 . It can be shown [32] that the number of independent loop equations is e l b − n + 1 n 2 0 − 2n 0 + 1. Thus, there are altogether e n + e l n 2 independent equations describing the electrical circuit. Each voltage measurement between the nodes of a contact pair provides an additional equation in the form U i R ci I ci .
There are n 0 unknown contact resistances and b unknown branch currents: these are the contact currents I ck , k 1, . . . , n 0 and the currents within the bodies I j , j 1, . . . , b − n 0 . Thus, if the voltage is measured at each of the n 0 contact pairs, the number of equations and unknowns is the same: b + n 0 .
To examine the solvability of the problem, more information is needed about the structure of the equations. Independently of the number of nodes, it is always possible to set up the set of equations based on an algorithm that is detailed in the Appendix (see Sect. 7).

The structure of the equations
As is shown in the Appendix (Sect. 7), there are 2n 0 −1 independent nodal equations with n 0 unknown currents in each of the nodal equations. In the example of Fig. 3, The symmetry of the equations is apparent: adding 2n 0 10 to the indices of non-contact currents and changing their signs in the first four equations, one obtains the last four equations. If one assumes that the upper and lower bodies can be characterized by the same internal resistances, the symmetry of the setting implies that the currents are of opposite sign in the corresponding branches. In this case, the number of independent nodal equations drops to n 0 . There are n 0 − 1 loop equations where the contact resistances appear. In the case of the example, these equations assume the form Each loop is obtained by closing a branch between two nodes of the so-called spanning tree that is shown by red line sections in Fig. 3. Since only one of the contact branches is included in the spanning tree (with resistance R c5 in the example), contact resistances occur in the loop equations only if the closing branch is one of the contact branches. Thus, only two contact resistances appear in each of these equations: the contact resistance of the tree contact branch (R c5 ) and one of the remaining contact resistances. There are additional n 2 0 −3n 0 +2 loop equations that contain only known resistances and unknown currents. In the case n 0 5: In addition to that, there are six more similar equations, where all the indices are larger by 10. Just as in the case of the nodal equations, the symmetry of the contacting bodies implies that these b 0 pieces of additional equations-where the currents (I 11 . . . I 20 ) of the lower body appear-can be dropped. It means that if the symmetry condition is fulfilled, there are only n 2 0 − b 0 b 0 + n 0 independent Kirchhoff's equations for the b 0 + 2n 0 unknowns. Thus, with n 0 measured voltages, the number of equations and unknowns is still the same.

Solvability and conditionality of the inverse problem
Equations (2) and (7) are linear, while Eqs. (3)-(6) comprise nonlinear terms: products of two unknowns-a contact current and a contact resistance-appear here.
If all the contact voltages are measured, one obtains n 0 additional equations in the form By substituting these values into Eqs. (3)-(6), the nonlinear terms disappear. The contact currents can be expressed by the other currents using n 0 nodal equations. Since the contact currents can be eliminated from all the other equations, these expressions can be discarded temporarily.
The remaining n 0 − 1 nodal equations and the n 2 0 − 2n 0 + 1 loop equations form an independent set of linear equations of size n 2 0 − n 0 for the same number of unknowns. This set of equations has a unique solution under quite general circumstances. As a result, the non-contact branch currents can be expressed by the measured voltages and known resistances. After that, the contact currents can be obtained using the previously discarded n 0 nodal equations, and the contact resistances can be determined using (8). One can show by means of a similar train of thoughts that the reduced set of equations-obtained by assuming symmetry-can be solved, as well.
The sensitivity of the results to variations or imperfections in the measured data can be characterized by the conditionality of the coefficient matrix of the derived equations. Since the contact resistances and the internal resistances can assume values in a broad interval, it is difficult to make general mathematical statements about the conditionality of the problem. However, our first experiments indicated that the conditionality improves if the resultant of the contact resistances is less than 30 times the maximal material resistance within a single body.
In certain realistic situations, a rather fine resolution might be necessary for the proper evaluation of the measurements, leading to very large sets of equations. Thus, to be able to solve such complex problems, we developed a neural network-based evaluation method that is trained by data generated by the finite element analysis of the contacting bodies.
Using this neural network, the expected reliability of the method can be checked in practical cases, when relevant information is available about the contacting bodies. Based on a simple, preliminary experiment, the order of magnitude of the resultant contact resistance can be estimated. This value can be used for the finite element analysis of the contacting bodies and the training of the neural network. The error histogram of the trained network-obtained with a previously unused, numerically generated data set-clearly characterizes the reliability of the calculation (see Sect. 4.7).

Practical considerations
As was shown in the previous section, the discretized problem can be solved if the current generator is attached to a pair of contact nodes and the voltages between all the contact pairs are measured. This measurement configuration is easily realizable if two thin metal sheets are in contact.
If the determination of contact properties is critical, it is also possible to make preparations for the measurement already during the design of the contacting parts of an assembled structure. For example, drilling thin tunnels into the parts may allow the attachment of special electrodes in the vicinity of the contact surface. Certainly, the access to all the nodes cannot be ensured in all of the practically relevant cases. In these cases, instead of performing voltage measurements at each node, one can specify the set of reachable nodes, where the current source and/or the voltmeter can be connected to the bodies. There are several possible ways to attach the current generator to the pair of contacting bodies. Each of these scenarios leads to a set of equations like the one analysed in the previous section. Certainly, the currents will be typically different, but the contact resistances remain unchanged. This approach leads to the increase of the number of equations, but opens the way to evaluate the contact resistances based on measurements performed on the perimeter of the contact region.
The nonlinear nature of the equations is not eliminated if not all the contact voltages are measured. Still, our numerical experiments based on the derived Kirchhoff's equations showed that the solutions remain unique if the problem is restricted to realistic positive resistance only.
Thus, we focused on the improvement of the method's precision. A trivial choice is the measurement of voltage in the immediate vicinity of the current generator's electrodes. This arrangement proved to be practical from the point of view of the measurement set-up, because it was easier to move the four electrodes together. Moreover, our experiences indicated that the increase of the current led to a more clearly recognizable contact pattern. Consequently, it is plausible to think that the voltages measured far from the current source's electrodes are typically less sensitive to the variation of contact resistances than the ones that are closer to them.
The characterization of the contact in areas far from the free surfaces requires quite refined measurement techniques. For the sake of simplicity in the present paper, we restrict ourselves to the one-dimensional problem, when the unknown contact resistance varies along a closed curve γ -illustrated in Fig. 4. The assumption of a line-like contact is an appropriate approach in several cases, like in case of the sealing contact between cylinder head and engine block in internal combustion engines, or the contact between flanges, etc. Certainly, insulated bolt connections must be used for the application of the proposed method in these cases. It is also worth to mention that preliminary studies showed that the majority of the contact force between certain conforming surfaces-e.g., the conical surfaces of tool holders and the main spindle in milling machines-is distributed near the edges. Therefore, the introduced model can be suitable for characterizing these sort of connections, too. One of the motivations of the present research work was the need for a method with which the repeatability of the tool holder-spindle connection can be monitored [33].

Measurement of electric potentials and the net external force
The first step of the measurement is the partition of the contact area that determines the resolution of the resistance and force distributions. In the most favourable arrangement, electric measurements take place along measurement lines, each belonging to a small region of the partition along the perimeter of the contact, as is shown in Fig. 5.
During the measurement, four electrodes touch the surfaces of the contacting bodies at points A, B, C and D. The distance of these points must be chosen such that singularities in the numerical model (see Sect. 4.6) can be avoided. The current I AD is generated by a current source between points A and D at each measurement line, while the voltage U BC is measured between points B and C on the same line.
Clearly, the voltage U BC divided by the current I AD is not equal to the contact resistance of the actual peripheral region, since current can flow through other contact spots, too, even rather far from the measurement line. Even in the hypothetical scenario when there is only one finite region with nonzero conductivity, there The partitioned contact surface is illustrated by the tessellated pattern. Measurement lines are assigned to each region along the perimeter of the contact. Current is generated between the outer points A and D, while voltage measurement is realized at the in-between points B and C at each measurement line will be a finite potential difference between each node-pairs on the two opposing bodies. Thus, we introduce a new notation: will be referred to as effective resistance. Since the effective resistance is measured at each of the N u measurement lines, one can introduce the following vector: To determine the conductivity distribution on the contact surface based on the measured effective resistances To gain information about the specific resistivity ρ 0 of the contacting materials, too, a slightly modified measurement procedure is conducted before or after the described contact resistance measurements, in such a way that all the four electrodes are shifted along the measurement lines so that they touch the surface of the same body. The value of ρ 0 can be determined by the finite element model of the contacting bodies. This information can be used during the processing of contact resistance measurements, and it is essential for the training of the neural networks.
Although there are published results about the connection between the contact resistance and the contact force [29], most of the publications deal with non-ferrous materials and the standard deviations of the published measurements are usually very large. Thus, it is worth to supplement the aforementioned voltage measurements with the measurement of the net force F that presses the bodies together, since it can be used for the calibration of the measured force distribution. Moreover, if the net moment of the external force system is zero (that can be assured in several cases using self-aligning bearings) or the components of the moment (M x and M y ) can be measured, additional constraining equations can be formulated for an even more accurate calibration.

The measurement set-up of contacting tubes
The considerations presented in the previous section led to the development of the measurement set-up (see Fig. 6), with which the contact between two tubes can be evaluated.
The contacting tubes (1) are made of structural steel (S355 JR), with diameter D 40 mm, height h 35 mm and wall thickness w 3 mm. Before the measurements, the contact surfaces had been turned or grounded flat, that led to a surface roughness of about R a ≈ 1.3−1.8 µm. The tubes are pressed together by a bolt (2) and a nut (3). Electrical isolation is realized by polyamide parts (4) that are equipped with self-aligning ball bearings (5). The net force and two components of the net torque in the contact plane can be measured by the force sensor (6). The generated current flows through electrodes marked by (7) in the picture. Between them, there are the electrodes (8) that are used for voltage measurement. All the four electrodes are placed onto a printed circuit board (PCB), and the tips of them can be impressed telescopically. These tips are supported by soft springs that ensure near-constant contact forces on the electrodes. Nevertheless, the contact resistances between the electrodes and the tube surfaces are eliminated from the voltage measurement, since the inner resistance of the voltmeter can be considered to be infinite compared to the contact resistance added by the electrodes. A computer controlled actuator (9) adjusts the angular position of the tubes, while another actuator (10) attaches the electrodes during voltage measurement and detaches them while the tubes are being positioned. The schematic wiring of the measurement is shown in Fig. 7a. A multimeter (see https://en.wikipedia.org/ wiki/Multimeter) V 1 (Keithley 2000) registers the voltage between the inner electrode pair with input resistance 10 GΩ. The device was used in the range 100 mV, where its resolution is 100 nV. Since the typical readings were on the order of only 3-4 mV, the absolute accuracy is approx. ± 30 ppm × range ± 3 µV. Note, however, that from the point of view of the proposed method, the systematic errors of the measurement do not play an important role. In order to evaluate the magnitude of the random error, we performed a measurement by shorting out the inputs of the multimeter. Setting the averaging and the bandwidth according to the expected speed of contact resistance measurements, the random error was found to be less than ± 5 ppm , i.e. ± 0.5 µV.
A current source (Matrix MPS-7033) adjusts the demanded current of I 0 ≈ 6 A. The actual value of current is determined by another multimeter V 2 (Keithley 2000) by measuring the voltage drop on a high precision shunt resistor of resistance R shunt 250 mΩ (type Y14680R25000D9L, 0.05 ppm/C o ). Although the shunt resistor's resistance somewhat depends on its temperature, the application of a practically constant current led to almost identical temperatures during the measurements. The 7V voltage provided by an ultra precision voltage reference (type LTZ1000, 0.05 ppm/C o ) was set to match the voltage drop on the shunt (approx. 1.5 V), using a digital-analog converter board (AD5791, 20 bit, 1 ppm resolution, 0.05 ppm/ • C) as a voltage divider. As the result of this set-up, the read voltage between the inputs of the voltmeter was very close to zero, practically eliminating the inaccuracy related to the reading. Multimeter V 2 was also used in the range 100 mV; thus, the absolute accuracy of the current measurement at I 0 ≈ 6 A and U ≈ 1.5 V is approx. ± 12 µA. Without the systematic error, this value corresponds to accuracy better than ± 2 µA.
The accuracy of the current measurement could have been improved even further using an instrumentation amplifier (see Fig. 7b) that reduces the noise of the measurement system considerably.
Once the voltage and current values have been detected, the electrodes are detached and actuator (9) rotates the tube assembly by the increment that corresponds to the partitioning N u along the circular curve γ of the apparent contact. With the introduced test rig, N u 600 voltage measurements were performed during each of the conducted experiments with the increment of 360 • /N u 0.6 o . The voltage and current data of the particular measurement series that is evaluated in Sect. 5.2 have been made available to download [34].

Data filtering and compensation
For the reliable processing of the measured voltages and currents, the data set had to be compensated and filtered. This procedure comprises several steps: • Due to small electrical noises, the multimeters display nonzero current (I ref ) and voltage (U ref ) even if no current is led through the assembled tubes. These reference values were detected at each measurement line.
After turning on the current source, the voltage U i and current I i values were also detected and the effective resistance values were determined: • To take the temporal variation of contact resistance into account, a linear compensation was applied. Assuming that the mechanical contact properties do not change in the same manner as the resistance, the compensation is based on the effective resistance measured at the first measurement line R eff 1 and the effective resistance R eff N u +1 measured at the same position, but at the end of the measurement: It is worth to note that the variation of contact resistance was found to be negligible after approx. 60-120 min.; thus, the compensation step can be skipped if the bodies had been in contact for a long time [28]. It can happen that the temporal change of the contact resistance during the measurement is significant or even larger than its spatial variation along the curve. In this case, the error of the applied linear compensation may worsen the resolution of the method. This problem can be solved by a quicker measurement procedure or the simultaneous detection of voltages at all of the measurement lines. • The noise level and the frequency content of the compensated resistance data set R eff,c i was reduced by three filters: • First, a moving median filter was applied-in the cases discussed in this paper, the window length was set to 5 on a data set of N u 600. • Then, a moving average filter was used-in our case, the window length was set to 10.
• Since the data are processed by a neural network, trained by FE method, the last filter was designed to fit to the resolution n c 30 of the FE model (see Sect. 4.6). According to the well-known Nyquist-Shannon sampling theorem [35], at least two sampled values are needed in each period for the reconstruction of a certain frequency component. This is why a low-pass filter was needed with a cut-off period set to n c /2 15. According to our previous experiences, a widely used Chebyshev IIR digital filter was designed with transition band (0.5n c , 0.6n c ) (15,18). In order to obtain zero-phase digital filtering, the input data was processed in both forward and reverse directions. As Fig. 18 illustrates, the application of this filter does not have significant effect on the R eff curve, but it makes the calculation of the second discrete derivative (see Sect. 4.5) 2 R eff more reliable.

A simple evaluation of measured resistances
In the 1D case, the contact properties can be characterized by the vector of contact conductivities To understand the physical connection between the effective resistance R eff and the contact conductivity c, one can imagine the scenario when there is a point along the perimeter of the contact with high contact conductivity (i.e. low contact resistance). As the electrodes approach this point, the effective resistance drops. Similarly, the effective resistance increases as the electrodes move away from this point. It means that the local minima of the effective resistance correspond to well-conducting spots.
In case of continuous functions, the second derivative of the function is positive at the local minima. Thus, we expect to find conducting spots at high values of the second derivative of R eff with respect to the position coordinate s along the perimeter of the contact. Since the measurement provides a set of discrete effective resistance values R eff i , a discrete version of the differentiation must be used for the contact evaluation by the proposed simple approach. Let us introduce the second discrete derivative of R eff where the elements of 2 R eff can be calculated as for i 2, . . . , N u − 1. Since γ is a closed curve, the second derivatives at the first and last positions are respectively. As it will be illustrated in Sect. 5.1, the variation of the second discrete derivatives provides information about the conductivity distribution and the contact pattern along the contact curve. However, in general, there is no closed-form solution for determining the connection between 2 R eff and the vector c of contact conductivities. Indeed, the shift of the contact resistance values by the same constant would clearly influence the conductivities, while the second derivatives would not change. Therefore, a numerical method is involved to quantify the assumed relationship.

Numerical determination of effective resistances based on contact conductivity distribution
We set up a finite element (FE) model of the two contacting steel tubes that is capable of the determination of effective resistances based on prescribed contact conductivity distributions.
The FE model was built by using 20-node quadratic solid electric elements with a single degree of freedom, that is, the electric potential under constant current. A commercial finite element software, ANSYS, was used for the simulations, and the applied element type was SOLID231. The finite element mesh is illustrated in Figs. 9 and 10, while Fig. 11 represents the electric potential distribution for a given input of contact resistance distribution c.
The model consists of the two bodies in contact and a thin boundary layer between them, as presented in Figs. 8 and 10. The thickness t of this boundary layer must be chosen carefully: if it is sufficiently small, it has negligible influence on the simulation results.
The boundary layer was divided into n c pieces (n c is not necessarily equal to the number of measured voltages N u ), and the specific resistivity of each piece was set individually in order to achieve the demanded conductivity for the corresponding region.  Fig. 9: thin, sectioned boundary layer (coloured elements) between the contacting bodies to model non-uniform contact resistance Fig. 11 Electric potential distribution for a given distribution of effective resistances on the boundary layer By using the finite element model, any discretized non-uniform resistivity distribution can be simulated, and the resulting electric potential field can be inquired. Clearly, the geometry of the bodies and the arrangement of the electrodes must be identical to the measurement set-up. Just as in the case of the measurements, the current source is swept around, and electric potentials are registered at the pairs of points shown in Fig. 8. By realizing the finite element simulations in a scripted form, several thousands of randomly generated resistivity distributions can be evaluated and analyzed automatically.
For the sake of the simpler evaluation of the measured voltage values, a dimensionless conductivity parameter was introduced. The electric resistance R c of the boundary layer can be expressed as where ρ is the specific resistivity of the contact layer in the FE model that can be varied along the circumference cell by cell. Parameter t is the thickness and A is the area of the virtually added layer. Let us introduce a new parameter where ρ 0 is the specific resistivity of the materials in contact, and λ represents an equivalent thickness of the mating materials that has the same resistivity as the virtually added layer. This latter parameter may help to visualize the magnitude of the contact resistance. As discussed in Sect. 2, we assume that the majority of the contact force is transmitted through a number of domains with high contact force, yet the greater part of the apparent contacting regions is subjected to relatively small or even no force. With respect to this assumption, the use of a dimensionless conductivity parameter c is beneficial instead of using parameter λ. The advantage of using a conductivity-like quantity instead of a resistivity-like one is that the no-pressure zones will be characterized by c 0 instead of λ ∞, hence providing more favourable range of numerical values for the analysis. By definition, where L is a characteristic linear dimension of the examined pair of bodies-it is in the order of the maximal possible length of a current line that passes through the contact. In case of the performed experiments, L was chosen to be approximately equal to the perimeter of the tubes, i.e. L 125 mm. The thickness of the boundary layer was 0.1 mm since-according to our experiences-this is sufficiently small not to influence the simulation results.
The voltage measured at a certain position is influenced by the distance of the well-conducting spots. The farther these spots are, the longer distance must be covered by the electrons in the material. During the measurement of voltages along the perimeter, these distances vary, implying that the effective resistances vary, as well-this is why the positions of the spots with various conductivity values can be determined. Thus, to be able to accurately register the variation of these distances, it is beneficial if the contribution of the material resistance is significant compared to the contact resistance. According to (20), parameter c characterizes the ratio of these two components of effective resistance, thus, high values of the maxima of c lead to a better resolution.
Consider the input of the FE simulation as a vector containing the prescribed dimensionless conductivity values at the contact regions. The output vector contains the effective resistances calculated by the simulation. Here N u > n u > n c in a typical case-see the next paragraph and Sect. 4.7. The finite element simulation can evaluate the function F(c) R eff , and the aim is to identify the inverse of the function F, such that F −1 (R eff ) c.
Excessive numerical studies revealed that the inverse problem can be solved with good accuracy if three effective resistance values are determined for each of the n c contact regions: two at the endpoints and one in the middle. Since the values at the endpoints belong to the neighbouring regions, too, altogether n u 2n c effective resistance values must be calculated numerically for a given conductivity distribution. Figure 12 shows the results of a finite element simulation. In the presented case, the angular resolution of the boundary layer was 6 degrees, i.e. n c 60, and voltage probes were placed above the centres and at the borders of each boundary layer division, resulting in n u 120.
As it is also depicted in Fig. 12, the second derivatives 2 R eff seem to be running together with the dimensionless conductivity, that supports the possibility to identify a simple qualitative rule of assignment between the two.
Indeed, based on the numerical results presented in Fig. 12 and on the measurements detailed in Sect. 5.1, the pattern of the contact can be reliably characterized by the second discrete derivative of the effective resistance. However, the magnitude of contact resistances cannot be determined using the second derivatives.

Training the neural network
According to the previous section, a more accurate-possibly nonlinear-relation must be established for the quantitative description of the contact.
This goal can be achieved by using a neural network realized in the Deep Learning Toolbox of MATLAB, trained and tested on the basis of FE simulations. By means of this procedure, the dimensionless conductivity distribution can be determined for any given effective resistance vector.
The training inputs of the network were the R eff vectors resulted by FE simulations, and the targets were the inputs of the FE simulations, i.e. the randomly generated vectors c with mean values in the range (10 −6 , 30). In case of the evaluation of a measurement, the mean value of the measured voltages is compared to the simulated values. Using the information provided by the measurements, the range of vectors c can be narrowed. Besides that, the feasibility of the method and uniqueness of the solution can be evaluated by a simplified preliminary measurement. For training the network, several training methods and training parameters-e.g., the number of hidden layers-were tested. Finally, the so-called Levenberg-Marquardt backpropagation method was chosen.
Using this nonlinear least-squares solver, the most beneficial number of effective resistance values n u for a given number of boundary layer elements n c was heuristically found: the optimal ratio is n u 2n c .
The resolution n u of the FE model is not necessarily equal to the desired resolution N u of the contact evaluation. The reason of this difference is that the processing of the numerical results becomes more difficult-and even unfeasible-over a certain value of n c , due to a combinatorial explosion. This is why we applied the following strategy: The number N u of measured effective resistances was set to a relatively high value (N u 600 in the performed experiments) to gain as much information about the contact properties as possible. However, the boundary layer was divided only to significantly fewer pieces (n c 30 and n u 60 led to reliable results in the examined case). Thus, 600/30 20 effective resistance values could be assigned to each of the n c 30 boundary layer elements of the FE model in 20 different ways, by shifting the data. The generated 20 smaller data sets could be processed by the neural network without numerical difficulties. Combining the results of these calculations (each providing n c 30 conductivity values), the resolution could be increased to N u 600.
The error histogram of the trained network is shown in Fig. 13, where c is the relative error between the neural network output and the reference conductivity for a given R eff input. The training data set consisted of 1,500 simulations, while the test data set consisted of another 3,000 simulations of those, that had not been utilized during the training.
Thus, we can conclude that the used neural network is able to find the accurate conductivity distribution. If this numerical tool is supplemented with the results of a quick, preliminary measurement of the resultant contact resistance, the reliability of the method can be checked.

Determination of contact pattern and contact force distribution
The further processing of the conductivity distribution data set c depends on the goals of the analysis and the availability of additional information about the contact.
• A possible practical goal can be the assurance of repeatable fixing between frequently exchanged machine parts. In this simplest case, even the pattern of the conducting areas may reveal the imperfections in the contact, like asymmetries or accumulation of contaminating particles. As is illustrated in Sect. 5.1, even tiny irregularities can be found using the second discrete derivatives of the effective resistances. • If the distribution of the contact force is to be determined, one has the following options (listed in the order of increasing accuracy): • There are published results about the relation between the normal force F and the contact resistance R c for a great number of contacting materials. The most common form of this relation is Thus, in principle, two parameters: a coefficientâ and an exponent n are needed from the literature to determine the force distribution. The value of exponent n is usually reported to be between 1/3 and 1.4, but n 0.5 is the most frequently used approximate value in the engineering practice-see [6,14,29]. The application of the appropriate empirical formula (23) for each of the previously introduced regions provides a rough estimate of the force distribution. The drawback of this approach is that the measurement results typically have a large deviation. • Since the force-resistivity relation depends on many factors (temperature, surface quality, relative humidity, etc.), the accuracy of this approach can be improved by conducting measurements under the required conditions, with pre-fabricated specimen and quasi uniform force distribution. The information gained by conducting such experiments can be utilized for the calibration of the measurements, performed with the real structure under investigation. • To obtain more reliable results, the magnitude of the net force, pressing together the surfaces, must also be known. Based on this information, the force distribution can be scaled such that its resultant becomes equal to the known force. • To achieve maximal accuracy, the components of the net moment of the external force system must be known, too. These pieces of information provide additional constraints that can be exploited to enhance the results.

Validation efforts
To be able to validate the results and make us assured that the contact pattern measurement provides repeatable output, our method was tested on special arrangements. Since it is difficult to identify a physical model with the output of a neural network, we used the second discrete derivative for the evaluation of the contact.

Artificially introduced line-like connections-quick evaluation
To test the method for finding areas of good conductivity, three thin (d ≈ 0.2 mm) copper wires were placed in the radial directions between the contacting tubes, hence realizing line-like mechanical and electrical contacts-that correspond to point-like contacts in the one-dimensional model. This arrangement must result in a pressure distribution with three definite peaks at the places of the interposed wires. Following the measurement process detailed in the previous sections, voltages and currents were measured at N u 600 circumferential pairs of points.
The linear and nonlinear filters, described in Sect. 4.4, were applied on the measured values, resulting in the vector that can be considered as a kind of sampling of an R eff (ϕ) effective resistance function along the circumference. The results are depicted in Fig. 14. The top chart presents the filtered effective resistance values, while the bottom chart shows the second discrete derivatives along the perimeter.
As it can be seen, the peaks of the second discrete derivative of R eff denote the places of the interposed wires unambiguously. Certainly, the heights of the peaks depend on the sizes of the contact areas: the more point-like the contacts are, the sharper and higher the peaks will be. Thus, we can conclude that in case of a few small contact points, the pattern of contacts can be identified based on the introduced second discrete derivative.

Experiment with photographic film
To examine how a well fitted contact can be characterized by the proposed method, the mating surfaces of the two tubes used for the experiments were turned flat by a turning machine and the other ends of the tubes were covered. A piece of photographic film was placed inside the tubes that were attached to each other with the bolt and nut shown in Fig. 6. Light could enter only along the microscopic gaps along the contact line. The assembly was lit by bright sunshine from various directions for a while, and then, the effective resistance measurement was performed.
The picture of the developed film, together with the second discrete derivatives is shown in Fig. 15. As it can be seen, three elongated black spots appeared along the contact line of the tubes (emphasized by three circles in the figure). These non-contacting gaps can be identified by the small values of 2 R eff .
Apparently, the three small gaps appear equidistantly in the developed film; however, one would expect a more stochastic pattern. The explanation of this surprising result was revealed only after the experiment: the three-segmented chuck of the used turning machine slightly distorted the tubes during turning. Thus, the turned surfaces had been almost planar while the tubes were in the chuck, but once the chuck was released, the tubes' surfaces became wavy. The height difference between the high and low spots was estimated by FE calculations: the calculation resulted in approximately 12 μm, consequently, the size of the gaps was less than 24 μm. We can conclude that the method is capable of the detection of discrete non-conducting spots.

The effect of current magnitude
In order to evaluate how the current magnitude influences the results, measurements at different current values had been carried out without the disassembly of the tubes. Since we found that the second discrete  Fig. 16 The effect of current magnitude. The tubes were not separated between individual measurements derivative 2 R eff of the effective resistances characterizes the variation of contact resistance reasonably well, the comparison of the results is based on this quantity. At first, electric potentials were measured twice around the circumference at 2.5 A. For the next measurement, the current was increased to 5 A. Near position ϕ ≈ 36 • , the current was gradually raised further to 6.5 A. Finally, two measurements were performed at the current value of 9 A. According to Fig. 16, we can conclude that the variation of the current has a negligible effect on the measurement results. This experiment also shows that despite the temporal variation of the contact resistance, the same contact pattern can be obtained by several subsequent measurements using the filtering and compensation procedure described in Sect. 4.4. Note that the application of a higher current has some benefits, since it enhances the signal/noise ratio. However, the increase of current is limited by practical issues. To fit to the characteristics of the measurement system (see Sect. 4.3), we used the current of 6 A for the other measurements presented in this paper.

Check of repeatability and application of machinist's blue
To test the method further, we marked the relative positions of two contacting tubes with ground mating surfaces and measured the effective resistance. One of the contacting surfaces was coated with machinist's blue. After that, the whole assembly was disassembled, the surface was cleaned and re-coated and the tubes The results are shown in Fig. 17.
As can be seen, the two curves are similar, reflecting that the measurements really characterize the contact pattern. The small shift in the results is the consequence of the unavoidable difference in the fitting of the two tubes. It is worth mentioning that the mean effective resistance was approximately twice as large during the second measurement than initially, due to the difference in the compressing force.
Although the contact pattern can be recognized based on the second discrete derivative, the application of machinist's blue did not lead to usable results. Only a small amount of paint appeared on the originally unpainted surface. Photographs were taken for the evaluation of the distribution of the appeared paint after both measurements, but the results were neither similar to the curves shown in Fig. 17, nor to each other. Indeed, for the proper application of machinist's blue, the two surfaces must be rubbed to each other. Since rubbing changes the surface pattern, the tubes were simply attached to each other.
It is worth to note that the contact pattern was tried to be revealed by the photoelastic coating technique, too. However, the resolution of this method was not enough for obtaining meaningful results.

Application of the method to a case intended to be in full contact
Several measurements were conducted with the two tubes without any pre-defined contact patterns. The contact surfaces had been turned flat, sanded with a scouring pad and cleaned with a rag. After that, the test rig was assembled and the compressing force was set by the bolt and nut assembly. The force and moment values were measured by a purpose-made force sensor. The obtained current and voltage values have been made available for download [34]. The obtained data sets were filtered and evaluated by the neural network, introduced in Sect. 4.7. Figure 18 depicts the values of the measured effective resistance and the calculated conductivities along the perimeter for five different values of the compressing force.
To check the neural network, its output was used to set the conductivity values in the FE model. The numerically obtained effective resistances agree with the measured ones with excellent accuracy. It is worth to emphasise that relatively good electric contact was assured during the measurement (c max ≈ 0.4). In case of smaller conductivity values, the reliability of the method deteriorates. In the examined case, the minimal value of c max at which the method could be used was approximately c max ≈ 0.04, but this can be dramatically improved in case of simultaneous voltage measurements around the circumference by eliminating the time dependence of contact resistances.
Note that the peaks of the last two measurements (t 4,5 ) are slightly shifted with respect to the first three curves (t 1,2,3 ) in Fig. 18. This deviation is the consequence of the small, but unavoidable frictional forces among the elements of the test rig. Namely, despite the application of ball bearings, some torque was acting on the contacting tubes when the already pre-stressed nut was tightened between measurements 3 and 4. This small perturbation changed the course of the curves. The determination of force distribution is based on the known formula [28,29] of the force-contact resistance relationship. According to (23), the dimensionless conductivity can be expressed as Since the parameters a and n are not known, we exploit that the contact curve lies on a plane (the xy plane) and the net compressing force F and the components of the net moment M x and M y can also be measured. If the tube assembly is supported by self-aligning bearings, these moments are almost zero. However, since we could measure M x and M y , the original bearings were replaced with non-aligning bearings for the measurements. If the force distribution was continuous in the form F(ϕ), such that the total force applied between ϕ and ϕ + dϕ is F(ϕ)dϕ, the net force and the components of the net moment could be expressed as where the mean radius of the tubes is r ≡ (D − w)/2. Since the contact line is a closed curve, the force distribution is periodic. Consequently, it can be expanded into a Fourier series: The number of terms was determined in accordance with the filtering of the measured effective resistances, where the cut-off period was set to 15. Substituting this expression into Eqs. (26) and exploiting that the sin and cos functions are orthogonal, one obtains Thus, three coefficients of the Fourier expansion can be obtained as The force distribution can be discretized along the perimeter of the tubes by defining a vector of angular positions ϕ such that Using these discrete angular values, the force acting between ϕ i and ϕ i+1 is expressed as Thus, applying the Fourier expansion and exploiting the relations (29), the force values can be formulated as Thus, there are altogether 30 unknown coefficients: a s 2 , a s 3 , ..., a s 16 , a c 2 , a c 3 , ..., a c 16 , besides the parameters a and n in (25). We assume that the contacting surfaces can be characterized by the same values of parameters a and n along the total length of the contact. We also suppose that this relationship can be applied independently on each of the N u pre-defined sections of the partitioned contact surface.
The data vector of conductivities c (c 1 , c 2 , . . . , c N u ) is known. Thus, using (25), one obtains N u equations in the form For the solution of this set of equations, a least square error minimizing procedure was used. As a result, an optimal set of the unknown parameters could be determined. To check the results, the values of F n i /a were compared with the measured dimensionless conductivity values c i .
The unconstrained optimization procedure resulted in unrealistic values of the exponent n, at about n ≈ 2.3. Since n 0.5 is an accepted value in the engineering practice, we prescribed n ≤ 0.5. This modification led to an increase in the range of force values with respect to the unconstrained calculation: the difference of minimal and maximal force was increased to approx. 30 N from 7 N. The obtained force distribution F i is shown in Fig. 19. The index j of t j in the legend refers to the temporal order of the measurements.
The effect of the performed correction procedure can be seen clearly if one compares Figs. 18 and 19: the three maximal peaks (at ∼120°, ∼180°, and ∼320°) of the conductivity diagram in Fig. 18 assume similar values, but the corresponding force values are rather different. Indeed-despite the careful preparation of the contacting surfaces, the local electric contact properties may vary from point to point due to surface contaminants or various physical or chemical reactions. This inhomogeneity is corrected by the correction step that is based on the measured net force and two components of the net torque. As a consequence of this correction, the relative heights of the peaks changed.  There is another effect that could have some influence on the determined conductivity curves: the temporal variation of contact resistance is accounted for by a linear compensation procedure, as described in Sect. 4.4. If the nonlinearity of this temporal change is significant, the linear compensation distorts the results by reducing the measured values in the middle part of the data set. Since the relative heights of the peaks at ∼20°and ∼210°a re almost the same in the conductivity and force diagrams, this effect can be neglected in the considered case.
To avoid this source of error, one can perform the voltage measurements simultaneously, or in a random-like pattern along the contact curve.
To improve the accuracy of the force measurement, it can be useful to identify those parts of the contact curve where the two bodies are not in actual contact. This way, a zero force level can be set. A possible solution can be the introduction of an artificial small dent or gap somewhere between the contacting surfaces. Certainly, this approach cannot be used in all of the possible practical applications, and further research work would be necessary to determine the optimal size of this non-contacting domain. In some cases, the positions of non-contacting areas can be revealed by the technique described in Sect. 5.1.2.
Since the contact pressure cannot assume negative values in reality, this restriction can be utilized to check the reliability of the applied procedure. For example, the conductivity values in Fig. 18 are strictly positive, rather far from zero. Still, the correction based on the net force and moments transformed the curves such that the force level dropped to near zero at several points in Fig. 19. Only a couple of points with small magnitude negative forces (F min ≈ −0.084 N) were found close to ϕ 194 • , but practically, the force values are still positive. Since the zero level of contact force is tuned only by means of the aforementioned correction, the found near-zero forces (that correspond to small, but appreciable conductivities) confirm the validity of the proposed method. Thus, we claim that the proposed procedure is able to estimate the order of magnitude of the force distribution. The reliability of the method is reflected by the consistency of the found values of parameter a, as listed in Table 1.
The proposed procedure may have several areas of application. Besides the quality control of electrical components and assembled machine parts, it can be utilized for the online monitoring of possibly varying contact force distributions and can provide reliable basis for contact and friction models in the future.

Conclusions
A novel experimental method was introduced for the evaluation of the actual mechanical contact between metallic bodies, without separating them. The proposed procedure is based on the measurement of voltage values between opposing points of the contacting bodies, while the electric current through the contact is kept constant.
The possibility of the determination of contact resistances was analysed using a lumped parameter model of the problem. Based on this model, it looks promising to apply the method to general, 2D cases. Still, to test the measurement technique in a simpler case, the present paper focuses on the contact of two bodies along a curve.
It was shown for the contact of two tubes that the second discrete derivative of the introduced effective resistance can be utilized for the characterization of the contact pattern. Both point-like contacts and small gaps could be clearly identified by this simple approach and the repeatability of the results was also checked.
In order to obtain a more accurate description of the contact properties, the electrical measurement was supplemented by a series of finite element calculations that was used to train a problem-specific neural network. It was shown that the electrical conductivity distribution can be obtained accurately by processing the properly filtered measurement data set by the neural network.
The most challenging part of the method is the determination of the mechanical force distribution along the contact line. Based on the known functional form of the relationship between the electrical contact resistance and the mechanical force for the given pair of mating materials, the conductivity distribution can be converted to the aimed force distribution. The accuracy of the results can be improved if the net force and the components of the net torque are also measured. A Fourier expansion-based procedure was introduced in order to properly utilize the results of force and torque measurements. Although the applied optimization technique had to be supplemented by a parameter value known from the literature, the uncertainty of the force distribution was kept below a half order of magnitude. While the force distribution cannot be determined as accurately as the electrical conductivity distribution, the authors believe that the developed method can outperform the photoelastic coating technique or the application of strain gauges. of the contact resistances. Then, the corresponding nodes on the other body are numbered from n 0 + 1 to 2n 0 . By connecting the nodes from 1 to n 0 and from n 0 + 1 to 2n 0 , and finally, by connecting the nodes n 0 and 2n 0 , one obtains the so-called spanning tree of the graph of the circuit. The spanning tree-that contains only one contact branch-is denoted by red lines in Fig. 3. In the example, the spanning tree from node 1 to node n 0 + 1 passes through nodes 2, 3, 4, 5, 2n 0 10, n 0 + 4 9, n 0 + 3 8, ..., and n 0 5. 2. The numbering of the nodes helps to set up the so-called incidence matrix that is used to formulate the nodal and loop equations. The rows of this matrix represent the nodes, while the columns correspond to the branches. A possible way to systematically construct the matrix is to consider the connections of the neighbouring nodes of one of the bodies first (i.e. nodes with consecutive indices), then the second neighbours, etc. To illustrate the procedure, we show a part A p of the incidence matrix, representing the connections within the upper body in Fig. 3. In the first column, the connection between the first and second nodes is shown by 1 (current flowing out from node 1) and − 1 (current flowing into node 2).
The same structure of connections can be found among the nodes of the other body. Taking into account the contact branches, too, one can systematically obtain the 2n 0 × b incidence matrixÃ 0 : where 0 denotes a zero matrix sized n 0 × 2n 0 , while U is a unit matrix of size n 0 . The matrices U and −U describe the outflowing and inflowing currents at the contact pairs, respectively. Since each column represents a branch, the column indices provide a numbering of the branches and the currents.
To obtain independent nodal equations, the last line-corresponding to node 2n 0 -can be dropped. The resulting 2n 0 − 1 × b matrix is referred to as fundamental incidence matrix A 0 . We found that the structure of equations is simpler if the current generator is connected to the contact branch of the spanning tree, just like in Fig. 3. In this case, the drop of the last row means that the node from where the current flows out is not considered in the equations. Consequently, the current of the generator appears only in one of the nodal equations. This choice does not imply any loss of generality, since a new spanning tree can be defined for each different attachment of the current generator. 3. For the algorithm-based compilation of the set of Kirchhoff's equations, one can construct the incidence matrix of the spanning tree, as well. In the example, it assumes the form The first 4-4 columns describe the connections within the upper and lower body, respectively, while the last column represents the contact branch in the tree. The size of this matrix is 2n 0 × 2n 0 − 1. By dropping the last row, one obtains the quadratic fundamental incidence matrix A T of the spanning tree. 4. The exchange of columns in the incidence matrix (36) corresponds to a renumbering of branches. Similarly, by a proper renumbering of the fundamental incidence matrix A 0 , an equivalent fundamental incidence matrix can be constructed in the form Here, the columns of the fundamental incidence matrix A T of the tree form the last 2n 0 − 1 columns of the fundamental incidence matrix A. Using matrix A, the so-called fundamental cycle matrix B can be calculated in the form where U is a unit matrix of size b − n + 1 and 5. In the next step, the currents of the branches are arranged into a vector I of length b. The 2n 0 −1 independent nodal equations are obtained as where the current I 0 of the current generator appears at position n 0 if it is connected to the contact branch of the tree (see Fig. 3).
In order to obtain the loop equations, the resistances characterizing the branches are arranged into a diagonal matrix-the so-called impedance matrix Z-according to the numbering defined by the columns of A. It can be shown that the loop equations assume the form (BZ)I 0.