*Article* **Effect of the Transverse Functional Gradient of the Thin Interfacial Inclusion Material on the Stress Distribution of the Bimaterial under Longitudinal Shear**

**Yosyf Piskozub 1,\* , Liubov Piskozub <sup>2</sup> and Heorhiy Sulym <sup>3</sup>**


**Abstract:** The effect of a functional gradient in the cross-section material (FGM) of a thin ribbon-like interfacial deformable inclusion on the stress–strain state of a piecewise homogeneous linear–elastic matrix under longitudinal shear conditions is considered. Based on the equations of elasticity theory, a mathematical model of such an FGM inclusion is constructed. An analytic–numerical analysis of the stress fields for some typical cases of the continuous functional gradient dependence of the mechanical properties of the inclusion material is performed. It is proposed to apply the constructed solutions to select the functional gradient properties of the inclusion material to optimize the stress–strain state in its vicinity under the given stresses. The derived equations are suitable with minor modifications for the description of micro-, meso- and nanoscale inclusions. Moreover, the conclusions and calculation results are easily transferable to similar problems of thermal conductivity and thermoelasticity with possible frictional heat dissipation.

**Keywords:** functionally graded material; thin inclusion; composites; nonperfect contact; frictional heating; crack; stress intensity factor

### **1. Introduction**

Many structural materials contain numerous thin inhomogeneities in the form of inclusions of different origins [1–7]. Quite often, these inclusions are used as elements to reinforce the structural parts of machines and structures or as fillers in composite materials or coatings [8–15]. The use of nanocomposites with specific properties in engineering and technology has significantly shifted the interest from the study of objects at the macro level (100–10−<sup>1</sup> m) and micro level (10<sup>−</sup>3–10−<sup>6</sup> m) to the nano level (10−<sup>9</sup> m) [16–20].

One of the typical examples of composite materials is the structure with thin ribbon inclusions. Structural elements made with the use of FGM have proven to be rather effective in practice [21–23]. In this way, it is possible to achieve a significant improvement in their mechanical, rheological, thermal, or other properties or the formation of protective thin layers [18,24–27].

The mathematical modeling of nanostructures requires the construction of more complex constitutional laws in comparison with the macro level [28,29]. Therefore, it is important to construct methods for studying the stress–strain state of such structures. To model a thin inclusion, there are mainly two basic approaches using analytical methods. The first one is based on the use of Eshelby's analytical solution [21,30,31] for an ellipsoidal inclusion, in which a limiting transition with a decrease in one of the characteristic dimensions of the inclusion is performed. However, its application to thin interphase inclusions

**Citation:** Piskozub, Y.; Piskozub, L.; Sulym, H. Effect of the Transverse Functional Gradient of the Thin Interfacial Inclusion Material on the Stress Distribution of the Bimaterial under Longitudinal Shear. *Materials* **2022**, *15*, 8591. https://doi.org/ 10.3390/ma15238591

Academic Editor: Francisco J. G. Silva

Received: 3 November 2022 Accepted: 28 November 2022 Published: 2 December 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/).

is impossible. The second one is based on the principle of the conjugation of continua of different measurability [32–37] and the method of jump functions [28,29,38–41]. According to this method, the inclusions are replaced by a certain surface (in the two-dimensional case, a line) of the discontinuity of the physical–mechanical fields, which describes the perturbing effect of a thin inclusion. Successful attempts were made in [17,33,42–45] to apply it to consider the influence of various physical and contact nonlinearities in the antiplane problem of elasticity theory for two compressed half-spaces with interfacial defects. The frictional slip with possible heat generation for contacting bodies [34,44,46–49] and the boundary element approach [50–52] were also considered here. Inhomogeneity of the mechanical properties of structural materials can be both designed for a specific purpose (FGM) and a consequence of technological processes of obtaining new materials and their processing (FSW, ball-burnishing process, etc.) [53,54]. Such factors cause additional complexity in the constitutive relations for the mathematical modeling of the behavior of such materials. However, the case of a thin inclusion of an inhomogeneous material has not been practically studied.

The process of improving the mathematical models of FGM [11,25,26,55–58] is complicated by the complex geometry of structural elements and the consideration of imperfections in the contact of their components. This is especially important to ensure their qualitative design, both in terms of mechanical strength [11,59–63] and in terms of the consideration of thermal, magnetic, and piezoelectric load factors [13,21,24,44,53,54,58,64,65].

The works [17,31,35,39] have been devoted to the consideration of surface energy and stresses in nanocomposites. Consideration of the heterogeneity of inclusions' properties at the micro and nano scale is particularly important because the heterogeneity density of matter (discrepancy variation in mechanical and other properties) of matter bodies with a decrease in their scale usually increases, and the impact of this heterogeneity increases even further.

This publication aims to develop the method of jump functions and to construct a convenient structured and highly versatile approach to study the longitudinal displacement and thermal heating of bodies with thin inclusions of arbitrary physical nature, including those made from FGM.

#### **2. Formulation of the Problem**

Consider an unbounded isotropic structure consisting of two half-spaces with elastic moduli *Gk* (*k* = 1, 2), which is subjected to an external longitudinal shear load determined by uniformly distributed infinity stresses *σ*<sup>∞</sup> *yz* and *σ*<sup>∞</sup> *yz*, concentrated forces of intensity *Qk* and screw dislocations with vector Burgers component *bk* at points *<sup>ς</sup>*∗*<sup>k</sup>* ∈ *Sk*(*<sup>k</sup>* = 1, 2). The stresses must satisfy the conditions *σ*<sup>∞</sup> *xz*2*G*<sup>1</sup> = *<sup>σ</sup>*<sup>∞</sup> *xz*1*G*<sup>2</sup> at infinity to ensure the straightness of the interface.

Let us investigate the stress–strain state (SSS) of the body section with a plane *xOy* perpendicular to the direction *Oz* of its longitudinal displacement (external problem). The plane sections of half-spaces perpendicular to this axis form two half-planes *Sk* (*k* = 1, 2) and the abscissa axis corresponds to the interface *L* ∼ *x* between them (Figure 1). At the interface of half-spaces (plane *xOz*), there is a tunnel section *L* = [−*a*; *a*] in the direction of the shear axis *z*, in which a certain object of general thickness 2*h* (*h a*) is inserted.

According to the paradigm of the method of jump functions [36], the presence of a thin inclusion in the bulk is modeled by jumps in the components of stress and displacement vectors in [38,40,41]:

$$\begin{array}{l} \left[ \sigma\_{yz} \right]\_h \cong \sigma\_{yz}^- - \sigma\_{yz}^+ = f\_3(\mathbf{x}), \ \mathbf{x} \in L'\\ \left[ \frac{\partial w}{\partial \mathbf{x}} \right]\_h \cong \frac{\partial w^-}{\partial \mathbf{x}} - \frac{\partial w^+}{\partial \mathbf{x}} = \left[ \frac{\sigma\_{\overline{\mathbf{x}}}}{G} \right]\_h \equiv \frac{\sigma\_{\overline{\mathbf{x}}}^-}{G\_1} - \frac{\sigma\_{\overline{\mathbf{x}}}^+}{G\_2} = f\_6(\mathbf{x}); \end{array} \tag{1}$$

$$f\_3(\mathbf{x}) = f\_6(\mathbf{x}) = 0,\text{ if } \mathbf{x} \notin L'. \tag{2}$$

**Figure 1.** The loading and geometric scheme of the problem.

It is hereinafter marked [•]*<sup>h</sup>* = •(*x*, −*h*) − •(*x*, +*h*), •*<sup>h</sup>* = •(*x*, −*h*) + •(*x*, +*h*); superscripts "+" and "−" correspond to the boundary values of the functions on the upper and lower banks of the line *L*.

The mathematical model of a thin inclusion is given as complicated conditions of imperfect contact between opposite matrix surfaces adjacent to the inclusion (internal problem) [28,38–40,65]. The general model of a thin, physically nonlinear inclusion is presented in [28,29,38], where the methods of modeling thin objects involve the integration of the defining relations describing the physical and mechanical properties of the material of the inclusion, with the subsequent consideration of the smallness of one of the linear dimensions of the inclusion.

Let us consider a similar model for a thin inclusion, assuming that the mechanical properties of the inclusion material are coordinate-dependent. This will allow us to model the inclusions of a functionally graded material:

$$\begin{cases} \left. G\_x^{in}(\mathbf{x}) \left< \frac{\partial w^{in}}{\partial \mathbf{x}} \right>\_{\hbar}(\mathbf{x}) - 2\sigma\_{\mathbf{x}\mathbf{2}}^{in}(-a) - \frac{1}{\hbar} \int\_{-a}^{\mathbf{x}} \left[ \sigma\_{\mathbf{y}\mathbf{2}}^{in} \right]\_{\hbar}(\mathbf{f}) d\mathbf{f} = 0, \\\ \left. G\_y^{in}(\mathbf{x}) \left[ w^{in} \right]\_{\hbar}(\mathbf{x}) + h \left< \sigma\_{\mathbf{y}\mathbf{2}}^{in} \right>\_{\hbar}(\mathbf{x}) = 0. \end{cases} \tag{3}$$

Here, *Gin <sup>x</sup>* (*x*), *Gin <sup>y</sup>* (*x*) are the variable shear moduli of the inclusion's material. As a special case, considering their values to be constant, we obtain Hooke's law. The upper index "*in*" denotes the terms describing the inclusion material's SSS components.

Contact between matrix components and the inclusion at *L* and between the bimaterial structure components along a line *L*\*L* is supposed to be mechanically perfect,

$$\begin{aligned} w(\mathbf{x}, \pm h) &= w^{in}(\mathbf{x}, \pm h), & \sigma^{in}\_{yz}(\mathbf{x}, \pm h) &= \sigma\_{yzk}(\mathbf{x}, \pm h), & \mathbf{x} \in L',\\ w(\mathbf{x}, +0) &= w(\mathbf{x}, -0), & \\ \sigma\_{yz2}(\mathbf{x}, +0) &= \sigma\_{yz1}(\mathbf{x}, -0), & \mathbf{x} \in L \backslash L', \end{aligned} \tag{4}$$

or frictional in some areas *x* ∈ *Lf* ⊂ *L* , as was considered in works [41,46–48],

$$
\sigma\_{yz}^{\rm in}(\mathbf{x}, \pm h) = \sigma\_{yzk}(\mathbf{x}, \pm h) = -\text{sgn}\left[w^{\rm in}\right]\_h \tau\_{yz}^{\rm max}.\tag{5}
$$

Here, *τ*max *yzK* is the limit value of shear stresses, at which the slippage begins. In this case, however, additional iterative methods should be applied to determine the area of the slip zones depending on the specific types of external loading of the composite [41].

### **3. Materials and Methods**

Expressions for the components of the stress tensor and the derivatives of displacements on the line *L* of the infinite plane *S* = *S*<sup>1</sup> ∪ *S*2, as well as inside the latter, can be obtained by applying the results of [37] to the solution of the external problem

$$\begin{array}{l} \sigma\_{yz}(z) + i\sigma\_{xz}(z) = \sigma\_{yz}^{0}(z) + i\sigma\_{xz}^{0}(z) + \\ + ip\_{k}\mathfrak{g}\_{3}(z) - \mathbb{C}\mathfrak{g}\_{6}(z) \quad (z \in \mathbb{S}\_{k}; \ r = \mathfrak{3}, \mathfrak{6}; \ k = 1, 2); \end{array} \tag{6}$$

$$\begin{array}{l} \sigma\_{yzk}^{\pm}(\mathbf{x}) = \mp p\_k f\_3(\mathbf{x}) - \mathsf{C}g\_6(\mathbf{x}) + \sigma\_{yz}^{0\pm}(\mathbf{x}),\\ \sigma\_{xzk}^{\pm}(\mathbf{x}) = \mp \mathsf{C}f\_6(\mathbf{x}) + p\_k g\_3(\mathbf{x}) + \sigma\_{xz}^{0\pm}(\mathbf{x}).\end{array} \tag{7}$$

where the notation [28,41] is introduced:

$$\begin{array}{ll} \mathcal{g}\_r(z) \equiv \frac{1}{\pi} \int \frac{f\_r(\mathbf{x})d\mathbf{x}}{\mathbf{x} - \mathbf{z}} \; , \quad \mathcal{s}\_r(\mathbf{x}) \equiv \int\_{-a}^{\mathbf{x}} f\_r(\mathbf{x})d\mathbf{x} \; \\\ p = \frac{1}{\mathcal{G}\_1 + \mathcal{G}\_2} \; \; p\_k = \mathcal{G}\_k p\_r \; \mathcal{C} = \mathcal{G}\_1 \mathcal{G}\_2 p. \end{array} \tag{8}$$

The superscript "+" corresponds to *k =* 2 and "–" corresponds to *k =* 1. Values marked with superscript "0" correspond to values in a continuous medium without inclusions under the same external load (homogeneous solution) [28,41].

Using (7), (8) and boundary conditions (4), it is easy to obtain from model (3) a system of singular integral equations:

$$\begin{cases} (p\_2 - p\_1)f\_6(\mathbf{x}) + 2p\_3\mathbf{x}(\mathbf{x}) - \frac{s\_3(\mathbf{x})}{hG\_x^{in}(\mathbf{x})} = F\_3(\mathbf{x}),\\ (p\_2 - p\_1)f\_3(\mathbf{x}) + 2Cg\_6(\mathbf{x}) - \frac{G\_y^{in}(\mathbf{x})s\_6(\mathbf{x})}{h} = F\_6(\mathbf{x}),\\ F\_3(\mathbf{x}) = \frac{2}{G\_x^{in}}\sigma\_{\mathbf{x}2}^{in}(-a) - \left(\sigma\_{\mathbf{x}22}^0(\mathbf{x})/G\_2 + \sigma\_{\mathbf{x}21}^0(\mathbf{x})/G\_1\right),\\ F\_6(\mathbf{x}) = \left\langle \sigma\_{yz}^0 \right\rangle(\mathbf{x}) - G\_y^{in}\left\langle \frac{\sigma\_{yz}^0(\mathbf{x})}{G\_k} \right\rangle - \frac{G\_y^{in}}{h}\left[\frac{0}{w}\right](-a). \end{cases} \tag{9}$$

Balance conditions on the power balance and unambiguity of displacements while moving around the thin defect must be added to the solution of the external problem:

$$\begin{aligned} \int\_{-a}^{a} f\_3(\xi) d\xi &= 2h \left( \sigma\_{xz}^{in}(a) - \sigma\_{xz}^{in}(-a) \right), \\ \int\_{-a}^{a} f\_6(\xi) d\xi &= [w](a) - [w](-a). \end{aligned} \tag{10}$$

Solving (9) and (10) using the methods in [29,38,41], it is easy to obtain a system of linear algebraic equations with unknown coefficients of the decomposition of the jump functions *fr*(*x*) into a series by orthogonal Jacobi or Chebyshev polynomials.

An important aspect of the study of the strength of such structures is the improvement of their strength criteria. In fracture mechanics, it is acceptable to use the stress intensity factor to describe the behavior of the SSS in the vicinity of the crack tip [42,45,61–63,66]. This is not sufficient for the case of a thin deformable inclusion. In [45], the authors obtained the two-term asymptotical expressions for the distribution of SSS in the vicinity of the thin inclusion tips using the introduced generalized stress intensity factors (GSIF):

$$K\_{31} + iK\_{32} = \lim\_{\substack{r \to 0 \ (\theta = 0, \pi)}} \sqrt{2\pi r} (\sigma\_{yz} + i\sigma\_{xz}) \,. \tag{11}$$

Here, (*r*, *θ*) is a system of polar coordinates with the origin near the right or the left tip of the inclusion *z* = ±*r* exp(*iθ*) ± *a*.

Considering the well-known mathematical analogy [67], the obtained solutions to the antiplane problem can be regarded as solutions to the accordant heat conduction problem, if we take into account the correspondence of the values

 $w \sim T$ ,  $\frac{\partial w}{\partial x} \sim \frac{\partial T}{\partial x}$ ,  $\frac{\partial w}{\partial y} \sim \frac{\partial T}{\partial y}$ ,  $q\_X \sim \sigma\_{xx}$ ,  $q\_Y \sim \sigma\_{yz}$ ,  $\frac{\partial T}{\partial x} \sim \lambda\_{xy}$ ,  $K\_{31} \sim k\_{qy}$ ,  $K\_{32} \sim k\_{qx}$ .

The terms are as follows: *T*—temperature, *qx*, *qy*—heat flows, *λx*, *λy*—thermal conductivity coefficients, *kqy*, *kqx*—heat flow intensity factors [40].

#### **4. Numerical Results and Discussion**

Since the main focus of this article is to investigate the effect of the functional gradient on the mechanical properties of the inclusion material, we will limit ourselves to one of the most representative variants of the structure loading: homogeneous longitudinal shear *σ*<sup>∞</sup> *yz* = *τ* and *σ*<sup>∞</sup> *xzk* = *τ<sup>k</sup>* (*k* = 1, 2) at infinity. However, the calculations for loading by concentrated force factors or dislocations do not make any fundamental difference except for the necessity to consider the locality of their application [29,41].

The dependence *Gx*(*x*), *Gy*(*x*) on coordinate *x* for mathematical modeling can be defined as an arbitrary function (linear, exponential, power, periodic [58], etc.), which adequately reflects the desired practical properties of the material. To illustrate the method, let us consider one of the illustrative variants of the functional gradient of the inclusion material—the piecewise linear one:

$$\mathbf{G}\_{\mathbf{x}}(\mathbf{x}) = \mathbf{G}\_{\mathcal{Y}}(\mathbf{x}) = \begin{cases} (\mathbf{G}\_{01} - \mathbf{G}\_{0})\frac{\mathbf{x}}{a} + \mathbf{G}\_{01}, & \mathbf{x} \in [-a, 0]; \\ (\mathbf{G}\_{02} - \mathbf{G}\_{01})\frac{\mathbf{x}}{a} + \mathbf{G}\_{01}, & \mathbf{x} \in [0, a]. \end{cases} \tag{12}$$

where *G*0, *G*01, *G*02—some given constants.

To significantly reduce the number of calculations without loss of generality, it is convenient to use the following dimensionless quantities, marked by symbol "~" (tilde) on top:

*<sup>x</sup>* <sup>=</sup> *<sup>x</sup>*/*a*, *<sup>h</sup>* <sup>=</sup> *<sup>h</sup>*/*a*, *<sup>y</sup>* <sup>=</sup> *<sup>y</sup>*/*a*, *Gin <sup>x</sup>* (*<sup>x</sup>*) = *<sup>G</sup>in <sup>x</sup>* (*x*)/*Ggav*, *<sup>G</sup>in <sup>y</sup>* (*<sup>x</sup>*) = *<sup>G</sup>in <sup>y</sup>* (*x*)/*Ggav*, *<sup>τ</sup><sup>k</sup>* <sup>=</sup> *<sup>τ</sup>k*/*Ggav*, *<sup>τ</sup>* <sup>=</sup> *<sup>τ</sup>*/*Ggav*, *Ggav* <sup>=</sup> <sup>√</sup>*G*1*G*2, *G*<sup>0</sup> = *G*0/*Ggav*, *G*<sup>01</sup> = *G*01/*Ggav*, *G*<sup>02</sup> = *G*02/*Ggav*, *<sup>σ</sup>xz*(*<sup>x</sup>*) = *<sup>σ</sup>xz*(*x*)/*Ggav*, *<sup>σ</sup>yz*(*<sup>x</sup>*) = *<sup>σ</sup>yz*(*x*)/*Ggav*. *K* <sup>31</sup> <sup>=</sup> *<sup>K</sup>*<sup>+</sup> 31 <sup>2</sup>*CG gav*√*π<sup>a</sup>* , *K* <sup>32</sup> <sup>=</sup> *<sup>K</sup>*<sup>+</sup> 32 <sup>2</sup>*p*2*Ggav*√*π<sup>a</sup>* ,

where *K*<sup>+</sup> <sup>31</sup> *<sup>K</sup>*<sup>+</sup> <sup>32</sup> are the GSIFs near the tip *x* = +*a* of the inclusion.

Figures 2–11 illustrate the dependence of the stress–strain behavior of the matrix in the inclusion vicinity on the variation in the parameters *G*0, *G*01, *G*02, the values of which were chosen to reveal a qualitative picture of the FGM effect on the stress–strain parameters. It can be immediately concluded from Figures 2 and 3 that under the load *<sup>τ</sup>*, the dimensionless *<sup>K</sup>*<sup>+</sup> <sup>31</sup> are expected to decrease with the increasing shear moduli of any part of the inclusion, while at *<sup>K</sup>*<sup>+</sup> <sup>32</sup> they appear to increase with increasing load *<sup>τ</sup><sup>k</sup>*.

The effect of changes in the moduli *<sup>G</sup><sup>x</sup>*(*x*), *<sup>G</sup><sup>y</sup>*(*x*) on the stresses *<sup>σ</sup>yz*, *<sup>σ</sup>xz* on the inclusion surface is more obvious if we choose a linear growth law for them along the inclusion axis (Figures 4–6). The magnitude of the surface stresses increases significantly in the stiffer part of the inclusion. The larger the stiffness gradient, the more significant the increase.

The choice of the piecewise linear law of moduli *G<sup>x</sup>*(*x*), *G<sup>y</sup>*(*x*) change in the Formulae (14) as *G*<sup>01</sup> = *G*<sup>02</sup> (variant 1) or *G*<sup>0</sup> = *G*<sup>01</sup> (variant 2) has a more contrasting effect on the surface stresses *<sup>σ</sup>yz*, *<sup>σ</sup>xz*, especially in the vicinity of the gradient breaking point *<sup>x</sup>* <sup>=</sup> <sup>0</sup> (Figures 7 and 8). Moreover, variant 2 of the functional dependence of the inclusion material moduli leads to partial unloading in the softer part of the inclusion near the breaking point *x* = 0 (Figure 8).

**Figure 2.** Influence of the parameters *<sup>G</sup>*0, *<sup>G</sup>*01, *<sup>G</sup>*<sup>02</sup> on the GSIF *<sup>K</sup>*<sup>+</sup> <sup>31</sup> under the load, uniformly distributed on infinity stress *σ*<sup>∞</sup> *yz* = *τ*.

**Figure 3.** Influence of the parameters *<sup>G</sup>*0, *<sup>G</sup>*01, *<sup>G</sup>*<sup>02</sup> on the GSIF *<sup>K</sup>*<sup>+</sup> <sup>32</sup> under the load, uniformly distributed on infinity stress *σ*<sup>∞</sup> *xzk* = *τ<sup>k</sup>* (*k* = 1, 2).

**Figure 4.** Stress distribution along with the upper interface (inclusion–matrix half-space *S*2) with a linear distribution of material stiffness.

**Figure 5.** Stress distribution along with the upper interface (inclusion–matrix half-space *S*2) with a linear distribution of material stiffness.

**Figure 6.** Stress distribution along with the upper interface (inclusion–matrix half-space *S*2) with a linear distribution of material stiffness.

**Figure 7.** Stress distribution along with the upper interface (inclusion—matrix half-space *S*2) with a piecewise linear distribution of material stiffness.

**Figure 8.** Stress distribution along with the upper interface (inclusion—matrix half-space *S*2) with a piecewise linear distribution of material stiffness.

**Figure 9.** Stress distribution in the matrix at the vicinity of the inclusion with a linear distribution of material stiffness.

**Figure 10.** Stress distribution in the matrix at the vicinity of the inclusion with the piecewise linear distribution of material stiffness.

**Figure 11.** Stress distribution in the matrix at the vicinity of the inclusion with the piecewise linear distribution of material stiffness.

Figures 9–11 illustrate the changes in the stress field in the matrix in the inclusion vicinity under different variants of the law of functional change of the inclusion material moduli. The trends towards a decrease in the stress magnitudes in the vicinity of the stiffer parts of the inclusion are visible.

#### **5. Conclusions**

The proposed sufficiently simple and mathematically correct methodology made it possible for us to construct, for the first time, a mathematical model of a deformable thin linear interfacial inclusion with essentially inhomogeneous linear mechanical properties. Such a model can be used to simulate a thin inclusion from a functionally graded material and to solve the corresponding problems of defining the stress–strain field of the corresponding micro- or nanostructures by efficient analytical–numerical methods (the jump function method and its modifications), without the need to involve purely numerical approaches (in particular, FEM).

The calculations of the stress–strain field components for simple test cases of the functional dependence of the shear moduli of inclusion material have demonstrated the expected qualitative picture of their effect on the variation in the FGM parameters. In particular, (1) the stress magnitude increases significantly in the vicinity of the inclusion regions with increased stiffness; (2) the combination of the inclusion materials from parts with piecewise linear mechanical characteristics may lead to partial unloading of the inclusion and matrix in their softer part in the vicinity of the breaking point of the gradient dependence of the inclusion material parameters; (3) the contrast of the stress field changes of the inclusion and matrix is proportional to the increase in the gradient dependence.

The conclusions and calculation results are easily transferable to analogous problems of thermal conductivity and thermoelasticity with possible frictional heat generation and can be used for recommendations on the optimal operating parameters of structures.

The discussed conclusions can be useful in designing the functionally gradient mechanical properties of the material of inclusions and in the optimization of engineering structures to increase their strength and service life. The proposed method is effective for solving a wide class of problems of deformation of solids with thin deformable inclusions of finite length and can be used for SSS calculation for different FGM inclusions.

**Author Contributions:** Conceptualization, Y.P. and H.S.; methodology, Y.P. and H.S.; software, Y.P. and L.P.; validation, Y.P. and H.S.; formal analysis, Y.P.; investigation, Y.P. and H.S.; resources, H.S.; data curation, Y.P.; writing—original draft preparation, Y.P.; writing—review and editing, Y.P.;

visualization, Y.P. and L.P.; supervision, H.S.; project administration, H.S.; funding acquisition, H.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This investigation was performed within the framework of research project No. 2017/27/B/ ST8/01249, funded by the National Science Centre, Poland, and with project financing through the program of the Minister of Education and Science of Poland named "Regional Initiative of Excellence" in 2019–2022, project No. 011/RID/2018/19; amount of financing: 12,000,000 PLN.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are openly available at https://ua1lib.org /book/665574/5c937e (accessed on 9 January 2022), reference number [37]; https://doi.org/10.3390/ma 15041435 (accessed on 9 January 2022), reference number [38]; and https://doi.org/10.3390/ma14174928 (accessed on 9 January 2022), reference number [41].

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

### **Nomenclature**


### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Materials* Editorial Office E-mail: materials@mdpi.com www.mdpi.com/journal/materials

Academic Open Access Publishing

www.mdpi.com ISBN 978-3-0365-8462-1