*Article* **Structural, Electrical, and Mechanical Properties Investigation of Open-Cell Aluminum Foams Obtained by Spark Plasma Sintering and Replication on Polyurethane Template**

**Alexandra Kosenko 1,\*, Konstantin Pushnitsa 1, Artem Kim 1, Pavel Novikov <sup>1</sup> and Anatoliy A. Popovich <sup>1</sup>**

Institute of Machinery, Materials, and Transport, Peter the Great Saint Petersburg Polytechnic University, Politechnicheskaya ul. 29, 195251 Saint Petersburg, Russia; pushnitsa.k@gmail.com (K.P.); artem\_7.kim@mail.ru (A.K.); novikov.p.a@gmail.com (P.N.); popovicha@mail.ru (A.A.P.) **\*** Correspondence: alxndra.kosenko@gmail.com; Tel.: +7-929-105-73-38

**Abstract:** The present paper illustrates a comparison of open-cell aluminum foams. The foams were fabricated by two different methods: spark plasma sintering and replication on a polyurethane template. The influence of pressure, temperature, and diameter of space holding material on foam obtained by the spark plasma sintering method was investigated. Additionally, the aluminum powder content in slurry and atmosphere during thermal processing of foam prepared by the replication technique were studied. The morphology and structure of obtained samples were analyzed by scanning electron microscopy and X-ray diffraction analysis. Supplementarily, mechanical properties and electrical conductivity were studied. The porosity of obtained samples was 83% for the SPS sample and 85% for the replication sample. The results of the studies carried out gave us an understanding that the SPS method is more promising for using the obtained foams as cathode current collectors in lithium-ion batteries due to excessive aluminum oxidation during sintering in the furnace.

**Keywords:** open-cell foams; porous materials; spark plasma sintering; replication technique; PU foam template

#### **1. Introduction**

Metal foams are currently a special class of lightweight materials that can be used in various fields of the industry due to their physical, chemical, and mechanical properties [1–3]. Using metal foam instead of solid metal is more efficient and reduces the weight and cost of the final product [2]. Metal foam might be able to substitute a variety of conventional materials. Currently it is used in a wide range of engineering applications. Therefore, this fact allows applying these materials in various fields such as automobiles, aerospace, and naval industries [3]. It is possible to manufacture metal foam from pure metals, for example, titanium, nickel, and aluminum [4,5], and in the same way from different alloys [6,7]. Due to the lightweight nature of this type of material, metal foam is an extremely promising material for the automotive industry; in particular, open-cell foams made of aluminum and Al alloys are of special interest. Remarkable properties of open-cell aluminum foam such as low density, high specific electrical and thermal conductivity, and high surface area allow to be applied as a component of the lithium-ion battery, specifically to use aluminum foam instead of aluminum foil as cathode current collectors [8–11].

A variety of fabrication processes have been developed to produce metal foams, in particular aluminum foams [12–20]. The common processes are either via the addition of a foaming agent into molten aluminum alloy or through bubbling gas into a molten aluminum alloy. The main drawbacks are that aluminum foams produced by these methods are limited in terms of control of the pore content and size and microstructure of cell walls. The methods of production for open-cell aluminum foam can be divided into three classes:

**Citation:** Kosenko, A.; Pushnitsa, K.; Kim, A.; Novikov, P.; Popovich, A.A. Structural, Electrical, and Mechanical Properties Investigation of Open-Cell Aluminum Foams Obtained by Spark Plasma Sintering and Replication on Polyurethane Template. *Materials* **2022**, *15*, 931. https://doi.org/ 10.3390/ma15030931

Academic Editors: Marc Cretin and Zhenghua Tang

Received: 13 December 2021 Accepted: 21 January 2022 Published: 26 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

solid-state methods, liquid metallurgy methods, and vapor state methods [21]. Solid-state methods include the sintering of metal powder. It is well known that spark plasma sintering (SPS) decreases the sintering time and temperature, compared with sintering in an electric furnace, resulting in suppression of grain growth during sintering [22–25]. Therefore, it is worthwhile to investigate the microstructure and mechanical properties of aluminum foam processed by SPS. In addition, one of the promising ways to manufacture open-cell metal foams is the sponge replication technique, and this method has been successfully applied for the production of open-cell titanium, copper, and aluminum foams [26–33]. The main disadvantage of the replication method is the high oxidation rates of powders during the sintering step [34,35].

In the present paper, open-cell aluminum foams, fabricated by spark plasma sintering and by sponge replication technique, are investigated. The necessity of the research performed was in changing the standard current collector in lithium-ion batteries, aluminum foil, to aluminum foam, due to aluminum foam being a more lightweight and effective type of material. The main difference between this research and previously published ones is that no additives, such as foaming and thickener agents, were used. This fact can lead to a low impurity content in composition of the resulting foams, which is important in the further use of the foam. The investigation focused on foam production with a highly open structure, with low impurity content, a lower amount of oxide film on the surface, high electrical conductivity, and sufficient mechanical strength. The last two points are high priority in this research. Mechanical strength is necessary since the aluminum foam must not collapse when the active cathode material is applied. High electrical conductivity is important for low internal resistance and high specific energy capacity of the battery.

#### **2. Materials and Methods**

#### *2.1. Materials*

Figure 1 represents an SEM image of the aluminum alloy powder used for SPS method (a), sodium chloride (b), and aluminum powder used for replication technique (c) obtained in the backscattered electron detection mode. Aluminum alloy powder particles are predominantly rounded, the particle size distribution is as follows: d10 = 20.9 μm, d50 = 51.6 μm, d90 = 95.9 μm. It means that 10% of particles have a size less than 20.9 μm, 50% of particles have a size less than 51.6 μm, and 90% of particles have a size less than 95.9 μm. Additionally, there are some submicron size particles. The powder has a homogeneous structure. In Figure 1b sodium chloride particles are presented. Particles have mostly spheroidal and round form, the powder was sieved and fraction 500–1000 μm was chosen. Figure 1c shows the aluminum powder particles, all of the particles have a spherical shape, a homogeneous structure, and the particle size distribution is d10 = 2.4 μm, d50 = 5.5 μm, d90 = 10.8 μm.

**Figure 1.** (**a**) Aluminum alloy powder, (**b**) sodium chloride, (**c**) aluminum powder.

#### *2.2. Spark Plasma Sintering Method*

The commercially available aluminum alloy powder (AMD5, supplier Normin LTD., Oxford, UK) with the content of Al 94.8 wt.%, Mg 4.8 wt.%, Ti 0.4 wt.%, with an average

particle size of 90 μm, was used as a starting material. It is well known that Al–Mg alloys are characterized by a combination of satisfactory strength, good ductility, very good weldability and corrosion resistance [36]. Titanium in Al alloys also enhances microstructural and mechanical properties, including specific strength, oxidation and corrosion resistance [37]. Sodium chloride (supplier LenReactiv JSC, Saint Petersburg, Russia) with average particle size 500–1000 μm, was used as a space holder for SPS. The weight ratio of the aluminum powder to NaCl was determined to be 10:90 and the volume ratio 8:92. To select the optimal weight ratio of aluminum powder content and porosity of the final Al foam, mixtures with a content of 5% and 10% were prepared and foams with the same sintering conditions were made. In the case of 5%, the sample retains its shape, but begins to crumple with a slight pressure with tweezers, in the case of samples of 10% this is not observed. The spark plasma sintering method is schematically presented in Figure 2. In general, the process consists of 4 steps: mixing the aluminum powder and sodium chloride, pressing, sintering, and leaching in water. Firstly, the aluminum powder and space-holding material were mixed in the vibrating mixer (C50.0 Vibrotechnik LLC, Saint Petersburg, Russia) for 90 min. After the ingredients were homogeneously mixed, the resulting powder was pressed in graphite mold at a pressure of 5 MPa. Then, sintering was carried out on the equipment for spark sintering FCT system HPD 25 with the following parameters: sintering temperature 550 ◦C, pressure 38 MPa, sintering time 5 min, pulse duration 2–4 ms, pause duration 2 ms. The sintered sample after the sintering was placed into hot water for 12 h to leach out the embedded sodium chloride particles, to obtain the aluminum foam with porous structure.

**Figure 2.** Schematic illustration of spark plasma sintering (SPS) method.

#### *2.3. Sponge Replication Technique*

The commercially available aluminum powder with purity 99% and with particle size <10 μm was used (ASD6, supplier Normin LTD). The aluminum powder was mixed with earlier-prepared 0.4 wt.% cornstarch solution in distilled water. Cornstarch is a binder in this solution. The content of Al powder in final dispersion was 80 wt.%. The dispersion was mixed by using an ultrasonic bath to make it homogeneous. This method consists of 3 steps: (1) coating a polyurethane (PU) sponge with a powder slurry; (2) burning the PU and binder out; (3) final sintering of the foam. The schematic illustration is below in Figure 3. An open-cell polyurethane foam cylinder (Polinazell PPI 20. supplier United

Service Company LLC, Chicago, IL, USA) with a height of 3 mm and a diameter of 19 mm was used as a template.

**Figure 3.** Schematic illustration of replication technique.

As a first step, the PU foam was drowned into previously mentioned dispersion, soaked for 5 min, and then the excess was removed by mechanical wringing. After this step, the dispersion-coated template was dried for 17 h in a vacuum at a temperature of 95 ◦C. Further, the sample was placed in a tube electric furnace to burn out the PU foam and binder. The final thermal processing step was carried out also in a tube electric furnace. All burning and thermal processes were conducted in the Argon atmosphere, the flow rate was 40 mL/min. The conditions of processing are listed below in Table 1. Due to the cooling effect of Argon flow, the sintering temperature is higher than the melting point of aluminum.



#### *2.4. Characterisation*

The morphology, microstructure, and chemical composition of the samples were studied using Mira 3 Tescan scanning electron microscope with an EDX Oxford Instruments X-max 80 energy dispersive detector for X-ray spectroscopy. X-ray phase analysis was carried out on a Bruker Advance D8 diffractometer in the range of angles from 25◦ to 85◦ with a step of 0.02◦ and an exposure of 1.5 s at each step. The wavelength was 1.5406 Å. The lattice parameters were measured by the Le Bail method using software TOPASS. The particle size (PS) distribution of the powders was studied using a Fritsch Analysette 22 NanoTec plus laser diffraction unit. To calculate the particle size distribution, the Fraunhofer model was used. The electrical conductivity was measured using an AKIP-2101 voltmeter, a GPD-74303S current source, and a FLUKE-289 multimeter. Mechanical compression tests of the specimens were carried out on a Zwick/Roell Z100 testing machine.

#### **3. Results and Discussion**

#### *3.1. Aluminum Foam Formation*

In Table 2 different conditions of sintering are demonstrated. First, it is worth noting the fact that at a temperature of 600 ◦C and a pressure of 19 MPa and above, aluminum powder melts and flows out of the mold. At 500 ◦C, and pressure 38 MPa and below, spherical non-deformed particles are observed, which indicates worse sintering. however, at pressure above 38 MPa the powder starts to melt. The best results of sintering were obtained at temperature 550 ◦C and pressure 38 MPa.

**Table 2.** Different conditions of spark plasma sintering method.


In Figure 4 the obtained samples of aluminum foam are presented. Figure 4a is the SPS sample, and Figure 4b is the replication sample. At first glance, it seems that the SPS sample has more open porosity than the sample fabricated by replication technique; in the next subsections the physicochemical and mechanical properties are considered in more detail. The shape and size of the samples were selected on the basis of further research as current collectors in lithium-ion batteries. The diameter and height of the samples should be no more than 20 mm and 2 mm, respectively, so that they can be placed in the CR2032 case for further testing.

**Figure 4.** Obtained open-cell aluminum foam by SPS (**a**), by replication method (**b**).

#### *3.2. Morphology and Microstructure Analysis*

Figure 5 shows the surface morphology of a sample obtained by the method of spark plasma sintering. The sample has a cellular structure with open porosity. It is apparent that the material has a homogeneous structure with a pore diameter of about 600 μm. It is possible to control the pore diameter of fabricating foam by choosing different fractions of space holding material—sodium chloride.

**Figure 5.** SEM image of foam obtained by SPS method: (**a**) general view, (**b**) pore size investigation.

Figure 6 shows images of the original polyurethane sponge and a sample of aluminum foam obtained by the replication method. The sample has open porosity. The pore diameter is in a wide range (300–1000 microns) and it is comparable to the pore size of the PU sponge. Small pores can be seen on the replication sample due to the incomplete repetition of the template relief, which is justified due to the high viscosity of the suspension.

**Figure 6.** SEM image of PU sponge (**a**) and obtained foam by replication technique (**b**).

From EDX spectra (Figure 7b) it can be concluded that spreading of alloy elements in the SPS sample after sintering is uniform, and the amount of sodium and chlorine is less than 1 wt.%; the presence of sulfur is due to its impurity in the composition of the salt. The results of the replication sample EDX analysis of Figure 7d represents that the chemical composition after thermal processing is characterized by 70.8 wt.% aluminum, and a large amount of the oxygen is recorded as well. However, this method gives a large error when measuring the number of light elements.

**Figure 7.** SEM image (**a**,**c**) and corresponding energy-dispersive X-ray spectroscopy spectra (**b**,**d**) of SPS sample (**a**,**b**) and replication sample (**c**,**d**).

Figure 8 shows the diffraction patterns of Al alloy powder (a), Al powder (b), the sample obtained by the SPS (c), and by the replication method (d). The main peaks of all diffraction patterns are indexed by the cubic Al phase (PDF 03-065-2869) with the space group Fm-3m. The unit cell parameters obtained by the Rietveld method are presented in Table 3. In the sample obtained by the replication method, the presence of the γ-Al2O3 phase is observed.

**Table 3.** Unit cell parameters were obtained by the Le Bail method.


**Figure 8.** Diffraction patterns of Al alloy powder (**a**), Al powder (**b**), SPS sample (**c**), and replication sample (**d**).

#### *3.3. Porosity*

The density of the aluminum foams was determined by hydrostatic weighing (Archimedes method). As a first step, the samples were weighed in air, then they were coated with paraffin to prevent the fluid from entering inside the foam pores. After this step, the samples were weighed in air and further in ethanol, the density of which amounts to 0.785 g/cm<sup>3</sup> at 25 ◦C. The relationship between the porosity *Pr* (%) of the aluminum foam and the relative density (*ρf*/*ρs*) of the foam are presented in Equation (1), where *ρ<sup>f</sup>* and *ρ<sup>b</sup>* represent densities of the foam sample and bulk aluminum (about 2.7 g/cm3), respectively. The density of the aluminum foam was determined by hydrostatic weighing. The measurements of samples densities and porosity are presented in Table 4.

$$P\_r = \left(1 - \frac{\rho\_f}{\rho\_b}\right) \times 100\tag{1}$$

**Table 4.** The measurements of density and porosity of aluminum foam.


#### *3.4. Electrical Conductivity*

The electrical conductivity of samples was determined by measuring the resistance using the 4-wire testing method. Electrodes for measuring the potential difference were installed directly on the sample with a constant distance between electrodes (10 mm) for each measurement. The accuracy of measuring the potential difference with the used device was 0.01%.

To determine the resistance of the samples, a direct current was applied. The accuracy of the applied current was ±0.5%. To clarify the measured results, the current was measured with an ammeter with an error of no more than ±0.2%. The potential of aluminum foil and samples was measured. Aluminum foil is a frequently used material as a cathode current collector in lithium-ion batteries. In this test, it was used as a reference to compare with. The diameter of the foil was 20 mm, the thickness, 14 μm. For each sample, three measurements were taken at different locations. The resistance of the samples for each measurement was calculated using Ohm's law equation. The results of the measurements are presented in Table 5.


**Table 5.** The resistance measurement of samples by using the 4-wire testing method.

The key factors affecting the electrical conductivity of the samples are the degree of sintering and the amount of aluminum oxide on the surface of the samples. The determination of the electrical conductivity of the samples showed results consistent with the X-ray analysis—the sample obtained by replication technique has the worst electrical conductivity, which can be explained by the formation of aluminum oxide film on the sample surface and the surface of grains. In addition, the calculated resistances of aluminum foil and SPS sample gave us an understanding that the SPS sample is the closest in electrical conductivity to aluminum foil. The obtained results demonstrate that the SPS sample obtained due to its properties is better for application as a three-dimensional current collector than the sample obtained by replication.

#### *3.5. Mechanical Properties*

The compression stress–strain curves of the obtained samples are shown in Figure 9. The dimensions of samples were: diameter—19 mm, height—2 mm. Three samples of each obtaining method were taken for the compressive test. For both samples, coefficient of variation was calculated. The average coefficient of variation for SPS samples was 12%, and for replication samples was 13%.

**Figure 9.** Compression strain curves for aluminum foams: SPS sample, replication sample.

The stress–strain curves of the SPS sample are characterized by three distinct regions (designate by solid lines). The first region is an elastic region characterized by a linear elastic region at very low deformation without peak stress. The second area is the plateau area on which pores collapse (destruction of the bridges), which may indicate the growing deformation without increasing pressure. Fluctuations of the curve in this area can be associated with the uneven distribution of pores in the sample. The third area is the area of compaction, where a sharp increase in compressive stress occurs.

The stress-strain curve of the replication sample is characterized by four regions (designate by dash lines). The first region is the region of elastic deformation; relative to the SPS sample, the region of elastic deformation ends at a lower pressure. In the second region, the pressure decreases with the continued growth of deformation, which may indicate not plastic compression of the bridges, but their brittle fracture. This behavior is demonstrated by all samples. The brittle fracture occurs due to the presence of aluminum oxide in the sample. The third region is less pronounced than for the SPS sample due to the smaller amount of aluminum in the sample. In the fourth region, the same process occurs as in region three for the SPS sample.

For the samples, the adsorbed energies were also compared during mechanical strength tests. Due to the small thickness of the sample, the tests were carried out without using an extensometer, which would not allow the calculation of the nominal adsorbed energy. However, it is admissible to compare the adsorbed energies of the two samples. The calculation of the adsorbed energies was carried out according to the method presented in the article [38]. The ratio of the adsorbed energies of the samples (a)/(b) was 5.24. This value shows that, upon compression, sample (a) absorbed 5.24 times more energy, which also indirectly indicates the metallic structure of sample (a) and the nonmetallic structure of sample (b).

#### *3.6. Properties Comparison of Samples*

Babcsán et al. have used the foam in their studies, manufactured by the batch casting process. It is possible to obtain aluminum foam with closed porosity by this method. The main difference between our method and theirs is that we pursued the goal of obtaining open-cell aluminum foam. Additionally, our research is distinguished from the Babcsán research since no additives were used in the foam production. In their method, calcium and titanium hydride were used as thickener and foaming agent, respectively. Using our method, as a result, we obtained a sample with a minimum content of impurities. Impurities, which were found by EDX analysis, also were present in the initial materials.

The comparison of samples is presented in Table 6. The low conductivity of the SPS sample can be explained by incomplete sintering of powder during the formation of the sample. The replication sample result was previously explained: the large amount of aluminum oxide has influence on the conductivity.


**Table 6.** Properties comparison of open-cell aluminum foams.

#### **4. Conclusions**

Open-cell aluminum foams were successfully obtained by two various methods: spark plasma sintering and replication on polyurethane template. The porosity of final samples was 83% and 85%, respectively. According to the investigations carried out on the microstructure and morphology, electrical conductivity, and mechanical properties of obtained samples, we can conclude that the sample obtained by the replication method has more aluminum oxide film on its surface than the SPS sample. The consequence of this fact

is the lower electric conductivity of the replication sample. The electrical conductivity of the SPS sample is comparable to the aluminum foil. This fact allows us to assume that the electrochemical behavior of the SPS sample will be similar to the foil as a part of lithium-ion batteries. Mechanical tests of the foams also demonstrated the presence of a brittle phase in the replication sample. The SPS sample showed metallic properties and determined sufficient mechanical strength. Due to investigated properties of the SPS sample, it is evident that the spark plasma sintering method is the most promising way to produce open-cell aluminum foam for further use in lithium-ion batteries.

The future research plan is to find the optimal balance of temperature, pressure and sintering time in order to improve the processing mode of the SPS method. Additionally, the further improvement of the samples' characteristics is possible by using fluxes, which limit the oxide film growth and improve the sintering of the aluminum particles. The next step is to assemble mock-ups of CR2032 batteries to study the influence of a new type of current collector on batteries' behavior using cyclic voltammetry and electrochemical impedance spectroscopy methods.

**Author Contributions:** Conceptualization, P.N. and K.P.; methodology, A.K. (Alexandra Kosenko), A.K. (Artem Kim), and K.P.; software, A.K. (Alexandra Kosenko) and K.P.; validation, A.K. (Alexandra Kosenko) and K.P.; formal analysis, A.K. (Alexandra Kosenko), A.K. (Artem Kim), and K.P.; investigation, A.K. (Alexandra Kosenko) and K.P.; resources, A.A.P. and P.N.; data curation, A.K. (Alexandra Kosenko); writing—original draft preparation, A.K. (Alexandra Kosenko); writing review and editing, A.K. (Alexandra Kosenko) and K.P.; visualization, A.K. (Alexandra Kosenko) and K.P.; supervision, K.P., P.N., and A.A.P.; project administration, P.N.; funding acquisition, A.A.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research is funded by the Ministry of Science and Higher Education of the Russian Federation as part of the World-Class Research Center program: Advanced Digital Technologies (contract No. 075-15-2020-934 dated 17 November 2020).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Multi-Steps Magnetic Flux Entrance/Exit at Thermomagnetic Avalanches in the Plates of Hard Superconductors**

**Viktor Chabanenko 1,\*, Adam Nabiałek <sup>2</sup> and Roman Pu ´zniak <sup>2</sup>**


**Abstract:** Avalanche cascades of magnetic flux have been detected at thermomagnetic instability of the critical state in the plates of Nb-Ti alloy. It was found that, the magnetic flux Φ enters conventional superconductor in screening regime and leaves in trapping regime in the form of a multistage "stairways", with the structure dependent on the magnetic field strength and magnetic history, with approximately equal successive portions ΔΦ in temporal Φ(*t*) dependence, and with the width depending almost linearly on the plate thickness. The steady generation of cascades was observed for the full remagnetization cycle in the field of 2–4 T. The structure of inductive signal becomes complex already in the field of 0–2 T and it was shown, on the base of Fourier analysis, that, the avalanche flux dynamic produces, in this field range, multiple harmonics of the electric field. The physical reason of complex spectrum of the low-field avalanche dynamics can be associated with rough structure of moving flux front and with inhomogeneous relief of induction. It was established that the initiation of cascades occurs mainly in the central part of the lateral surface. The mechanism of cascades generation seems to be connected to the "resonator's properties" of the plates.

**Citation:** Chabanenko, V.; Nabiałek, A.; Pu ´zniak, R. Multi-Steps Magnetic Flux Entrance/Exit at Thermomagnetic Avalanches in the Plates of Hard Superconductors. *Materials* **2022**, *15*, 2037. https:// doi.org/10.3390/ma15062037

Academic Editors: Sophie Tingry, Marc Cretin and Zhenghua Tang

Received: 12 February 2022 Accepted: 8 March 2022 Published: 10 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Keywords:** hard superconductors type II; thermomagnetic instability; avalanche flux dynamics; screening and trapping regimes; cascades of flux

### **1. Introduction**

Magnetic properties of hard type-II superconductors are usually described in the frames of Bean critical state model [1,2]. In this model, the superconductor volume is screened by a current of the critical density, *j*c. In increasing external magnetic field, magnetic flux enters into the superconductor's volume as a series of tiny flux jumps of different scales [3]. These jumps allow to relax the critical state and were studied experimentally [4,5].

In the terms of thermodynamics, the critical state is metastable. At certain conditions, small fluctuation of an external magnetic field, temperature, or tiny magnetization fluctuations in superconductor may lead to appearance of a catastrophic thermo-magnetic avalanche [6,7]. Magnetic flux jumps and accompanying heat release are the phenomena commonly observed both in conventional and in high temperature superconductors [8]. During the thermo-magnetic avalanche, the magnetic flux enters/exits into/from superconductors abruptly and magnetic moment decreases sharply. The temperature range, in which the critical state is unstable, is determined by the material parameters and can reach the temperatures as high as ten Kelvin. Low temperatures, where the critical current density increases significantly, are especially interesting for practical application. However, in this case, the instability of the critical state increases strongly too, leading to a giant flux avalanches.

An experimental study of avalanche flux dynamics, in particular, visualization of magnetic field penetration patterns, began more than fifty years ago [9–13]. Direct observation

with high speed cinematography of the flux distribution and its change in time over the entire surface of a superconducting sample made it possible to study the kinetic of flux jumps. Particular attention was paid to the study of the flux front velocity, both in pure metals and in alloys [10–14]. It was found that the flux front velocity for conventional superconductors is of an order of tens of meters per second [11,14]. At the same time, the structure of magnetic flux avalanche spot in disks was found to be very diverse: from circular spots [12] to flux fingers and "dendrite-like" structures ("irregular" jumps [13]) in inhomogeneous bulk superconductors. Here, "irregular" jumps are the jumps leading to a complex, irregular distribution of the magnetic flux in the superconductor. The analysis of the visualization patterns of the flux dynamics, presented in the above papers, allowed, within the framework of Bean's concept, to understand the mechanism of local inversion of the surface self-field in the volume of conventional NbTi superconductor [15,16] and the complete inversion of the magnetic moment in the YBa2Cu3O7−<sup>δ</sup> single crystal [17].

Avalanche flux dynamics in thin films, with penetration field depth greater than their thickness, attracts special attention because of several phenomena, recently discovered by means of magneto-optic, including dendritic flux patterns [18–24]. Such avalanches consist of sharp impulses of the magnetic flux rushing into the specimen. The "branches" of dendritic structures repel each other as they grow. Numerical simulations of the vortex avalanches in superconducting films confirmed a phenomenon of magnetic flux fragmentation [25,26], occurring when the hot spot of the thermomagnetic avalanche reflects from inhomogeneities or the boundary, on which magnetization currents either vanish or change direction. As a result, thermomagnetic avalanche produce complex flux patterns similar to dendrite.

Variation of Nb film geometry extended diversity of magnetic flux structures that penetrate superconductor during avalanches. One of such structures is referred to as "huge compact avalanches" [27] and it is different from dendrites. A recent study found an unambiguous relationship between the size and shape distributions of avalanches and thermomagnetic conditions of instability development [28]. It was established that the transition from a curved "finger-like" structure of avalanche at low fields to a dendritic one at higher fields allows dividing the regimes of dynamical process, where either thermal or magnetic diffusivity prevail.

In addition to the variety of avalanche structures, the study of the dynamic properties of the flux led to an unexpected result: the dendrite front propagation reaches the velocity up to 160 km/s [18,22]. An impressive experimental fact, discovered recently, is the rate of nanoportions of the flux, the Abrikosov vortex, penetrating into the superconducting Pb film [29,30] or into the Nb-C superconductor [31], being up to tens of kilometers per second and exceeding the theoretical estimates of pair breaking speed limit for superconducting condensate. Moreover, as it was shown, calculating the electromagnetic response of a nonequilibrium superconducting state, much lower velocities of the superconducting condensate lead to the appearance of multiple harmonics of the electric field [32], which means that the condensate dynamics is a highly nonlinear process. These results indicate that the study of the dynamic properties of magnetic flux helps to uncover unknown properties of the superconducting state. However, limited attention only was paid to the effective mass of Abrikosov's vortex [33–35], which can manifest itself, for example, in the experiments with circular dichroism of far-infrared light [36].

The purpose of current research was to establish the regularities of dynamic processes, occurring in a superconductor, when the path length of the finger-shaped avalanche is significantly limited by the size of the superconductor, i.e., it was supposed to implement the conditions of the "resonator" for a directed electromagnetic shot in the plates of the most popular hard superconductor. We present, here, an evidence of a multi-step "stairway" structure of magnetic flux dynamics at the thermomagnetic avalanches in the volume of the plates of conventional NbTi superconductor, being the number one material in the world in the terms of wide application. The width of the step, depending on the thickness of the plate, was determined in inductive measurements. The phenomenon of multi-steps

magnetic flux dynamics was observed, both at the flux avalanche entry (screening mode) and at its exit (trapping mode). Spectral analysis of the structure of inductive signal allowed us to characterize the avalanche flux dynamics in the entire region of instability of the critical state, and to establish the region of magnetic fields where steady state generation of cascades is realized. The experimental observation of multi-step "stairway" structure of magnetic flux dynamic at avalanches is reported here for the first time.

The experimental results presented here demonstrate the new dynamical properties, appearing during thermomagnetic avalanches, in the plates of type II superconductors. It was found that, the magnetic flux enters conventional superconductor in screening regime and leaves in trapping regime in the form of a multistage "stairways" (cascades), with the structure dependent on the magnetic field strength and magnetic history, with approximately equal successive portions in temporal dependence, and with the width depending almost linearly on the plate thickness. The mechanism of cascades generation seems to be connected to the "resonator's properties" of the plates.

#### **2. Experiments and Materials**

Time resolved measurements were performed with the aid of different type of inductive sensors, where a voltage (*U*coil (*t*)~dΦ/d*t*, *t*—time) caused by the avalanche flux jump was induced. Figure 1a shows flux jump at *T* = 1.7 K in Nb disc, frozen in the form of "fingers". Figure 1b–d shows the arrangement of pickup coils on superconducting plate for the avalanche registration: in the whole plate's volume (b), on the three height levels (c), and in the magnetic stray field (*L*ext) (d). Time-dependent voltage *U*coil (*t*) was registered with NI DAQ 6115S (National Instruments, Austin, TX, USA) (data acquisition board) simultaneously in four channels with time resolution of ~10−<sup>7</sup> s. External magnetic field (*H*ext) applied in direction parallel to the plate was swept with the rate of 0.6 T/min and local magnetic induction (*B*surf) at the center of the sample surface was measured with Hall sensor. The magnetization *M* of superconducting plate was determined by the difference between local and external magnetic field induction, μ0*M* = *B*surf − μ0*H*ext. The specimen was immersed in liquid helium bath and the temperature down to 2 K was reached in pumped helium.

**Figure 1.** (**a**) Flux jump frozen in the form of "fingers"; Nb disc, *T* = 1.7 K (Figure 1a from Ref. [13]). The arrangement of pickup coils on superconducting plate for the avalanche flux registration: (**b**) in the whole plate's volume, (**c**) on three height levels, and (**d**) in magnetic stray field (*L*ext). Hall sensor line and probe for magnetic induction (*B*surf) measurement are presented in (**b**,**c**); *U*coil (*t*)—voltage on the pickup coil; CP is the central part, where the avalanches are mainly triggered.

Initial plate with the size of 20 × <sup>14</sup> × 7 mm<sup>3</sup> has been cut from extruded cylindrical rod of NbTi 50 at% alloy with diameter of 15 mm. Hot extrusion of Nb-Ti 50% alloy was carried out according to the standard technology [37] through a deformation of matrix at a temperature of 750 ◦C along the route from Ø 50 mm → Ø 15 mm (Ø—diameter of the rod) with a draw value *R* = *S*before/*S*after ≈ 11, where *S*before and *S*after are the areas of the sample sections before and after deformation, respectively. The side surfaces of the plate were ground with diamond (corundum) powder in order to remove the layers of the material strongly deformed during cutting. Such procedure was repeated at each stepwise thinning of the plate to the required thickness *d*. The pickup coils were close-wound with 0.1–0.2 mm diameter of copper wire and consisted of 15–100 turns.

#### **3. Results and Discussion**

The salient features of the manifestation of thermomagnetic avalanches in the plates at 2 K are their realization in the form of a multistage "stairway" manner of magnetic flux dynamics (Figure 2a,b). One example of the voltage *U*coil(*t*) signal inherent to a cascade induced on a pick up coil at the avalanche in trapping regime is presented in Figure 2a (2nd quadrant of hysteresis loop on Figure 2e, Jump 34). The signal consists of a series of regularly spaced impulses (portioned avalanches) and shows a multi-step "stairway" structure of magnetic flux dynamic [Φ(*t*)], received as a result of integration of voltage *U*coil(*t*) curve (Figure 2b below Figure 2a). Thus, the magnetic flux enters (in screening regime) or leaves superconductor (in trapping regime) in "steps-like" manner (Figure 2b). Such well-structured cascades are observed at the temperature of 2 K for plates of different thicknesses with *d* = 2.7–6 mm in a certain region of the magnetic field, marked with rectangles and arrows on the hysteresis loops *M*(*H*ext) (Figure 2d–f, right column). Remagnetization loop at 4.2 K for the plate with thickness of 6 mm—for comparison with *M*(*H*ext) curve at 2 K—is shown in Figure 2c.

**Figure 2.** (**a**) The voltage *U*coil(*t*) at the flux cascade in a trapping regime [Jump 34 on hysteresis loop (**e**)]; (**b**) time-integrated voltage *U*coil(*t*)—a multi-step "stairway" structure of magnetic flux Φ(*t*)); (**d**–**f**) remagnetization curves *M*(*H*ext) for the plates with the different thickness: (**d**) *d* = 2.7 mm, (**e**) 4 mm, (**f**) 6 mm at the temperature of 2 K; the ranges of cascades are shown by the rectangles and arrows on hysteresis loops. (**c**) Hysteresis loop at 4.2 K, *d* = 6 mm.

Figure 3 presents a panorama of the voltage impulses structures *U*coil(*t*) at 2 K, induced on a pick up coil at the flux cascades, originating for the plates with various thickness, i.e., with thickness of 2.7 mm (a), 3.1 mm (b), 4 mm (c), and 6 mm (d). Left side of each figure contains cascades in screening mode (1st quadrant) and the right one in trapped mode (2nd quadrant). The most regular step-like structures were found for the plates with diameter of 2.7 mm, 3.1 mm, and 4 mm. Nevertheless, the step-like structure is also visible in a narrow range of magnetic fields for the plate with a diameter of 6 mm (Figure 3d). The maximal number of steps in the cascades was found for magnetic field in the range of 2–4 T.

**Figure 3.** A panorama of the voltage impulses *U*coil(*t*) at the flux cascades for the plates with various thickness: 2.7 mm (**a**), 3.1 mm (**b**), 4 mm (**c**), and 6 mm (**d**). Each panel presents cascades recorded in four various fields and in the left column there are cascades in screening mode (1st quadrant of hysteresis loops) while in the right one—in trapped mode (2nd quadrant).

#### *3.1. Evolution of the Structure of Cascades at the Change of Magnetic Field*

*Screening regime.* Time evolution of the voltage impulses *U*coil(*t*) at *T* = 2 K at the avalanches with increasing magnetic field in screening regime, measured after cooling in zero magnetic field (ZFC—zero field cooling), for the 1st quadrant in the plate with *d* = 2.7 mm, is shown in Figure 4a. Amplitude of the signals presented was normalized to the maximal peak value *U*maxcoil for each cascade. The absolute value of voltage impulses decreases exponentially with increasing external magnetic field, as it is presented in Figure 4b. The structure of the avalanche is relatively simple for the first and for the last jump: only one single peak of voltage was found, as it is seen in Figure 4a.

**Figure 4.** (**a**) Evolution of the normalized voltage impulses structure *U*coil(*t*)/*U*maxcoil at the avalanche during the increase of magnetic field; (**b**) the amplitude of voltage impulse vs magnetic field and exponential fitting to the experimental data: *d* = 2.7 mm, *T* = 2 K; (**c**) the voltage impulse at first flux jump: *T* = 4.2 K, μ0*H*ext = 0.57 T, *d* = 4 mm. Screening regime; ZFC.

However, for the first avalanche already, the process of excitation or triggering of cascades of avalanches is observed (Figure 5a), despite that, the "oscillations" *U*(*t*) associated with this process are still weakly expressed. With increasing magnetic field, the voltage impulse expands into a shape of a crown with numerous peaks at the fourth avalanche (Figure 5b), and the step structure of the flux becomes apparent. In stronger fields, the "crown" is transformed into a structure of clearly separated, isolated avalanche impulses (Figures 3a and 4a). It means that, the whole avalanche, occurring in the period of several tens of milliseconds, exhibits a multi-step structure, i.e., the magnetic flux enters (leaves) in approximately equal, successive portions ΔΦ, the number of which can reach the value of ten (Figures 2b and 3). The range of magnetic field between 2 and 4 T may be defined as a range of steady state of cascades generation. It should be noted that duration of entering of each individual portion of the flow ΔΦ in the cascade is the same as in the single avalanches that are observed at the beginning and at the end of the region of instability of the critical state (Figure 4a). The nucleation of cascades is shifted into stronger magnetic fields with thickening the plate to *d* = 3.1 mm (Figure 5c).

It should be noted that with an increase in temperature to 4.2 K, the signal from the excitations of the cascades is practically not observed on the avalanche pulse (Figure 4c, *d* = 4 mm). At the same time, one can see, on the trailing part of the signal, a certain sequence of small peaks of the flux entry.

**Figure 5.** (**a**–**c**) The voltage impulses *U*coil(*t*) (left ordinate) and magnetic flux Φ(*t*) steps (right ordinate) at the avalanches: (**a**) first jump (μ0*H*ext = 0.54 T), (**b**) fourth jump (μ0*H*ext = 1.44 T) for the plate with thickness *d* = 2.7 mm; (**c**) fourth jump (μ0*H*ext = 1.56 T), *d* = 3.1 mm; ZFC, screening regime (1st quadrant), *T* = 2 K.

#### *3.2. The Flux Step-Like Structure in Avalanche Cascades*

*Periodicity analysis.* To study the patterns of stepped structure of the flux Φ(*t*) (Figure 2b) under the influence of external factors, its derivative *U*coil(*t*) should be analyzed, which allows to consider changes in the repetition period *T* of individual avalanche pulses (flux portions) in a cascade, or to analyze their characteristic frequencies *F* (pulses) of their repetition. For simple harmonic signals, this is, in principle the same, since the relationship between these quantities is trivial: *T* = 1/*F*. In the case of the studies performed here, it is not always easy to visually determine the period from the signal structure. This is especially true for the region of relatively weak magnetic fields −1.5 T–0–+1.5 T in the second and third quadrants.

Let us first consider the periods *T*<sup>c</sup> between peaks of voltage *U*coil(*t*) in cascades for different forms of studied signal. In Figure 6a,c and Figure 7a,c, some characteristic and their detailed analysis, for 2.7 mm and 3.1 mm samples, respectively, are shown, with the voltage impulse *U*coil(*t*) (left ordinate) inherent to a cascade at the avalanche, and a multi-step "stairway" structure of magnetic flux Φ(*t*) (right ordinate) in a screening (a) and a trapping (c) regime. Two types of signals have been selected here: in the first case (a) avalanche impulses overlap slightly (the end of the "crown", Figure 4a), in the second—Figure 6c—avalanche impulse are spaced (separated) in time. The latter case is typical for avalanches in strong magnetic field, near the boundaries of the region of instability of the critical state.

The periods *T*<sup>c</sup> between jumps in cascade vs number of impulses *n* (or steps) are shown in the main frame of Figure 6b,d, left ordinate. As follows from presented data, the screening and trapping regimes *T*c increase over time in both cases. The experimental points are scattered around the straight lines, constructed with least squares method. The values of the frequency pulse *F*<sup>T</sup> = 1/*T*av in cascades (*T*av is the average period) are given below the lines. Importantly, a similar period is observed in three quadrants for the plates of different thicknesses, in the range of steady state of cascades generation.

The increase in the *T*<sup>c</sup> period with time can be associated with an increase in the dissipation of the induction-current system in the critical state as the flux enters the superconductor. Such increase is characteristic for oscillatory with an increase in the attenuation. Structured cascades with a uniquely visually identified frequency rarely appear during magnetization reversal in the range of 0–±1.5 T, in the second and third quadrants. Here, avalanches exhibit, mainly, a complex structure. An example of the signals for a plate with a thickness of *d* = 3.1 mm, where it is difficult to visually assess the main period of *T*c and the corresponding frequency, is presented in Figure 7c. Here, spectral analysis helps to establish the frequency components of the spectrum and, if necessary, to determine the appropriate periods.

**Figure 6.** (**a**,**c**) The voltage impulse *U*coil(*t*) at the avalanche in screening (**a**) and in trapping (**c**) regime (left ordinate) and multi-step "stairway" structure of magnetic flux Φ(*t*) (right ordinate); (**b**,**d**) main frame—periods *T*c between jumps in cascade vs number of jumps *n*. The duration of separate flux impulses in the cascade is given by bottom line in (**d**). Data for the magnitude in fast Fourier transformation of signal *U*coil(*t*) are presented in the inserts (**b**,**d**); *d* = 2.7 mm, *T* = 2 K.

*Spectral analysis.* With an increase in the magnetic field, the structure of avalanche pulses passes smoothly from a single pulse with small oscillations (Figure 5a) to almost periodic oscillations limited in time (Figure 3). Proper analysis of the data requires information on the amplitude and frequency of the signal spectrum components, i.e., on the amplitude-frequency spectrum and thus, fast Fourier transformation (FFT) from the Origin program (OriginLab Corporation, Northampton, MA, USA) was utilized in the analysis of the *U*coil(*t*) spectra in the entire region of instability of the critical state. The characteristic *F*FFT frequencies in the signal spectrum, determined with this package and compared with the values of the pulse repetition frequency *F*T, obtained directly from the average period *T*av, exhibited its good efficiency (see, the inset in Figure 6b,d and Figure 7b,d—the amplitude spectra of avalanche signals with characteristic *F*FFT frequencies). The values of the pulse frequency gained from the average period *F*<sup>T</sup> for plates with a thickness of *d* = 2.7 mm and 3.1 mm, are presented in Figures 6 and 7, respectively. The spectrum component with the maximum amplitude presented in Figure 6b corresponds to the frequency *F*ma = 952 Hz, which is in good agreement with the frequency value estimated from the average period: *F*<sup>T</sup> = 927 Hz. The spectrum shown in the inset in Figure 6d is more complex. Here, there is a fundamental frequency *F*ma = 470 Hz with a maximum amplitude, which practically coincides with the frequency *F*<sup>T</sup> = 487 Hz. In addition, the spectrum contains two more components, the amplitude of which exceeds the amplitude level of 50% of the fundamental harmonic. A good agreement of spectral analysis is presented in Figure 7: *F*ma = 724 Hz and *F*<sup>T</sup> = 730 Hz (b) and *F*ma = 551 Hz and *F*<sup>T</sup> = 564 Hz (d).

**Figure 7.** (**a**,**c**) The voltage impulse *U*coil(*t*) at the avalanche in trapping (**a**) and screening (**c**) regime (left ordinate) and multi-step "stairway" structure of magnetic flux Φ(*t*) [right ordinate in (**a**)]; (**b**,**d**) magnitude data from fast Fourier transformation of signal *U*coil(*t*); *d* = 3.1 mm, *T* = 2 K.

It should be noted that the signal *U*coil(*t*) presented in Figure 7c exhibits a spectrum with three component frequencies (Figure 7d), with similar amplitude values. This is definitely due to the strong beats in the signal. Apparently, the presence of such frequency components, with close amplitude values, leads to "jumps" in the experimental dependences of the *T*<sup>c</sup> periods vs number *n* (Figure 6b,d) around the approximating straight lines. The data indicate that the harmonic *F*ma, with the maximum amplitude in the signal spectrum, is suitable for analyzing the effect of external influences and parameters of a superconducting plate on the phenomenon of multi-steps magnetic flux dynamics in the entire region of instability of the critical state.

*The results of studying signals from avalanches* for three quadrants in a plate with a thickness of 3.1 mm are shown in Figure 8 and the features of the hysteresis loop *M*(*H*ext) associated with the excitation of cascades are shown in Figure 8a. A stepwise decrease in the amplitudes of magnetization jumps Δ*M* is clearly expressed in the region of cascades for three quadrants of hysteresis loop and for all thicknesses (Figure 2d,e), which is a result of the excitation of cascades. This is evidenced by magnetic field dependence of the avalanche duration Δ*t*(*H*), shown in Figure 8b (main frame) and by their structures in the region of the amplitude step in field dependence of magnetization. The avalanches in this place [inserts into (b), jump 5 and 6] undergo threshold changes. During the sixth avalanche, the origin of multistage nature of the flux dynamics sharply increases, while the duration of the avalanche suddenly increases by three times, and the incoming flux is divided into six portions spaced in time (jump no. 6). The explanation of the step in magnetization is quite simple. Step-like behavior is associated with a change in the localization of the avalanche spot of the incoming flux and with Hall sensor method of induction measurements, where a local field change at the center of the surface is registered. During avalanche no. 5, the magnetic spot was located mainly in the central part, under the sensor. In the case of cascades (avalanche 6), the magnetic flux was dispersed like a comb of six teeth along

the entire side of the plate. In this case, the level of induction under the sensor becomes significantly smaller.

**Figure 8.** (**a**) Remagnetization curve *M*(*H*ext) of superconducting plate; the ranges of steady state cascades generation are shown by the rectangles and arrows. (**b**) Main frame—the duration of avalanches Δ*t*(*H*ext) in screening mode; insert—the structure of jump 5 and jump 6. (**c**) Main frame frequencies *F*ma with the maximal amplitude in specters of signals vs. magnetic field for different quadrants; inserts—the voltage impulses induced at avalanches in fields μ0*H*ext = 1.56 T and 0.73 T (2nd quadrant) and −0.32 T (3rd quadrant). (**d**) Main frame: the frequencies *F*ma marked by asterisk (✩) and frequencies of components with the amplitude closest in the value vs magnetic field marked by full circles (•); insert—the structure of signal in field μ0*H*ext = 2.01 T with two characteristic pronounced frequencies, 1040 Hz and 674 Hz, marked by crossed open circles (⊕), 2nd quadrant. *d* = 3.1 mm, *T* = 2 K.

Here, it is appropriate to emphasize practical implication of the importance of phenomenon of cascades at avalanche in the volume of superconductors: the incoming avalanche flux is dispersed in time and becomes non-localized in high field region. One of the main problem challenges in superconducting applications is to avoid thermomagnetic breakdown of critical state (quench) and as a result, the sudden disappearance of beneficial properties. To avoid creation of conditions for non-localized place of quench is one of the important practical issues solved by technologists today [38]. The discovery of a multi-steps magnetic flux entrance/exit at thermomagnetic avalanches in the plates of hard superconductors could make a difference in the case of applications of electric motors or generators built with superconducting elements in the form of plates. In the plates with certain thicknesses, this phenomenon reduces the risk of failure at maximum loads.

The evolution of the avalanche structure in the first quadrant (shielding mode) was previously analyzed in detail, using the example of a plate with *d* = 2.7 mm (Figures 4a and 5). Now, we will consider the properties of cascades in the second and third quadrants. The instability of the critical state is observed here in a field range from +4T to −4T (Figure 8a,c). The frequency values *F*ma for the components with the maximum amplitude in the signal

spectra are shown in Figure 8c (main frame). The *F*ma frequencies, in the range of high magnetic field ±2–±4 T (steady state cascades generation), are marked by full squares (-). Here, the frequency decreases linearly with an increase of magnetic field for all quadrants. Such behavior of the *F*ma may be interpreted as a result of increasing of dynamical dissipation of critical state and by decreasing of the avalanche start velocity as a result of the critical current density *j*c(*H*ext) decreasing.

The stable cascades with a small number of flux steps ΔΦ appear in the field region −2 T–+2 T, as it is shown, for example, in the left and in the center inserts at the top of Figure 8c. The data in the main panel, described by asterisks (✩), indicate the frequencies *F*ma. The signal *U*coil(*t*) exhibits here a complex structure with a wide spectrum (Figure 7c,d) and it seems that, in such a case, determination of the fundamental harmonic *F*ma (asterisks in Figure 8b) is very difficult, if it is possible at all. Therefore, an attempt was made to present, for each signal *U*coil(*t*), the frequencies of all spectrum components with amplitudes higher than 50% of the maximal value. The result for the second quadrant is shown in Figure 8d. Here, the frequencies *F*ma, described by asterisks (✩), as in Figure 8c, are shown, together with frequencies of the harmonics with slightly lower amplitudes, described by the full circles (•). In the field region *H* < 2 T, one can visually distinguish between certain lines, designated with the numbers 1, 2, and 3, on which most of the experimental points lie on the array. Field dependence for lines 1 and 3 is similar to that of line 2, which is a continuation of the straight line *F*c(*H*ext) from the field range of steady state generation of cascades. In terminology taken from other branches of physics, the lines 1, 2, and 3 can be described as follows: Frequency "level" 2 (line 2) for *F*ma(*H*ext) corresponds to the dynamic state of the flux with steady state generation of cascades and it splits in a magnetic fields below 2 T into two "levels", corresponding to the states of lines 1 and 3. Such behavior can be expressed as a certain degeneration of the "frequency levels" on lines 1 and 3, which occurs in a magnetic field μ0*H*ext > 2 T.

We managed to observe, for one of the cascades in the transition region (inset in Figure 8d), a signal, occurring at the boundary in *H*ext = 2.01 T, where two frequencies marked by crossed open circles (⊕) are clearly visually present: the upper *F* = 1040 Hz belongs to the straight line 1 and the lower frequency *F* = 674 Hz, which falls on line 2 for frequencies corresponding to stable cascades. Analysis of the avalanche flux dynamics in the second quadrant indicates the presence, in the region of weak fields (μ0*H*ext < 2 T), of some effective mechanism, leading to the appearance of electromagnetic radiation of a complex spectral composition. In a thinner plate with *d* = 2.7 mm, the field dependences of the frequency *F*ma, which characterize the avalanche flux dynamics look similar—Figure 9a. In the range of steady state generation of cascades, the frequency dependences *F*ma(*H*ext) for the first and third quadrants practically coincide. The data for the second quadrant also exhibit a significant "scatter-rocking" of the experimental values of *F*ma(*H*ext) and are lowered, in the same field range, in comparison with those in the first quadrant. Two measurements are presented for the second quadrant [marked by two different kind of stars (✩)] in order to show the repeatability and rocking range of experimental data in complex spectra.

The amount of magnetic fluxes in avalanches significantly differs between the first and the second quadrant (Figure 9b): in the same magnetic field much more flux enters in the shielding mode than exits in the trapping mode. Moreover, in the latter case, experimental values for the exiting fluxes are extensively scattered, as do the values of the frequencies in the dependence *F*ma(*H*ext). It should be noted that the amount of magnetic fluxes decreases exponentially with increasing external magnetic field for the 1st quadrant (Figure 9b) as well as the absolute values of voltage impulse (Figure 4b). Figure 9c shows the period *T*ma = 1/*F*ma between jumps in cascade as a function of magnetic field for various quadrants. It becomes clear here, the pulse frequency in the inductive signal depends linearly on magnetic field.

**Figure 9.** (**a**) The frequencies *F*ma with the maximal amplitude in specters of signals vs magnetic field for different quadrants. (**b**) Quantity of the magnetic flux enters (1st quadrant) and exits (2nd quadrant) at avalanches vs magnetic field; (**c**) the period *T*ma = 1/*F*ma between jumps in cascade vs magnetic field for different quadrants; *d* = 2.7 mm, *T* = 2 K. Symbol (•)—the time during which the avalanche flux front travels a distance of 3 mm in Nb and NbZr at 1.6 K and at magnetic field μ0*H*ext = 0.2 T [14].

*The role of plate thickness.* Magnetic field dependence of the jump frequency *F*ma for cascades in screening regime (1st quadrant) for the plates with various thickness *d* is shown in Figure 10. It can be seen from these data that the region of the magnetic field, where cascades are stably realized, rises up along the frequency axis when the plate becomes thinner. At the same time, the range of magnetic fields and frequencies of cascades is significantly expanded. Insert in Figure 10 shows the dependence of period 1/*F*ma between jumps in cascades as a function of plate thickness along vertical lines at μ0*H*ext = 2.5 and 3 T in the main frame. The period between avalanches inside cascades increases almost linearly with increasing plate thickness. It may indicate on the functional role of the plate thickness in triggering of cascades.

**Figure 10.** Frequency *F*ma vs. *H*ext for plates with different thickness *d*; screening regime (1st quadrant). Insert—the value of period *T*ma = 1/*F*ma along vertical lines at the magnetic field μ0*H*ext = 2.5 and 3 T on main frame vs. thickness *d* of plate.

#### *3.3. Triggering, Propagation and Localization of the Avalanche Flux Cascades*

*The place of avalanche origin.* An avalanche appears, usually, in a small volume, almost at a point compared to the size of the plate, and a bundle of vortices trapped at pinning centers can stimulate an avalanche. The process of magnetic field penetration constantly accompanies jumping of the bundles [3,4]. Simultaneous registration of the signal at different levels along the height of the plate (Figure 1c) opened possibility to find the level where the signal appears first and to estimate the propagation velocity of the disturbance along the magnetic induction line. The analysis of the time delay Δ*t*delay, between signals (more than 60 signals) from the avalanche beginning, allowed us to fix the most probable place of their origin. Statistics showed that the avalanches originate in 68% near the central region (CP, Figure 1c) of the lateral surface of the plate, and in 32% only—at its edges. This result agrees with calculation of field penetration into the plate [39]. In the central part, the screening is maximal, the magnetic field pressure is the highest, and, accordingly, the initiation of instability is most likely here. An attempt, to find side surface, where an avalanche originates (and to feel from which side a separate jump of the flux enters), using two external inductive sensors on different sides of the plate (*L*ext, Figure 1d), while recording signals, was unsuccessful. It means that, the delay (or accelerates) of the avalanche onset on one of the lateral surfaces and the asymmetry of scattering field dynamics on opposite sides of the plate, during the stepwise entry of the avalanches, were not detectable. A beautiful shaped cascade when the signals from two sensors were added and its almost complete compensation to the noise level, when they were subtracted, was observed. It should be noted, here, that, in wide coils, covering a significant part of the side surface of the plate, the stages looked to be more structured. This may be due to the fact that the origin of the cascades travels along the vertical or along the lateral surface and their appearance is somewhat different in different sections. Most probable speed of propagation of the disturbance along the direction of the magnetic field, *V*||, was determined on the basis of statistical data analysis and it was found that, it corresponds to the maximum velocity distribution function. The velocity reaches *V*||max = 0.5 *h*/Δ*t*delay = 286 ± 20 m/s, where *h* is the plate height.

*Spatial-temporal cascades behavior*. In order to establish the location of the parts of the incoming flux on the central section of the plate, the pickup coils for the avalanche flux registration were placed on superconducting plate (Figure 11a) in the whole plate's width (*W*) − (*L*), in the 1/3*W* (*L*1/3) and 2/3*W* (*L*2/3). In order to cover individual sections of the plate, a hole with a diameter of 0.7 mm was made for the placement of inductive sensors *L*1/3 and *L*2/3. Self-consistent cascade for the plate with *d* = 3.1 mm, recorded at *T* = 2 K in μ0*H*ext = 3.35 T simultaneously by three coils *L*, *L*1/3, and *L*2/3, is shown in Figure 11b. Full signal registered by coil *L* is divided proportionally between the coils *L*1/3 and *L*2/3, confirming spatial separation of impulses in the cascades along the lateral surface of the central section. Schematic view of possible spatial localization of successive flux avalanches (steps) in the form of fingers is shown on the cross-section of a plate orthogonally to magnetic field (Figure 11c). The numbers given in Figure 11b,c indicate the sequence of avalanches occurrence. The individual flux pulses are shown here in the form of "fingers" since it is known [13] that the shape of avalanche spots transforms at low temperatures from round at *T* = 4.2 K to "finger-like" at 1.8 K. "Fingers" of avalanches in the scheme (Figure 11c) cross the plate all the way to the opposite side, rather than stopping at its middle, which is in agreement with the results of experimental studies of the distribution of surface induction—signal registered by an array of the Hall probes (Figure 1b) before and after appearance of thermomagnetic avalanches in increasing and decreasing external magnetic field [15]. Here, local induction inversion as a result of a thermomagnetic avalanche over the entire thickness (*d* = 4 mm) of the studied NbTi plate occurred already at *T* = 4.2 K (Figure 3b,c from Ref. [15]), which may indicate that the avalanche front crosses the entire plate. However, an unambiguous answer can only be given by direct observation of the flux front using magneto-optics.

**Figure 11.** (**a**) The arrangement of pickup coils for the avalanche flux registration in the whole plate's width (*W*) − (*L*), in the 1/3*W* (*L*1/3) and in 2/3*W* (*L*2/3); (**b**) one of the cascades recorded simultaneously by three coils *L*, *L*1/3, and *L*2/3; *d* = 3.1 mm, *H* = 3.35 T, *T* = 2 K. Full signal registered by coil *L* is divided proportionally between the coils *L*1/3 and *L*2/3, confirming spatial separation of impulses in cascades. (**c**) A schematic cross-section of a plate oriented orthogonally to magnetic field and spatial localization of successive flux avalanches. Weak spots appear in the concave corner, near the finger's gate.

Simulation of avalanches movement in superconducting films [25,26] indicates that the avalanche front tends to displace the superconducting current from the surface to the opposite side of the sample. It means that the superconducting current is pushed out of the hot zone of the avalanche finger into the cold zone of the superconductor [40]. In shielding mode (diamagnetic induction), an appearance of "paramagnetic" profile is expected if the flux front manages to cross the middle of the sample, which takes place in the simulation. This is exactly what was observed in the experiment (Figure 3b from Ref. [15]).

#### *3.4. Possible Mechanism of Cascades Triggering*

The linear dependence of the cascades period on plate thickness suggests that the period of step structure of magnetic flux avalanche dynamics Φ(*t*) in cascades may be controlled by plate thickness as well as the limitation of path length of the front of an individual avalanche finger in the cascade, introduced by plate thickness, can lead to the appearance of a certain sequence of excitation in the induction-current system of the critical state of the superconductor.

Let us consider a possible scenario of cascades triggering. When the critical induction drop Δ*B*FJ is reached, the initial flux step (finger) in cascades can be stimulated in the superconductor critical state layer by small jumps of the flux bundles. Penetrated first finger's flux redistributes current in the sample, and thus creates weak spots for penetration of the consequent avalanche finger, etc., [40]. Weak spots appear in the concave corner (Figure 26 from Ref. [41]), near the finger's gate (Figure 11c), where the current lines looming (become thicken) occur. Subsequent avalanches in the cascade can trigger a weak spot by a locally amplified pulse of a magnetic or electric field [42] and, accordingly, a current Δ*J* due to deceleration of the previous finger of the flux reaches the opposite side of the plate. The plate in this way can show its "resonator" properties.

The expansion of the range of magnetic fields, where cascades stably arise, with decreasing plate thickness, as well as the disappearance of this phenomenon for plates with *d* > 6 mm (Figure 10), testifies in favor of the important role of finger's-type flux reflection from the opposite side (resonator properties of the plate) in the process of self-excitation.

Let us try to estimate the time it takes for the flux front to cross our plate in order to compare it with the period of multistage "stairway" structures. The dynamics of the avalanche front in conventional superconductors Nb and NbZr at a temperature of *T* = 1.4 K and in a magnetic field of 0.2 T was studied by the magneto-optical method [14] and time dependence of path traveled by the avalanche front is shown here. The magnetic flux front runs a distance of 3 mm in 0.4 ms for NbZr and in 0.6 ms for Nb, respectively, which is plotted in Figure 9c, presenting the data on the pulse periodicity in cascades for a 2.7-mm thick plate, corresponding well to the values of the periods of cascades. Hence, we can conclude that the plate, showing its "resonator" properties, can serve as a metronome that sets the period of magnetic flux dynamics in plates and lead to a multi-step "stairway" structure of flux Φ(*t*). Existence of avalanches in the form of a cascade of almost equidistant events can be unambiguously confirmed by simulating it only. However, it is quite plausible since penetrated flux redistributes current in the sample, and thus creates weak spots for penetration of the consequent avalanches, etc., [40].

#### *3.5. Impact of Magnetic Prehistory*

The change in the profile of magnetic field in the plate, with the change from increasing to decreasing field in the second quadrant, leads to a displacement of the boundary of critical state *H*<sup>2</sup> tfj stability region into weaker fields (Figure 8b), as compared to the first quadrant (*H*<sup>1</sup> tfj). The physical reason for this is a decrease in the critical current density due to the trapped flux, since the trapped flux increases, in a given external magnetic field, the average effective field in superconductor, as compared to that in shielding regime [43]. The profile of magnetic field induction can affect the spectral characteristics of the avalanche dynamics of the front flux. The magnitude of induction gradient in superconductor sets the impulse to the flux in the avalanche at the "start". The subsequent movement of the avalanche front in the plate passes along the induction relief formed before this, with the exception of the first avalanche in the screening mode. A schematic representation of Bean's magnetic induction profile *B*(*x*) in a plate corresponding to situation before first flux jumps in screening regime and for trapping regime in strong magnetic field (field ramping backward) is shown in Figure 12 [left side of (a) and (b), respectively]. The result of avalanche flux entry/exits process is presented on the right side of Figure 12a,b. Here, a cross-section of a plate perpendicular to magnetic field is shown.

*Screening regime.* Avalanches appear above the field of the first instability *H*ext > *H*1fj (Figure 12a) and the regime of cascades generation appears at the first avalanche, as it was depicted in Figure 5a. The critical current is maximal here, the pressure is high, and the flux front speed is high. Thus, a large amount of magnetic flux enters the region of the Meissner state of the superconductor (Figure 12a, right side). The central strong avalanche appeared, apparently, in the form of a finger. However, when its propagation was limited by the second side of the plate, the width of the finger under the pressure of the flow of vortices began to increase in the orthogonal direction, occupying a significant part of the plate. According to the signal structure (Figure 5a), subsequent flux fingers in the cascade were weaker. As the magnetic field increases, the critical current density decreases. Therefore, in the next cascades, the amount of flux that entered at the first pulse decreases, and to fill the plate with the flux, the number of subsequent avalanches in the cascade increases.

**Figure 12.** (**a**,**b**) Left side—a magnetic induction *B*(*x*) of Bean's profile into plate before the 1st flux jumps in screening regime (**a**) and before the cascade in trapping regime appears (**b**); right side is cross section of plate perpendicular to magnetic field direction: a schematic representation of the flux avalanche entry process, corresponding to the 1st flux jump (**a**) and to the cascade creation (**b**).

*Trapping regime.* The induction-current structure of the critical state in the flux trapping mode is fundamentally different from that in the screening mode. With shielding, only one critical current direction is present (Figure 12a, right side). In the second quadrant, due to a change from increase to decrease of external field, a second current loop in the opposite direction appears in the current structure (Figure 12b, right side). In this case, a boundary arises that separates these counter currents. In Bean model it is a straight line. In reality, due to the inhomogeneity of the pinning in the superconductor or the stochasticity of the vortex dynamics, this boundary is a jagged line [44–47]. The dynamics of this jagged (zigzag) current boundary can lead to turbulent formations along the front and circular currents (vortices). Moreover, these vortex formations can also keep the opposite direction of the circulation of currents. Therefore, the destruction of the critical state and the breaking the flux off a dome-like profile of induction can be accompanied by a dynamic "mixture" of oppositely directed electromagnetic fields, which leads to a more complex spectral composition of the induction signal than in the screening mode. This electromagnetic mixture, superimposed on the "resonator manifestations" of the plate, can strongly influence the triggering of first finger and steady state of cascades generation. As a consequence, a scatter-swing of experimental points in field dependences of the *T*c(*H*ext) periods (Figure 6b,d) and frequencies in the *F*FFT(*H*ext) spectra (Figure 8c,d) may appear.

The possible physical reasons of the complicated signal spectrum, accompanying avalanche flux dynamics, can be associated with the structure roughness of the moving flux front and inhomogeneous relief of induction. The electromagnetic radiation, accompanying the avalanche, is formed as a result of the acceleration and deceleration of individual sections of the front of moving flux, during interaction with fixed vortex bundles, followed by their separation from the pinning centers. A non-monotonic increase in the amount of flux entering the superconductor may occur and a moving flux can sweep away pinned flux hills (vortex bundles), like a tsunami, as simulation of thermomagnetic avalanches demonstrates [40]. In this case, the magnitude of the pinning inhomogeneity determines the level of electromagnetic "noise", accompanying the flux dynamics. The magnitude of inhomogeneities can decrease, with an increase in magnetic field, due to a decrease in the critical current density. In the third quadrant, when the direction of the external field changes, regions of the superconductor arise, where the magnetic field induction has the opposite direction of the field lines compared to the trapped flux. In the fields of 0–−2T, this can be strongly reflected in the excitation conditions and, accordingly, in the spectrum of avalanche pulses. In stronger fields μ0*H*ext ≥ −2 T, the induction of external field direction only remains in the plate and its distribution is similar to that in the first quadrant [43], leading to the coincidence of frequency dependences *F*ma(*H*ext) of cascades (Figure 9a) for two quadrants.

#### **4. Conclusions**

The cascades of magnetic flux avalanches at the remagnetization of superconducting bulk plates, with the magnetic flux Φ entering (screening mode) or leaving (trapping regimes) in the form of a multistage "stairway" manner with approximately equal successive portions ΔΦ, were discovered. The width of steps in temporal dependence Φ(*t*) was found to increase almost linearly with an increase of the plate thickness. Cascades in the *screening mode* after ZFC are realized in the form of fairly stable flux steps in the entire range of the instability of critical state. In this case, the number of flux steps gradually increases, with an increase in the field, reaching ten in the field range of 2–4 T. When the magnetic field is ramping backward (flux *trapping mode*) and then its orientation is changed to the opposite, the structure of the cascades becomes much more complicated.

Spectral analysis of the inductive signal allowed us to characterize the avalanche flux dynamics in the second and third quadrants in the entire region of critical state instability. It was found that, in the magnetic field range of ±2–±4 T, an avalanche is realized in a form of multistage "stairway" manner (steady state of cascades generation). In the 0–2 T region, stable cascades with small quantity of steps were rarely observed. Here, the signal spectrum consists of several spectral components with approximately equal amplitudes. Such a spectrum can arise as a result of an avalanche rearrangement of the current configuration of the critical state under the action of magnetic flux front. In the second quadrant, there are two powerful current circuits with the opposite direction of circulation. These are the shielding current and flux trapping current circuits.

The possible physical reasons of the complicated signal spectrum accompanying avalanche flux dynamics can be associated with rough structure of moving flux front and inhomogeneous relief of induction, formed as a result of previous avalanche processes. In the third quadrant, when the direction of the external field changes, the regions of superconductor arise, where the induction has the opposite direction of the field lines compared to the orientation of the trapped flux. This, in the fields of 0–−2 T, is strongly reflected in the conditions of excitation of cascades and, accordingly, in the spectrum of avalanche pulses. In stronger fields μ0*H*ext ≥ 2 T, the induction of external field direction only remains in the plate. In this region of fields, the distributions of induction in the plate in the third and first quadrants behave in the same way. Similarity in behavior is characteristic of the periodicity of cascades.

The use of various sensors allowed us to establish that, the initiation of cascades occurs mainly in the central part of the lateral surface of the plate. The speed of propagation of an avalanche front along the magnetic field line in a superconductor was experimentally determined and possible mechanism of cascades generation, related to the "resonator's properties" of plate, was suggested.

**Author Contributions:** V.C. conceived an experiment. V.C. and A.N. performed the magnetic measurements and analyzed the data. V.C. proposed a mechanism for cascades triggering. V.C., R.P. and A.N. wrote the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Acknowledgments:** We thank I. S. Aronson for fruitful discussions of the results, for the proposed possible mechanism for the occurrence of self-consistent cascades, and for the simulation of some avalanche processes. Thanks are due to S. Vasylyev for his valuable help in taking part of the data acquisition that are reported here.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **References**


### *Article* **Computation of the Electrical Resistance of a Low Current Multi-Spot Contact**

**Gideon Gwanzuwang Dankat and Laurentiu Marius Dumitran \***

Laboratory of Electrical Materials, Faculty of Electrical Engineering, University Politehnica of Bucharest, Splaiul Independentei 313, 060042 Bucharest, Romania; gdankat@elmat.pub.ro **\*** Correspondence: dumitran@elmat.pub.ro

**Abstract:** In high complexity electrical systems such as those used in the automotive industries, electric connectors play an important role. The automotive industry is gradually shifting its attention to electric cars, which means more electrical connectors for sensors and data collection. A fault in connectors for sensors used in a vehicle can cause drastic damage to capital equipment and, in the worst case, the loss of life. The studies of faults or degradation of electrical contacts are essential for safety in vehicles and various industries. Although such faults can be due to numerous factors (such as dust, humidity, mechanical vibration, etc.) and some yet to be discovered, high contact resistance is the main factor causing erratic behavior of electrical contacts. This paper presents a study on the computation of electrical contact resistance of two metal conductors (in the form of a disk) with analytical relations and a numerical computation model based on the finite element method (FEM) in COMSOL Multiphysics. The contact spots were considered to have a higher electrical resistivity value (ρcs) than those of the two metal conductors (ρCu). Studies such as the one in view that is carried out on a microscopic level are often difficult to investigate experimentally. Therefore, with the help of a simplified numerical model, the consequences of the degradation of electrical contacts are investigated. To validate the FEM model, the numerical results were compared to those obtained from analytical models.

**Keywords:** electrical contact; contact resistance; electrical connector; numerical analysis; FEM; COMSOL Multiphysics

#### **1. Introduction**

Electric contacts are one of the essential components in electronics and automotive systems. Estimating the contact resistance determines the accuracy of connectors. In machinery assembly, evaluating the reliability and tightness of contact can be done using the contact resistance of metals contact surface. For better conductivity and lower resistance, the point of contact between two bodies must have a large area with fewer impurities [1].

The growing number of electronic devices that equip electric vehicles calls for an increase in the number of electric connectors to be adapted in the vehicles' electrical system. The automotive market demands optimum safety requirements and software-enabled features at a low price; to attain this, automotive manufacturers are trying to miniaturize electrical contacts to reduce weight by manufacturing them with high-performing alloys such as CuNiSi (copper: balance; 0.8–1.8% nickel; 0.01–0.05% phosphorus; and 0.15–0.35% silicon). Additionally, aluminum or copper alloys (CuMg, CuAg) are used for smaller wire cross-sections to reduce weight [2]. An important aspect of electrical contact technology is contact housing. They are usually made of polymer and perform functions such as electrical shielding, insulation, and protection of all connector components. Contact between two conductors occurs at discrete spots owing to the rough nature of the surfaces (asperities). Therefore, the real point of contact is a smaller portion of the apparent contact area [3].

Holm [4] in the 1930s presented a well-known theory describing contact systems and contact resistance. Based on several hypotheses, Holm's model was established for the

**Citation:** Dankat, G.G.; Dumitran, L.M. Computation of the Electrical Resistance of a Low Current Multi-Spot Contact. *Materials* **2022**, *15*, 2056. https://doi.org/ 10.3390/ma15062056

Academic Editors: Marc Cretin, Sophie Tingry and Zhenghua Tang

Received: 2 February 2022 Accepted: 8 March 2022 Published: 10 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

contact resistance of a single contact spot without considering the effect of surface films. Greenwood [5] presented a more elaborate interpretation of Holm's model; he derived a formula for calculating contact resistance based on the number of a set of circular contact spots and the distance between them. Similar to Holm's model, the influence of interface films was ignored.

Over time, Holm and Greenwood's theory has been further developed and serves as a starting point for researchers in designing quality and reliable electrical contacts. Boyer [6] generalized the Greenwood formula in his analytical solution; he introduced the influence of interface films. Nakamura and Minowa [7–10] computed numerically the contact resistance using the boundary element method (BEM) and finite element method (FEM) of a contact model consisting of two cubic electrodes that communicate through square spots [8,9]. Nakamura also analyzed the contact resistance of regular and irregular conducting spots and their behavior [10]. Timsit [11] examined the electrical conduction through small constrictions and the dependence of electrical resistance on the shape and dimensions of the contact spots [12]. Shujuan et al. [13] built a contact resistance model based on a rough surface contact model by surface profile measurement and statistical analysis. Ren et al. [14] and Lim et al. [15] both analyzed contact resistance by applying a combination of experimental, numerical, and analytical methods. Malucci [16,17] established several models, considering the influence of interface films to predict the performance of degraded contacts. Furthermore, Fukuyama et al. [18] evaluated the contact resistance of a simplified contact sample and compared the result to the contact resistance calculated using Holm and Greenwood equations.

Although few studies of contact resistance focus on comparing numerical simulation and analytical solution, this paper presents a simplified contact model to evaluate the contact resistance by analytical approaches (using Holm equation and Greenwood equation) and by numerical simulation (using FEM in COMSOL Multiphysics). The model consists of two metallic disks coming in contact through multiple regular circular spots (a-spots) that have a higher electrical resistivity compared to the two metallic disks (ρcs < ρCu) in view to simulate the contact aging. Though analytical models, in theory, are considered to be the exact solution, they are derived based on multiple assumptions. Despite all the assumptions considered in the analytical models, the present study aims to explore the limits of numerical calculation for the contact resistance in the case of a simplified model. Moreover, this paper examines the influence of the size of the contact spots on the value of contact resistance.

#### *Aging of Electrical Contacts*

Electrical contacts undergo different stresses during operation (i.e., mechanical wear due to induced vibrations, atmospheric conditions, corrosion, etc.). Swingler et al. [19] classified these stresses on automotive connectors into two significant groups. One is external stresses due to the environmental conditions, while the other is internal stresses created by the vehicle. Studies on contact resistance calculations and degradation of electrical contacts can be very tasking and complex because of the contact interface's coarse nature and asperities that comprise numerous spots constricting the flow of electrons and the numerous degrading factors manifesting at the interface.

Environmental factors such as temperature, air pollution, humidity, condensation, etc., affect electrical connectors because they promote the rate of corrosion (Figure 1); Corrosion is much more severe at the interface of connectors; it gradually reduces the area of contact at the interface, which increases resistance over time and eventually leads to electrical failure of the connectors.

**Figure 1.** Examples of corroded connectors: (**a**) corroded Ac blower wire harness connector in an Acura TL [20]; (**b**) corroded wiper control relays in a BMW 7 series [21].

On the other hand, fretting due to vibration, dust, the passage of current, temperature, etc., also contributes and accelerates the degradation process of electrical contacts. Swingler et al. [22] reported that in vehicles, the connectors situated close to the exhaust system can possess high-temperature regimes such as 85, 105, 125, and 155 ◦C. Additionally, Jedrzejczyk [23] experimented on CuSn/Ni/Sn contact material system (Figure 2); he reported that during the process of degradation, electrical resistance depends on the nonoxidized area of contact at the interface. i.e., when all the areas of contact are fully covered by oxidized films, then the electrical resistance value increases above a normal threshold value.

**Figure 2.** Fretting scars as seen in an SEM micrograph obtained after 2600 fretting cycles in a CuSn/Ni/Sn contact material system: (**a**) lower sample (**b**) upper sample [23].

#### **2. Analytical Models for Contact Resistance**

Numerous studies on contact resistance and electrical contacts have made mention of Greenwood's and Holm's analytical model for contact resistance calculation [6–18]. The research reported in this paper used Holm's and Greenwood's models to calculate the contact resistance.

#### *2.1. Holm's Analytical Model for Contact Resistance*

In Holm's analytical model for the calculation of contact resistance [4], two cylinders (C1 and C2) in contact were considered (Figure 3).

**Figure 3.** Holm contact resistance model: cylinders in contact showing the apparent area of contact (*Aa*), the real area of contact (*Ac*), and the radius *R* with a circular constriction of the current lines of radius *a*.

Though both cylinders were assumed to be clean and free of any impurity, they come in contact through a smaller portion of the apparent area (*Aa*). This smaller portion is the actual point of contact (*Ac*). Constriction resistance occurs when the current is restricted from its normal flow to pass through (*Ac*). The voltage between both cylinders can be measured as the current flows through and subsequently the resistance between both surfaces. Holm's analytical model calculates the total constriction resistance of one circular spot between two electrically conducting cylinders. It has the expression:

$$R\_{\mathfrak{c}} = \mathfrak{p} \cdot (2a)^{-1} \, , \tag{1}$$

where ρ is the resistivity of the conductor and *a* is the radius of the constriction.

In situations where there is a difference in material properties (when the bodies in contact have different resistivities), Equation (1) becomes:

$$R\_{\mathbb{C}} = \left(\wp\_1 + \wp\_2\right) \cdot \left(4a\right)^{-1} \tag{2}$$

The assumptions in Holm's analytical model are:


#### *2.2. Greenwood's Analytical Model for Contact Resistance*

In 1966, J.A. Greenwood published a paper [5] in which he made a detailed interpretation of Holm's analytical model. He considered a set of circular spots (multiple a-spots) within a single cluster located at the interface between two electrodes (Figure 4). The metallic electrodes come in contact through the circular spots with no interface film between them. He obtained the formula for calculating the constriction resistance based on multiple spots within a single cluster (3) by treating the current flow between the metallic electrodes similar to that of an electrostatic charge distribution problem.

$$R\_{\rm cG} = \mathfrak{p} \cdot (2\sum a\_i)^{-1} + \mathfrak{p} \cdot \pi^{-1} \left(\sum\_{i \neq j} (a\_i a\_j) \cdot (s\_{\bar{i}\bar{j}})^{-1}\right) / (\sum a\_i)^2,\tag{3}$$

where ρ is the resistivity, s*ij* is the distance between the centers of the *i*th and *j*th spot, *ai* and *aj* is the radius of the *i*th and *j*th spot. The first term ρ·(2∑*ai*) <sup>−</sup><sup>1</sup> is the resistance of all the spots in parallel, the second term <sup>ρ</sup>·π−<sup>1</sup> (∑*<sup>i</sup>* <sup>=</sup>*<sup>j</sup>* (*aiaj*)·(*sij*) <sup>−</sup>1)/(∑*ai*) <sup>2</sup> is the resistance due to the interaction between all the spots.

**Figure 4.** Greenwood's contact resistance model showing multiple spots (*ai*—radius of the *i*th spot; *aj*—radius of the *j*th spot; s*ij* the distance between the *i*th and *j*th spot).

In cases where all the contact spots have the same size, Greenwood further approximated (3) and derived:

$$R\_{cG1} = \wp \cdot (2\sum a\_i)^{-1} + \wp \cdot (\pi n^2)^{-1} \sum\_{i \neq j} (s\_{ij})^{-1} \tag{4}$$

Furthermore, the hypothesis used in Holm's model also holds in Greenwood's analytical model.

#### **3. Numerical Model**

As shown above, the calculation of contact resistance in the case of a discontinuous interface (spot contact) can be done analytically. However, extensive multiphysics studies in which the electrical problem is coupled with thermal and/or mechanical (contact pressure) problems can only be solved numerically. Thus, in this paper, a numerical model is proposed that allows the calculation of the contact spot resistance, but can be further improved and completed later for the study of the thermal regime or the consideration of other parameters that can influence the value of the contact resistance. Therefore, to validate the proposed numerical model, a simple geometry very similar to that treated by Greenwood and consisting of two copper disks communicating through multiple contact spots (Figures 5 and 6) was considered. At the same time, such a model is indicative of the contact between two metal electrodes through a set of circular spots. To attain the simplification of the model, we considered circular spots that are equal in terms of area and relatively evenly distributed across the contact interface. At the terminals, a low direct current (200 mA) of density *J* is injected and the electric field distribution is calculated using the finite element method (FEM) in COMSOL Multiphysics software 5.6.

**Figure 5.** Schematic view of the contact configuration: (**a**) copper disks in contact; (**b**) constriction of the current lines flowing through the multiple contact spots.

**Figure 6.** Geometrical model of the problem investigated: (**a**) two copper disks in contact; (**b**) contact spots and the insulating layer of polyethylene (PE) at the interface.

#### *3.1. Geometrical Model*

The two copper disks model investigated come in contact through multiple spots (28 identical circular spots). Both copper disks have a radius *r* = 5 mm and a thickness *h* = 1 mm, while the homogenous contact spots have a radius *a* = 0.2 mm. Apart from the contact spots, the apparent contact area has an insulating layer of polyethylene 30 μm thick; it restricts the current to flow only through the multiple contact spots (Figure 6).

#### *3.2. Mathematical Model*

In the COMSOL Multiphysics software, the numerical analysis of the investigated problem was calculated in the stationary electrokinetic regime. It consists of an electromagnetic problem that calculates the distribution of a constant electric current of density *J* flowing through the multiple contact spots. The essential equations dictating the problem are the electromagnetic induction law (5), the electric conduction law (6), the electric charge conservation law (7), and the electric field strength, which was evaluated as a function of the electric potential *V* (Equation (8)).

$$\text{rot } E = 0,\tag{5}$$

$$
\mathbf{J} = \boldsymbol{\sigma} \cdot \mathbf{E} \,, \tag{6}
$$

$$\text{div}\,\mathbf{J} = \mathbf{0},\tag{7}$$

$$E = -\text{grad}\, V,\tag{8}$$

where *E* is the electric field (*V*/m) and σ = 1/ρ is the electric conductivity (S/m).

Boundary Conditions

The boundary conditions applied in this study are:

• Continuity (9): this signifies that the normal components of the injected current flowing through the copper disks are continuous and conserved across the interior boundaries of both disks:

$$m \cdot (\mathbf{J}\_1 - \mathbf{J}\_2) = 0\tag{9}$$

• Insulation (10): this specifies that no current flows across the boundary. It applies to all surfaces except the contact areas:

$$
\boldsymbol{m} \cdot \mathbf{J} = \mathbf{0} \tag{10}
$$

The discretization of the computational domain shown in Figure 7 has a total of 2,938,081 domain elements.

**Figure 7.** Discretization of the computational domain: (**a**) two copper disks in contact; (**b**) contact spots and the insulating layer of polyethylene (PE) at the interface.

#### **4. Results**

For an injected current of 200 mA, Figure 8 shows the current lines as they become constricted to flow through the multiple contact spots from one medium of the copper disk to the other.

**Figure 8.** The current line constricted at the contact spots: (**a**) two-dimensional view and (**b**) threedimensional view.

Figure 9 shows the voltage computed for an injected current of 200 mA, the resistivity of the copper disks (ρCu = 1.72·10−<sup>8</sup> <sup>Ω</sup>·m), the resistivity of the contact spots (ρCs = 1.72·10−<sup>2</sup> <sup>Ω</sup>·m), and the resistivity of the polyethylene (ρPE = 10+17 <sup>Ω</sup>·m). The result shows that for a contact spot radius of 0.1 mm, the maximum voltage is 0.12 V.

After conducting a parametric study on the radius of the contact spots from *a* = 0.1 mm to *a* = 0.5 mm, Figure 10 shows the voltage drop computed for each contact spot radius. As expected, the voltage decrease as the contact spot radius increases. The results show that when the contact spot radius is 0.1 mm, the voltage is 0.12 V; for a radius of 0.3 mm, the voltage is 12.8 mV; and when the radius is 0.5 mm, the voltage becomes 4.6 mV.

**Figure 9.** Computed voltage for contact spot radius of 0.1 mm: (**a**) two-dimensional view and (**b**) three-dimensional view.

**Figure 10.** Voltage-drop for different values of contact spot radius *a* = (0.1 mm/0.5 mm).

Figure 11 shows the contact resistance evaluated in the case where for resistivity of the copper disks is the same as that of the contact spots <sup>ρ</sup>Cu <sup>=</sup> <sup>ρ</sup>Cs = 1.72·10−<sup>8</sup> <sup>Ω</sup>·m. The result shows that the contact resistance decreases from 3.25·10−<sup>6</sup> to 4.34·10−<sup>7</sup> <sup>Ω</sup> (numerical), 3.07·10−<sup>6</sup> to 1.02·10−<sup>6</sup> <sup>Ω</sup> (Holms), and 4.76·10−<sup>6</sup> to 2.71·10−<sup>6</sup> <sup>Ω</sup> (Greenwood) when the contact spot radius increases from 0.1 to 0.3 mm. When the contact spot radius increases to 0.5 mm, the contact resistance decreases to 3.47·10−<sup>7</sup> <sup>Ω</sup> (numerical), 6.14·10−<sup>7</sup> <sup>Ω</sup> (Holm), and, respectively, 2.3·10−<sup>6</sup> <sup>Ω</sup> (Greenwood).

Figure 12 shows the contact resistance calculated for resistivity of the copper disks <sup>ρ</sup>Cu = 1.72·10−<sup>8</sup> <sup>Ω</sup>·m, with the resistivity of the contact spots <sup>ρ</sup>Cs = 1.72·10−<sup>6</sup> <sup>Ω</sup>·m (a) and <sup>ρ</sup>Cs = 1.72·10−<sup>4</sup> <sup>Ω</sup>·m (b). The result from Figure 12a shows that the contact resistance decreases from 6.03·10−<sup>5</sup> to 7.03·10−<sup>6</sup> <sup>Ω</sup> (numerical), 3.07·10−<sup>4</sup> to 1.02·10−<sup>4</sup> <sup>Ω</sup> (Holms), and 4.76·10−<sup>4</sup> to 2.71·10−<sup>4</sup> <sup>Ω</sup> (Greenwood) when the contact spot radius increases from 0.1 to 0.3 mm. When the contact spot radius increases to 0.5 mm, the contact resistance decreases to 2.58·10−<sup>6</sup> <sup>Ω</sup> (numerical), 6.14·10−<sup>5</sup> <sup>Ω</sup> (Holm), and 2.30·10−<sup>4</sup> <sup>Ω</sup> (Greenwood).

**Figure 11.** Contact resistance showing values for Holm and Greenwood (analytical) and numerical analysis calculated for (ρCu <sup>=</sup> <sup>ρ</sup>Cs = 1.72·10−<sup>8</sup> <sup>Ω</sup>·m).

**Figure 12.** Contact resistance showing values for Holm and Greenwood (analytical) and numerical analysis calculated for <sup>ρ</sup>Cu = 1.72·10−<sup>8</sup> <sup>Ω</sup>·m with (**a**) <sup>ρ</sup>Cs = 1.72·10−<sup>6</sup> <sup>Ω</sup>·m and (**b**) <sup>ρ</sup>Cs = 1.72·10−<sup>4</sup> <sup>Ω</sup>·m.

In Figure 12b, when the contact spot increases from 0.1 to 0.3 mm, the contact resistance decreases from 6·10−<sup>3</sup> to 6.4·10−<sup>4</sup> <sup>Ω</sup> (numerical), 3.07·10−<sup>2</sup> to 1.02·10−<sup>2</sup> <sup>Ω</sup> (Holm), and 4.76·10−<sup>2</sup> to 2.71·10−<sup>2</sup> <sup>Ω</sup> (Greenwood). At a contact spot radius of 0.5 mm, the contact resistance becomes 2.30·10−<sup>4</sup> <sup>Ω</sup> (numerical), 6.14·10−<sup>3</sup> <sup>Ω</sup> (Holm), and, respectively, 2.30·10−<sup>2</sup> <sup>Ω</sup> (Greenwood).

Figure 13 shows the contact resistance obtained for resistivity of the copper disks <sup>ρ</sup>Cu = 1.72·10−<sup>8</sup> <sup>Ω</sup>·m with a resistivity of the contact spots <sup>ρ</sup>Cs = 1.72·10+2 <sup>Ω</sup>·m (a) and <sup>ρ</sup>Cs = 1.72·10+4 <sup>Ω</sup>·m (b). In Figure 13a, the result shows that for an increase in contact radius from 0.1 to 0.3 mm, there is a decrease in the value of the contact resistance from 5.25·10+3 to 6.39·10+2 <sup>Ω</sup> (numerical), 3.07·10+4 to 1.02·10+4 <sup>Ω</sup> (Holms), and 4.76·10+4 to 2.71·10+4 <sup>Ω</sup> (Greenwood). For a contact radius of 0.5 mm, the contact resistance becomes 2.30·10+2 <sup>Ω</sup> (numerical), 6.14·10+3 <sup>Ω</sup> (Holm), and, respectively, 2.30·10+4 <sup>Ω</sup> (Greenwood).

**Figure 13.** Contact resistance showing values for Holm and Greenwood models and numerical analysis calculated for <sup>ρ</sup>Cu = 1.72·10−<sup>8</sup> <sup>Ω</sup>·m, <sup>ρ</sup>Cs = 1.72·10+2 <sup>Ω</sup>·m (**a**), and <sup>ρ</sup>Cs = 1.72·10+4 <sup>Ω</sup>·m (**b**).

The result in Figure 13b indicates that when the contact radius increases from 0.1 to 0.3 mm, the contact resistance value decreases from 5.75·10+5 to 6.4·10+4 <sup>Ω</sup> (numerical), 3.07·10+6 to 1.02·10+6 <sup>Ω</sup> (Holms), and 4.76·10+6 to 2.71·10+6 <sup>Ω</sup> (Greenwood). For a contact radius of 0.5 mm, the contact resistance becomes 2.3·10+4 <sup>Ω</sup> (numerical), 6.14·10+5 <sup>Ω</sup> (Holm), and, respectively, 2.30·10+6 <sup>Ω</sup> (Greenwood).

Figure 14 presents the contact resistance values obtained only by numerical analysis in a case where the resistivity of the copper disks remains <sup>ρ</sup>Cu = 1.72·10−<sup>8</sup> <sup>Ω</sup>·m and the 28 contact spots all have resistivity of <sup>ρ</sup>Cs = 1.72·10−<sup>6</sup> <sup>Ω</sup>·m (-), <sup>ρ</sup>Cs = 1.72·10−<sup>4</sup> <sup>Ω</sup>·m (•), <sup>ρ</sup>Cs = 1.72·10−<sup>2</sup> <sup>Ω</sup>·m (), and when the 28 contact spots have different resistivities in the range from 1.72·10−<sup>6</sup> to 1.72·10+6 <sup>Ω</sup>·m (). The resistivity was imposed evenly, i.e., (4 contact spots <sup>ρ</sup>Cs1 = 1.72·10−<sup>6</sup> <sup>Ω</sup>·m, 4 contact spots <sup>ρ</sup>Cs2 = 1.72·10−<sup>4</sup> <sup>Ω</sup>·m, 4 contact spots <sup>ρ</sup>Cs3 = 1.72·10−<sup>2</sup> <sup>Ω</sup>·m, 4 contact spots <sup>ρ</sup>Cs4 = 1.72 <sup>Ω</sup>·m, 4 contact spots <sup>ρ</sup>Cs5 = 1.72·10+2 <sup>Ω</sup>·m, 4 contact spots <sup>ρ</sup>Cs6 = 1.72·10+4 <sup>Ω</sup>·m, and, respectively, 4 contact spots <sup>ρ</sup>Cs7 = 1.72·10+6 <sup>Ω</sup>·m). The results show that the contact resistance value varies between 0.6 <sup>Ω</sup> (ρCs = 1.72·10−<sup>2</sup> <sup>Ω</sup>·<sup>m</sup> at *<sup>a</sup>* = 0.1 mm) and 2.58·10−<sup>6</sup> <sup>Ω</sup> (ρCs = 1.72·10−<sup>2</sup> <sup>Ω</sup>·m at *<sup>a</sup>* = 0.5 mm).

**Figure 14.** Contact resistance showing values for numerical model calculated for <sup>ρ</sup>Cu = 1.72·10−<sup>8</sup> <sup>Ω</sup>·<sup>m</sup> with <sup>ρ</sup>Cs = 1.72·10−<sup>6</sup> <sup>Ω</sup>·m; <sup>ρ</sup>Cs = 1.72·10−<sup>4</sup> <sup>Ω</sup>·m; <sup>ρ</sup>Cs = 1.72·10−<sup>2</sup> <sup>Ω</sup>·m, and, respectively, <sup>ρ</sup>Cs = (1.72·10−<sup>6</sup> <sup>Ω</sup>·m/1.72·10+6 <sup>Ω</sup>·m).

#### **5. Discussion**

As seen in Figures 11–13, it is clear that the contact resistance decreases as the contact spot radius increases in all cases; a similar outcome was reported by Fukuyama et al. [18]. He experimented on a configuration consisting of samples with SiO2 insulator layers that have gold (Au) contact spots distributed on the circumference of the electrodes. The numerical and analytical results obtained in this study are consistent with the results presented in [18] in terms of the variation of the contact resistance as a function of the contact spots area. However, comparing the numerical calculation of contact resistance to that of the analytical calculation results (Equations (1) and (4)), it is clear that the numerical results are below the two analytical results. Concerning the data presented in [18], it can be noted that the experimental results were between the two analytical calculations (Holm and Greenwood). Additionally, from the point of view of the relationship between contact resistance and the area of contact, similar variation results were reported by Shujuan et al. [13], Ren et al. [14], Lim et al. [15], and Su et al. [24].

In Figure 14, the contact spots with a high resistivity showed a high value of contact resistance and vice versa, as expected. In practice, this is typical of aged electrical contacts where on the one hand, the parts at the interface with a large concentration of contaminants possess high resistivity to the flow of current and thus, have high contact resistance. On the other hand, the parts with fewer/or no contaminants have low contact resistance. Shibata et al. [25] performed a detailed analysis of the contact distribution of fretting corrosion on samples with low phosphorus–bronze alloy covered with tin plating. He reported that the contact resistance distribution results of fretting corrosion wear as observed by scanning electron microscopy (SEM) and energy dispersive X-ray (EDX) analysis indicated that the parts with few oxide deposit layers showed low contact resistance, while those with a large oxide deposit layer showed high contact resistance. It must be noted that both Holm and Greenwood have some drawbacks. In reality, the actual connector consists of multiple small contact areas (a large number of a-spots), but Holm assumed these a-spots to be equivalent to one solid circular contact. On the other hand, Greenwood's model is suitable for multiple spots within a cluster, and it depends on the distances between the set of the circular spots within that cluster. Though this is an upgrade from Holm's model, it will be difficult to calculate the constriction resistance when the contact consists of more than one cluster [26].

#### **6. Conclusions**

This paper presents a study on contact resistance of a model of two metallic disks made of copper that communicate via 28 homogenous circular spots. In the first aspect of the study, the contact resistance calculation was done using Holm's and Greenwood's analytical model, and in the second aspect, using a numerical model in COMSOL Multiphysics. The study was carried out keeping the number of contact spots and the electrical resistivity of the copper disks constant and changing the radius and resistivity of the contact spots as a possible aging consequence, between 0.1 and 0.5 mm and 1.72·10−<sup>6</sup> to 1.72·10+4 <sup>Ω</sup>·m, respectively. As expected, the resulting values of contact resistance in all cases indicate that the contact resistance decreases as the contact spot radius increases; it validates that the contact resistance decreases as the effective contact area increases. This means that the interface contains few weaknesses and the contact pressure between the two conductors takes appropriate values. Additionally, the results show that the value of contact resistance obtained for both analytical models is higher with at least an order of magnitude than that of the value obtained by numerical analysis in all cases.

Furthermore, without considering the model used or the contact spot radius, as the resistivity of the contact spot increases in all cases (Holm, Greenwood, and numerical), the contact resistance value obtained increases by two orders of magnitude. In the future, detailed work on the electrothermal coupling of the model will be done to help understand the thermal regime and its effect on contacts.

**Author Contributions:** Conceptualization, L.M.D.; methodology, G.G.D.; software, G.G.D.; validation, L.M.D. and G.G.D.; formal analysis, G.G.D. and L.M.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Molecular Dynamics Simulation of Sintering Densification of Multi-Scale Silver Layer**

**Peijie Liang 1, Zhiliang Pan 1, Liang Tang 1, Guoqi Zhang 1,2, Daoguo Yang 1, Siliang He 1,3,\* and Haidong Yan 1,4,\***

	- <sup>4</sup> Hang Zhou Global Scientific and Technological Innovation Center, Zhejiang University, Hangzhou 311200, China
	- **\*** Correspondence: siliang\_he@guet.edu.cn (S.H.); haidong\_yan@zju.edu.cn (H.Y.)

**Abstract:** Based on molecular dynamics (MD), in this study, a model was established to simulate the initial coating morphology of silver paste by using a random algorithm, and the effects of different sizes of particles on sintering porosity were also analyzed. The MD result reveals that compared with the sintering process using large-scale silver particles, the sintering process using multi-scale silver particles would enhance the densification under the same sintering conditions, which authenticates the feasibility of adding small silver particles to large-scale silver particles in theory. In addition, to further verify the feasibility of the multi-scale sintering, a semi in-situ observation was prepared for a sintering experiment using micro-nano multi-scale silver paste. The feasibility of multi-scale silver sintering is proved by theoretical and experimental means, which can provide a meaningful reference for optimizing the sintering process and the preparation of silver paste for die-attach in powering electronics industry. In addition, it is hoped that better progress can be made on this basis in the future.

**Keywords:** molecular dynamics; initial coating morphology; random algorithm; semi in-situ observation

#### **1. Introduction**

Sintering using silver particles is considered to be one of the ideal bonding techniques for the wide bandgap semiconductors, e.g., silicon carbide (SiC) and gallium nitride (GaN) devices, because this technique is capble of forming a high-melting-point sintered silver joint with sintering temperatures lower than 300 °C [1–4].

Nanosilver paste is considered an ideal bonding material since it can sinter at a temperature lower than 250 °C without pressure and operate at a temperature higher than 400 °C [5–8]. However, the difficulty of preparation and preservation of nanosilver paste results in a high cost in the sintering process. With the continuous improvement of power performance requirements of power electronic equipment, realizing the maximum commercial value of sintered silver is no longer limited to optimizing the sintering process, sintered joint performance, and energy consumption. Rather, optimizing the formula of silver paste is also very important There has been a rapid increase in the demand for power electronics devices because of the strong growth of electric vehicles. In order to maximize commercial value, while considering the sintering behaviors, mechanical properties of joint, energy consumption, and safety of silver particles, researchers are focusing on optimizing the silver particle formula in silver particles paste to achieve compromise solutions. Figure 1a shows the scanning electron microscope (SEM) image of

**Citation:** Liang, P.; Pan, Z.; Tang, L.; Zhang, G.; Yang, D.; He, S.; Yan, H. Molecular Dynamics Simulation of Sintering Densification of Multi-Scale Silver Layer. *Materials* **2022**, *15*, 2232. https://doi.org/10.3390/ ma15062232

Academic Editors: Marc Cretin, Zhenghua Tang and Sophie Tingry

Received: 28 January 2022 Accepted: 15 March 2022 Published: 17 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

silver nanoparticles after sintering without pressure. Sintering neckings formed between the nano particles can be clearly observed due to the easy diffusion between the adjacent particles [9]. However, the rapid diffusion of silver nanoparticles formed a sintered layer with high porosity during the sintering process. Figure 1b shows the SEM image of silver micron particles sintered under the same sintering conditions as Figure 1a. Compared with Figure 1a, the necking between micron particles is few and inconspicuous. It is difficult to form diffusion necks in the pressureless sintering process due to the low surface energy of large particles. Therefore, applying 10–30 MPa or even higher assist-pressure during the sintering process was used to obtain a reliable micron silver sintered joint [10,11].

**Figure 1.** The SEM image of Ag sintered layer with (**a**) nano-silver and (**b**) micro-silver by Kazohiko Suganuma et al., adapted with permission from Ref. [9], 2022 Chuantong Chen; (**c**,**d**) are the sintering morphology of multi-scale silver at different sintering temperatures by H D. Yan et al., adapted with permission from Ref. [12]. 2022 Haidong Yan; (**e**) super large silver particles and micron silver particles formed silver neck after sintering.

Figure 1a,b showed the SEM of different particles after pressureless sintering at 250 °C for 60 min by Suganuma's team from Osaka University. It demonstrates the advantages and disadvantages of different scale granular silver sintered solder joints, which are verified by using the experimental method of controlling variables. Based on this study, a pressureless sintering process was proposed by adding the silver nano particles in the silver micro paste to improve the sintering driving force and reduce the sintering temperature, as shown in Figure 1c,d [12]. In our previous study, we were surprised to find that, after pressureless sintering at 250 °C for 60 min, a particle with a diameter of 3.6 μm can form a sintering necking with a silver particle with a diameter of 27.6 μm, as shown in Figure 1e. These results prove the feasibility of low cost nano-micron multi-scale silver paste which prepared by appropriate addition of nanoparticles in low-cost micron silver paste as the potential bonding materials for sintering. The reliability of multi-scale silver sintering as a bonding layer has been clarified in many previous research works [12–15]. However, in past demonstrations of multi-scale silver sintering, most of the experiments are based on the densification analysis at the macro level, and the densification analysis at the micro molecular level is not explained. The initial dynamic diffusion evolution of sintering using multi-scale silver particles, which is important to reveal the sintering mechanism, is not clear and has thus far been difficult to observe. In order to further reveal the micro mechanism of mixed scale particle sintering, we use molecular dynamics simulation to reveal the difference of diffusion mechanism between mixed particles and large particles.

MD is a multidisciplinary technology that relies on Newtonian mechanics and integrates physics, chemistry, and mathematics to simulate the molecular system's motion, which is also widely used in studying the diffusion mechanism of silver sintering [16]. In this study to simulate more authenticity, the sintering of different sizes of particles for MD analysis was prepared with a random algorithm. It is interesting to note that using the random algorithm could well simulate the randomness of particles stacking and arrangement in the silver paste layer, which will be described in the next chapter. After simulation, to further verify the feasibility of the multi-scale sintering, a semi in-situ observation was prepared for a sintering experiment using micro-nano multi-scale silver paste.

#### **2. Method**

#### *2.1. Methodology of Simulation*

#### 2.1.1. Molecular Dynamics Interatomic Potential

Sintering technology forms a dense joint based on the particles' diffusion below the melting temperature, according to the embedded atom method (EAM) [17]. Using the Large-Scale Atomic Molecular Massively Parallel Simulator (LAMMPS), this study combined the EAM to Simulate the diffusion between particles with different scales at different temperatures. The simulation atom-style was selected as atomic because sintering bonding is a diffusion process between particles. The timestep unit is related to the selected atom-style; In this study, the timestep size is set to the default values 0.001 ps to steadily solve Newton's equation of motion.

The molecular dynamics model is composed of atoms arranged in a specific order, and it could cause geometric distortion if the simulation's model size is too small. In order to reduce the model distortion, we had established the minimum particle radius of 10 Å in this study. The three style particles with radii of 10 Å, 30 Å, and 50 Å were designed by the "equal interval" method. As in past studies [18,19], the model's crystal orientation could affect the simulation results. To make the simulation more closely fit the actual situation, the particles in the simulation box as described in the next section were set at different crystal orientations by random placement.

#### 2.1.2. Random Algorithm

Matlab is a matrix-based language that can provide a most natural expression of computational mathematics and provides researchers with the required algorithm programming for designing and analyzing the data [20].

In order to simulate a random mix like the particles of the silver paste coating on the substrate in the real application, Matlab used mathematical analysis and computation to design a simulation box with the random algorithm. The circles with different radius are randomly arranged in proportion according to the random function in the box. However, it is difficult to avoid intersecting and overlapping circles because generating random circles generates random dots and draws circles. Thus, in this step, we design the circles with a proposed ratio by the random function and remove the circles which intersect or overlap, and, finally, output the random circle's position and size which had been retained.

In this study, we had designed a Half-two-dimensional (H2D) box of 500 × 6 × 300 Å along *X*, *Y*, and *Z* directions with the periodic boundary to simulation the densification process of one cross-section of the silver layer to reduce simulation redundancy, shorten the running time, and achieve a more precise understanding of the densification process In order to understand the effect of multiscale on densification, the following three groups of different mass ratio groups are set in this study as shown in Table 1.

**Table 1.** Grouping with different mass ratios.


However, the mass ratio in the grouping must transform into the corresponding particle ratio before modeling. In the H2D simulation in this study, the conversion of mass ratio and particle ratio is associated with cross-sectional particle area, and the around ratios of the corresponding particles are shown in Table 2.


**Table 2.** Particle ratio corresponding to different mass ratio (base on H2D).

The particle ratio was substituted into the random equation in MatLab, and we obtained the random arrangement of various circles shown in Figure 2. Position coordinates of the circle generated in the random function are extracted and substituted into the MD model. In order to calculate the porosity easily and make the simulation the silver paste coating on the substrate, we designed the upper and lower silver plates in the *Z* direction.

**Figure 2.** The circles' random placement by MatLab with diameter ratio from Table 2 with (**a**) group A, (**c**) group B, and (**e**) group C. (**b**,**d**,**f**) were the MD model based on randomly placed data to (**a**,**c**,**e**), respectively.

#### *2.2. Semi In-Situ Observation Sintering of Multi-Scale Silver*

In-situ sintering observation is a method to observe the diffusion bonding process of silver particles by scanning electron microscope without resampling during the sintering, which observes the morphology change of the particles better than the traditional subsection sampling method [21,22].

In this study, we used a scanning electron microscope equipped with a protochips single tilt heating holder with a maximum heating temperature of 1200 °C to observe the multi-scale silver paste sintering semi in-situ.

During the sintering, in order to avoid electron beam radiation, only open the gun valve when sampling the sample [23]. We referred to the micron silver sintering profile by K. S. Siow [24], and we set the sintering temperature to 250 °C for 10 min on the pressureless setting.

#### **3. Result and Discussion**

#### *3.1. Sintering Densification of Silver Particles with Different Proportions*

Figure 3 shows the morphology of the sintering densification process of groups A, B, and C sintering at 300 °C. Comparing the sintering processes of A and B, it is obvious to find that group B with small particles is better than group A in the reduction of porosity, the reduction of which is more significant for group C with smaller silver particles based on group B. In the initial stage, driven by thermal energy, the particles in the sintering layer begin to "move" and "attract" the particles around them, and begin to gather to prepare for the start of sintering; which also called the first stage of sintering. With the continuous sintering, the grain boundary energy of particles begins to weaken and the particles begin to diffuse each other, which is a longer process (10–1000 ps). However, with the subsequent sintering, the change of sintering pore decrease becomes weaker, which is more clearly reflected in the change curve of porosity of groups in Figure 4. Multiscale group C not only shows lower porosity than group A and B after sintering but also reduces the area of pore due to its complementarity in the gap of large particles. In the actual production process, the large pore will lead to uneven heat dissipation, resulting in heat concentration and local thermal resistance increasing, which will also affect the reliability of packaging.

**Figure 3.** Morphology of sintering densification process of (**a**) group A, (**b**) B, and (**c**) C at 0 ps, 10 ps, 500 ps, and 1000 ps.

Figure 4a shows the changing profiles of different particle ratios' porosity and Figure 4b shows the change rate of porosity. The initial stage is the early stage of sintering, so this stage has the maximum driving force in the whole pressureless sintering process. It shows that the porosity decreases rapidly which is shown in Figure 4a, and the diffusion rate is the fastest in the whole sintering, as shown in Figure 4b. With the continuous sintering, the porosity changes gradually from the initial sharp change to stability. We also found that the larger the particle, the earlier it tends to be stable; see Figure 4b. It is also found

from Figure 3 that the larger pores which are found at 500 ps hardly change after 1000 ps, while the smaller pores will be eliminated or further reduced. This is because the bonding completed early greatly weakens the surface energy in the system, which leads to the application temperature being insufficient to provide sufficient diffusion driving force. Porosity slight shrinkage leads to the diffusion rate slowing down.

Combined with the curves from Figure 4, we intuitively find that the porosity and porosity change rate of the simulation group with small particles is better than without the small particles. The relative surface energy is inversely proportional to the particle size, and adding small particles to large particles can improve the sintering driving force. The relationship between particle size and corresponding drive force follows [16]:

$$
\sigma = \gamma \left( \frac{1}{R\_1} + \frac{1}{R\_2} \right) \tag{1}
$$

where *σ<sup>i</sup>* is the sintering driving force, *σ<sup>i</sup>* is the material surface energy, and *R*<sup>1</sup> and *R*<sup>2</sup> are the radii of two particles that diffuse each other. It can be seen from the graph in Figure 4b that the maximum sintering diffusion change rate in the initial stage of mixed particle sintering group C can reach 0.67, which is twice that of group A with only large particles reaching 0.32.

**Figure 4.** The profile of (**a**) porosity change with sintering time and corresponding (**b**) change rate of porosity.

From the above results, with the same sintering temperature and time, multi-scale silver is better than large particle sintering in sintering porosity. Moreover, the small particles in multi-scale silver act as a "filler", which also greatly reduces the average volume of pores. In order to further understand the diffusion mechanism of silver particles with different sizes in the sintering process, different silver particles' change of necking was simulated and analyzed. This will be described in the next section.

#### *3.2. Evolution of Sintering Silver Neck Size*

As shown in Figure 5 the bonding of the following three different size particles pair was studied. The particles pair in the model set different crystal orientations to avoid the contingency of the same crystal orientation. In this chapter, we studied the necking of different particles' radii. We cut 474, 13,240 and 61,370 silver atoms from a large silver FCC supercell to form silver particle pairs with radii of 10 Å, 30 Å and 50 Å respectively as shown in Figure 5, the initial distance of 2 Å between particle pair. As shown in Figure 4, the porosity in the sintering process decreases rapidly within the first 200 ps, and then the densification slows down. Therefore, in this chapter, the simulation mainly takes the first 300 ps for study.

**Figure 5.** The Silver particles pair with different radius: (**a**) 10 Å, (**b**) 30 Å, and (**c**) 50 Å.

Figure 6 shows the simulated value and the simulated value's fitting curve of *x*/*r* (*x* is the radius of diffusion neck and *r* in the radius diffusion particle.) of different particles' radii during the sintering. As shown in Figure 6, the necking of particles with radius of 10 Å grows and diffuses rapidly, and the *x*/*r* reaches 1 at about 8 ps (the necking was completed). Figure 7 shows the profile of total potential energy during the necking of particles with a radius of 10 Å. Two particles approach and begin to diffuse under the driving force provided by thermal energy. The relative potential energy of the system decreases to provide diffusion kinetic energy. When the diffusion necking is completed, the system's kinetic energy begins to be constant, the relative potential energy stops decreasing and remains stable. As the thermal energy in the system continues to be maintained, some atoms will cross the potential barrier and further diffuse, and potential energy will also change, as shown in the potential energy fluctuation at 250 ps in Figure 7.

**Figure 6.** The measurement data and fitting curve of *x*/*r* with the time of silver particles pair with different radius sintered at 300 °C. The fitted curve conforms to the growth of the logarithmic function, which is consistent with the rapid change to slow the change of measured data.

However, under the same simulation conditions, the growth rate of particles with radii of 30 Å and 50 Å is slower, and the *x*/*r* with the radius of 30 Å tends to be stable when it increases to about 0.7, and the radius of 50 Å tends to be stable when it increases to about 0.52. This difference further indicates that small particles can diffuse more easily than large ones at the same temperature This difference further indicated that small particles are more easily diffused than large particles at the same temperature. The process of necking size accords with the diffusion formula [25]:

$$(\frac{\chi}{r})^n = \frac{Bt}{(2r)^m} \tag{2}$$

where *t* is the sintering time, *n* and *m* are constants dependent on the specific transport mechanism, and *B* is a term made of material and geometric constants. Combined with Equation (2) and the trend of measured data in Figure 6, the radii of 30 Å and 50 Å will continue to increase until the maximum values are reached if the simulation time or sintering temperature increases.

**Figure 7.** Potential energy with time during sintering of particle pair with a radius of 10 Å.

Small silver particles were doped into large silver particles to shorten the sintering bonding time and increase the sintering neck of large particles in this study, as shown in Figure 8. Figure 8a shows the sectional of 50 Å particles dope with 10 Å particles. In order to observe the diffusion process, the particles were divided into different groups and colored, respectively. As driven by thermal energy, the small particles approach two large particles and diffuses at 1 ps. With the sintering continuing, small particles diffuse into and fill the bonding gap of two large particles, and the diameter of the bonding neck reached 67.8 Å at 10 ps. However, under the same conditions, the diffusion diameter of only large particles is only 28.69 Å. With continued sintering, small particles continue to "climb" along with two large particles, further forming bonds with large particles. At the end of sintering, the neck of the sintering group with small particles shrinks to 66.46 Å, which is larger than 53.66 Å sintering with only large pair particles.

**Figure 8.** MD profile of particle pair with the radius of 50 Å; (**a**) doping the silver particles with radius of 10 Å, and (**b**) control group without any doping.

Small particles can "close" two large particles far away to form necking to increase the density of sintering, as shown in Figure 9. The initial distance of the two large particles in Figure 9a,b is 10 Å and the difference between them is that Figure 9a had added some small particles in the middle and the lower part between two large particles. Figure 9a shows that large particles pair get close together and form necking due to the pulling of small particles' diffusion at the bottom. This is because small particles provide a certain pulling force for large particles to bond. Driven by this pulling force, the particles can close

together and form necking. This can also be confirmed from the particle pair in Figure 9b which, without small particles, always keep the distance during the simulation process.

From the above results, whether from the simulation results of silver layer sintering or silver neck diffusion analysis, the evidence points to the superior characteristics of multi-scale silver sintering. From the perspective of molecular dynamics, the small silver particles of multi-scale silver do have a good densification effect on the sintering of large particles. However, this micro analysis does lack some support in practical application. Considering this, in order to further verify this feasibility, we will carry out semi in-situ sintering of multi-scale silver sintering in the next chapter.

#### *3.3. Result of Semi In-Situ Observation*

From Figure 10, the diffusion bonding process of Ag particles during the sintering by semi in-situ SEM is clearly shown. From the in-situ sintering diagram, particle pair B which is closed in the initial stage begins to form an obvious diffusion bonding neck when the sintering is carried out for 6 mi, as shown in Figure 10b. In the subsequent sintering, the diffusion silver neck of the B particle pair continues to grow. At the same time, nanosilver particles around the particle pair continue to diffuse to the large particles, also resulting in the also growth of the large particle pair while bonding. The bonding reaction between micron particles and nanoparticles in the A particle pair increase the feasibility of multi-scale sintered silver. For the A particle pair, we can see that the two particles are not in contact before sintering. However, after 8 min sintering, the particle pairs are driven by the nanoparticles and form a necking connection. Also as the particle pair which initial distance was 2.06 μm is shows in Figure 10a, and shorten to 0.6 μm after sintered 8 min. At the same time, the radius of the particle pair also increases from the initial 0.53 μm to 1.14 μm. This sintering change is consistent with the simulation results in Figure 9, which shows that the small silver particles of multi-scale silver paste can provide good sintering driving force in large particles, which can provide better density.

In this study, the micro silver paste with nano-silver is sintered for 8 min without pressure to form a densification bonding, which is much shorter than the traditional pure micro silver paste. Multiscale silver sintering is shortened not only the sintering time but also reduces the sintering temperature, which is also one of our future research directions.

**Figure 10.** The SEM image of In-Situ multiscale silver sintering at different time steps (**a**) dried, (**b**) after sintering 6 min, (**c**) after sintering 7 min, and (**d**) after sintering 8 min.

#### **4. Conclusions**

Using molecular dynamics, combined with MatLab random arrangement,he sintering densification mechanism of silver with different mixed particle ratios has been studied, and the neck growth mechanism of a silver particle pair in the sintering process was analyzed. Combined with the validation experiment, the following conclusions can be summarized:


**Author Contributions:** Conceptualization, P.L. and H.Y.; methodology, G.Z., H.Y., D.Y., S.H. and Z.P.; software, P.L., Z.P. and L.T.; formal analysis, P.L., H.Y. and Z.P.; experiment, H.Y., L.T. and S.H.; data curation, P.L., H.Y., S.H., L.T. and; writing, P.L., S.H. and H.Y.; supervision, H.Y., S.H. and D.Y.; project administration, H.Y., S.H. and D.Y.; funding acquisition, H.Y., S.H. and Z.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** Project supported by National Natural Science Foundation of China (51967005, and 12172096); Science and Technology Planning Project of Guangxi Province (GuiKe AD20297022, and Guike AD20159081); Guangxi Key Laboratory of Manufacturing System & Advanced Manufacturing Technology (Grant No. 20-065-40-003Z, and Grant No. 19-050-44-006Z).

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare that there is no conflict of interest regarding the publication of this article.

#### **References**


### *Article* **Improving Energy Harvesting from Bridge Vibration Excited by Moving Vehicles with a Bi-Stable Harvester**

**Zhiyong Zhou 1,2, Haiwei Zhang 1, Weiyang Qin 3, Pei Zhu 1,\* and Wenfeng Du <sup>1</sup>**


**Abstract:** To monitor the health status of the bridge, many sensors are needed to be mounted on it. Converting bridge vibration energy to electrical energy is considered as a potential solution to the problem of providing reliable electric power to these sensors. The objective of this work is to present an operable strategy for improving the electric energy output of a piezoelectric energy harvester installed on a bridge by introducing bi-stable characteristics. A bi-stable harvester is proposed. By adjusting the tip and fixed magnets, different types of potential energy can be generated, and then the harvester can exhibit the linear, mono-stable and bi-stable characteristics. In the bi-stable state, the harvester triggers snap-through motions easily and generates large outputs. The corresponding prototype was fabricated, and the experiment was carried out to validate the advantage of the bistable energy harvester. The experiment results show that the bi-stable energy harvester outperforms the classical linear harvester over the whole range of vehicle speed. As the vehicle speed exceeds a critical one, the snap-through motion will happen, which is beneficial to enhancing the electricity output. This conceptual design may provide guidance for promoting the performance of bridge energy harvesting.

**Keywords:** energy harvesting; bridge vibration; bi-stable energy harvester; linear energy harvester; snap-through motion

### **1. Introduction**

Bridges are becoming a more and more important part of the modern transportation system [1]. It is necessary to monitor the health of bridges at all times to ensure their safety. Although deploying wireless sensors in the detection region is an ideal scenario to monitor the state of the bridge, providing cheap and reliable power to the sensors is a well-known difficulty [2]. At present, the conventional electrochemical battery is still the primary power source of wireless sensors, but replacing the battery costs a lot, especially in dangerous locations [3]. Moreover, the heavy metals in wasted batteries pose a threat to the health and living environment of human beings. Energy harvesting from vibrations is considered a great potential energy solution to powering wireless sensors. In practice, there exists abundant vibration energy when vehicles pass through the bridge. This vibration energy results from passing vehicles and is not affected by the weather, geographical conditions and other factors, so it emerges as a promising solution for the power of bridge health monitoring systems [4]. Moreover, the piezoelectric energy harvester is usually simple in structure and easy to be manufactured. Thus, it can be produced in a large amount, so the vibration-based energy piezoelectric harvester has received more and more attention [5].

In recent decades, there are many studies focusing on converting the energy of bridge vibrations into electric energy. Ali et al. [6] investigated a linear energy harvester for the highway bridge under the excitation of a moving vehicle. The highest efficiency

**Citation:** Zhou, Z.; Zhang, H.; Qin, W.; Zhu, P.; Du, W. Improving Energy Harvesting from Bridge Vibration Excited by Moving Vehicles with a Bi-Stable Harvester. *Materials* **2022**, *15*, 2237. https://doi.org/10.3390/ ma15062237

Academic Editors: Marc Cretin, Sophie Tingry and Zhenghua Tang

Received: 16 February 2022 Accepted: 16 March 2022 Published: 17 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of energy harvesting was sensitively affected by the resonance frequency. Peigney and Siegert [7] investigated energy harvesting from a concrete highway bridge by a cantilever energy harvester. The results showed that the electric energy output could reach 0.03 mW. Zhang et al. [8] used a cantilever beam energy harvester to harvest energy from four bridges excited by truck models. The results showed that poorer road conditions on bridges were helpful for harvesting vibration energy. Karimi et al. [9] considered inertial effects on energy harvesting from a bridge with a heavy vehicle traveling on it. Zhang et al. [10] analyzed the output voltages of the piezoelectric energy harvester for nine different conditions. The vehicle speed had a significant effect on the energy output. Since the bridge response is the excitation source of the linear piezoelectric energy harvester (LPEH), the resonant frequency of LPEH must be tuned to the fundamental frequency of the bridge to obtain a high harvesting efficiency. Thus, some significant efforts have been devoted to studying the vibration response of a bridge and estimating the fundamental frequency and stiffness parameters. Aloisio et al. [11] conducted a modal analysis for a multi-span concrete bridge to estimate the compressive strength of concrete specimens with seven spans. Bonopera et al. [12] predicted the short-term dynamics of the uncracked prestressed concrete members, which was also valuable for designing new concrete bridges. However, it should be noted that the vibration energy of a bridge has a broadband frequency range since the speeds and weights of vehicles crossing the bridge change with time. Thus, the LPEH is inefficient in harvesting the broadband random vibration energy.

In recent years, some researchers tried utilizing nonlinearity to improve the operating bandwidth. The investigations on nonlinear energy harvesters with mono-stability [13], bi-stability [14], tri-stability [15] and quad-stability [16] have been carried out. For example, Stanton et al. [17] studied the hardening and softening response of the nonlinear monostable piezoelectric energy harvester. Sebald et al. [18] conducted a theoretical analysis of the nonlinear Duffing oscillator. The simulation results showed that the highest solution could obtain a huge gain in output power. Erturk and Inman [19] introduced a bi-stable piezoelectric energy harvester (BPEH) to improve the performance of power generation, which had a significant advantage over the LPEH. Arrieta et al. [20] introduced bi-stable characteristics in the cantilever composite to enhance the energy harvesting effectiveness under low-amplitude input levels. Masana et al. [21] designed an axially loaded BPEH, which could produce large output voltages. Lan et al. [22] added a small magnet in BPEH to improve the energy harvesting ability. Pellegrini [23] introduced a guideline to classify the BPEHs according to the feasibility of fabrication. Zhou et al. [24] proposed a broadband tri-stable piezoelectric energy harvester to further improve the energy harvesting ability. The results demonstrated that the tri-stable piezoelectric energy harvester could be superior to the bi-stable one, which was more suitable to the low input acceleration level. Litak et al. [25] investigated multiple solutions of a tri-stable energy harvester. It was found that the initial condition played an important role in the formation of the basins of attraction. Zhou et al. [26] designed a quad-stable energy harvester to improve the energy harvesting efficiency. The results proved the energy harvester designed in a quad-stable state made it much easier to attain snap-through motion, and its energy harvesting efficiency was significantly improved. Therefore, introducing nonlinear multi-stability is an effective way for efficient piezoelectric energy harvesting from the excitation in the shaking table test. Recently, Zhou et al. [27] exploited the mono-stable characteristics generated by the magnetic interaction to harvest bridge vibration energy, which brought about a noticeable improvement in the harvesting performance. However, to the authors' best knowledge, the concept of harvesting energy from bridge vibrations by bi-stable characteristics has not been reported so far. In bridge vibration energy harvesting, the bi-stable characteristics may help improve the harvesting performance significantly.

In this paper, the bi-stable characteristics generated by the magnetic interaction were exploited to enhance the performance of energy harvesting from a bridge under the excitation of moving vehicles. When the excitation triggers the BPEH to execute snap-through motion and oscillates in a high-energy orbit, it can generate a large output. In order to

study the dynamic characteristics of the BPEH, the potential energy, restoring force and stiffness are analyzed. Then, the output voltages and dynamic strains of LPEH and BPEH are measured in experiments and compared, as a vehicle model is designed to run on the bridge with different speeds. Finally, some important findings and conclusions are drawn.

#### **2. Energy Harvester for Bridge Vibration Induced by a Moving Vehicle**

The schematic sketch of the BPEH is shown in Figure 1. The BPEH is composed of a substrate partly covered by a piezoelectric patch or two, a tip magnet and two fixed magnets. A tip magnet is placed at the free end of the substrate, while two fixed magnets are mounted on the frame and positioned near the tip magnet. The tip magnet and fixed magnets have opposite polarities and face each other to produce a nonlinear magnetic attraction force. For clarity, we define two distances: one is the gap distance between two fixed magnets (*dg*); the other is the separation distance between the tip and the fixed magnets (*d*). Both of them can be adjusted to make the BPEH own two stable equilibrium positions. When the BPEH is excited by bridge vibrations, the piezoelectric patch will deflect and generate voltage. There is load resistance in the circuit. The BPEH can jump between the two stable equilibrium positions and execute snap-through motions, and then it will generate high voltage outputs. If the two fixed magnets are removed, the system degenerates to a classical LPEH. For the LPEH, as a linear system, it usually owns a narrow resonant bandwidth. Because the acceleration spectrum of bridge vibration usually exhibits a random characteristic and distributes over a wide range, the energy harvesting efficiency of LPEH under this excitation will be far from satisfactory.

**Figure 1.** Schematic of the BPEH.

Figure 2 shows the schematic of harvesting vibration energy from a bridge with a moving vehicle. As the vehicle runs on the bridge, the bridge is under the excitation of a moving load, so it will oscillate. In analysis, for simplicity, the vehicle can be modeled as a mass traveling on an elastic beam with length *Lb* and thickness *Tb* [9]. When the vehicle moves at a speed of *v* and passes the bridge, the bridge's oscillation will make the BPEH vibrate and produce snap-through motions.

**Figure 2.** Energy harvesting from a bridge traversed by a vehicle.

#### **3. Modeling Potential Energy, Restoring Force and Stiffness**

As we know, the LPEH has a linear restoring force and constant stiffness. It could exhibit a satisfactory work efficiency only when the frequency of the excitation source coincides with the resonance frequency. In fact, the vibration energy in the realistic environment usually distributes over a broadband frequency range, so the LPEH is inefficient in harvesting vibration energy. In contrast, the BPEH may overcome this defect.

To show the characteristics of the BPEH, its potential energy, restoring force and stiffness are presented. The total potential energy of the system includes three parts: the elastic potential energy of the substrate and piezoelectric patch, the gravity potential energy and the magnetic potential energy.

The substrate can be modeled as a Euler–Bernoulli beam, whose strain is proportional to the second spatial derivative of deflection. Thus, the elastic potential energy of the piezoelectric beam can be given by [28]

$$\mathrm{d}L\_{b} = \frac{1}{2} \mathrm{E}\_{\mathrm{s}} I\_{\mathrm{s}} \int\_{0}^{L\_{\mathrm{s}}} \left( \frac{\partial^{2} y(\mathbf{x}, t)}{\partial \mathbf{x}^{2}} \right)^{2} d\mathbf{x} + E\_{p} I\_{p} \int\_{0}^{L\_{p}} \left( \frac{\partial^{2} y(\mathbf{x}, t)}{\partial \mathbf{x}^{2}} \right)^{2} d\mathbf{x} \tag{1}$$

where *Es* and *Ep* are the Young's modulus of the substrate and that of the piezoelectric patch, respectively; *Is* and *Ip* are the inertial moments of the substrate and the piezoelectric patch, respectively; *Ls* and *Lp* are the lengths of the substrate and the piezoelectric patch, respectively; *y*(*x*, *t*) is the displacement of the piezoelectric beam.

The magnetic potential energy of the system comes from the magnetic attraction between the tip and fixed magnets, and the point dipole approximation is employed to model the magnets. Then, the magnetic potential energy considered in the BPEH can be given by [29]:

$$\mathcal{U}L\_{\rm H} = -\frac{\mu\_{0}a\_{1}a\_{2}}{2\pi} \left\{ \left( y(L\_{\rm s},t) - \frac{d\_{\rm \mathcal{S}}}{2} \right)^{2} + d^{2} \right\}^{-\frac{3}{2}} - \frac{\mu\_{0}a\_{1}a\_{3}}{2\pi} \left\{ \left( y(L\_{\rm s},t) + \frac{d\_{\rm \mathcal{S}}}{2} \right)^{2} + d^{2} \right\}^{-\frac{3}{2}} \tag{2}$$

where *μ*<sup>0</sup> is the magnetic permeability constant; *a*<sup>1</sup> and *a*<sup>2</sup> (*a*3) are the effective magnetic moments of the tip and fixed magnets, respectively; *d* is the distance between the tip magnet and the fixed magnets in the horizontal direction (as shown in Figure 1); *dg* is the vertical distance between the two fixed magnets.

The gravitation potential energies of the substrate, piezoelectric patch and tip magnet can be given by the following equation:

$$\mathcal{U}L\_{\mathcal{S}} = -m\_l \text{gg}(L\_s, t) - \rho\_s h\_s b\_s \text{g } \int\_0^{L\_s} y(\mathbf{x}, t) d\mathbf{x} - 2\rho\_p h\_p b\_p \text{g } \int\_0^{L\_p} y(\mathbf{x}, t) d\mathbf{x} \tag{3}$$

where *mt* is the mass of the tip magnet; *g* is the gravitational acceleration; *ρ<sup>s</sup>* and *ρ<sup>p</sup>* are the densities of the substrate and piezoelectric patch, respectively; *hs* and *hp* are the thicknesses of the substrate and piezoelectric patch, respectively; *bs* and *bp* are the widths of the substrate and piezoelectric patch, respectively.

Thus, the total potential energy of BPEH can be given as follows

$$
\mathcal{U}\_t = \mathcal{U}\_b + \mathcal{U}\_m + \mathcal{U}\_\mathcal{\mathcal{S}} \tag{4}
$$

The dynamic characteristics of the energy harvester can be inferred from its potential energy shape. Figure 3 shows the potential energy shapes of the harvester for different separation distances (*d*). The corresponding parameters are listed in Table 1. As is clear from Figure 3, when *d* is quite large (e.g., *d* = 90 mm), the harvester acts like a linear system since the magnetic force generated by the magnetic attractive effect is weak. Then, as *d* decreases (e.g., *d* = 17 mm), it turns into a mono-stable one due to the introduction of the magnetic force. It can be found that the curve of the potential energy function has a potential well and is asymmetric due to the gravitation effect. When *d* is less than a critical value, the attractive magnetic force will be greater than the resultant force of the elastic force and gravitational force. Thus, a potential barrier and two asymmetric potential energy wells are formed in the potential energy diagram. The asymmetry will become more obvious with the decrease in *d*. Subsequently, as *d* decreases further (e.g., *d* = 14 mm), the depths of two potential energy wells will increase, forming two deeper potential wells. For the bi-stable system, its snap-through motion may be excited to happen, which can make the piezoelectric beam have a large deflection and thus generate a relatively large voltage output.

**Figure 3.** Potential energy shapes of the harvester for various separation distances.


To further study the characteristics of the BPEH, Figure 4 shows the restoring force and stiffness of the BPEH with respect to different separation distances (*d* = 14, 15, 16, 17, 18 or 90 mm). As *d* decreases, the equilibrium position corresponding to the zero value of the restoring force gradually shifts toward the negative direction due to the gravitation effect, as shown in Figure 4a. When the attractive magnetic and elastic forces get to the same order of magnitude, the coupling effect becomes stronger. Now, the gravity effect plays an increasingly important role. As the separation distance decreases, it will influence the equilibrium position of tip deflection, which corresponds to the point of zero restoring force. When there appear three zero points in the restoring force curve, the system exhibits bi-stable behavior. In particular, by adjusting *d*, the restoring force and stiffness shown in Figure 4 can be reduced to a very small value and even become negative. The restoring force in the bi-stable state is extremely low and exhibits a negative value near zero tip deflection, which is expected to give rise to a large deflection and generate a high output voltage under weak vibration of the bridge.

**Figure 4.** (**a**) Restoring force and (**b**) stiffness of the BPEH for various separation distances.

#### **4. Experimental Setup**

To validate the energy harvesting efficiency of the BPEH, corresponding experimental tests were carried out. The experimental setup is shown in Figure 5. The vehicle is represented by a steel ball (1003 kg). The model bridge used for experimental tests is made from an acrylic sheet (1300 mm × 80 mm × 8 mm), as shown in Figure 5a. Two guardrails are attached to the bridge to prevent the ball from falling off the bridge. Because the first mode plays an important role in bridge vibration, the first modal displacement is a determinant element for the installation of the harvester. For collecting high electric energy, the BPEH is installed under the midpoint of the bridge to harvest the vibration energy. As for the excitation, a steel ball is released from different heights and rolls down the acceleration section (inclined track). Then, it passes the bridge, as shown in Figure 5b. The prototype of the BPEH is fabricated, as shown in Figure 6. The two stable equilibrium positions of the BPEH are shown in Figure 7. The substrate is made of stainless steel and has dimensions of 190 mm × 10 mm × 1 mm. The piezoelectric patch has the dimensions of 5 mm × 5 mm × 0.25 mm, which is connected to a high-precision resistance box. A strain sensor (120–5 AA) is bonded on the substrate to measure the dynamic strain. As for the magnets, three permanent magnets have the same dimensions of 10 mm × 10 mm × 5 mm. If the two fixed magnets are removed, the BPEH will degenerate to a classical LPEH (as shown in Figure 8). In all experimental tests, the LPEH and BPEH are put in the same experimental conditions to make a direct comparison. A data acquisition device (DH5922N, Donghua, JingJiang, China) is adopted to measure the output voltage and dynamic strain. The schematic diagram of the experimental procedure is illustrated in Figure 9.

**Figure 5.** The experiment setup: (**a**) test section and (**b**) acceleration section.

**Figure 6.** The prototype of the BPEH.

**Figure 7.** The two stable equilibrium positions of the BPEH: (**a**) stable equilibrium position 1 and (**b**) stable equilibrium position 2.

**Figure 8.** The prototype of the LPEH.

**Figure 9.** Schematic diagram of experimental procedure.

#### **5. Results and Discussions**

The dynamic behavior and electrical output are selected as the key indicators in evaluating the dynamic behavior and output performance. The steel ball is released from a titled track, which is named the acceleration section. In the acceleration section, for the ball, different positions could produce different initial speeds to pass the bridge. In the experiment, eleven positions are chosen to obtain different initial speeds. The running interval of the steel ball for measurement is set to be from 0.52 s to 1.18 s, which is consistent with the actual interval of a vehicle traveling through a bridge. Figure 10 shows the variance of strain and the power density (the load resistance is about 0.9 MΩ) of LPEH and BPEH for different moving speeds. It is clear from Figure 10 that, for BPEH, the variances of strain increase greatly at the speed of 1.68 m/s, resulting in a significant increase in the output power density. This phenomenon indicates that the BPEH is triggered to execute the snap-through motion and oscillates in a high-energy orbit. Thus, the BPEH can significantly improve the energy harvesting performance compared to the LPEH. For example, at the speed of 2.07 m/s, the BPEH's maximum power density can reach 430 W/m3, whereas the LPEH's maximum RMS voltage is only 261 W/m3.

**Figure 10.** (**a**) Variance of strain and (**b**) power density versus moving speed for LPEH and BPEH.

To reveal the vibration response of the LPEH and BPEH, Figure 11 illustrates the time histories of strain for six moving speeds (*v* = 1.10, 1.52, 1.82, 2.07, 2.30 and 2.50 m/s). In the time response, the tag "On bridge" indicates the instant when the vehicle enters the bridge, while the tag "Leaving bridge" indicates the instant when the vehicle gets off the bridge. At a very low moving speed *v* = 1.10 m/s or 1.52 m/s, the BPEH's oscillation cannot cross the potential barrier and thus is trapped in a potential well. Thus, both BPEH and LPEH have a small vibration amplitude now. When the moving speed increases to 1.82 m/s (Figure 11c), the BPEH produces a large-amplitude vibration, for it begins to realize jumps between two stable equilibrium positions. Now the maximum strain of the BPEH can reach 4.5 × <sup>10</sup>−4, whereas the maximum strain of the LPEH is only 3.4 × <sup>10</sup>−4. The amplitude of BPEH is larger than that of LPEH, and the vibration period of BPEH is longer as well. During the running period on the bridge, for the BPEH, the snap-through motion keeps happening. Finally, as the vehicle gets off the bridge, the response will fall into a potential well and oscillate in it. As the moving speed increases to 2.07 m/s, the vibration amplitude of BPEH will have a significant increase, whose maximum strain will increase up to 9.5 × <sup>10</sup>−<sup>4</sup> (Figure 11c). If the moving speed increases further to a high one (e.g., 2.30 m/s or 2.50 m/s), the frequent snap-through motion will happen in the BPEH, resulting in a much larger amplitude than the LPEH, as shown in Figure 11e,f. For example, at 2.30 m/s, the maximum strain of BPEH can reach 1.16 × <sup>10</sup>−3, nearly two times that of the LPEH (6.8 × <sup>10</sup><sup>−</sup>4). Thus, through comparing the time histories of the strain of two systems, it can be concluded that the BPEH is preferable to the LPEH in producing large deformations and giving large outputs for the same running vehicle excitations.

**Figure 11.** Strain responses of the LPEH and BPEH for the vehicle traveling on and leaving the bridge: (**a**) *v* = 1.10 m/s; (**b**) *v* = 1.52 m/s; (**c**) *v* = 1.82 m/s; (**d**) *v* = 2.07 m/s; (**e**) *v* = 2.30 m/s; and (**f**) *v* = 2.50 m/s.

Figure 12 shows the open-circuit output voltages of LPEH and BPEH for different moving speeds. It is well known that the PZT output is proportional to its deformation or strain. As is clear from Figure 12, the BPEH is superior to the LPEH in open-circuit voltage output at every speed. Specifically, at the low moving speed, e.g., *v* = 1.10 m/s, as Figure 12a shows, although the oscillation of BPEH is trapped in a potential well, the open-circuit voltage output of BPEH is still larger than that of LPEH due to the tailored potential energy shape. Then, when the moving speed increases to 1.82 m/s (Figure 12c), the BPEH's maximum open-circuit output voltage has a large increase and reaches up to 4.09 V, for the snap-through motion happens now, whereas the LPEH's maximum opencircuit voltage is only 2.44 V. Next, when the moving vehicle runs at a speed of 2.07 m/s, the open-circuit output voltage is further improved and is up to 6.72 V. However, if the moving speed increases continuously, the open-circuit output voltage of BPEH will not increase greatly, even will decrease. For example, as the moving speed reaches 2.30 m/s, the peak voltage of BPEH is 7.73 V. Then, as the moving speed increases a little to 2.50 m/s, the peak voltage is 7.56 V, a little smaller than that at 2.30 m/s. This is due to the protective effect of the fixed magnets. If the cantilever is excited to execute a very large vibration, the cantilever is likely to be damaged rapidly. Therefore, the attractive force between the tip and fixed magnets will take effect and constrain the excessively large vibration to protect the BPEH.

**Figure 12.** Open-circuit output voltages of the LPEH and BPEH for the vehicle traveling on and leaving bridge: (**a**) *v* = 1.10 m/s; (**b**) *v* = 1.52 m/s; (**c**) *v* = 1.82 m/s; (**d**) *v* = 2.07 m/s; (**e**) *v* = 2.30 m/s; and (**f**) *v* = 2.50 m/s.

In order to show the characteristics of the open-circuit voltage response of LPEH and BPEH, corresponding frequency spectra for different moving speeds are shown in Figure 13. It can be found that the peak of the LPEH is located near 14.3 Hz for any speed of the moving vehicle, implying that it is the system's resonance frequency. However, for the BPEH, the distribution of frequency spectra is different depending on the moving speed. At a relatively low moving speed (e.g., *v* = 1.10 m/s), the peak of the spectrum of BPEH output is located at 8.5 Hz, the BPEH's amplitude is much larger than that of the LPEH. When the moving speed increases to 1.82 m/s, the main peak shifts to 5.2 Hz, and the frequency energy of BPEH distributes over a wide range of 0–11.4 Hz. This is because the BPEH executes snap-through motions and exhibits a strong nonlinear characteristic. As the moving speed continues to increase, the spectrum of BPEH distributes more widely in a frequency range of 0~17.2 Hz (as shown in Figure 13d,f). Therefore, it can be concluded that the BPEH has a wider bandwidth and exhibits strong nonlinearity, which could help harvest more energy from bridge vibration.

**Figure 13.** Frequency spectra of the open-circuit output voltages for the vehicle traveling on and leaving bridge: (**a**) *v* = 1.10 m/s; (**b**) *v* = 1.52 m/s; (**c**) *v* = 1.82 m/s; (**d**) *v* = 2.07 m/s; (**e**) *v* = 2.30 m/s; and (**f**) *v* = 2.50 m/s.

#### **6. Conclusions**

In this work, a bi-stable PEH was developed by introducing a tip and two fixed magnets. The bi-stable characteristic of this BPEH is generated by the magnetic attractive interaction, which helps enhance the energy harvesting from bridge vibrations. For the BPEH, corresponding theoretical analysis and simulations were carried out, and the validation experiment was conducted. The results show that the BPEH can significantly improve the energy harvesting performance by realizing snap-through motions. Compared to the LPEH, the BPEH has a wide frequency spectrum, which varies with the moving speed. The voltage response of the BPEH outperforms that of the LPEH over the whole range of moving speeds. The BPEH could significantly increase the deflection under the same excitation, which is beneficial for harvesting more energy from bridge vibration induced by traveling vehicles. Compared to the classical designs, the BPEH is simple in structure and can be easily adjusted in size to obtain a satisfactory harvesting performance. This design may lead to the realization of an efficient and reliable energy harvester. However, there remain some aspects to be studied in future work. For example, the separation and gap distances can be optimized so as to make the BPEH realize snap-through motions more easily. In practice, multiple vehicles will run on the bridge simultaneously, the BPEH's response and output under this excitation should be studied.

**Author Contributions:** Conceptualization, P.Z.; methodology, W.Q.; software, H.Z.; formal analysis, Z.Z. and W.D.; data curation, H.Z.; writing—original draft preparation, P.Z.; writing—review and editing, Z.Z.; supervision, W.Q.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (grant # 52005155), the China Postdoctoral Science Foundation (grant # 2020M673470), the Key Scientific Research Project of Colleges and Universities in Henan Province (grant # 20A130001) and the Key Research Development and Promotion Project in Henan Province (grant # 202102310249, 212102310952).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

