*Article* **Cement Composites with Graphene Nanoplatelets and Recycled Milled Carbon Fibers Dispersed in Air Nanobubble Water**

**Anastasia I. Patrinou 1, Eirini Tziviloglou 2, Athanasios Varoutoglou 1, Evangelos P. Favvas <sup>3</sup> , Athanasios C. Mitropoulos 1, George Z. Kyzas <sup>1</sup> and Zoi S. Metaxa 1,\***


**Abstract:** The individual effect of nano- and micro-carbon-based fillers on the mechanical and the electrical properties of cement paste were experimentally examined in this study. The objective of the study was to separately examine the effects of size and morphology (platelets and fibers) of nano- and micro-reinforcement. Three different sizes of Graphene Nanoplatelets (GNPs), at contents of 0.05% and 0.20% and recycled milled carbon fibers (rCFs), at various dosages from 0.1–2.5% by weight of cement, were incorporated into the cementitious matrix. GNPs and rCFs were dispersed in water with air nanobubbles (NBs), an innovative method that, compared to common practice, does not require the use of chemicals or high ultrasonic energy. Compressive and bending tests were performed on GNPs- and rCFs-composites. The four-wire-method was used to evaluate the effect of the conductive fillers on the electrical resistivity of cement paste. The compressive and flexural strength of all the cementitious composites demonstrated a considerable increase compared to the reference specimens. Improvement of 269.5% and of 169% was observed at the compressive and flexural strength, respectively, at the GNPs–cement composites incorporating the largest lateral size GNPs at a concentration of 0.2% by weight of cement. Moreover, the rCFs–cement composites increased their compressive and flexural strength by 186% and 210%, respectively, compared to the reference specimens. The electrical resistivity of GNPs- and rCFs-composite specimens reduced up to 59% and 48%, respectively, compared to the reference specimens, which proves that the incorporation of GNPs and rCFs can create a conductive network within the cementitious matrix.

**Keywords:** graphene nanoplatelets; recycle carbon fibers; air nanobubbles; cement-based composites and nanocomposites; mechanical properties; electrical properties

#### **1. Introduction**

Recognizing cement as one of the dominant materials of the construction industry, researchers from around the world focus on optimizing its structural health. Cement-based materials have been extensively used due to their high availability and low construction and maintenance costs. However, their brittleness and low tensile strength may significantly decrease their functionality [1]. Cement-based materials are prone to cracking, which might affect the integrity of structures over time [2]. The presence of cracks can cause cement degradation by increasing its permeability, and might make cement vulnerable to corrosion, chloride penetration, etc., leading to accelerated deterioration [3]. Therefore, it is important to enhance the cracking resistance of cement-based materials. Reinforcement of cementitious composites in micro- and macro- scale using fibers and other fillers is a common method for the improvement of the intrinsic brittle behavior of cement-based

**Citation:** Patrinou, A.I.; Tziviloglou, E.; Varoutoglou, A.; Favvas, E.P.; Mitropoulos, A.C.; Kyzas, G.Z.; Metaxa, Z.S. Cement Composites with Graphene Nanoplatelets and Recycled Milled Carbon Fibers Dispersed in Air Nanobubble Water. *Nanomaterials* **2022**, *12*, 2786. https:// doi.org/10.3390/nano12162786

Academic Editors: Mohammad Malikan and Shahriar Dastjerdi

Received: 22 June 2022 Accepted: 10 August 2022 Published: 14 August 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/).

materials [4]. Nano- and micro-reinforcement can increase the tensile strength of the cement matrix by enhancing the load-transfer capacity and the crack bridging mechanism.

Recent advances in materials science and nanotechnology enable the use of various types of carbon-based nanomaterials, such as nanofibers, nanotubes, and nanoparticles, as reinforcing agents in the cement to control the creation of cracks at the nanoscale. Carbon nanomaterials demonstrate outstanding mechanical chemical and physical properties, such as low density, high specific surface area, hardness, and excellent chemical and thermal stability. Carbon nanomaterials incorporated in cement-based materials at relatively low proportions provide important improvement in the compressive, tensile, and flexural strength [5–7], as well as new functions such as stress- and strain-sensing, temperature monitoring, and energy harvesting [8–11].

Lately, graphene-based materials have captured the attention of scientists. Graphene has extraordinary thermal, electrical, and mechanical properties, and it has been described as the most powerful, tough, thin, and dense material in the world [12]. Both its electrical and thermal conductivity are considered superior to the materials used to date. Recently, a number of graphene composites have been developed utilizing the excellent properties of graphene [13]. Graphene Nanoplatelets (GNPs) are of particular interest for combining superior properties, such as two-dimensional planar structure, great aspect ratio, remarkable strength and durability, thermal, electrical, and chemical stability, with low cost and weight, and large-scale production potential. They are also electrical and thermal conductors [14]. Compared to pristine graphene, GNPs are a promising alternative as nanofillers in material science, due to their availability in the market and their affordable prices [15].

By incorporating GNPs in cement to develop innovative and advanced composites, improvements in the stability and longevity of structures are possible [16]. Research has shown that the use of GNPs accelerates the hydration reactions of cement, improves the porosity, and creates a denser microstructure [17,18]. It has also been reported that GNPs are able to increase the compressive and flexural strength of the matrix. Silva et al. [19] found that the incorporation of GNPs in mortar, in a small dosage 0.021% by weight, enhanced the compressive strength by 95.7% after 28 days. Tong et al. [20] used 0.1% by weight GNPs and recorded 19.9% increase in the compressive strength. Other studies refer that cement composites modified with GNPs 0.05% and 0.06% by weight, show flexural strength increase of 15–24% and 27.8%, respectively [21,22]. GNPs have also been used to create conductive cement-based composites, which can monitor their own strain by detecting alterations in the electrical resistivity values [23,24].

Natural and synthetic fibers are widely used in fiber-reinforced cement-based materials [25]. The material, the type, and the amount of the fibers used determine the properties of the composite materials [26]. The implementation of microfibers can greatly improve the mechanical strength of cement and restrain the crack growth into the matrix, however, they cannot prevent crack formation [27]. Due to their unique properties, carbon fibers stand out compared to other types of fibers. Carbon fibers (CFs) have many advantages, including excellent strength-to-weight ratio, high tensile strength and stiffness, low weight, chemical and thermal stability, and low thermal expansion [28,29]. They have high electrical conductivity and can be used to reduce the electrical resistance and provide self-detecting capabilities, as well as to create a shield of electromagnetic protection in building materials [30]. An increase of 85% in flexural strength and 22% in compressive strength has been reported, along with the introduction of a small amount of CFs (0.189% by volume) into concrete [31]. Cholker et al. [32] studied the behavior of "smart" carbon reinforced concrete with microfibers. The specimens were assessed under loading-unloading cycles to examine the stress and damage sensing ability. It was derived that the presence of CFs in 1.5% and 2% by weight of cement provides a sufficient link between stress and strain of concrete and electrical resistivity. In other research, it was found that cement-based composites reinforced with short CFs are capable of sensing their own damage, stress, and temperature [33].

CFs, despite their numerous advantages in modifying the performance of cement composites, have excessive cost, which is a deterrent to their use. Nowadays, the use of recycled CFs (rCFs) as reinforcing material has gained increasing attention owing to their durability and outstanding properties with decreased cost [34]. rCFs reveal superior mechanical performance, attributed to their rugged surface, which creates a stronger bond with the surrounding cement matrix [35]. Faneca et al. [36] used rCFs to fabricate low cost electrically conductive cement composites that can be employed not only in laboratory scale, but also in the civil engineering industry. Akbar and Liew [37] investigated the effect of elevated temperatures on the reinforcement mechanism of cement composites containing rCFs. The amount of fibers and elevated temperatures remarkably affect the mechanical properties and the mass loss of the rCFs–cement composites. Mastali and Dalvand incorporated rCFs 1.25% by volume in plain concrete, which resulted in an increase of 65.10% and of 66.93% on the compressive and the flexural strength of specimens, respectively. Moreover, the impacted resistance increased 6.48 times with rCFs concentration of 1.25% [38]. rCFs represent an alternative to obtain superior performance and low-cost cementitious composites, and they are commercially available by several companies. However, the scientific research on rCFs–cement composites has been limited. The effect of rCFs on the mechanical properties, as well as the interaction between the rCFs and the surrounding cement matrix, has not been investigated yet.

The homogeneous dispersion of carbon materials in cementitious composites is another critical issue that strongly influences the reinforcing efficiency and the ultimate performance of the cement composites. GNPs and rCFs, due to strong Van der Waals forces that are developed among them, have the tendency to agglomerate. Their hydrophobic nature makes their homogeneous dispersion in aqueous suspensions challenging [39]. However, homogeneous dispersion is usually succeeded either mechanically, by using high ultrasonic energy, or chemically, by using various water-reducing agents, such as superplasticizers as surfactants [40–42]. A new promising technology called 'nanobubble technology' is utilized by many important fields for applications, including flotation, nanomaterials dispersion, and crystal growth [43–45]. Nanobubble nucleation is favored at hydrophobic particle surfaces, enhancing the stability of micro-bubble suspensions, and thus, facilitating particle flotation [46]. The use of water enriched with air nanobubbles improves the morphological characteristics of the dispersed nanoparticle clusters without the presence of aggregates [45]. To our knowledge, the application of this technology on dispersing nano- and micro-carbon materials for use in cement pastes has not been previously studied.

Recently, extensive research has been conducted on the improvement of the mechanical and electrical properties of cementitious composites using carbon-based materials. Their size and morphology strongly affect their enhancing action. However, the study on the reinforcing mechanism of different fillers based on their dimensions or size scales is narrow. The influence of the different type and scale of the reinforcing materials on the mechanical and the electrical properties of cement-based composites has not been explicitly clarified yet.

In this research, two commercially available carbon-based materials, GNPs and milled rCFs, were incorporated in the cementitious matrix to serve as reinforcement in nano- and micro-scale. Composite mixtures with three different sized GNPs, at contents of 0.05% and 0.20% and rCFs at proportions from 0.1–2.5% by weight of cement, were fabricated. The prepared specimens were assessed for their compressive strength, flexural strength and electrical resistivity. Overall, recycled milled carbon fibers with improved carbon footprint and lower cost were investigated in the present study, as they differ significantly in size from short carbon fibers reported in the literature (rCFs length is at the μm scale while typical CFs have length at the mm scale) and the research of their effect on the mechanical properties of cement-based materials is limited. Another innovative aspect of this work is the usage of water with air nanobubbles (NBs) for the proper dispersion of GNPs and rCFs. It is noteworthy that compared to frequently used practices, this method does not require the use of chemicals or high ultrasonic energy, and it can be easily applied on the construction site.

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

#### *2.1. Materials and NBs Method*

A commercial Portland cement CEM II 32,5 N (TITAN Cement Company S.A Athens, Greece) was used for this study. Three types of GNPs were used, namely N008-100-N, N008-100-P10, N008-100-P40 (Angstron Materials Inc. Dayton, OH, USA), exhibiting the same thickness but different lateral size (diameter). The physical properties of the GNPs are shown in Table 1. Moreover, milled rCFs, produced by Haufler Composites GmbH & Co. KG, Blaubeuren, Germany, were utilized in this study. According to the manufacturer CF−milled100 is a mixture of all origins carbon and graphite ex-PAN fibers, obtained from virgin carbon fibers, chosen and milled for compounding. Their mechanical and physical properties are shown in Table 2.


**Table 1.** Physical properties of GNPs.

**Table 2.** Properties of the carbon fibers (given by the manufacturer).


All specimens were prepared with tap water (low salinity, S < 350 mg/L) enhanced with air nanobubbles (NBs). Air nanobubbles were created using a generator which utilized the counterflow hydrodynamic cavitation. Details about the used NBs generator can be found in our previous work [47]. The mean size and the concentration of the produced NBs are shown in Figure 1. The optimal processing time, regarding the size of the produced NBs, was 30 min [47]. The NBs methodology was chosen for the adequate dispersion of the GNPs and rCFs. Due to the efficacy of the method, no additives, such as superplasticizers, water reducing agents, or surfactants, were employed in this study. Another advantage of this method is that zero energy consumption is required, as the water feed pressure must be at 2.5 bar, equal to the water supply network pressure, and of course it can be easily applied in the cement industry. The difference in the dispersion quality between aqueous solutions with and without NBs is easily noticeable (Figure 2).

**Figure 1.** Nanobubbles concentration versus size.

**Figure 2.** GNP suspensions in (**a**) plain tap water and (**b**) tap water with nanobubbles.

#### *2.2. Preparation of Composites Specimens*

Two concentrations of each GNPs-type were used, namely 0.05% and 0.2% by weight (wt.) of cement. Furthermore, six mixtures of milled rCFs were prepared at different concentrations, namely 0.0%, 0.1%, 0.5%, 1.0%, 1.5%, 2.0%, 2.5% by weight of cement. The reinforcing materials were added into water with NBs, and the aqueous suspensions were then poured in the mixer (Figure 3). The pastes were prepared in a stainless-steel standard mixer, according to ASTM C305, at low speed (140 rpm) for 30 s and then at high speed (285 rpm) for 60 s. The water-to-cement ratio (w/c) was constant at 0.50 for all the mixtures. The proportions of the materials used for reference, GNPs- and rCFs-composites are presented in Tables 3 and 4, respectively.

**Figure 3.** GNP–cement composite preparation.

**Table 3.** Mix proportions for reference and GNP-composites.


**Table 4.** Mix proportions for reference and rCFs-composites.


The prepared GNPs- and rCFs–cement composites were cast in 30 mm × 60 mm cylindrical steel molds and 20 mm × 20 mm × 80 mm prismatic plastic molds, to investigate

the compressive and flexural strength, respectively. To measure the electrical resistance of the composite specimens, right after molding, four stainless steel electrodes were embedded into the specimens covering their entire cross-section (Figure 4a,b). The stainless steel mesh used had openings of 1.74 mm × 1.74 mm, wire thickness of 0.8 mm, and dimensions of 20 mm × 50 mm (Figure 4c). After 48 h, the specimens were de-molded and placed in tap water tanks with calcium hydroxide, where they were kept for 28 days.

(**a**) (**b**) (**c**)

**Figure 4.** Specimens for the electrical resistance tests (**a**) inside the mold write after casting, (**b**) after demolding, and (**c**) stainless steel mesh embedded inside the samples.

#### *2.3. Experimental Measurements*

The effect of GNPs and milled rCFs on the compressive strength pastes was evaluated after 28 days of curing. Three specimens of each testing series were assessed. The compressive strength was measured using the Instron 8801 testing machine (Instron, Norwood, MA, USA) with maximum load capacity of 100 kN at a constant displacement rate of 0.2 mm/min. Loading was applied monotonically until the failure of the specimen (Figure 5).

**Figure 5.** Experimental apparatus for compressive stress test (close up shows a fractured sample).

The effect of GNPs and milled rCFs on the flexural strength of the cement paste was evaluated after 28 days of curing. Flexural strength measurements were carried in a 4-pointbending configuration (Figure 6) using an MTS Insight testing machine (MTS Systems Corporation, Eden Prairie, MN, USA) with maximum load capacity of 10 kN applying a constant displacement rate of 0.002 mm/min (Figure 4b). The average flexural strength of three specimens was recorded and presented as representative flexural strength of the cement composites.

**Figure 6.** 4-point bending test experimental set up.

The electrical resistance was measured using the four-wire-method (Figure 7). Three specimens of each testing series were tested. After 28 days of curing, the specimens were placed in a laboratory drying oven where they remained for 72 h, at a temperature of 80 ◦C, to evaporate the water trapped in the pores. A 34450A Keysight Laboratory Digital Bench Multimeter, Keysight Technologies, Santa Rosa, CA, USA )was used for the electrical measurements. The internal electrodes were used to measure the voltage, while the external electrodes were used to supply the current as shown in Figure 7.

**Figure 7.** Four-wire method set-up.

During the measurements, a rubber material was placed under the specimens to insulate the specimen. The electrical resistance was recorded every 2 s over a period of 30 min to avoid potential deviations caused by the effects of electric polarization. Data from the last five minutes of the measurements were used to calculate the average resistance values. The resistivity, *ρ*, of the nanocomposites was calculated using Ohm's law as

$$
\rho = \mathcal{R} \frac{S}{L}
$$

where *R* is the electrical resistance, *S* is the cross-section of the sample, and *L* is the distance between internal electrodes.

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

#### *3.1. GNPs Strenthening Mechanism*

The average measured values of compressive and flexural strength of cement composites with different sizes and contents of GNPs are illustrated in Figure 8. The abbreviations of N, P-10, and P-40 shown in the graph refer to the N008-100-N, N008-100-P10, and N008- 100-P40 GNPs types followed by their concentration 0.05 and 0.2% by weight of cement, i.e., P10\_0.05 is the N008-100-P10 GNPs-type at an amount of 0.05%.

**Figure 8.** Average (**a**) compressive strength and (**b**) flexural strength of cement composites with GNPs at 28 d.

The incorporation of GNPs into the cement paste is beneficial for its compressive and flexural strength. The compressive strength of the plain cement paste was found to be 9.5 MPa and raised to 30.8 MPa and 35 MPa for the P10\_005 and P40\_02 specimens, corresponding to an increase of 224.8% and of 269.5%, respectively, which is more than reported in the literature [19,21,48–50]. The GNPs-composites also demonstrated an important increase in the flexural strength compared to the plain cement specimens. The highest flexural strength values varied between 6.14 MPa and 6.25 MPa, corresponding to an increase of 164.6% and 169.4%, for the P10\_005 and P40\_02 specimens, respectively, compared to the reference specimens (2.32 MPa).

This high improvement in the mechanical properties of GNPs–cement composites can be attributed to the interaction between the GNPs and the paste. GNPs possibly have developed good affinity/bond with the cement matrix, which can enhance the loadtransfer capacity by absorbing energy more efficiently and offer a better resistance to crack spreading [17,22,51,52]. Upon load application, GNPs effectively transfer the stress by enhancing the crack bridging mechanism and restraining the crack formation [50,53]. Cracks are blocked and diverted or branched (Figure 9) [19].

**Figure 9.** GNPs action under 4-point bending load.

According to the results, there were two reinforcing mechanisms that were developed and associated with the GNPs size. When GNPs with a low average diameter of 5 and 10 μm are incorporated, a greater increase in the compressive strength of cement composites is observed at low concentration (0.05 wt.%). This agrees with other studies claiming that GNPs demonstrate better performance at lower concentrations, since a more homogeneous GNPs network can be developed at lower concentrations [49,53]. On the contrary, GNPs type N008-100-P40 demonstrated the maximum compressive stress improvement at the higher concentration (0.2 wt.%). The P-40 type demonstrates up to 10 times larger diameter compared to the N and P-10 GNPs-types. This finding is supported by the fact that larger sized nanomaterials are more easy to disperse [54,55]. Therefore, for the P-40 GNPs, the development of a homogeneous network at a larger concentration is feasible.

#### *3.2. rCFs Strenthening Mechanism*

The average compressive and flexural strength values of cement composites reinforced with milled rCFs and their deviation from the reference specimens are shown in Figure 10.

**Figure 10.** Average (**a**) compressive strength and (**b**) flexural strength of cement composites with milled rCFs at 28 d.

The mechanical strength of cement paste filled with different dosages of rCFs was significantly enhanced. This improvement seems to depend on the rCFs content. The compressive strength of all the specimens was enhanced gradually, with the increase in fiber content ranging from 0.1% to 2.5 wt.%. When the fiber content was 1% and 2.5% by weight of cement, an increase of 143.7% and of 185.9% was recorded, respectively. Similarly, Akbar et al. observed a 47% increase in compressive strength of cement composites by the addition of 1% by volume of milled rCFs [35]. Moreover, Mastali and Dalvand reported 65.10% improvement on the compressive strength of concrete specimens filled with 1.25% by volume rCFs [38].

The rCFs-composites showed considerable improvement in flexural strength compared to the plain cement specimens. A maximum flexural stress of 7.2 MPa was recorded for CF\_2.0 mixture (210.3% higher than the reference mixture). This improvement is related to the crack bridging action of CFs. According to the results, the flexural performance of the specimens improves as the fiber content increases. This trend is in line with the findings of Yoo et al. [56], Donnini [30], and Han et al. [57]. As the fiber concentration in the matrix increases, the crack resistance mechanism and the fracture energy absorption improve. However, after a certain concentration, the CFs agglomerate and the mechanical strength of the composite reduces [57].

#### *3.3. Electrical Resistivity of GNPs–Cement Composites*

Typical curves of the electrical resistivity over testing time for the plain cement paste for the three different GNPs-types in both concentrations are shown in Figure 11. High electrical resistivity values (1.29 MOhm·cm) were recorded for the reference mixture. At the beginning of the test, increased values of approximately 1,42 MOhm·cm were observed. This can be explained by the electrical polarization phenomenon. Electric polarization occurs when a dielectric material is exposed to a DC electric field. Next, positive and negative charges are pushed in opposite directions, resulting in separation of positive and negative charges throughout the mass of the material. Consequently, the resistance of the material, especially at the beginning of measurement, is constantly changing. The resistance stabilizes after approximately 3 min; however, marginal scatter of about 5% may be observed in the measurements.

**Figure 11.** (**a**) Typical electrical resistivity over testing time curves of the GNPs–cement composites; (**b**) Average electrical resistivity results of the GNPs–cement composites.

It is evident that GNPs-composites, regardless of the GNPs' type and amount, exhibit lower resistivity values compared to the plain cement paste [24,49]. The values of the resistivity were lower when the concentration of GNPs was 0.05 wt.%, which may be related to the GNPs dispersion state, i.e., the more uniform the dispersion of GNPs, the lower the electrical resistance in the cement matrix. The lowest *ρ*, 59% reduction compared to the plain cement paste, was measured for the specimens with GNPs-type N008-100-P40. Although the concentration of GNPs remained the same in all specimens, the results indicate that N008-100-P40 created the most efficient conductive network within the cementitious matrix. This is probably associated with its diameter, which was larger (44μm) than the other two types of GNPs. Therefore, the distance between the GNPs is small enough to form conductive paths in the matrix, favoring the passage of the current [17]. A significant decrease in the resistivity of nanocomposites modified with low contents of GNPs has been also reported [51,53].

When the concentration of GNPs increased to 0.2 wt.%, all the composite specimens demonstrated similar reduction in the resistivity close to 35%. This finding can be attributed to the inadequate dispersion of the nanomaterials. According to the literature [9,58,59], the resistivity of a material depends on the quantity and sufficient dispersion of conductive additives (Figure 12). Inadequate dispersion leads to GNPs-agglomeration, which increases the electrical resistance of the composite.

**Figure 12.** Electrical conductivity of GNPs–cement composites.

#### *3.4. Electrical Resistivity of rCFs–Cement Composites*

The typical electrical resistivity curves over testing time for the reference sample, and for rCFs–cement composites, are shown in Figure 13. The rCFs-composites demonstrated reduced resistivity compared to the reference specimens. The rate of resistance reduction is affected by the concentration of rCFs within the matrix. The highest reduction rate of 48% was measured for the specimens with fiber concentration of 1 wt.%. When the fiber content is low, the current flow is difficult, since there might be a large space among the fibers. However, by increasing the concentration of the rCFs, the fibers approach each other, and a stable conductive network is established in the matrix [57]. The value of the resistivity decreases until the percentage of rCFs reaches 1 wt.%, indicating that the percolation threshold, i.e., when the fibers form a continuous electrical network, has been reached. More rCFs cause an increase in the value of the electrical resistance.

**Figure 13.** (**a**) Typical electrical resistivity over testing time curves of the investigated rCFs–cement composites; (**b**) Average electrical resistivity results of the GNPs–cement composites.

#### **4. Conclusions**

The mechanical and electrical behavior of cement paste reinforced with GNPs and milled rCFs were experimentally examined in this research. The addition of GNPs and rCFs into cement paste had a positive effect on the mechanical and electrical properties.

Among the three types of GNPs used in this study, the most effective proved to be the N008-100-P40 with the larger lateral size. Specimens with 0.2 wt.%. GNPs concentration increased both the compressive and the flexural strength of the cement composites by 269% and 169%, respectively. Furthermore, N008-100-P40 at a content of 0.05 wt.%. created the most effective conductive network within the cementitious matrix, attributed to their larger diameter.

Increasing the content of milled rCFs enhanced the mechanical properties of cement paste. The compressive and flexural strength of the rCFs-composites increased by 186% and 210%, with the incorporation of fibers at contents of 2.5 wt.% and 2 wt.%, respectively. The rate of reduction in electrical resistivity was affected by the concentration of rCFs in the composite. Percolation threshold was reached at a fiber content close to 1 wt.% of cement, as an increase of the fiber dosage was not able to further reduce the electrical resistivity of the cement paste.

**Author Contributions:** Conceptualization, Z.S.M.; methodology, A.I.P., E.T. and Z.S.M.; investigation, A.I.P., E.T., A.V., E.P.F. and G.Z.K.; resources, G.Z.K. and A.C.M.; writing—original draft preparation, A.I.P., E.T.; writing—review and editing, E.P.F. and Z.S.M.; supervision, Z.S.M.; project administration, A.C.M.; funding acquisition, Z.S.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was co-financed by the European Regional Development Fund of the European Union and Greek national funds through the Operational Program Competitiveness, Entrepreneurship, and Innovation. Projects: (i) "Nanoreinforced concrete for pavement deicing" (acronym: NEA ODOS; project code: T1EΔK-02692; MIS: 5030194) under the call Research– Create–Innovate.

**Acknowledgments:** E.P.F., A.M., G.Z.K. and Z.S.M. acknowledges funding from Research–Create– Innovate call under the framework of project NEA ODOS.

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

#### **References**


## *Article* **Study on Nanoporous Graphene-Based Hybrid Architecture for Surface Bonding**

**Xiaohui Song 1,\* , Mingxiang Chen 2, Jingshuang Zhang 3, Rui Zhang <sup>3</sup> and Wei Zhang <sup>1</sup>**

	- Wuhan 430074, China; chimish@163.com

**Abstract:** Graphene-copper nanolayered composites have received research interest as promising packaging materials in developing next-generation electronic and optoelectronic devices. The weak van der Waal (vdW) contact between graphene and metal matrix significantly reduces the mechanical performance of such composites. The current study describes a new Cu-nanoporous graphene-Cu based bonding method with a low bonding temperature and good dependability. The deposition of copper atoms onto nanoporous graphene can help to generate nanoislands on the graphene surface, facilitating atomic diffusion bonding to bulk copper bonding surfaces at low temperatures, according to our extensive molecular dynamics (MD) simulations on the bonding process and pull-out verification using the canonical ensemble (NVT). Furthermore, the interfacial mechanical characteristics of graphene/Cu nanocomposites can be greatly improved by the resistance of nanostructure in nanoporous graphene. These findings are useful in designing advanced metallic surface bonding processes and graphene-based composites with tenable performance.

**Keywords:** surface bonding; nanoporous graphene

#### **1. Introduction**

Surface bonding is an essential step in three-dimensional electronic packaging, as it provides mechanical support, heat transfer, and electrical integration [1]. A low-temperature bonding technology and reliable connection interface greatly influence electronic systems' performance and service life, especially when packaging density dramatically increases with device scaling [2,3]. Traditional surface bonding techniques in electronic assembly rely on high-temperature processes such as reflow soldering [4,5] and thermo-compression bonding [6], which can lead to undesirable thermal damage, toxic solder materials pollution, and a thermal mismatch at the bonding interface. Recently various nanometal materials such as metal nanowires, nanoparticles, and nanocones-based surface bonding are being studied to lower the bonding temperature and pressure [7–10]. However, these bonding technologies introducing low deformation resistance joints have been entangled in thermomechanical stresses and aging degradation issues, limiting their reliability. Therefore, it is necessary to develop high stretch and shear deformation resistance of interconnection materials architecture and a low-temperature bonding process for a complex interface of three-dimensional packaging structures.

Since its isolation as a two-dimensional (2D) system, graphene-copper nanolayered composites have been widely acclaimed and proved to be the effective alternatives to pure metal materials owing to their outstanding mechanical and heat transfer properties [11–15], showing great promise for next-generation electronic and optoelectronic device packaging. However, the reinforcing effect of graphene is limited by vdW non-bonded interfacial interaction between graphene and Cu, resulting in a low load transfer rate from Cu matrix

**Citation:** Song, X.; Chen, M.; Zhang, J.; Zhang, R.; Zhang, W. Study on Nanoporous Graphene-Based Hybrid Architecture for Surface Bonding. *Nanomaterials* **2022**, *12*, 2483. https:// doi.org/10.3390/nano12142483

Academic Editors: Mohammad Malikan and Shahriar Dastjerdi

Received: 28 June 2022 Accepted: 10 July 2022 Published: 20 July 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/).

to graphene during mechanical deformation [16]. To pave the way toward practical applications, a variety of approaches have been developed, such as chemical functionalization and metal coating [17,18]. The controllable and precise solution that uses graphene-metal composites as bonding layers remains a major challenge, despite considerable progress.

Here, we report a simple new Cu-nanoporous graphene-Cu based bonding technology that is of low bonding temperature and high reliability. We transfer nanoporous graphene to the copper surface and form copper nanoparticles on nanopores by depositing copper atoms before thermocompression bonding. Since the copper nanoparticles deposited in the nanopores form a metal bond connection with the copper bonding surface, they effectively enhance the interfacial connection and, consequently, lead to higher interfacial shear strength between graphene and Cu surface. At the same time, the metal nanoparticles contribute to solid-state diffusion and intermixing of surface atoms between the interfaces at low temperatures. The present study aims to quantitatively investigate the abovementioned process. Comprehensive molecular dynamic simulations are carried out, and the results reveal that the copper-nanoporous graphene-copper sandwich structure can be bonded at low temperature by the thermocompression method, and the interfacial shear strength of the interface can be significantly increased by nanoporous-nanometal particles hybrid structure.

#### **2. Computational Methods**

For molecular dynamics simulation of Cu-nanoporous graphene-Cu bonding process and the reinforcing effect evaluation, the setup of the simulations performed in the present work is shown in Figure 1. Initially, the composite model is constructed by covering singlelayer graphene with nine nanopores onto the single crystal copper surface that contains 34,924 atoms with the size of 161.78 Å × 99.65 Å × 23.35 Å, as shown in Figure 1a. The boundary condition for the simulation box is periodic in the *x* and *y* directions while the nonperiodic boundary condition is applied along the *z*-direction. The equilibrium structures are achieved using the canonical ensemble (NVT) [19] of the constant volume and temperature for 20 ps with a time step of 0.5 fs.

**Figure 1.** Atomic configurations. (**a**) Cu-nanoporous graphene composite model; (**b**) Cu atoms deposition onto the nanoporous graphene surface; (**c**) Cu-nanoporous graphene-Cu thermocompression bonding; (**d**) Pull-out simulation model.

Next, the bottommost two layers are fixed. The other atoms above the fixed layers are defined as temperature control layers that control the system temperature using the canonical ensemble (NVT), and the time step is set to 0.5 fs. In the condition of copper sputtering deposition simulation, the deposition rate of incident Cu atoms is 0.25 atoms/ps. The coordinates of incident Cu atoms were distributed randomly within the defined insertion volume, which is 50 nm above the graphene surface, as shown in Figure 1b. In each simulation, after all Cu atoms are released, a relaxation process of 2500 ps is conducted to enable the deposited atoms to reach full thermal equilibrium.

Then the bonding model is constructed by introducing a single crystal copper cell of the same size as the previous composite model as shown in Figure 1c. The distance between the two surfaces is set to nm. The uppermost two atomic layers are set as a fixed layer to produce the pressure of the system, and five atomic layers connected to the fixed layers are set up as temperature-controlled layers using the canonical ensemble (NVT). The other free copper atoms and nanoporous graphene at the bonding interface use the micro-canonical ensemble (NVE) of the constant volume and energy with a time step of 0.5 fs. During the first stage of simulation, no pressure is loaded on the fixed layer in 2 ns to reduce internal stress and make the model be in a steady state at a certain temperature. The pressure is loaded on the upper fix layer for 4 ns. At last, the pressure is unloaded and the model holds for 5 ns.

Finally, the pull-out simulation is used to investigate the interfacial shear strength of Cu nanocomposite. During the pull-out simulations, the upmost two layers and the bottommost two layers of Cu atoms are fully fixed. The pull-out simulations using the canonical ensemble (NVT) are performed by applying a constant velocity to the graphene along the *x*-direction until it is completely pulled out from the bonding interface as shown in Figure 1d. To ensure that the model is fully relaxed in each step and eliminate the effect of velocity on the pull-out force, a relatively low velocity of 0.8 × <sup>10</sup>−<sup>5</sup> Å/fs is adopted to 245 atoms treated as a rigid body on the rightmost side of the graphene. A temperature of 300 K is employed for all simulations, and a time step of 0.5 fs is adopted throughout the whole simulation process. The boundary condition for the simulation box is periodic in the *y* directions while the nonperiodic boundary condition is applied along the *x* and *z* directions.

The embedded atomic method (EAM) [20] models successfully used in modeling various bulk metals are used to describe the Cu-Cu interactions. For the vdW interactions between the graphene surface and Cu atoms, the Lennard–Jones (LJ) pair potential involving nonbonded long-range interactions is used with parameters σ = 3.0825 Å, ε = 0.02578 eV [21], and a cutoff radius = 4 σ. Our simulations are based on the massively parallel LAMMPS code [22]. The visualization is based on OVITO [23].

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

First, the Cu atoms deposition onto the nanoporous graphene surface is investigated. As shown in Figure 2a, in the initial state, the Cu atoms prefer to deposit in the nanoporous graphene and form metallic connections with the attached copper surface until the porous are fully filled. The reason for this is that Cu-Cu interactions are stronger than Cu-graphene. More copper atoms are then deposited onto the copper structure in the nanopores and gradually grow to form nanoislands. The growth pattern observed in this simulation approximates the insular growth shown in Figure 2b. When the nanoislands reach a certain size, they join with neighboring nanoislands to form larger ones as shown in Figure 2c,d.

Then thermocompression bonding simulations are performed to exploit the bonding mechanisms of nanoporous graphene-based hybrid architecture. A constant temperature of 300 K and pressure of 0.5 MPa are applied to the model. In the initial state, two contact interfaces form between nanoislands and the bulk copper surface as shown in Figure 3a. Compression deformation of the nanoislands increases until the bulk copper surface and the graphene are joined together as shown in Figure 3b,c. In this process, the nanoislands are gradually dispersed on the graphene surface and fill the interface

while atomic diffusion occurs with the bulk copper surface. During this state, nanoislands undergo a typical process from elastic deformation to plastic deformation, and the crystal structure of nanoislands has been mostly destroyed. As a result, more and more active copper atoms on the nanoislands surface make contact with the bonding surface to form metallic bonding connections with time increasing.

**Figure 2.** Atomic configurations of Cu atoms deposition onto the graphene at various stages. (**a**) Filling of nanopores at the initial stage; (**b**) Insular growth to nanoislands; (**c**) Further growing of nanoislands; (**d**) Joining with neighboring nanoislands.

**Figure 3.** Cross-section configurations of Cu-nanoporous graphene-Cu during the bonding process at various stages. (**a**) Contacting of bonding surface; (**b**) Compression deformation of the nanoislands; (**c**) Final bonding structure.

Pull-out simulations of the nanoporous graphene from the metal matrix are used to investigate the reinforcing mechanisms of Cu-nanoporous graphene-Cu nanocomposites. Along the *x*-axis, a constant velocity is delivered to the graphene. Figure 4 shows the pull-out force for graphene/Cu composites with and without nanoporous surfaces in terms of sliding distance. The atomic configurations during the nanoporous pull-out process are depicted in the insets of Figure 4. The pull-out force without nanoporous grows fast until the sliding distance reaches roughly 8 Å as illustrated in Figure 4. After this, it swings at around 13 nN for the next 140 Å. The pull-out force then steadily declines until graphene is totally pulled out in about the last 10 Å of sliding distance. In this condition, the van der Waals interaction between the graphene surface and the copper atoms at the interface dominates the pull-out force. The pull-out force with nanoporous, on the other hand, rises quickly to about 60 nN at the start of the pull-out process, nearly six times higher than without nanoporous. Copper atoms in nanopores and their van der Waals contact with the graphene surface are credited with resisting the pull-out motion of nanoporous graphene. In particular, the graphene pull-out process requires breaking through the resistance of the nanostructures in the pores, which greatly increases the pull-out force. Therefore, much larger pull-out force is generated in graphene. As a result, a larger pull-out force is needed for the nanoporous graphene, giving rise to better interfacial properties than that of pure graphene.

**Figure 4.** The pull-out force changes as a function of the sliding distance for the graphene/Cu composite with and without nanoporous.

When the sliding distance is larger than 4 Å, the force decreases to about 18 nN when the first column of nanoporous is fully pulled out. Figure 5 presents the displacement of copper atoms at the interface of graphene/Cu nanocomposites for this process. As the sliding distance increases, the atoms in the nanoporous are largely deformed until they move along with the graphene. The fluctuation of the pull-out force corresponds to the deformation-movement process of copper structures in nanoporous. Consequently, van der Waals interactions play an increasingly important role, and the copper atoms in nanopores are becoming less resistant to graphene. Next, as the sliding distance increases, there are two processes of increasing and decreasing pull-out force until it becomes 0. It is noted that these two processes occur after the first and second columns of nanoporous are pulled out, respectively, as shown in Figure 3. Since the graphene is regarded as a flexible body during the pull-out process, during the constant velocity pull-out process, the outer porous are subjected to more deformation, which makes the unpulled part move relatively slowly, so the atoms in the nanoporous move more slowly, and the deformation effect of the structure makes it have more resistance. Therefore, the pull-out force appears to increase. ˄ ˅

**Figure 5.** The displacement of copper atoms at the interface of graphene/Cu nanocomposites with various sliding distances.

As mentioned above, the deformation resistance of the nanostructures obtained by depositing copper atoms in nanoporous is the key to enhancing the interfacial mechanical properties of composites. For this purpose, the effect of the number of deposited copper atoms on the pull-out force is investigated in the present work. Figure 6 shows the maximum pull-out force change as a function of the number of deposited copper atoms. A larger maximum pull-out force can be observed as-deposited atoms onto the graphene increase. The reason for this is that more deposited atoms help fill the nanopores and form nanoislands to form diffusion bonds with the bulk copper bonding surface. The nanostructures in nanoporous can increase the deformation resistance effect, consequently leading to improved interfacial mechanical properties. However, when the deposited atoms exceed 1200, the maximum pull-out force decreases instead, mainly because the oversized nanoisland structure is not sufficiently deformed in the bonding, which makes the bonding surface not fully in contact with the graphene surface and reduces the van der Waals force interaction between copper and graphene surface.

**Figure 6.** The maximum pull-out force changes as a function of the number of deposited copper atoms for the graphene/Cu composite.

The influence of nanoporous diameters on the reinforcing effect of graphene/Cu nanocomposites is also further studied. We have considered three different nanoporous diameters: small (3.45 Å), medium (8.07 Å), and large (11.93 Å). The number of deposited atoms is 500. As shown in Figure 7, the larger the nanoporous, the greater the maximum pull-out force. The reason is that larger nanoporous helps to form bigger and stronger nanostructures in the nanoporous, which are not deformable and require greater pullout forces to destroy the nanostructures. However, larger holes make the effective area of graphene reduced, which may affect the ability of graphene to perform thermal and electrical transport functions.

**Figure 7.** The maximum pull-out force changes vs. nanoporous diameters.

#### **4. Conclusions**

In this study, we report a novel Cu-nanoporous graphene-Cu based bonding technology that is of low bonding temperature and high reliability. The process of depositing copper atoms and thermocompression bonding has been numerically investigated by employing MD simulations. Numerical results show that deposition of copper atoms onto nanoporous graphene can help to generate nanoislands on the graphene surface, facilitating atomic diffusion bonding to bulk copper bonding surfaces at low temperatures. Moreover, the resistance of nanostructure in the nanoporous can dramatically improve the interfacial mechanical properties of graphene/Cu nanocomposites. It is worth mentioning that an increase in the number of deposited atoms enhances the maximum pull-out force, but too many atoms deposited will reduce the interfacial strength. While larger nanoporous helps to form bigger and stronger nanostructures in the nanoporous, which are not deformable and require greater pull-out forces to destroy the nanostructures. Decidedly, the research findings of this study provide a feasible and facile route for low-temperature metal surface bonding with a high-performance metal nanocomposites interface layer reinforced by nanoporous graphene.

**Author Contributions:** Conceptualization, X.S. and M.C.; methodology, X.S.; software, J.Z.; validation, R.Z., W.Z. and X.S.; formal analysis, J.Z.; investigation, X.S.; resources, R.Z.; data curation, W.Z.; writing—original draft preparation, X.S.; writing—review and editing, M.C.; visualization, W.Z.; supervision, X.S.; project administration, X.S.; funding acquisition, X.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by [Zhongyuan Science and Technology Innovation Leadership Program of China] grant number [214200510014]. And The APC was funded by [Zhongyuan Science and Technology Innovation Leadership Program of China].

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

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

#### **References**


## *Review* **Multiscale Mechanical Performance of Wood: From Nano- to Macro-Scale across Structure Hierarchy and Size Effects**

**Yuri I. Golovin 1,2, Alexander A. Gusev 1,3,4,\* , Dmitry Yu. Golovin 1, Sergey M. Matveev <sup>3</sup> and Inna A. Vasyukova <sup>1</sup>**


**Abstract:** This review describes methods and results of studying the mechanical properties of wood at all scales: from nano- to macro-scale. The connection between the mechanical properties of material and its structure at all these levels is explored. It is shown that the existing size effects in the mechanical properties of wood, in a range of the characteristic sizes of the structure of about six orders of magnitude, correspond to the empirical Hall-Petch relation. This "law" was revealed more than 60 years ago in metals and alloys and later in other materials. The nature, as well as the particular type of the size dependences in different classes of materials can vary, but the general trend, "the smaller the stronger", remains true both for wood and for other cellulose-containing materials. The possible mechanisms of the size effects in wood are being discussed. The correlations between the mechanical and thermophysical properties of wood are described. Several examples are used to demonstrate the possibility to forecast the macromechanical properties of wood by means of contactless thermographic express methods based on measuring temperature diffusivity. The research technique for dendrochronological and dendroclimatological studies by means of the analysis of microhardness and Young's modulus radial dependences in annual growth rings is described.

**Keywords:** wood; nano-, micro-, meso-, and macro-structure; multiscale mechanical properties; size effects; Hall-Petch law; dendrochronology

#### **1. Introduction**

Interest in wood and other cellulose-containing materials, composites in particular, had considerably increased by the beginning of the 21st century. The studies of wood nanoand micro-structures have especially intensified in the last decade (Figure 1) [1–5]. Several reasons can be named. Mineral resources (especially various metallic and nonmetallic materials, coal, oil, and natural gas) are being extracted at continually rising rates, and opencycle processing technologies create ever-growing volumes of industrial and household waste. This poses a threat to the biosphere because of the environmental pollution and increased carbon dioxide concentration in the atmosphere, while only a fraction of the manufactured materials are recycled and reused. The situation is aggravated by a sharp increase in polymer material manufacturing for packaging, which are seldom recycled and mostly non-biodegradable. The surging pressure on the environment requires more and more efforts for its neutralization.

In this regard, a wider use of biogenic materials as well as substituting them for traditional ones seems a promising step. Such cellulose-containing substances as modified

**Citation:** Golovin, Y.I.; Gusev, A.A.; Golovin, D.Y.; Matveev, S.M.; Vasyukova, I.A. Multiscale Mechanical Performance of Wood: From Nano- to Macro-Scale across Structure Hierarchy and Size Effects. *Nanomaterials* **2022**, *12*, 1139. https:// doi.org/10.3390/nano12071139

Academic Editors: Mohammad Malikan, Shahriar Dastjerdi and Takuya Kitaoka

Received: 2 March 2022 Accepted: 24 March 2022 Published: 29 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/).

wood and various agricultural vegetable wastes, and especially the nano- and microcellulose they contain, offer the best potential for numerous applications.

Cellulose is the most common natural linear polymer polysaccharide (C6H10O5)*<sup>n</sup>* in the biosphere. The materials formed on its basis provide vast advantages:


**Figure 1.** Growth dynamics of the number of research papers on the structure, properties, and applications of the materials comprising natural nano- and micro-fibers. Data obtained from Scopus using the following search parameters "TITLE-ABS-KEY" with keywords «Nanocellulose or microcellulose» on 28 February 2022.

Certainly, wood and other cellulose-containing materials have several disadvantages. They require special treatment, as they are flammable and hygroscopic. High humidity makes them lose some of their strength properties, while low humidity causes deformation. They succumb to rot and unwanted biodegradation. Besides, in their original, state the mechanical properties in every sort of wood strongly depend on the conditions of its growth, usage, and testing humidity, the structure of cell walls and annual growth rings, proportion of young and mature wood, stress condition, size of the sample or stressed area, and also the direction, rate, and duration of load application. The aforementioned considerations have obstructed identifying the universal patterns that form the mechanical properties of different species of wood. Nevertheless, some generalizations can be derived from the literature and from accumulated experience, as outlined below.

In the present review, we explore the methods and results of a multiscale study of the mechanical properties of various wood species, in connection with peculiarities of their nano-, micro-, and meso-scale structural levels of material organization. The analysis of literature data shows that, in a huge range of characteristic sizes of the structural units (about six orders of magnitude), mechanical properties of wood generally follow the Hall-Petch relation, which is well known in material science. The possibilities for non-destructive assessment of the mechanical properties of wood by means of contactless measurement of the temperature diffusivity tensor components are discussed, as well as using the scanning nanoindentation method for evaluating woods' micromechanical characteristics, in order to obtain dendrochronological and dendroclimatological data. The main scopes of the review are presented in Figure 2.

**Figure 2.** The main scopes of the review.

#### **2. The Hierarchical Structure of Wood**

From the point of view of material science, wood is a hierarchically organized natural composite with a complicated structure and a clear heterogeneity and anisotropy of all its properties, as well as an ability to regenerate [5,31–35]. In the wood architecture, one can distinguish, though only tenuously, several size and hierarchical levels (Figure 3), namely atomic–molecular, nano- (nanocrystals, nanofibrils), micro- (microfibers, cell walls), meso- (cells, large vessels, radial rays), and macro-level (annual growth rings, macroscopic

structural defects, cracks, etc.) [31–33,35]. They all contribute to forming the complex of physical, chemical, and mechanical properties [31,36]. A large range of the characteristic size of the structural components of wood (about six to eight orders of magnitude) and a wide scope of tasks and questions emerging from the study of this material all require a varied arsenal of research techniques and means to implement them. They will be briefly analyzed in the following section.

**Figure 3.** Scale hierarchy in wood structure and its main components. *R\** is the characteristic size.

Identification of patterns in the formation of the macro-properties of wood, as derived from its nano-, micro, meso-, and macrostructure, is the most important task in wood science. There are many reasons for the interest in the relations between macromechanical properties of wood and its nano- and microstructures, as well as physical characteristics, thermal characteristics in particular. Let us enumerate the most important ones. Firstly, the relevant patterns help to elucidate the nature and mechanisms of formation of the parameters most significant for practical applications of wood in the macroscale, i.e., its mechanical and thermal properties. Secondly, nanomechanical strength properties, being much higher than those at the micro- and macro-scale, indicate the potential for strengthening, which may approach the ultimate tensile strength of nanocrystalline cellulose (~10 GPa). Thirdly, the increased use of composite and nanocomposite materials in different spheres of engineering, construction, biochemical technologies, and medicine paves the way to replacing traditional metals and alloys with more lightweight and ecologically friendly composites. For example, the bodies of the most recent Boeing and Airbus airplane models consist, by weight, of more than 50% of fiber-reinforced composites. Their popularity in the auto industry, shipbuilding, sports equipment manufacturing, etc. is growing fast. However, glass and basalt fibers used for composite reinforcement, not to mention carbon micro- and nano-fibers, have some adverse properties from the ecological point of view; they are quite expensive and still unable to conquer the wide market for consumer goods. Cellulose fiber is by about an order of magnitude less expensive than fibrous glass while having almost

the same mechanical characteristics. Therefore, it is important to understand the nature of strength and damage mechanisms in microcellulose fibers, and to find approaches to improve their strength, thus enhancing the properties of textiles, non-woven materials, and the composites they are used to reinforce. When correlations are revealed between the mechanical characteristics and other physical properties, for example thermal properties, this information will be of great use for developing non-destructive contactless thermophysical methods for evaluating mechanical characteristics, instead of applying labor-intensive destructive techniques. Fourthly and finally, many tree species have a lifespan of several hundred or even thousands of years, with sequoia as an example. In their nano-, micro-, meso-, and macro-structure they accumulate a vast amount of information about the climatic conditions during their growth and about ecological catastrophes they have witnessed. The variations in composition and structure are inevitably reflected in the local physico-mechanical properties of wood. This natural archive can serve as a valuable source of information for climatology and for dating various events in earth's history (dendrochronology).

#### **3. Methods of Studying the Structural and Mechanical Properties of Wood at Various Levels of Scale**

The aim of the classic wood science is to discover and describe the dependence of woods' macromechanical, physico-chemical, and service properties on its inner structural characteristics, humidity, and external thermodynamic factors [6,37,38]. Wood type classification and its grading, according to mechanical properties, is an important pragmatic task [39]. Since the end of the last century, more and more attention is being paid to the fine structure of wood at the nanoscale. This interest was brought forward, on the one hand, by the growth in nanotechnology and nanometrology, and on the other hand by realization of what untapped resources are hidden at the nanoscale.

In the recent 15–20 years, numerous modern methods and means traditionally used in solid state physics and material science are being applied for studying the micro-structure and physico-chemical properties of wood [2,6,33,35,40].

Micro-structure is studied by means of transmission and scanning electron microscopy, scanning probes (mostly atomic force), confocal laser, and optical microscopy in various modes. Numerous X-ray methods are used to determine the composition and the parameters of atomic- and micro-structures. The character and degree of order of cellulose molecules in nano-fibers, the angle between the micro-fibers and the long axis of the cell, are determined by X-ray diffractometry and microtomography, as well as small-angle (SAXS) and wide-angle (WAXS) X-ray scattering. Elemental and molecular composition is revealed by spectroscopic methods, such as X-ray fluorescence, various types of spectroscopy such as infrared (IR), Fourier transform IR (FTIR), Raman, Brillouin, nuclear magnetic resonance (NMR), and other analytical methods. Together, they cover a huge spatiotemporal range of structures and events in them, namely more than twelve orders of magnitude in time and about eight orders of magnitude in length (Figure 4) [32]. The comparative analysis of the most widely employed physical methods for studying the molecular, sub-cellular, and cellular structures of wood can be found in most recent reviews [32–34,41].

To study mechanical properties at the nano- and micro-scale, a number of nano-/micromechanical testing (SSMT—small scale mechanical testing) methods [42–45] are employed. Atomic force microscopy (AFM) [46–49] and nano-indentometry (NI) [50–59] can be named as the most widely used ones.

They have similar structure flowcharts (see Figure 5) and capabilities [45,50,52]. In both cases, a high precision driver brings a probe, with the radius of its curve being from a few (in AFM) to a few tens (in NI) of nanometers, close to the studied surface and the probe starts interacting with it. The force *P* and penetration depth *h* of the probe are measured continuously, and their alteration kinetics are registered throughout the testing cycle (Figure 6a). Most commonly a *P*–*h* diagram (similar to a σ–ε diagram created during macro-testing) is built using the obtained data (Figure 6b), and standardized algorithms are applied to calculate about ten various mechanical characteristics, such as Young's modulus, contact stiffness, hardness, fracture toughness, creep rate, etc., at the nano-/micro-scale. In NI, a three-sided diamond Berkovich tip is used, as it is better calibrated from the point of view of the real tip geometry than the one used in AFM, thus providing more accurate and reliable quantitative data. Among the variety of proposed techniques for mechanical characteristics extraction from raw data, the method proposed and developed by Oliver and Pharr [60–62] has become the most widely used and has been incorporated into ISO standard [63], so that this method has been used for processing all NI experiments described in this review.

SSMT methods were used to examine the mechanical properties of individual cellulose nanofibrils and microfibers [64–68], cell walls [32,69–75], layers of early and late wood in annual growth rings [74], and to obtain plenty of other interesting data. However, there are very few papers that compare and analyze several scale levels at once [76–79]. Thus, connection between the properties of individual nano- and micro-structural elements of wood and their influence upon macro-mechanical characteristics cannot yet be traced.

The analysis of the structure and role of annual growth rings in shaping wood macroproperties requires, at the very least, one-dimensional, or better yet, two-dimensional scanning of certain physical characteristics. Three-dimensional imaging can be applied as well. Dating archaeological finds, works of art, climate changes and events, and ecological catastrophes based on the changes in growth ring structure and width is a separate issue. These approaches are known as dendrochronology, dendroecology, and dendroclimatology, respectively. The width of rings, proportion of early wood (EW) to late wood (LW), and changes in their morphology reflect the specific growing conditions each long-lived plant witnessed during each vegetation season.

**Figure 4.** Schematic of the spatiotemporal ranges of the most popular physical methods for studying wood structures. Adapted with permission from Ref. [32]. Copyright 2021, Wiley-VCH.

**Figure 5.** Schematic diagrams of (**a**) the nano-indentometer and (**b**) the atomic-force microscope.

While microstructure and physico-chemical properties of wood are studied with elaborate modern equipment, examination and analysis of annual growth rings for dendrochronological and dendroclimatological applications is carried out using simple optical methods, where primary information is derived from the difference in reflectivity between EW and LW. Quite often, the same approach is employed while assessing wood strength and other service properties. A detailed description of dendrochronological methods developed by the beginning of this century is given in [80]. These methods reveal only geometrical and morphological characteristics of the studied object (annual growth rings width, proportion of EW and LW in them, variations from ring to ring, etc.) and allow for comparison between the data obtained by different methods.

Numerous attempts have been made to improve traditional dendrochronological methods, mostly by modifications introduced to the sample preparation techniques, staining, use of blue light instead of white, application of computer vision technologies, and mathematical data processing (see e.g., [81–88]). However, despite this, the capabilities of the approach based on the analysis of transversal section images and photodensitometry

remain severely limited, as the reflective optical properties of wood are variable and their connection with other wood characteristics, such as mechanical and thermal, are either ambiguous or very weak.

**Figure 6.** Two methods of representing the data obtained by means of normal nanoindentation: (**a**) as kinetic curves *P*(*t*) and *h*(*t*); (**b**) as a *P*−*h*-diagram (*P*—force, *h*—indenter displacement). Adapted with permission from Ref. [45]. Copyright 2021, Springer Nature. The circled letters (from A to E) mark the characteristic points on the loading curves and the indenter position relative to the sample surface. Inset F shows the vector diagram depicting correlation in the complex plane between the vectors of oscillating force and the resulting indenter displacement in the CSM method. Five loading regimes are marked from 1 to 5. The indices at *P* and *h* mean as follows: *up*—increase; *cr*—creep; *down*—drop; *e*—elasticity; *v*−*e*—viscoelasticity; *max*—maximum value; *We*—elastic energy; *Wpl*—energy absorbed and dissipated by the sample in a single load–unload cycle.

In order to expand the capabilities of the analysis of mechanical properties in their connection with the architecture of wood ring structures, the following methods have been used: two-dimensional mapping of properties on cross-sections of tree trunks by AFM methods [89–92] and NI [93] scanning, 3D X-ray [94,95] and NMR tomography [96], and synchrotron-based X-ray microscopy [97]. However, these methods are complicated, labor-intensive and require expensive or unique equipment; therefore, they are used only sporadically. The method of X-ray densitometry [98] presents fewer difficulties, but it requires access from both sides to a perfectly flat sample cut exactly perpendicular to the long axis of wood cells.

It should be noted that the mechanical properties of wood and cellulose-containing materials show a significant dependence on the rate of monotonous loading, oscillating load frequency, and duration of load application. They can vary between samples and change over time in a significantly greater range than in similarly structured technogenic composites (e.g., in glass and carbon fiber-reinforced plastics) [2,5,40,41]. Such variability of properties makes identification of common regularities in their formation even more challenging.

To sum up, we should mention that, thus far, the overall links between the properties of all the scales and hierarchical levels—from cellulose nanocrystals (CN) to macro-samples—require additional study. However, there is a considerable volume of information on every individual level. The following sections present the examples of the most representative data from the lowest to the highest scale levels of the structure.

#### **4. Nanocellulose and Elementary Nanofibrils**

Cellulose is the most common natural polymer and the major structural component that provides strength to wood and stems of grass, reed, bamboo, etc. [5,40]. Cellulose is a macro-molecular polysaccharide (C6H10O5)*n*, consisting of linear chains of several tens to many hundred *n* of β-(1→4) linked glucose molecules (Figure 2). In its origin, cellulose can belong to three types: plant, regenerated, and bacterial [5]. The current state of affairs in extraction and functionalization of cellulose nanofibers is described in Handbook [5] and in the most recent reviews [99–102].

Cellulose molecules easily form nanocrystals with a lateral size of 3–10 nm and being 100–300 nm long. These nanocrystals form nanofibrils 5–20 nm in diameter and up to many hundreds of nanometers long. Inside, nanofibril cellulose is present in an amorphous– crystalline state as a series of alternating domains. The amorphous phase, to some extent, reduces the strength of the nanofibril, but makes it more supple and elastic. The most typical structural characteristics of nanocellulose-based formations are presented in Table 1 [67].

**Table 1.** Structural characteristics of cellulose nanocrystals and individual nanofibrils (CNC and CNF, respectively) [67]. Copyright 2017, Wiley books.


Mechanical properties of nanocrystalline cellulose have been characterized by various methods, including calculations of bond strength inside macro-molecules and between them, computer generated simulations, processing data from IR and Raman spectroscopy, AFM, WAXS, and others [5,40,65,68]. A brief overview of the mechanical characteristics of nanocellulose is given in Table 2 [67]. The variability of data is explained by specific characteristics of the calculation schemes, models, raw data processing algorithms, and also by the difficulty of carrying out direct measurements at the nanoscale. The differences in age, structure, and origin of wood affect the experimental results as well. Besides, the mechanical properties of nanocellulose samples depend on their size significantly. For instance, transversal Young modulus reduction by a factor of 1.6 has been reported in [78] for increasing NC diameter from 2.5 nm to 6.5 nm.

**Table 2.** Mechanical characteristics of cellulose nanocrystals and individual nanofibrils (CNC and CNF, respectively) [67]. Copyright 2017, Wiley books.


#### **5. Cellulose Microfibers**

The typical hierarchy of wood structure at higher levels continues with nano- and micro-fibers. They are formed by elementary nanofibrils, mainly due to hydrogen bonds. Nanofibrils form strands surrounded by a matrix composed of lignin (an aromatic polymer polyphenol), hemicellulose (low molecular weight branched polysaccharide), pectin (gelforming polysaccharide), and water [5,40,68]. Cellulose content in the fibers can vary in a wide range. For example, it is 40–60% in the wood fibers of various species and can exceed 96% in cotton fibers [5,23,24,31,68].

Nano- and micro-structures of cellulose materials and their properties strongly depend upon the specifics of interaction between nanocrystals in elementary fibrils and the ordering and binding of the latter in nano- and micro-fibers [103–107]. Mechanical, strength in particular, properties of cellulose nano- and micro-structures are structure sensitive, just as those of most other other organic and non-organic materials. In turn, their morphology and inner structure depend upon plant species, their growth conditions, and cellulose extraction technology. The dominant role in determining fiber properties belongs to the cellulose content in nanofiber, the degree of its crystallinity, and specifics of nanofiber binding at the material. The angle between the nanofibril axis and nanofiber or cell axis has significant impact too. A comprehensive review [104] contains various data concerning morphology, microstructure, and mechanical properties of micro-fibers of various origin (Figure 7) and examples of their application for polymer composite reinforcement.

**Figure 7.** Dependence of cellulose fiber tensile strength upon Young modulus for various plant materials. Adapted with permission from Ref. [104]. Copyright 2018, Elsevier.

The strongest of the studied micro-cellulose fibrils have demonstrated Young modulus *E* = 75–85 GPa and tensile strength σ*<sup>b</sup>* = 1.6–1.7 GPa, so that the ratio *E*/σ*<sup>b</sup>* ≈ 50. One of the possible techniques allowing researchers to reach such high mechanical properties is described in [66]. The authors have used the efficient technique of double hydrodynamic ordering of nanocrystals and nanofibrils to produce the fibers with diameter 6–8 μm. Their tensile strength reached 1.1 GPa. Nanofibril cross-linking has increased the fibers' strength up to σ*<sup>b</sup>* = 1.57 GPa.

As follows from fundamental considerations, the theoretical strength of any defectfree material can reach 0.1 *E*, while the strongest micro-cellulose fibers mentioned above have values around 0.015–0.020 *E*. Hence, even the strongest studied micro-fibers have the potential of increasing their strength by 3–5 times.

It should be mentioned that the data concerning cellulose micro-fibers mechanical properties differ significantly depending on the measurement technique (see Table 3) [104]. Results obtained using AFM and NI are in agreement with each other regarding the measurement accuracy, despite using different probes and measurement techniques, so

that they are just as reliable as the ones obtained using the undebatable method of uniaxial tension. Usually, its tensile strength is two to three times higher than compression strength or hardness [5,40,104,106,107], unlike void-free materials, where the reversed value is quite typical. For example, the Tabor ratio for metals is well known, where hardness is three times higher than the yield stress. We suppose that this difference is due to specific nanofibril behavior when subjected to tensile and compressive stress or in hardness measurements. In the first case, the molecular chains are strained and partially oriented along the fiber axis, which increases their strength. Indentation, on the other hand, used both in NI and AFM, promotes the arising of compression stress that causes micro-fibril bending, and micro-fibril buckling failure, which occurs earlier than their failure in uniaxial tension.


**Table 3.** Microfibers Young moduli *E* obtained using different methods for three materials [104].

For many applications such as aviation, space aeronautics, the automotive industry, sport equipment etc., the most important mechanical characteristics are not absolute but specific ones, i.e., normalized on material density *ρ*. Figure 8 shows the absolute σ*b*, specific strength σ*b/ρ*, and Young modulus *E* for highly oriented cellulose nano- and micro-fibers when compared to those for macroscopic wood and other materials. As it can be seen, the specific strength of defect-free nanocellulose can be manifold higher than that of aluminum or titanium alloys or constructional steels. However, nano- and micro-cellulose materials are inferior to metals at thermal and crack resistance, failure deformation, and other related energy characteristics. Only some polymers such as kevlar and carbon microfibers can contest nanocellulose at specific strengths. Single-wall carbon nanotubes and graphene are manifold superior at specific strength to any other known material.

**Figure 8.** Mechanical properties of nanocellulose (NC) and cellulose microfibers (MC) in comparison to common and perspective constructional materials. Crosshatched areas denote absolute values of tensile strength *σb*, while non-crosshatched ones are the strength normalized over material density *ρ*.

Lastly, let us mention one more advantage of natural cellulose fibers produced from the wood. It is several times cheaper than the flax fibers and an order of magnitude cheaper than ecologically unsafe glass fibers widely used in composite reinforcement applications [104]. Highly ordered cellulose microfibers have nearly the same strength as glass fibers already and have a good prospects of further strength increases.

#### **6. Cells and Cell Walls**

While a tree grows, cellulose microfibers integrating with other components such as lignin, hemicellulose, pectin, water, etc., form walls of cells that are highly elongated in the direction of tree trunk axis. Several layers are discerned within cell wall including primary wall and multilayered secondary walls that usually consists of three layers, named S1, S2, and S3, which differ in the angle μ between cellulose microfibers and cell long axis. The secondary wall provides the main contribution to cell stiffness and mechanical strength. Cell size diminishes, cell wall width increases, and the cross-section of internal capillary reduces while going from early wood (EW) that is a part of the annual ring formed at the first stage of vegetation, to late wood (LW), formed at the second stage.

To study mechanical properties of wood cells, various SSMT methods are used [45], and the most widespread ones are AFM [89,92] and NI [70,71,74,75,93]. Let us present some typical examples of NI application with load *Pmax* = 0.1–1 mN to this problem. The authors of [74] studied radial dependence of cell wall longitude Young modulus *ENI* and nanohardness *HNI* in two annual rings of common pine (*Pinus sylvestris* L.) wood, corresponding to the ages of 7 and 74 years. As could be seen at Figure 9a, in going from EW to LW, *ENI* increases by nearly 50%, while *HNI* (Figure 9b) increases by just 5–7%.

**Figure 9.** Dependencies of cell wall longitude Young modulus *ENI* (**a**) and nanohardness *HNI* (**b**) upon cell sequential number in the annual growth ring for two rings with ages of 7 and 74 years. Adapted with permission from Ref. [74]. Copyright 2020, PAN.

A number of other papers report similar data supporting that cell wall nanohardness varies not too much at different layers, rings, or even tree species. For instance, the following results are reported: *HNI* = 0.35–0.42 GPa for *Pinus massoniana* Lamb. In [107], *HNI* = 0.41–0.53 GPa for *Masson pine*, coinciding within the measurement accuracy for EW and LW in [108] and *HNI* = 0.34–0.54 GPa for *Pinus taeda*, not discerning EW and LW, in [109]. Similar results are reported for *HNI* in cell walls junction through the middle lamella for *Norway spruce*. Nanohardness in the cell corner middle lamella was estimated to be 0.34 ± 0.16 GPa [110].

The NI technique allows more detailed measuring cell wall elastic properties and determining the main components of elasticity tensor. So, the measured value of Young modulus of secondary wall S2 has been reported to be 26.3 GPa in the longitudinal direction and 4.5 GPa in the lateral one [111].

The most informative experiments are carried out in situ in an electron microscope column using a sharp indenter or flat piston [112]. Simultaneous recording of loading diagrams and obtaining visual information concerning the nano-/micro-structure evolution allow the studying of the micro-mechanisms of deformation and failure [113].

There are a number of papers studying the mechanical properties of micro-pillars cut from cell walls by focused ion beams (FIB). The pillars were subjected to uniaxial compression in a scanning electron microscope column (SEM) [114,115]. It allows simultaneous recording of σ-ε diagram and specific features of micro-pillar deformation.

Such works are not numerous due to high labor content and rather complex and expensive equipment required for sample preparation. However, they provide direct confirmation of quantitative information obtained using NI, allow for obtaining unique information concerning various deformation modes, buckling, and failure mechanisms of cell walls, and verification of various behavior models of wood hierarchic structures under load.

#### **7. Annual Growth Rings**

As follows from the data presented in the section above, cell wall Young modulus *ENI* does not vary significantly with the indent location, be it late or early wood layer, secondary wall, or cell conjugation region with middle lamella. The ring age, weather conditions during its formation, or particular tree species does not change it more than by a factor of 1.5–2. Nanohardness *HNI* dependence upon these factors is even weaker. Nevertheless, macromechanical properties woods of differing origin can differ manifold reaching up to an order of magnitude or even more. Evidently, this weak correlation between nano- and macro-properties is due to the difference in cell wall thickness, the relative share of late wood, and the number of large tracheides and other wood structure elements which reduce wood macroscopic strength. To close the gap between nano- and macro- scale mechanical properties, the nanoindentation tests were carried out at loads ranging from 5 to 500 mN and reported in the set of papers [116–119], unlike the 0.1–1 mN range usually used in studying cell walls. It extended the deformed region over the whole cell or several cells up to 50–150 μm, as opposed to precise targeting at the cell wall.

The values of Young modulus *Eeff* and microhardness *Heff* obtained this way can be considered as effective, as long as, being the result of averaging over the cells structure, they incorporate not only mechanical properties of cell walls but also their thickness, material porosity, and microdefects, just as in macroscopic mechanical tests. However, indentation size was at least an order of magnitude less than annual ring width. It allowed obtaining *Eeff* and *Heff* spatial distributions across several annual rings.

Typical *Eeff* and *Heff* radial dependencies measured at cross-sections of common pine (*Pinus sylvestris* L.), which represent coniferous trees, are shown at Figure 10, and pedunculate oak (*Quercus robur* L.), which represent hardwood trees, are shown at Figure 11 [116,117]. As could be seen, both species manifest pronounced periodicity of local mechanical properties. Positions of abrupt changes in *Eeff* and *Heff* coincide with annual growth ring boundaries determined by wood color change using the standard optical method. Changes of *Eeff* and *Heff* in transition from EW to LW within each annual ring are gradual in pine and abrupt in oak. The linear dependence between *Eeff* and *Heff* with almost the same slope *m* = 0.017 ± 0.002 has been observed (see Figures 10b and 11b). In other words, the ductility ratio DR = *Eeff*/*Heff* is found to be around 60 for both tree species. This value is quite typical for many other species, for example DR in gum-tree lies within the 54 to 68 range, with an average of 61 [120], and in beech it is around 55 [121].

As could be seen at Figures 10a and 11a, *Eeff* and *Heff* vary not too much within each EW layer and across several EW layers (relative variance of both values is around 10–15%), despite that the weather conditions during layer growth may differ significantly. For instance, the year 2010 was very dry. It affected the ring width *wo*, so that it is less than half of an average one, however, *Eeff* and *Heff* values for EW are almost the same as appropriate values for other years (Figure 10a). The lateral size of EW cells in different rings does not

differ significantly either, but cell wall width does. Thus, ring width variation is mainly due to the difference in cell morphology and its quantity in the layer, while the mechanical properties of cell wall material are almost the same.

**Figure 10.** Micromechanical properties of common pine annual growth rings measured at *Pmax* = 500 mN [116]. (**a**) Spatial dependencies of *Heff* and *Eeff* over radial distance *r* for six successive rings. Ring boundaries are shown using dashed lines. The extraordinarily draughty 2010 year is highlighted by red color. (**b**) Dependence of hardness *Heff* upon Young modulus *Eeff* for six successive rings.

**Figure 11.** Micromechanical properties of pedunculate oak annual growth rings measured at *Pmax* = 500 mN [116]. (**a**) Spatial dependencies of *Heff* and *Eeff* over radial distance *r* for six successive rings. Ring boundaries are shown using dashed lines. (**b**) Dependence of hardness *Heff* upon Young modulus *Eeff* for six successive rings.

As follows from the data presented, the ratio of averaged Young moduli for LW and EW is around 3.1 for pine and 3.5 for oak. The ratio of averaged hardness for LW and EW is close to those measurements, and equal 3.7 for pine and 3.0 for oak. These values stay for average values calculated for individual rings, but in some years, they can differ from the mean substantially. Obviously, such variations are due to anomalous weather conditions during these years, and narrower growth rings corresponding to these years evidence the same. However, year-to-year variation in mechanical properties is much higher than in annual ring widths, which is usually used for climate reconstruction. Thus, *Eeff* and *Heff* measurements can be a much more sensitive and accurate dendroclimatological method than annual growth rings width measurement.

The values of *Eeff* and *Heff* measured as described above are two to three times lower than those of cell walls. This is expected as long as the former are affected by material porosity. However, *Heff* is two to three times higher than macroscopic Brinell hardness [8,76]. It can be formally qualified as size effect (SE) in wood mechanical properties. However, determining the relative contributions of void-free material properties, porosity, and the related differences in deformation mechanisms in such SEs requires separate research.

#### **8. Mechanical Properties of Wood at Macroscopic Scale**

The largest part of the literature discussing mechanical properties of the wood concerns properties at the macroscopic scale [4,8,32,35,122]. As long as mechanical properties of the wood demonstrate prominent anisotropy, reference data typically comprise Young modulus, hardness, strength, and other properties for directions along and across fibers separately [35]. The methodological review of generally used approaches and experimental techniques of mechanical testing of the wood can be found in [123]. Standard methods of timber mechanical testing are described in [124]. Methods of timber strength classification can be found in [35,125].

The most general relationships between the mechanical properties and structure of the wood are discussed below. Going from juvenile wood (JW) that is formed during the first 5–20 years of tree growth to mature wood (MW), specific gravity, cell length and wall thickness, percentage of late wood, and strength increase, while fibril angle μ, moisture content *W*, and annual ring width *w* decrease. Moisture content *W* increases by 1% in the range between 10–12% and 50–60% results in a decrease in uniaxial compression strength by 5% and in uniaxial tensile strength by 2–2.5% both along and across fibers [126]. Young modulus across the fibers diminishes with growing *W* too but at the slower rate of around 1.5% per 1% of *W*. Typical tensile strength of softwoods (fir, pine, spruce, cedar, etc.) along the fibers is between 45 and 112 GPa, while that of hardwoods (beech, oak, maple, elm, etc.) is between 70 and 120 GPa [8]. Cooling from room temperature to –195 ◦C results in increase in compressive strength of dry wood with *W* = 12% by a factor of 2–2.5, while heating to 50 ◦C rises it by 10–20% [8]. Holding the wood under load for a long time diminishes its strength. So, it drops by 10–15% for an hour, by 20–25% for a month, and 30–35% for a year when compared to 1 min load [8]. Fracture toughness of macroscopic wood samples range from 250 kPa m1/2 for Western white pine to 517 kPa m1/2 for Yellow poplar [35]. In addition to the species, growth conditions and moisture content, structure defects significantly affect wood macromechanical properties [8,122,123]. Their variation can reach 15–35% from sample to sample due to such sensitivity. More detailed data concerning mechanical properties for various sample size and testing conditions can be found in [4,8,32,35,122,125,126] and Table 4 in Section 10.

#### **9. Modification and Hardening of Wood and Cellulose**

Drastic reduction of wood mechanical characteristics at increasing characteristic size of the samples stimulates development of various modification techniques for both macroscopic wood products [30,127,128] and nano- and micro-cellulose [5,129,130]. Several classes of wood modification techniques could be distinguished, including (a) chemical processing (acetylation, furfurylation, resin impregnation etc.); (b) thermally-based processing; (c) thermo–hydro–mechanical processing (surface densification); (d) microwaves, plasma, and laser light treatment; (e) mineralization; (f) biological treatment [30]. Besides improvement of mechanical properties, wood modification can be aimed at reducing water absorption or susceptibility to rotting and biodegradation, enhancing fire resistance or antiseptic properties, improving dimensional stability or resistance to acids or bases, ultraviolet radiation etc. Getting back to mechanical properties, let us mention some examples of wood modification leading to significant improvement of such properties. So, the widely used Compreg technique, consisting of wood compression before the resin is cured within the material, leads to considerable compressive strength increases that are even higher than wood density increases; tensile strength, and flexural strength increase less than its density increases [30]. As long as wood density can be increased up to a factor of 2–2.5 during this processing, its strength can be increased nearly twofold. At the same time, material toughness drops by 25–50%. The hardness can be raised by Compreg more substantially, and the factor of 10 to 20 has been reported [30]. Another modification technique is described in review [31]. The two-step process involves the partial removal of lignin and hemicellulose from the natural wood via a boiling process in an aqueous mixture of NaOH and Na2SO3, followed by hot-pressing, leading to the total collapse of cell walls and the complete densification of the natural wood with highly aligned cellulose nanofibers. The processed wood has a specific strength higher than that of most structural metals and alloys, making it low-cost and high-performance. More detailed information concerning techniques and results of wood modifications aimed at changing the mechanical and other service properties of the wood can be found in books [8,30,127], reviews [128], and original papers [131–136].

#### **10. Size Effects in Wood**

All the size effects in materials are usually divided into internal ones that depend on nano- and micro-structures of the object, and the external ones, depending on the size and shape of the sample, loading method, and size of the area under load. There is not much systematized information available regarding both types of SEs at different scale hierarchical levels of wood structure. In our survey, we present the most interesting and typical data. They are summed up in Table 4.

As follows from Table 4, the strength of cellulose nanocrystals, assessed both by calculations and by experimental techniques, is 4.9–10 GPa, while the strength of nanofibrils with a diameter 3–15 nm is close to the lower end of this range [5,40,65–68]. These values exceed the strength of cellulose microfibers 8–12 μm in diameter, which is 0.5–1.65 GPa, by about one order of magnitude [5,40,66,104,106,107]. The typical nanohardness values *HNI* for cell walls with a thickness of 2–5 μm are about 0.3–0.5 GPa [70–75,107–110], which is 2–3 times less than the strength of cellulose microfibers. As it is shown in Figures 10a and 11a, the effective values of microhardness, *Heff*, which take into account the porosity in layers EW and LW, are several times lower (from two to four times) than *HNI*. However, the *Heff* value is several times higher than Brinell macrohardness HB [76], and bending and uniaxial tensile strength obtained in macrotests [8].

It is evident that the values of effective Young's modulus and hardness at mesoand macro-levels fall so dramatically not only because of the internal reasons defined by molecular and supramolecular structures, but also due to a great influence of nano- /micro-porosity of the material that should be attributed to the number of pores, capillaries, and larger tracheides with a high aspect ratio, ubiquitous in in the wood structure. Their presence results in several considerable differences between the mechanical behavior of wood and non-porous bodies. Firstly, Tabor's rule, according to which the hardness of soft materials exceeds their yield stress or strength by about three times, is almost never met. On the contrary, in most wood species, their macroscopic hardness is several times lower than the yield stress and ultimate tensile strength determined for macroscopic samples. Apparently, in all similar events, the reason is that during indentation and uniaxial compression, the tested wood cell structure loses its stability long before the manifestation of plastic deformation and tensile rupture. These events strongly depend on the direction of the load application in relation to the long axis of the cell due to anisotropy of the mechanical properties of wood. In the longitudinal direction, they are higher by about an order of magnitude than in the transverse direction. It is difficult to provide an explanation for the various possible modes and mechanisms of wood deformation at the nano- and micro- scale because of insufficient experimental data. However, the quantitative information quoted above leads to the conclusion that there exist highly pronounced SEs in wood which cause a sharp decrease in strength/hardness from ~10 GPa in nanocrystalline cellulose to ~0.1 GPa or less in macrovolumes of wood. This means that all cellulosecontaining materials have a large strengthening potential that can be realized through optimally organized nano- and micro-structures, and employment of relevant technologies.

*Nanomaterials* **2022**, *12*, 1139


**4.**Somemechanicalpropertiesofwoodstructurecomponents.


**Table4.***Cont.*


44

Representation of the strength characteristics as a function of characteristic dimensions of the structure *R\** in double-logarithmic coordinates provides a distinctive hockey stickshaped diagram (Figure 12). The descending part of the curve in the nano- and microdomain has a slope close to −0.5, a feature it shares with Hall-Petch relation, which is well-known in materials science and described for the first time for polycrystalline metals more than 60 years ago [137–139]:

$$
\sigma\_y = \sigma\_0 + A(d)^{-0.5} \tag{1}
$$

where *σ<sup>y</sup>* is the yield stress, *d* is the crystallite size, *σ*<sup>0</sup> and *A* are material constants. In most cases, the value of *A* that used to be called the Hall-Petch constant turned out to be deformation-dependent [139]. Later it was clarified that external dimensions affect the mechanical properties in the similar way, though the index of power may differ from 0.5 [45].

**Figure 12.** The size dependence of cellulose-containing materials strength. (Data collected by the authors). *R\** is the characteristic size.

Relations, similar to (1) were discovered for hardness [140]

$$H = H\_0 + A(R^\*)^{-0.5} \tag{2}$$

where *R\** is the characteristic dimension of the locally deformed area, which, in the process of indentation with a Berkovich tip, is usually taken equal to the indentation depth *hmax*.

The specific internal and external SEs have been studied not only in void-free poreless materials (metals, alloys, rocks, composites, etc.), but also in such porous materials as ceramics, solidified foams [141–144], and organic gels [145].

Due to a discrete character of damage accumulation in porous materials under load, the loading diagrams obtained in NI and in microsample deformation contain deformation

jumps. They carry information about elementary events, their rate and statistics as a function of size of the area under load, deformation rate and other experiment conditions [143].

Obviously, the causes for drop in strength/hardness with increase in *R\** may differ in different groups of materials. Nevertheless, some similarities can be observed, for example, the principle "the smaller the stronger" works in wood as well, though the index of power for *R\** may vary within quite a wide range, namely from 0.2 to 1. Therefore, many authors suggest other dependences to account for SEs, such as 1/*R\*,* ln(*R\**)/*R\** and others [146,147].

The softening mechanisms turning on with the increase in both internal and external characteristic size of the system require additional study of interrelations between the multilevel nano-/micro-structure of wood and its physico-mechanical properties. However, there are reasons to suppose that micromechanics of thin filaments, walls, partitions, as well as the macromechanical behavior of wood may have a lot in common with those in other highly porous materials [148–154]. Therefore, the general approaches and models developed for the analysis of the latter can be applied to wood as well. Plausibly, such softening depends upon cell wall width, cell morphology, aspect ratio and adhesion, and percentage of the lumens in wood crosscut area.

Some SEs can also be observed at the macro-scale and, to some extent, at mesoscale, though they are less pronounced than that at the nano- or micro-scale. They can be attributed to a growing possibility of large defect (cracks, delaminations, knots, and other wood defects) occurrence in larger objects, and can be described using Weibull statistics [155,156]. However, such a study lies outside the scope of our present review.

#### **11. Nanomechanics in Dendrochronology**

Jumps in *Eeff* and *Heff* at the edges of annual growth rings made it possible to measure their width *wNI* using scanning NI. Then, the comparison was carried out with the *w*<sup>0</sup> value determined by an optical method (analysis of image contrast). The image processing method was similar to the one used in the widely used LINTAB equipment. The comparison of the data obtained by these two techniques for measuring tree-ring width is presented in Figure 13; one can see that differences between them do not exceed 2–3% for pine and 4–5% for oak and lime. Mean average deviation for six to seven rings was about 2%. In effect, this means that the scanning indentation method can be used as an alternative to the optical one, or can complement it, providing some additional data on the local mechanical properties.

**Figure 13.** The results of annual growth ring width measurements obtained by nanoindentation *wNI* and by the optical method *w*<sup>0</sup> (**a**), and the discrepancies between these methods (**b**) [117].

#### **12. Correlation between Thermal Diffusivity and Mechanical Properties of Wood**

Kinetic thermophysical characteristics (thermal conductivity λ and diffusivity *a*) and mechanical properties of wood (Young's modulus, strength, hardness) depend on the same factors, namely on composition, structure, density, porosity, humidity, and specifics of interconnections between microstructural units [5,8,10]. Moreover, in both groups, the same pattern can be observed: higher wood density is accompanied by higher values of the characteristics mentioned above. So, it seems reasonable to suggest that there is a certain correlation between these two groups of properties. Once revealed, such relations could allow switching from laborious and material extensive destructive mechanical tests to non-destructive contactless measurements of λ or *a*. Such approaches could be used to estimate relative mechanical properties, and to sort and grade materials and products made of wood, fiber-reinforced composites, etc.

It should be noted that despite a wealth of information concerning the mechanical and thermal properties of natural and modified wood, wood-based layered materials, as well as composites reinforced with artificial and natural organic fibers, they are measured in separate tests executed on different samples. Exceptions are several papers that employ thermal non-destructive testing to estimate structural damage of wood [157,158], defects [159,160], porosity [161], and anisotropy [162,163] of composites and their possible impact on material mechanical properties. A micromechanical model for the prediction of effective thermal conductivity in two- and three-phase composites is proposed in [164].

The dependence of thermal diffusivity tensor *aij* on Brinelle hardness HB in common pine (*Pinus sylvestris* L.), pedunculate oak (*Quercus robur* L.), and small-leaf lime (*Tilia cordata Mill.*) wood at various humidity levels was studied in [116,165,166].

The thermal diffusivity tensor components *aij* = λij/*ρCp*—where λij represents the thermal conductivity tensor components, *ρ* is the material density, and *Cp* is the specific thermal capacitance—were measured using an original non-destructive thermal imaging technique described in detail in [167–169]. The method is based upon local stepped heating at small spots on the sample surface by a focused laser beam while continuously monitoring the surface temperature distribution with a thermal camera. Heat propagation in such a setup is close to spherical symmetry in isotropic materials, while in orthotropic materials, the isothermal surfaces are close to three-axis ellipsoids with the axes fully determined by the main components of *aij* tensor and the time elapsed since heating onset, provided that the distance to the heating center is at least several times higher than the heating beam radius. Therefore, the *aij* values were determined by processing dynamic thermal images obtained on lateral, radial, and transverse crosscuts of wood samples, as described in [167–169].

The values of hardness HB*<sup>W</sup>* at current humidity values *W* normalized to the hardness HB*W10* at *W* = 10% for the lateral and radial faces were statistically indistinguishable (Figure 14); therefore, they were approximated as a single set by the following linear function (HB*W*/HB*W*10) = 1.59(*a*l/*an*) − 1.34 with the coefficient of determination *<sup>R</sup>*<sup>2</sup> = 0.75, where *a*<sup>l</sup> and *an* are longitudinal and lateral components of *aij* tensor, accordingly. The hardness measured on the face perpendicular to the fibers was found to be independent of tensor *aii* components, so that it could be used for method and equipment calibration.

Similar results were obtained for other porous materials. For example, a linear relationship between the coefficient of thermal conductivity λ and compressive strength σ*<sup>b</sup>* in lightweight porous cement composites containing aerogel has been reported in [170]. A decrease in porosity resulted in λ growth from 0.39 to 0.47 W/m·K, and simultaneous σ*<sup>b</sup>* increase from 7.5 to 15 MPa.

**Figure 14.** Dependence of the relative macrohardness HB*W*/HB*W*<sup>10</sup> of pine wood upon thermal diffusivity anisotropy—*a*1/*an*. HB—Brinell hardness (12.7 mm sphere diameter, 1 mm indentation depth), *a*1, *an*—thermal diffusivity coefficients along and across the fibers. HB*W*/HB*W*<sup>10</sup> is HB*<sup>W</sup>* values on radial (A) and tangential (B) faces, as shown in Figure 2, normalized to HB*<sup>W</sup>* at humidity *W* = 10% (HB*W*10) [116].

#### **13. Discussion**

Despite the great variety of plant domain and trees in particular, their structural composition is similar, and this similarity stays across different hierarchical scale levels. They are composed of the same organic molecules that form nanocrystals and nanofibrils, which in turn form microfibers and then cell walls. The tree trunk contains the cells organized in annual growth rings. It is interesting that not only such a structure is more or less universal, but its individual elements in various species have similar morphology and close mechanical properties. Additionally, the higher the scale level in this hierarchy the lower its mechanical characteristics. This decrease is roughly in accordance with the Hall-Petch relation, and this fact has not been noticed before. Reconciling mechanical strength and physiological functions, nature had to find some compromise. Capillaries, large tracheids, and gum canals weaken wood structure inevitably, so that their specifics determine in large part the macromechanical properties of timber products. Relations between wood mechanical properties and the architecture and features of its structure are studied insufficiently yet. In particular, quantitative relations between the polymerization degree of polysaccharides, microfibers and cell walls structure, morphology and shape of the cells, wood microarchitecture, turgor pressure etc., on one hand, and macroscopic mechanical characteristics such as stiffness, hardness, strength, strain at break, and fracture toughness on the other hand, are unknown. Meanwhile, the turgor of a plant cell, for example, averages around 0.44 MPa and can be as high as 2 MPa [35]. It is essential in providing mechanical stability of any plant. It is known that higher cellulose percentage and crystallinity generally provide higher strength and stiffness to microfibers [8,36,40,67,68], as long as other structural components have much lower mechanical characteristics. In particular, hemicellulose stiffness and strength are several tens times lower, and around two orders of magnitude lower for lignin than that of cellulose [36]. Besides, an increase in moisture content leads to much faster deterioration of their mechanical properties than that in cellulose. However, reliable systematic research into cellulose, hemicellulose, lignin, and pectin percentages impacts on wood mechanical properties at different structural levels are lacking. Nevertheless, the role of chemical composition is very significant, especially accounting for the high variance of the above percentages, even in plants of the same species, and more so for plants of differing species [104].

A better understanding of formation regularities of wood hierarchical structure and their relations with mechanical properties could allow for controlling the latter, as well as manufacturing materials with desired properties. In theory, it could be achieved by means of controlled growth of the plant in an artificial environment, modification of wood structure using more efficient techniques, or reconstruction of the material that could provide it with new design and properties. All these approaches are used to some extent already, but without physical substantiation their efficiency is low. The best results have been achieved at nano- and micro-scales by growing nanofibers, and aligning and crosslinking them to form high strength microfibers. A scientifically valid and optimized design of nanocellulose-based materials can fulfill natural potential of this polymer and allow manufacturing ecologically friendly, high-strength, economically efficient materials.

#### **14. Conclusions**

In spite of the fact that people have been using different kinds of wood and wood derivatives for many thousands of years, the variety of cellulose-containing materials is so vast, and their multilevel architecture is so complicated, that the nature of their properties has not been completely revealed yet. As a result, the potential of these natural materials is not fully used. In order to make accurate predictions about wood characteristics, to be able to control them during growth, or to modify them with various treatments, and to preserve them in service, we need a deep comprehension of which elements of woods' structure, their condition, orientation, interconnection, and evolution, as well as which defects and flaws in the nano-, micro-, meso-, and macro-structure, influence particular macroscopic properties. First of all, these properties should be understood as the most common and important key elements of wood microstructure, such as nanocrystalline and amorphous phases of cellulose, microfibrils, cell walls, and annual growth rings.

Various modern research methods provide multiscale data on wood mechanical properties at different structural levels, from nano- to macro-scale. The interrelated information from these levels can offer new approaches to the cultivation of wood with predetermined mechanical properties, for example, with high strength and elastic properties, necessary acoustic characteristics, or low creep rate, and would help to create effective physically substantiated methods of wood strengthening. High specific strength of nanocellulose, exceeding that of almost all the other construction materials, except nanocarbons, encourages designing new extra-strong, ecologically friendly materials that would take advantage of this potential.

The results of scanning nano- and micro-mechanical properties of wood across several successive annual growth rings provides the basis for innovative methods and techniques for dendroclimatological and dendrochronological applications that would complement the existing approaches. Since the average cell size in a wood cross-section is about 30–50 μm, and the average ring width is 1–3 mm, there are approximately 50–100 cells in each ring. Technically, nanoindentation provides a means to measure the mechanical characteristics of each individual cell. Therefore, the temporal resolution limit for NI in dendrochronological applications is close to one week, making it a much more precise method with better temporal resolution than the traditional optical ones.

**Author Contributions:** Methodology—Y.I.G.; project administration—A.A.G., I.A.V.; supervision—Y.I.G.; data processing—D.Y.G.; writing—draft, review & editing Y.I.G., D.Y.G., S.M.M., A.A.G.; visualization— S.M.M., I.A.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** The results were obtained using the equipment of the Center for Collective Use of Scientific Equipment of TSU named after G.R. Derzhavin. This research was funded by the Russian Scientific Foundation grant 20-19-00602 (measurement of thermal properties), and grant 21-14-00233 (macro- and micromechanical testing), and by the Ministry of Science and Higher Education of the Russian Federation, the contract 075-15-2021-709, unique identifier of the project RF-2296.61321X0037 (equipment maintenance).

**Data Availability Statement:** All the data is available within the manuscript.

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

#### **References**


## *Article* **Wear of Functionally Graded Coatings under Frictional Heating Conditions**

**Vladimir B. Zelentsov, Polina A. Lapina \* and Boris I. Mitrin**

Research and Education Center "Materials", Don State Technical University, 344000 Rostov-on-Don, Russia; vbzelen@gmail.com (V.B.Z.); boris.mitrin@gmail.com (B.I.M.)

**\*** Correspondence: polina\_azarova86@mail.ru

**Abstract:** Multilayered and functionally graded coatings are extensively used for protection against wear of the working surfaces of mechanisms and machines subjected to sliding contact. The paper considers the problem of wear of a strip made of a functionally graded material, taking into account the heating of the sliding contact from friction. Wear is modeled by a moving strip along the surface of a hard abrasive in the form of a half-plane. With the help of the integral Laplace transform with respect to time, the solutions are constructed as convolutions from the law of the introduction of an abrasive into the strip and the original in the form of a contour integral of the inverse Laplace transform. The study of the integrands of contour quadratures in the complex plane allowed determination of the regions of stable solutions to the problem. Unstable solutions of the problem lead to the concept of thermoelastic instability of the contact with friction and formed regions of unstable solutions. The solutions obtained made it possible to determine a formula for the coefficient of functionally graded inhomogeneity of the coating material and to study its effect on the occurrence of thermoelastic instability of the contact taking friction into account, as well as on its main characteristics: temperature, displacement, stress and wear of the functionally graded material of the coating. The effects of the abrasive speed, contact stresses and temperature on wear of the coating with the functionally graded inhomogeneity of the material by the depth were investigated.

**Keywords:** functionally graded material; thermoelasticity; sliding contact; wear; heating from friction; thermoelastic instability

#### **1. Introduction**

Design of multilayered materials increased rapidly in the 1980–1990s, thus enabling the introduction of such materials with properties variable by depth in a number of industries such as construction, transport and other fields. Overall, their necessity in such fields was due to the increasing requirements for strength, wear and tear of details and parts. Recently, methods and technologies for the manufacturing of protective coatings from functionally graded materials (FGM) to micro- and nanometer technological levels have been developed: the centrifugal method [1–7], the technique of pulsed-laser deposition [8,9], the technique of magnetron sputtering [9], the plasma-spray technique [10,11], electrophoretic deposition and anodizing [12,13], micro-arc oxidation [14] and others [15].

FGM coatings are widely used to protect friction surfaces from wear. Generally, the tribological properties of FGM coatings are investigated by experimental methods, such as disc–pad [1,5], ball–disc [3,9,11], pin-on-disc [4,13,15], ring–pad [8], pin–plate [10] and other tests [7,14]. The coating wear as a result of testing is determined by examining the wear track with a microscope or a profilometer [1–3,9,10,13] or by weighing the sample before and after the test [4–8,11,15]. In this case, tests are carried out at one value of the experimental load or for a limited set (3–5 options).

The need to optimize the design of FGM coatings, to predict wear of the working surfaces of mechanisms and machines, to diagnose and prevent abnormal situations ne-

**Citation:** Zelentsov, V.B.; Lapina, P.A.; Mitrin, B.I. Wear of Functionally Graded Coatings under Frictional Heating Conditions. *Nanomaterials* **2022**, *12*, 142. https://doi.org/ 10.3390/nano12010142

Academic Editor: Mohammad Malikan

Received: 13 December 2021 Accepted: 28 December 2021 Published: 31 December 2021

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

**Copyright:** © 2021 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/).

cessitates mathematical modeling of the wear process. In modern microelectronics, semiconductor materials are used in the form of thin plates with a thickness not exceeding 20–30 μm down to 100 nm, and the presence of a functionally graded inhomogeneity of the material (FGIM) by depth. The process of thinning such plates in practice is carried out by grinding, polishing and washing. Mathematical modeling of the grinding process of the workpiece material is associated with the elucidation of the degree of influence of the properties of the FGIM of the workpiece material on the grinding process; on the nature of the workpiece heating; on the possibility of the occurrence of thermoelastic instability of the sliding contact; on acceleration/deceleration the grinding process, etc.

A lot of research is devoted to mathematical modeling of the contact of bodies with FGM coatings [16–26]. However, insufficient attention has been paid to the problem of modeling of wear or grinding of FGM coatings. There has been little research conducted in this direction [27].

Today, when solving problems of material wear, Archard's relations are most often used [28]. Dow and Burton were among the first investigators who used the Archard relation in the study of wear under conditions of heat release from friction [29], where the conditions for the occurrence of thermoelastic instability of blade sliding along the surface of a half-space were investigated using the method of small perturbations. Aleksandrov and Annakulova considered contact problems taking into account heat release from friction and wear of the coating [30], as well as the problem of mutual wear of coatings [31]. An attempt to develop a thermodynamic model to describe thermomechanical phenomena at a contact, taking into account friction and wear, was undertaken in [32,33]. Beginning with [34,35] a new direction of the development of the model of contact of two elastic bodies considering friction, wear and heat release, based on the principle of virtual energy and the basic laws of thermodynamics, has emerged. The finite element implementation of such a model in a two-dimensional statement was conducted in [36]. In [37–41] the integral Laplace transform was used with the solution in the form of functional series over the poles of the integrands of the contour quadratures of the inverse Laplace transform to solve contact thermoelasticity problems regarding wear. The solution method allows one to establish the parameters of the boundary of thermoelastic instability of a contact with a friction to study the properties of the obtained solutions. In [42], using the method of integral transforms, the contact problem of sliding an elastic coating over the surface of another coating with friction, wear and heat release from friction was reduced to solving a differential equation, and the conditions of thermoelastic stability of such a system were considered. Quasi-static and dynamic uncoupled contact problems of thermoelasticity regarding friction and wear of a rod were considered in [43]. The conditions for the thermoelastic instability emergence during mutual wear of surfaces made of different materials were considered in [44,45]. Due to the large number of parameters in the problems of thermofrictional contact and wear, one-dimensional quasi-static problems were often considered. The connectivity of the deformation fields and temperature in the listed works was neglected and the problems of uncoupled thermoelasticity were considered. The coupled problem of thermoelasticity on wear of a coating taking into account frictional heat release in a quasi-static formulation was considered in [46].

In the present paper, we consider a problem of wear of an elastic coating in the form of a strip made of FGM by a hard abrasive in the form of a half-plane, which slides along the strip at a constant speed, heating it due to friction. The influence of the FGIM of the strip on the process of wear (grinding), heating from friction, the conditions for the occurrence of thermoelastic instability of the sliding contact are investigated.

#### **2. Problem Statement**

To investigate the effect of the FGIM of the coating on its wear we consider a quasistatic contact problem of sliding of a rigid heat-insulated abrasive, half-plane *I* (*h* ≤ *x* < ∞), sliding with constant velocity *V*, over the upper surface (*x* = *h*) of the elastic thermally conductive coating with thickness *h*(0 ≤ *x* ≤ *h*). The lower surface of the coating is perfectly

bonded to a rigid substrate, half-plane *II* (−∞ < *x* < 0). The coating shear modulus μ(*x*) varies with depth 0 ≤ *x* ≤ *h*. During abrasive (half-plane *I*) sliding, the coating wear takes place, which can also be thought of as the abrasive grinding of the coating surface. The frictional heat generated at the contact interface flows into the coating. From the initial time moment, the abrasive (half-plane *I*) slides along the *y* axis and deforms the upper (*x* = *h*) surface of the coating in the negative direction of the *x* axis according to the indentation law Δ(*t*). Until the initial time moment, the coating was at rest, and its temperature was equal to zero. The scheme of the contact problem is shown in Figure 1.

**Figure 1.** Scheme of the contact problem.

In the described problem formulation, the temperature, stresses and displacements distributions in the coating do not depend on the horizontal coordinate *y* and are functions only of the vertical coordinate *x* and time *t* [37–41]. In the case of a quasi-static formulation without body forces, the stress state of the coating is described by the differential equations of equilibrium

$$\frac{\partial \sigma\_{\text{xx}}}{\partial \mathbf{x}} = 0, \; \frac{\partial \sigma\_{\text{xy}}}{\partial \mathbf{x}} = 0, \; \; 0 \le \mathbf{x} \le h, \; \; t > 0 \tag{1}$$

where σ*xx*(*x*, *t*), σ*xy*(*x*, *t*) are the normal and shear components of the stresses in the coating. The heat equation has the form

$$\frac{\partial^2 T}{\partial \mathbf{x}^2} - \frac{1}{\kappa} \frac{\partial T}{\partial t} = 0 \quad 0 \le \mathbf{x} \le h, \ t > 0 \tag{2}$$

where *T*(*x*, *t*) is the coating temperature, κ is the thermal diffusivity.

The Duhamell–Neumann law takes place [47]

$$
\sigma\_{\rm XY} = \frac{2(1-\mathbf{v})}{1-2\mathbf{v}} \mu(\mathbf{x}) \left( \frac{\partial \mu}{\partial \mathbf{x}} - \frac{1+\mathbf{v}}{1-\mathbf{v}} \alpha T \right), \quad \sigma\_{\rm XY} = \mu(\mathbf{x}) \frac{\partial w}{\partial \mathbf{x}'} \tag{3}
$$

where *u*(*x*, *t*), *w*(*x*, *t*) are the vertical and horizontal components of the displacements in the coating, μ(*x*), ν, α are the shear modulus, Poisson's ratio, coefficient of linear heat expansion of the coating material.

The differential equations of linear uncoupled elasticity are represented by the system of equilibrium Equation (1) and heat Equation (2), which together describe thermoelastic state of the coating.

Boundary conditions for Equation (1) are as follows (*t* > 0):

$$
\mu(h, t) = -\Delta(t) + \mu\_w(t) \tag{4}
$$

$$
\sigma\_{\rm xy}(h, t) = -f \sigma\_{\rm xx}(h, t) \tag{5}
$$

$$u(0,t) = 0\tag{6}$$

$$w(0,t) = 0\tag{7}$$

where *f* is the coefficient of friction, *uw*(*t*) is the abrasive (half-plane *I*) displacement due to the coating wear. To find the solution, we use the abrasive wear model [23] in integral form

$$u\_{\rm av}(t) = -fVK^\* \int\_0^t \sigma\_{\rm xx}(h, \tau)d\tau \qquad t > 0\tag{8}$$

where σ*xx*(*h*, *t*) is the compressive normal stress on the contact interface, *K*<sup>∗</sup> is the proportionality coefficient between the work of friction forces and the volume of removed material. The boundary conditions for the heat Equation (2) are (*t* > 0)

$$K \frac{\partial T(h, t)}{\partial x} = Q(t) \tag{9}$$

$$K\frac{\partial T(0,t)}{\partial \boldsymbol{x}} = k(T(0,t) - T(0,0))\tag{10}$$

where *K* is the thermal conductivity of the coating material, *k* is the heat transfer coefficient through the coating–substrate interface, *Q*(*t*) = *f V*(−σ*xx*(*h*, *t*)) is the amount of frictional heat originatingfrom the contact interface [48]. From (9) it follows that all the heat at the contact is due to friction.

Initial conditions for displacements and temperature are equal to zero:

$$w(\mathbf{x},0) = w(\mathbf{x},0) = T(\mathbf{x},0) = 0, \quad 0 \le \mathbf{x} \le h,\tag{11}$$

Thus, the solution of the considered quasi-static thermoelastic contact problem for the elastic FGM coating wear (or grinding) by a certain depth with a hard abrasive in the form of a half-plane, taking into account heating from friction, is derived through the solution of the initial boundary value problem, including the system of the differential equations of elasticity (1) and heat conduction (2) with boundary (4)–(10) and initial (11) conditions. Note that vertical displacements *u*(*x*, *t*), normal stresses σ*xx*(*x*, *t*) and temperature *T*(*x*, *t*) are found separately from the horizontal displacements *w*(*x*, *t*). The horizontal displacements *w*(*x*, *t*) are determined from (1), (5), (7) knowing the normal stresses' distribution σ*xx*(*h*, *t*).

#### **3. Exact Solution for Arbitrary** μ**(x)**

We proceed to the solution of the formulated problem introducing the Laplace integral transform [49]

$$T^{L}(\mathbf{x}, p) = \int\_{0}^{\infty} T(\mathbf{x}, t) e^{-pt} dt,\ T(\mathbf{x}, t) = \frac{1}{2\pi\mathrm{i}} \int\_{-i\infty + c}^{i\infty + c} T^{L}(\mathbf{x}, p) e^{pt} dp \quad \mathrm{Rep} < c, \ c > 0 \tag{12}$$

To obtain the temperature distribution *T*(*x*, *t*) in the coating, we apply the Laplace transform to heat Equation (2). As a result, the Laplace image *TL*(*x*, *p*) of the coating temperature has the form

$$T^L(\mathbf{x}, p) = A\_1 \text{sh} \sqrt{\frac{p}{\kappa}} \mathbf{x} + A\_2 \text{ch} \sqrt{\frac{p}{\kappa}} \mathbf{x} \tag{13}$$

where *A*<sup>1</sup> and *A*<sup>2</sup> are arbitrary constants.

The vertical displacements *u*(*x*, *t*) are determined from the first equilibrium Equation (1) and taking into account for the first relation in (3). We use the Laplace transform (12) and obtain the Laplace image *uL*(*x*, *p*) of the vertical displacements

$$u^L(\mathbf{x}, p) = \frac{1+v}{1-v} \mathfrak{a} \frac{1}{\sqrt{\frac{p}{\kappa}}} \left( A\_1 \text{ch} \sqrt{\frac{p}{\kappa}} \mathbf{x} + A\_2 \text{sh} \sqrt{\frac{p}{\kappa}} \mathbf{x} \right) - A\_3 B(\mathbf{x}) + A\_4 \tag{14}$$

where *A*1, *A*<sup>2</sup> are from (13) and *A*3, *A*<sup>4</sup> are additional arbitrary constants. The function *B*(*x*) is defined through μ(*x*) as follows

$$B(\mathbf{x}) = \int\_0^\mathbf{x} \frac{d\mathbf{\tilde{y}}}{\mu(\mathbf{\tilde{x}})}, \quad 0 \le \mathbf{x} \le h \tag{15}$$

where μ(*x*) is a continuous function and does not vanish μ(*x*) = 0 for any *x* ∈ [0, *h*].

To determine constants *Ak*, *k* = 1 − 4 in (13), (14) we apply the Laplace transform to boundary conditions (4), (6), (9) and (10), which results in

$$u^L(\hbar, p) = -\Delta^L(p) + u\_w^L(p) \tag{16}$$

$$u^L(0, p) = 0\tag{17}$$

$$K\frac{dT^L(h,p)}{dx} = -fV\sigma\_{\text{xx}}^L(h,p)\tag{18}$$

$$K\frac{dT^L(0,p)}{dx} = kT^L(0,p)\tag{19}$$

where

$$u\_w^L(p) = -fVK^\* \frac{\sigma\_{\text{xx}}^L(h, p)}{p} \tag{20}$$

$$\sigma\_{xx}^{L}(\mathbf{x}, p) = \frac{2(1 - v)}{1 - 2v} \mu\_1 \left( \frac{du^L(\mathbf{x}, p)}{d\mathbf{x}} - \frac{1 + v}{1 - v} \alpha T^L(\mathbf{x}, p) \right) \tag{21}$$

where μ<sup>1</sup> = μ(*h*) is the value of the shear modulus at the upper coating boundary and Δ*L*(*p*) is the Laplace image of the function Δ(*t*), representing half-plane *I* displacement towards the coating.

Substituting (13), (14), (20), (21) into boundary conditions (16)–(19), we obtain the linear system for the determination of the constants *Ak*, *k* = 1 − 4. Solving this system, we obtain the Laplace images of the temperature, displacement and stress distribution in the coating as follows

$$T^L(\mathbf{x}, p) = \frac{1 - \mathbf{v}}{1 + \mathbf{v}} \frac{\hat{\mathcal{V}}}{\alpha h} \Delta^L(p) \frac{h \mathcal{B}'(h)}{B(h)} \frac{\mathcal{N}\_\Gamma(\mathbf{x}, z)}{\mathcal{R}(z)} \tag{22}$$

$$N\_T(\mathbf{x}, z) = \sqrt{z} \left( \text{Bish}\sqrt{z}\frac{\mathbf{x}}{\hbar} + \sqrt{z}\text{ch}\sqrt{z}\frac{\mathbf{x}}{\hbar} \right) \tag{23}$$

$$u^L(\mathbf{x}, p) = -\Delta^L(p) \cdot \frac{N\_\mu^0(\mathbf{x}, z)}{R(z)}\tag{24}$$

$$N\_{\rm ll}^{0}(\mathbf{x}, z) = zr(h, z)\frac{B(\mathbf{x})}{B(h)} - \dot{\mathcal{V}}\frac{hB'(h)}{B(h)}(r(\mathbf{x}, z) - \text{Bi})\tag{25}$$

$$\sigma\_{\rm xx}^{L}(\mathbf{x}, p) = \frac{2(1 - v)}{1 - 2v} \mu(\mathbf{x}) \Delta^{L}(p) \frac{hB'(\mathbf{x})}{B(h)} \cdot \frac{N\_{\sigma}^{0}(z)}{R(z)} \tag{26}$$

$$N\_{\sigma}^{0}(z) = zr(\hbar, z) \tag{27}$$

$$R(z) = zr(h, z) - \mathcal{V}\eta((1 - k\_w)r(h, z) - \text{Bi})\tag{28}$$

$$r(x, z) = \text{Bi ch}\sqrt{z}\frac{\mathcal{X}}{\hbar} + \sqrt{z}\text{sh}\sqrt{z}\frac{\mathcal{X}}{\hbar} \tag{29}$$

$$z = \frac{p}{\kappa}h^2,\ \text{Bi} = \frac{kh}{K},\ \ k\_{\text{\textquotedblleft}v} = \frac{1-v}{1+v}\frac{\text{KK}^\*}{\alpha\kappa},\ \ \ \hat{V} = \frac{fV\alpha}{K}\frac{2\mu\_1(1+v)h}{1-2v}.$$

For the image *u<sup>L</sup> <sup>w</sup>*(*p*) of the coating wear, we obtain

$$
\mu\_w^L(p) = k\_w \hat{\mathcal{V}} \frac{hB'(h)}{B(h)} \Delta^L(p) \frac{r(h, z)}{R(z)}\tag{30}
$$

We apply the inverse Laplace transform to the Laplace images *TL*(*x*, *p*), *uL*(*x*, *p*), σ*L xx*(*x*, *p*) and obtain the problem solution in the convolution form (*t* > 0)

$$T(\mathbf{x},t) = \frac{1-v}{1+v} \frac{\mathcal{V}}{\alpha h} \cdot \frac{hB'(h)}{B(h)} \int\_0^t \Delta(\mathbf{r}) f\_T(\mathbf{x}, t-\mathbf{r}) d\mathbf{r} \tag{31}$$

$$f\_T(\mathbf{x}, t) = \frac{1}{2\pi i} \int \frac{N\_T(\mathbf{x}, z)}{t\_\mathbf{x} R(z)} e^{z\tilde{t}} dz \tag{32}$$

$$u(\mathbf{x},t) = -\int\_0^t \Delta(\mathbf{r}) f\_{\mu}^0(\mathbf{x}, t-\mathbf{r}) d\mathbf{r} \tag{33}$$

$$f\_{\boldsymbol{\mu}}^{0}(\boldsymbol{x},t) = \frac{1}{2\pi i} \int\_{\Gamma} \frac{N\_{\boldsymbol{\mu}}^{0}(\boldsymbol{x},z)}{t\_{\boldsymbol{\kappa}}R(z)} e^{z\tilde{t}} dz \tag{34}$$

$$\sigma\_{\rm xx}(\mathbf{x},t) = -\frac{2(1-v)}{(1-2v)B(h)} \int\_0^t \Delta(\tau) f\_\sigma^0(\mathbf{x}, t-\tau)d\tau \tag{35}$$

$$f\_{\sigma}^{0}(\mathbf{x},t) = \frac{1}{2\pi i} \int \frac{N\_{\sigma}(\mathbf{x},z)}{t\_{\kappa}R(z)} e^{z\tilde{t}}dz\tag{36}$$

where *<sup>t</sup>* <sup>=</sup> *<sup>t</sup> <sup>t</sup>*<sup>κ</sup> , *<sup>t</sup>*<sup>κ</sup> <sup>=</sup> *<sup>h</sup>*<sup>2</sup> κ .

To obtain (35) from (26) we take into account that μ(*x*)*B* (*x*) ≡ 1. Coating surface wear *u*w(*t*) is determined by inverting *u<sup>L</sup>* <sup>w</sup>(*p*) in (30)

$$
\mu\_{\rm W}(t) = k\_{\rm W} \mathcal{V} \frac{hB'(h)}{B(h)} \int\_0^t \Delta(\tau) f\_{\rm w}(\mathbf{x}, t - \tau) d\tau \tag{37}
$$

$$f\_{\mathcal{W}}(\mathbf{x},t) = \frac{1}{2\pi i} \int\_{\Gamma} \frac{r(\mathbf{h},z)}{t\_{\mathbf{k}}\mathcal{R}(z)} e^{z\tilde{t}} dz \tag{38}$$

To find out the existence conditions of integrals in (32), (34), (36), (38) we analyze integrands for 0 ≤ *x* ≤ *h* for large values of the integration variable *z* (arg*z* = π/2):

$$\begin{aligned} N\_T(\mathbf{x}, z)R^{-1}(z) &= O\left(z^{-1/2}\right)|z| \to \infty \\ N\_\mu^0(\mathbf{x}, z)R^{-1}(z) &= \frac{B(\mathbf{x})}{B(h)} + O\left(z^{-1/2}\right)|z| \to \infty \\ N\_\sigma^0(\mathbf{x}, z)R^{-1}(z) &= 1 + O\left(z^{-1/2}\right)|z| \to \infty \\ r(h, z)R^{-1}(z) &= O\left(z^{-1}\right)|z| \to \infty \end{aligned} \tag{39}$$

From (39) it follows that the integrands in (34), (36) do not decay at infinity (at |*z*|→ ∞), and the corresponding integrals are divergent and understood in the generalized sense [50]. After regularization of integrals (34), (36) and separation of the generalized part, expressions for the displacements *u*(*x*, *t*) and stresses σ*xx*(*x*, *t*) can be written as follows (*t* > 0)

$$u(\mathbf{x},t) = -\frac{B(\mathbf{x})}{B(h)}\Delta(t) - \int\_0^t \Delta(\tau)f\_{\mu}(\mathbf{x}, t-\tau)d\tau \qquad 0 \le \mathbf{x} \le h, \ t > 0 \tag{40}$$

$$f\_{\mathfrak{u}}(\mathbf{x},t) = \frac{1}{2\pi i} \int \frac{N\_{\mathfrak{u}}(\mathbf{x},z)}{t\_{\mathfrak{k}}R(z)} e^{z\tilde{t}} dz \tag{41}$$

$$N\_{\rm ll}(\mathbf{x}, z) = N\_{\rm ll}^{0}(\mathbf{x}, z) - \frac{B(\mathbf{x})}{B(h)} \mathcal{R}(z) \tag{42}$$

$$\sigma\_{\rm xx}(\mathbf{x},t) = -\frac{2(1-v)}{(1-2v)B(h)} \left( \Delta(t) - \int\_0^t \Delta(\tau) f\_\sigma(\mathbf{x}, t-\tau) d\tau \right) \\ 0 \le \mathbf{x} \le h, t > 0 \tag{43}$$

$$f\_{\sigma}(\mathbf{x},t) = \frac{1}{2\pi i} \int \frac{N\_{\sigma}(\mathbf{x},z)}{t\_{\mathbf{x}}R(z)} e^{z\tilde{t}}dz\tag{44}$$

$$N\_{\sigma}(\mathbf{x}, z) = N\_{\sigma}^{0}(\mathbf{x}, z) - R(z) \tag{45}$$

where *<sup>t</sup>* and *<sup>t</sup>*<sup>κ</sup> are previously identified, and the contour of integration Γ = {*z* : −*i*∞ + *dt*κ, +*i*∞ + *dt*κ} is a straight line in the complex plane of integration variable *z*, which is parallel to the imaginary axis and a distance from equal to the value *dt*κ. The value of *d* is chosen so that the contour of integration passes to the right of all of the isolated singular points of integrands.

Integrands in Equations (40) and (43) at 0 ≤ *x* ≤ *h* are meromorphic functions and decay at large *z* (−π < arg*z* < π)

$$\begin{aligned} N\_{\rm ll}(\mathbf{x}, z) \mathbf{R}^{-1}(z) &= O\left(z^{-1/2}\right)|z| \to \infty \\ N\_{\sigma}(\mathbf{x}, z) \mathbf{R}^{-1}(z) &= O\left(z^{-1/2}\right)|z| \to \infty \end{aligned} \tag{46}$$

These properties of integrands allow us to apply the methods of complex analysis for calculation of the integrals and stability analysis. In quadratures (32), (38), a regularization is carried out to obtain integrands decreasing at infinity, after which they are investigated by the same methods as in (40), (43).

#### **4. Poles of the Integrands**

To investigate the stability of the solutions obtained in the previous subsection, we need to study effect of the problem parameters on the integrand poles in integrals (32), (38), (41), (44). The poles of the integrands are zeros of the equation

$$R(z) = zr(h, z) - \eta \hat{V}((1 - k\_w)r(h, z) - \text{Bi}) = 0, \ |\text{arg}z| < \pi, \ |z| < \infty \tag{47}$$

in the complex plane *z* = *ξ* + *i*η, where *R*(*z*) and *r*(*x*, *z*) at *x* = *h* are from (28) and (29), η = *hB* (*h*)/*B*(*h*) is the coefficient of the FGIM strip. Equation (47) contains four dimensionless parameters of the problem (*kw*, *V*ˆ , Bi, η), which itself contains dimensional parameters, described after (29). In this case, the dimensionless parameter η, the numerical value of which characterizes the FGM coating, consists of *B* (*h*) = μ−1(*h*) and *B*(*h*) from (15). From the mean value theorem, it follows for *B*(*h*) that there is a point *c* ∈ [0, *h*] satisfying the equation

$$B(h) = \mu^{-1}(c)h, \quad c \in [0, h] \tag{48}$$

Then, the parameter η can be represented as the averaged value of the shear modulus in the coating μ(*c*)*c* ∈ [0, *h*] divided by the value of μ(*h*) at the upper boundary of the coating

$$\eta = \frac{\mu(c)}{\mu(h)'}, \quad c \in [0, h] \tag{49}$$

which is its mechanical meaning.

Locating zeros *ζ<sup>k</sup>* of Equation (47) *R*(*ζk*) = 0, *k* = 0, 1, 2, ..., their movements in the complex plane *ζ* = *ξ* + *i*η depending on the change in the dimensionless parameters of the problem is the main goal of solving Equation (47).

Equation (47) is solved using numerical methods and complex analysis methods [51]. Analysis of *<sup>R</sup>*(*z*) (47) zeros is similar to [52–54] and is performed for *<sup>V</sup>*<sup>ˆ</sup> <sup>∈</sup> [0, <sup>∞</sup>) at fixed values of dimensionless parameters *k*w, η, Bi. Assuming *V*ˆ = 0, we find the following equation

$$zr(\mathbf{h}, z) = z(\mathbf{B} \text{ch}\sqrt{z} + \sqrt{z} \text{sh}\sqrt{z}) = 0\tag{50}$$

to determine the initial estimates *ζ*<sup>0</sup> *<sup>k</sup>* <sup>=</sup> *<sup>ζ</sup>k*(0), *<sup>k</sup>* <sup>=</sup> 0, 1, 2, ... of zeros *<sup>ζ</sup>k*(*V*ˆ), *<sup>k</sup>* <sup>=</sup> 0, 1, 2, ... of Equation (47). In general case, Equation (50) does not allow convenient analytic solutions. However, at Bi = 0 Equation (50) possesses analytic solution, and the initial estimations *ζ*<sup>0</sup> *k* , *k* = 0, 1, 2, . . . of zeros *ζk*, *k* = 0, 1, 2, . . . of Equation (47) are found from

$$\zeta\_k^0 = -\left(\pi k\right)^2, \quad k = 0, 1, 2, \dots \tag{51}$$

And do not depend on *kw*, *T*ˆ. At Bi = ∞, from (50) we obtain another initial estimations *ζ*0 *<sup>k</sup>* , *k* = 0, 1, 2, . . . as follows

$$\mathcal{L}\_k^0 = -\pi^2 (k + 1/2)^2, \quad k = 0, 1, 2, \dots \tag{52}$$

Asymptotic estimation of *ζ*<sup>0</sup> *<sup>k</sup>* for large numbers of *k* has the form of (51). From Formulas (51) and (52) it follows that the initial estimations *ζ*<sup>0</sup> *<sup>k</sup>* , *k* = 0, 1, 2, ... of function *R*(*z*) of zeros from (47) are located on the negative part of the real axis or at the origin. As *V*ˆ changes from 0 to ∞ at fixed *kw*, Bi, η the first two poles *ζ*<sup>0</sup> and *ζ*<sup>1</sup> can be located: I—on the real axis Re(*ζ*0, *ζ*1) < 0, Im(*ζ*0, *ζ*1) = 0 at 0 < *V*ˆ < *V*ˆ I; II—in the vertical strip *<sup>l</sup>*<sup>−</sup> <sup>&</sup>lt; Re(*ζ*0, *<sup>ζ</sup>*1) <sup>&</sup>lt; 0, <sup>|</sup>Im(*ζ*0, *<sup>ζ</sup>*1)<sup>|</sup> <sup>&</sup>lt; <sup>∞</sup> at *<sup>V</sup>*<sup>ˆ</sup> <sup>I</sup> < *V*ˆ < *V*ˆ II; III—in the vertical strip <sup>0</sup> <sup>&</sup>lt; Re(*ζ*0, *<sup>ζ</sup>*1) <sup>&</sup>lt; *<sup>l</sup>*+, <sup>|</sup>Im(*ζ*0, *<sup>ζ</sup>*1)<sup>|</sup> <sup>&</sup>lt; <sup>∞</sup> at *<sup>V</sup>*<sup>ˆ</sup> II < *V*ˆ < *V*ˆ III; IV—on the positive part of the real axis Re(*ζ*0, *ζ*1) > 0, Im(*ζ*0, *ζ*1) = 0 at *V*ˆ III < *V*ˆ < ∞.

The numbers I, II denote regions (we will call them domains of stability), where Re(*ζ*0, *ζ*1) < 0 at 0 < *V*ˆ < *V*ˆ II, while III, IV denote the regions (we will call them domains of instability), where Re(*ζ*0, *ζ*1) > 0 at *V*ˆ II < *V*ˆ < ∞. Figure 2 shows some examples of the poles *ζ*<sup>0</sup> *V*ˆ and *ζ*<sup>1</sup> *V*ˆ trajectories when *V*ˆ changes from 0 to ∞ at fixed Bi = 100 for three different values of η = 0.5, 1.0, 2.0 (green, red, blue colors) and following values of *kw*= 0.5 (curve *1*), 0.9 (*2*), 1.0 (*3*), 1.1 (*4*), 1.35 (*5*), 5.0 (*6*). Solid diamonds indicate locations of *ζ*<sup>0</sup> *V*ˆ and *ζ*<sup>1</sup> *V*ˆ at *<sup>V</sup>*<sup>ˆ</sup> = 0, while empty ones correspond to *<sup>V</sup>*<sup>ˆ</sup> <sup>→</sup> <sup>∞</sup>. The point of the trajectory marked as a crossed-out square is a point, after crossing which, with increasing *V*ˆ , the real poles *ζ*<sup>0</sup> *V*ˆ and *ζ*<sup>1</sup> *V*ˆ become a pair of complex conjugated poles, and vice versa. Note that even small change in the parameter *kw*, which is proportional to the ratio of the wear factor *K*<sup>∗</sup> and the thermal expansion coefficient α, leads to significant variation of *ζ*<sup>0</sup> *V*ˆ and *ζ*<sup>1</sup> *V*ˆ trajectories and, to a lesser extent, trajectories of other *ζ<sup>k</sup> V*ˆ , *k* = 2, 3, 4, .... If wear is extensive (*kw* > 1), both *ζ*0, *ζ*<sup>1</sup> and all others *ζk*, *k* = 2, 3, 4, ··· are located in regions I, II (curves *4*–*6* on Figure 2). When thermal expansion of the coating material prevails over wear (0 < *kw* < 1), the poles *ζ*<sup>0</sup> and *ζ*<sup>1</sup> move to the right complex half-plane to regions III, IV (curves *1*–*3* on Figure 2).

One significant property of the poles from regions II, III is noted; they are complex conjugated, i.e., *ζ*<sup>1</sup> = *ζ*<sup>0</sup> and *ζ*<sup>0</sup> = *ζ*1.

Wear causes a dramatic effect on the poles' behavior. In quasi-static sliding contact problems with frictional heating but without wear, the poles *ζ*<sup>0</sup> and *ζ*<sup>1</sup> remain on the real axis for any *<sup>V</sup>*<sup>ˆ</sup> <sup>∈</sup> [0, <sup>∞</sup>) [52,53]. In corresponding problems accounting for wear, the poles *ζ*<sup>0</sup> and *ζ*<sup>1</sup> can have a non-zero imaginary part, and in this case become complex conjugates *ζ*<sup>1</sup> = *ζ*0, *ζ*<sup>0</sup> = *ζ*1.

**Figure 2.** Movement of zeros of *R*(*z*) (47) in the complex plane when *V*ˆ changes from 0 to ∞.

#### **5. Formulas for Exact Solution**

When the poles *ζk*, *k* = 0, 1, 2, ... are known, to calculate the integrals in (32), (41), (44) we calculate the sum of the residues at the points *z* = *ζk*. For all simple poles we obtain

$$\frac{1}{2\pi i} \int \frac{N\_d(\varkappa, z)}{t\_\varkappa R(z)} e^{z\tilde{l}} dz = \sum\_{k=0}^\infty \mathbb{C}\_d(\varkappa, \zeta\_k) e^{\zeta\_k \tilde{l}} \,\tag{53}$$

$$\mathbf{C}\_{a}(\mathbf{x},z) = \frac{\mathbf{N}\_{a}(\mathbf{x},z)}{t\_{\mathbf{x}}R'(z)}\tag{54}$$

By replacing symbolic index *a* in (53), (54) with *T*, *u* or σ, Equation (53) can be used to calculate integral in (32), (41), (44), respectively. If *ζ<sup>k</sup>* and *ζk*<sup>+</sup>1, *k* = 0, 1, 2, ... are complex conjugates (*ζk*+<sup>1</sup> = *ζk*), then

$$\mathbb{C}\_{a}(\chi, z)e^{z\tilde{t}} = 2\text{Re}\frac{\mathrm{N}\_{a}(\chi, z)}{t\_{\mathsf{K}}R'(z)}e^{z\tilde{t}}\tag{55}$$

and summation in (53) can be entered by even numbers *k* = 2*n*, *n* = 0, 1, 2, ... for complex conjugate *ζk*, *k* = 0, 1, 2, . . .. Using (54) for (32), (41), (44) we obtain

$$f\_d(\mathbf{x}, t) = \frac{1}{2\pi i} \int\_{\Gamma} \frac{\mathcal{N}\_d(\mathbf{x}, z)}{t\_\mathbf{x} \mathcal{R}(z)} e^{z\tilde{I}} dz = \sum\_{k=0}^\infty \mathcal{C}\_d(\mathbf{x}, \zeta\_k) e^{\zeta\_k \tilde{I}} a = T\_\prime \, u\_\prime \, \mathbf{c} \tag{56}$$

and write out the problem solution in form of series

$$T(\mathbf{x},t) = \frac{1-\upsilon}{1+\upsilon} \frac{\hat{V}}{\propto} \frac{B'(h)}{B(h)} \sum\_{k=0}^{\infty} C\_T(\mathbf{x}, \zeta\_k) D(\zeta\_k, t) \quad 0 \le \mathbf{x} \le h, \ t > 0 \tag{57}$$

$$u(\mathbf{x},t) = -\frac{B(\mathbf{x})}{B(h)}\Delta(t) + \sum\_{k=0}^{\infty} \mathcal{C}\_{\mathbf{u}}(\mathbf{x}, \mathcal{J}\_{k}) D(\mathcal{J}\_{k}, t) \quad 0 \le \mathbf{x} \le h, \ t > 0 \tag{58}$$

$$\sigma\_{\text{xx}}(\mathbf{x},t) = -\frac{\Im(1-v)}{(1-2v)B(h)} \left( \Delta(t) - \sum\_{k=0}^{\infty} \mathbb{C}\_{\sigma}(\mathbf{x}, \mathbb{zeta}\_k) D(\mathbb{f}\_{k\prime}, t) \right) \quad 0 \le \mathbf{x} \le h, \ t > 0 \tag{59}$$

where *Ca*(*x*, *z*) is calculated either by (54) or (55), and *D*(*z*, *t*) is calculated as follows

$$D(z,t) = \int\_0^t \Delta(\tau) \exp(z(t-\tau)/t\_\kappa) d\tau \qquad t>0\tag{60}$$

Calculating *f* <sup>0</sup> *<sup>w</sup>*(*t*) in (38) by relation

$$f\_w^0(t) = \frac{1}{2\pi i} \int\_{\Gamma} \frac{r(\hbar, z)}{t\_\kappa \mathcal{R}(z)} e^{z\tilde{t}} dz = \sum\_{k=0}^\infty B\_w(\zeta\_k) e^{\zeta\_k \tilde{t}} \quad \mathcal{C}\_w = \frac{r(\hbar, z)}{t\_\kappa \mathcal{R}'(z)}\tag{61}$$

and substituting it into (37), we obtain a formula for calculation of the coating wear

$$u\_W(t) = k\_W \hat{\mathcal{V}} \frac{hB'(h)}{B(h)} \sum\_{k=0}^{\infty} \mathbb{C}\_{\mathbf{w}}(\zeta\_k) D(\zeta\_{k'} t), \quad t > 0 \tag{62}$$

The horizontal displacements *w*(*x*, *t*) are determined from (1), (5), (7) and after integrating (1) take the form

$$w(\mathbf{x}, t) = -fB(\mathbf{x})\sigma\_{\mathbf{x}\mathbf{x}}(h, t) \quad 0 \le \mathbf{x} \le h, \ t > 0 \tag{63}$$

#### **6. Domains of Stability and Instability**

Analysis of Formulas (57)–(59) for *T*(*x*, *t*), *u*(*x*, *t*), σ*xx*(*x*, *t*) shows that, in case of Re(*ζk*) < 0, *k* = 0, 1, 2, ... the solution is stable and tends towards a stationary state with increasing *t*. However, if at least one of *ζk*, *k* = 0, 1, 2, ... has Re(*ζk*) > 0, then the magnitude of the solution grows indefinitely when *t* → ∞ and oscillates with the frequency Im(*ζk*) = 0, which indicates instability of the problem. If the indentation law Δ(*t*) is a bounded function

$$m < \Delta(t) < M \quad \text{ } m, M > 0, \quad 0 < t < \infty$$

then the following estimate for the integral in (60) takes place

$$|D(\mathcal{J}\_{k'}t)| \ge m \left| \frac{1 - e^{\mathcal{J}\_k \tilde{t}}}{\mathcal{J}\_k} \right| \text{ at } \text{Re}(\mathcal{J}\_k) > 0 \quad k = 0, 1, 2, \dots$$

Trajectories of the poles *ζ<sup>k</sup> V*ˆ , *<sup>k</sup>* <sup>=</sup> 0, 1, 2, ...*V*<sup>ˆ</sup> <sup>∈</sup> [0, <sup>∞</sup>), lying in the left complex halfplane (Re(*ζk*) < 0), correspond to the stable solution, and therefore we refer to domains I, II as the domains of stable solutions. Domains III, IV, lying in the right complex half-plane (Re(*ζk*) > 0, *k* = 0, 1), can be referred to as the domains of unstable solutions, because in domain III lim *<sup>t</sup>*→∞*T*(*h*, *<sup>t</sup>*) and lim *t*→∞ σ*xx*(*h*, *t*) do not exist (because Im(*ζk*) = 0, *k* = 0, 1), and in domain IV we have lim *<sup>t</sup>*→∞*T*(*h*, *<sup>t</sup>*) = lim *t*→∞ σ(*h*, *t*) = ∞ (because Im(*ζk*) = 0, *k* = 0, 1).

Figure 3 presents domains of stability (I, II) and instability (III, IV) in plane (*V*ˆ , *kw*) and boundaries between them for Bi = 100 and different values of η = 0.1; 0.25; 0.5; 1.0; 5.0; 10.0. These plots show that, at fixed *k*w, *V*ˆ , η the parameter Bi affects the position of the boundaries of stability and instability. The point of intersection of boundaries I–IV, lying on axis *V*ˆ (Figure 3), shifts depending on Bi and has a coordinate *V*ˆ <sup>∗</sup> = 2Bi(2 + Bi) −1 η−<sup>1</sup> (*k*<sup>w</sup> = 0). At η = 1 this result coincides with [41]. The formula for *V*ˆ <sup>∗</sup> shows the effect of the parameter η (Figure 3) on the boundaries of domains I–IV.

**Figure 3.** Domains of stability (I, II) and instability (III, IV) of the problem solution at Bi = 100 for different values of η: (**a**) 0.1, (**b**) 0.25, (**c**) 0.5, (**d**) 1.0, (**e**) 5.0, (**f**) 10.0.

Note that the dimensionless parameter η, which characterizes inhomogeneity of the FGM of the coating, significantly affects the boundary between domains of stable (II) and unstable (III) solutions for any *k*w. These boundaries are presented in Figure 4 in detail. One can observe that with decrease in η, the stable solution region increases.

**Figure 4.** Boundaries between stability and instability domains of the problem solution in the (*V*ˆ , *k*w) plane at Bi = 100 for different values of η: *1*—0.25, *2*—0.5, *3*—1.0, *4*—2.0, *5*—5.0.

#### **7. Features of Wear of FGM Coating**

In Section 5 we obtained the exact formulas for the main characteristics of the problem: temperature *T*(*x*, *t*) (57), displacements *u*(*x*, *t*) (58), stresses σ*xx*(*x*, *t*) (59) and coating wear *u*w(*t*) (62) on sliding contact. These functions depend on the shear modulus μ(*x*) variation by the coating depth 0 ≤ *x* ≤ *h*. To illustrate the solution, we assume that the shear modulus μ(*x*) varies according to a power (parabolic) law

$$
\mu(x) = \mu\_0 \left( a \left( \frac{x}{h} \right)^2 + b \frac{x}{h} + c \right) \tag{64}
$$

$$a = 2\left(\frac{\mu\_1}{\mu\_0} - 2\frac{\mu\_{1/2}}{\mu\_0} + 1\right), \\ b = -\left(\frac{\mu\_1}{\mu\_0} - 4\frac{\mu\_{1/2}}{\mu\_0} + 3\right), \\ c = 1$$

where μ<sup>0</sup> = μ(0), μ<sup>1</sup> = μ(*h*), μ1/2 = μ(*h*/2). This means that the shear modulus is equal to μ<sup>0</sup> at the coating–substrate interface (*x* = 0), to μ<sup>1</sup> at the contact interface (*x* = *h*), and to μ1/2 in the middle of the coating. If μ1/2 = (μ<sup>1</sup> + μ0)/2 then *a* = 0 and the law (64) becomes linear.

Calculating integral (15) of the function in (64), we obtain the expression for *B*(*x*)

$$B(\boldsymbol{x}) = \frac{h}{\mu\_1} \begin{cases} \frac{2\chi}{\sqrt{-D}} \text{arctg}\frac{\sqrt{-D}}{\Phi(\boldsymbol{x})} & D < 0\\ \frac{\chi}{\sqrt{D}} \ln \left| \frac{1 + \Theta\_-(\boldsymbol{x})}{1 - \Theta\_+(\boldsymbol{x})} \right| & D > 0\\ -4\frac{\epsilon}{b} \frac{\chi}{\Phi(\boldsymbol{x})} & D = 0 \end{cases} \tag{65}$$

where *<sup>D</sup>* <sup>=</sup> *<sup>b</sup>*<sup>2</sup> <sup>−</sup> <sup>4</sup>*ac*, *<sup>χ</sup>* <sup>=</sup> <sup>μ</sup><sup>1</sup> μ0 , θ(*x*) = 2*c* + *b <sup>x</sup> <sup>h</sup>* , <sup>θ</sup>±(*x*) = <sup>2</sup>*<sup>a</sup> b*± <sup>√</sup>*<sup>D</sup>* · *<sup>x</sup> h* .

Equation (65) allows one to determine other characteristics depending on μ(*x*)

$$B'(\mathbf{x}) = \frac{1}{\mu(\mathbf{x})} = \frac{1}{\mu\_0 \left(a\left(\frac{\mathbf{x}}{h}\right)^2 + b\frac{\mathbf{x}}{h} + c\right)}\tag{66}$$

$$B'(h) = \frac{1}{\mu\_0 \chi} \tag{67}$$

$$B(h) = \frac{h}{\mu\_1} \begin{cases} \frac{2\chi}{\sqrt{-D}} \text{arctg}\frac{\sqrt{-D}}{2c+b} & D < 0\\ \frac{\chi}{\sqrt{D}} \ln \left| \frac{1+b^1\_-}{1-b^1\_+} \right| & D > 0\\ -4\chi \frac{c}{b(2c+b)} & D = 0 \end{cases} \tag{68}$$

where *<sup>D</sup>* = *<sup>b</sup>*<sup>2</sup> − <sup>4</sup>*ac*, *<sup>b</sup>*<sup>1</sup> <sup>±</sup> <sup>=</sup> <sup>2</sup>*<sup>a</sup> b*± <sup>√</sup>*<sup>D</sup>* .

Then, from (67), (68) we obtain an expression for the parameter η, which characterizes the inhomogeneity of the FGM coating. In case of the parabolic variation of μ(*x*) at 0 ≤ *x* ≤ *h*, according to (64), it has the form

$$\eta = \frac{hB'(h)}{B(h)} = \frac{h}{\mu\_1 B(h)}\tag{69}$$

The law Δ(*t*) of indentation of a hard abrasive (half-plane *I*) into the coating is given by the formula

$$
\Delta(t) = \Delta\_0 \left( (e^{\varepsilon t} - 1)H(t\_0 - t) + H(t - t\_0) \right), \ t > 0 \tag{70}
$$

where *H*(*t*) is the Heaviside step function. Formula (70) assumes that the time section of indentation 0 < *t* ≤ *t*<sup>0</sup> is active and, when *t* > *t*<sup>0</sup> it is passive, since Δ(*t*) = Δ<sup>0</sup> at *t* > *t*0.

The nature of the loss of thermoelastic stability of the main parameters of the contact (temperature *T*(*h*, *t*), contact stresses σ*xx*(*h*, *t*), coating wear *u*w(*t*)) is studied in detail in Section 6 depending on the value of the parameter η. The boundaries of the thermoelastic stability region on the set of parameter values (*V*ˆ , *k*w) at fixed values η and Bi are also indicated there.

Let us study the effect of the parameter η of the considered thermoelastic problem on wear by a hard abrasive (half-plane *I*) of an elastic strip made from aluminum with a graded content of alumina (Al2O3) on the main characteristics of the contact: temperature *T*(*h*, *t*) from (57), contact stresses *p*(*t*) = −σ*xx*(*h*, *t*) from (59), coating wear *u*w(*t*) from (62) and wear rate . *u*w(*t*). This FGM strip is characterized by an increased shear modulus μ<sup>1</sup> = μ(*h*) = 125.0 GPa at the contact and shear modulus at the interface with the substrate <sup>μ</sup><sup>0</sup> = <sup>μ</sup>(0) = 25.0 GPa, <sup>ν</sup> = 0.34, <sup>κ</sup> = 88.1 × <sup>10</sup>−<sup>6</sup> mm2/s, <sup>α</sup> = 22.9 × <sup>10</sup>−<sup>6</sup> 1/K, *K* = 209.3 BT/(M·K), *f* = 0.47, *h* = 20 mm, Δ<sup>0</sup> = 5 mm, *t*<sup>0</sup> = 5 s, *V* = 2.5 mm/s, *ε* = ln 2. We consider three different values μ1/2 = μ(*h*/2), which, together with their corresponding values η, are presented in Table 1.


**Table 1.** Parameters of the shear modulus variation with the coating depth.

Variation of the shear modulus μ(*x*) by the *x*-coordinate is illustrated in Figure 5.

**Figure 5.** Variation of the shear modulus μ(*x*) by the coating depth at the considered values of μ1/2: *1*—50, *2*—75, *3*—100 GPa.

The wear of the coating surface at the depth Δ<sup>0</sup> ends up at *t* = *tw*, when the coating wear *u*w(*t*) equals Δ0, and the contact stress turns to zero (*p*(*t*) = −σ*xx*(*h*,*t*) = 0). We call *t*<sup>w</sup> the time of the coating wear by amount <sup>Δ</sup>0. Assuming the wear factor *<sup>K</sup>*<sup>∗</sup> = 1.0 × <sup>10</sup>−<sup>11</sup> m2/N, we obtain the values of dimensionless parameters *k*<sup>w</sup> = 0.511 and Bi =105 using Formula (29).

Table 2 gives the coating wear time, together with maximum values of contact pressure *p*(*t*) and temperature *T*(*h*, *t*) depending on the shear modulus in the middle of the coating μ1/2 and sliding velocity *V*.

**Table 2.** Values of main contact characteristics depending on shear modulus in the middle of coating μ1/2 and sliding velocity *V* in case of parabolic variation of the shear modulus μ(*x*).


The effect of the parameter η on the main contact characteristics is illustrated in Figure 6a–d, presenting plots of *<sup>T</sup>*(*h*, *<sup>t</sup>*), *<sup>p</sup>*(*t*), *<sup>u</sup>*w(*t*), . *u*w(*t*). Curves denoted *1*, *2*, *3* are plotted for μ1/2 = 50, 75, 100 GPa, respectively.

**Figure 6.** Dependence of the main contact characteristics by time: (**a**) contact pressure *<sup>p</sup>*(*t*) = <sup>−</sup>σ*xx*(*h*, *<sup>t</sup>*), (**b**) contact temperature *<sup>T</sup>*(*h*, *<sup>t</sup>*), (**c**) coating wear *<sup>u</sup>*w(*t*), (**d**) coating wear rate . *u*w(*t*), for different values of shear modulus in the middle of the coating μ1/2: *1*—50, *2*—75, *3*—100 GPa.

Figure 6 shows that the coating wear is accompanied by an increase in contact stress *p*(*t*) and temperature *T*(*h*, *t*). The growth of the contact stress *p*(*t*) is explained by the proportionality of the stresses σ*xx*(*x*, *t*) to the shear modulus μ(*x*) according to (3). The temperature *T*(*h*, *t*) is proportional to the contact stresses σ*xx*(*x*, *t*), according to (9), and, hence, to the shear modulus μ(*x*). Since curve *3* indicates a greater stiffness of FGM than curve 1 (Figure 5), the values of *p*(*t*) and *T*(*h*, *t*) for curve *3* will exceed the values of these characteristics for curve 1 (Figure 6). The wear of the coating demonstrating higher stiffness (curve *3*) will occur faster than the wear of the less stiff coating (curve *1*), for the same reason according to (8), and also because the vertical movement of abrasive I carried out forcibly according to the law for Δ(*t*). Hence, it follows that it is possible to experimentally determine the depth-variable shear modulus of the FGM coating under wear conditions in accordance with the kinematic law for Δ(*t*) of abrasive I, using pressure and temperature sensors at the contact.

For the preparation of materials with grinding, the power supply of the abrasive is often used due to the contact pressure *p*(*t*), which must be kept at a certain level. This restriction makes it possible to calculate the kinematic law for Δ(*t*) with the condition of the limitation of *p*(*t*). Let the formula for the contact stresses be given by the following expression

$$p(t) = p\_0 \left( \left( e^{\delta t} - 1 \right) H(t\_0 - t) + H(t - t\_0) \right) \tag{71}$$

where δ = *t* −1 <sup>0</sup> ln 2, and *t*<sup>0</sup> characterizes the time of reaching the stationary mode at *t* > *t*0, at which *p*(*t*) = *p*0. To express Δ(*t*) in terms of *p*(*t*), we use formula (26) in the form

$$p^L(p) = -\frac{2(1-v)}{1-2v} \Delta^L(p) \frac{h}{B(h)} \cdot \frac{N\_\sigma^0(z)}{R(z)}\tag{72}$$

Δ*L*(*p*) is determined from (72) by the following formula

$$\Delta^L(p) = -\frac{1-2v}{2(1-v)} p^L(p) \frac{B(h)}{h} \cdot \frac{R(z)}{N\_\sigma^0(z)}\tag{73}$$

After inverting the previous formula, taking into account (27), we obtain the following formula for determining Δ(*t*) from *p*(*t*)

$$
\Delta(t) = -\frac{1-2v}{2(1-v)} \frac{B(h)}{h} \int\_0^t p(\tau) \cdot f\_p(t-\tau)d\tau \tag{74}
$$

$$f\_p(t) = \frac{1}{2\pi i} \int\_{\Gamma} \frac{\mathcal{R}(z)}{t\_\kappa z \, r(h, z)} e^{z\tilde{t}} dz \tag{75}$$

where *<sup>t</sup>* <sup>=</sup> *<sup>t</sup> <sup>t</sup>*<sup>κ</sup> , *<sup>t</sup>*<sup>κ</sup> <sup>=</sup> *<sup>h</sup>*<sup>2</sup> <sup>κ</sup> , *R*(*z*) and *r*(*h*, *z*) are taken from (28), (29), and *B*(*h*) is from the expression (15).

The temperature at the contact *T*(*h*, *t*) with limited *p*(*t*) will also be limited, and wear is determined according to (62).

Similarly, it is possible to determine Δ(*t*) so that the temperature is limited, for example, by the following formula

$$T(t) = T\_0 \left( \left( e^{\Theta t} - 1 \right) H(t\_0 - t) + H(t - t\_0) \right) \tag{76}$$

where θ = *t* −1 <sup>0</sup> ln 2, and *t*<sup>0</sup> is the time, when *T*(*t*) becomes *T*<sup>0</sup> at *t* > *t*0. Using formula (22) for *TL*(*x*, *p*) according to the above scheme, we obtain

$$T^L(p) = \frac{1-\text{v}}{1+\text{v}} \frac{\text{\textdegree V}}{\alpha h} \Delta^L(p) \frac{hB'(h)}{B(h)} \frac{\text{N}\_T(h, z)}{R(z)}\tag{77}$$

Δ*L*(*p*) is determined from (77) by the following formula

$$
\Delta^L(p) = \frac{1+\text{v}}{1-\text{v}} \frac{\alpha h}{\hat{V}} T^L(p) \frac{B(h)}{hB'(h)} \frac{R(z)}{N\_T(h, z)}\tag{78}
$$

After inverting (78), we obtain

$$
\Delta(t) = \frac{1+\mathbf{v}}{1-\mathbf{v}} \frac{\alpha h}{\mathcal{V}} \frac{B(h)}{hB'(h)} \int\_0^t T(\tau) \cdot f\_T(t-\tau)d\tau \tag{79}
$$

$$f\_{\mathcal{T}}(t) = \frac{1}{2\pi i} \int\_{\Gamma} \frac{R(z)}{t\_{\kappa} N\_{\mathcal{T}}(h, z)} e^{z\tilde{t}} dz \tag{80}$$

where *<sup>t</sup>* <sup>=</sup> *<sup>t</sup> <sup>t</sup>*<sup>κ</sup> , *<sup>t</sup>*<sup>κ</sup> <sup>=</sup> *<sup>h</sup>*<sup>2</sup> κ .

When Δ(*t*) is taken from the expression (79), not only will the temperature on the contact be limited, but also the contact stress *p*(*t*). The wear is determined through Δ(*t*) from (79) according to (62).

#### **8. Discussion**

The present study was aimed at solving the problem of wear-grinding of a FGM coating with a depth-varying shear modulus. According to the results of the studies presented in the paper, it becomes feasible to solve the inverse problem of determining the law of distribution of the depth-varying shear modulus of the FGM l of the coating. It will be necessary to find out what the minimum set of information is for the unambiguous determination of the shear modulus varying by depth. The solution to this problem is important for the design of FGM coatings with special properties (anti-friction, anticorrosion, wear-resistant, etc.).

#### **9. Conclusions**

The considered thermoelastic contact problem of sliding wear by a hard abrasive material at a constant speed over the surface of an elastic coating in the form of a strip made of FGM with an arbitrary shear modulus varying with the depth of the coating made it possible to establish:


**Author Contributions:** Problem statement and development of the solution method, V.B.Z.; construction of the solution, V.B.Z. and P.A.L.; calculation, B.I.M.; analysis of the solution, B.I.M. and P.A.L.; conclusions, V.B.Z.; writing—original draft preparation, V.B.Z. and B.I.M.; writing—translation and editing, B.I.M. and P.A.L.; supervision, V.B.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Government of the Russian Federation, grant number 14.Z50.31.0046.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are available on request.

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

#### **Nomenclature**


#### **References**


#### *Article*
