Development of a Model of Crack Propagation in Multilayer Hard Coatings under Conditions of Stochastic Force Impact

The article deals with the problems of cracking in the structure of multilayered coatings under the conditions of stochastic loading process. A mathematical model has been proposed in order to predict the crack propagation velocity in the coating while taking the influence of interlayer interfaces into account. A technique for calculating the probability density distribution of the coating fracture (failure rate) has been developed. The probability of a change in the crack growth direction is compared with the experimental data that were obtained as a result of the studies focused on the pattern of cracking in the Zr,Nb-(Zr,Nb)N-(Zr,Nb,Al)N and Ti-TiN-(Ti,Cr,Al)N coatings under the conditions of the real stochastic loading of cutting tools during the turning. The influence of the crystalline structure of the coating on the cracking pattern has been studied. The investigation has found the significant effect of the crystalline structure of the coating layers on the cracking pattern.


Introduction
Brittle fracture is a crucial and often dominant factor that reduces the efficiency of wear-resistant coatings and cutting tools [1]. In general, the coatings are very hard and rather brittle coverings with a high tendency to crack formation [2,3]. Accordingly, the prediction of the crack formation in coatings and creation of the coatings with enhanced resistance to cracking is an important scientific and practical challenge. It is necessary to consider the fact that the cutting of materials is a complex stochastic process, which combines the influence of force and temperature factors, as well as oxidation, diffusion, and other processes. To date, there is no model that is able to take the effect of all the above factors into account, and it is hardly possible to create such a global model. Meanwhile, some models have been developed, and researchers continue to develop models that consider the process of coating functioning from different points of view, while using various modeling techniques. A finite element model (FEM) is quite often applied in order to simulate crack formation. In particular, in [4], brittle coating on an elastic substrate in four-point bending is investigated. The studies considered the influence of such factors, as random strength distribution, coating thickness, residual stresses, and coating modulus. The influence on the crack formation of such parameters, such as the number of load The conditions of stochastic loading typical for the zone of elastic contact between the chips and the tool rake face are considered [38][39][40]. During the cutting, the tool material is being affected by stochastically influencing force factors. As is known, the external sliding friction acts in the zone of elastic contact, and the stochastic friction force between the chips and the tool rake face determines the resistance to the movement of the chips [39,40]. The tangential contact stresses in the section of the elastic contact gradually decrease, while the distance from the cutting edge increases. During the cutting process, the friction process in the zone of elastic contact is determined by overcoming intermolecular bonds, that is, the adhesive interaction. The considered area is characterized by the influence of the stochastic tangential contact stresses in combination with the effect of high temperature, and both of these factors gradually decrease while the distance from the cutting edge grows.
The analysis of the papers published earlier finds that the proposed models of crack propagation in the coatings cannot be considered to be fully consistent with the real operation conditions, since the models do not take the influence of the following key factors into account: • stochastic nature of the loading process, and • multilayered structure of the coating, in which adjacent layers have different physical and mechanical properties, and interlayer interfaces have significant influence on the cracking pattern.
Thus, the task of this paper was to develop a mathematical model for predicting the development of a crack in a multilayer coating under the conditions of a stochastic load process. Because of the complex challenge to take into account all of the conditions of the cutting process in a single mathematical model, this paper only considers the effect of force factors. No factor influencing temperature, oxidation, and diffusion processes were taken into account. To make the modeling process easier, the task was divided into two stages, as follows: • predicting the time during which a crack passes a layer of a multilayer coating, and • predicting a change in the direction of crack development, while a crack passes an interface between two coating layers with different physical and mechanical properties.
As is known, successive alternation of nanolayers with identical properties forms the nanolayer structure of the coating. Thus, it will be possible to easily predict the crack propagation in a nanolayer structure with any number of layers after overcoming the challenge of predicting the development of a crack in a two-layer or three-layer system.
An analysis of the microstructure and nanostructure of the samples was studied while using JEM 2100 (JEOL, Akishima, Japan) high-resolution (HR) transmission electron microscope (TEM) with an accelerating voltage of 200 kV. An FEI Quanta 600 FEG (Thermo Fisher Scientific, Waltham, MA, USA) scanning electron microscope (SEM) was also used to study the coating microstructure.
The mechanical and tribological properties of the Ti-TiN-(Ti,Cr,Al)N coating were investigated in our previous works [54,60].

Passage of a Crack through a Coating Layer
A coating layer on the 0xy plane, with thickness of l* (Figure 1) was considered. microscope (TEM) with an accelerating voltage of 200 kV. An FEI Quanta 600 FEG (Thermo Fisher Scientific, Waltham, MA, USA) scanning electron microscope (SEM) was also used to study the coating microstructure.
The mechanical and tribological properties of the Ti-TiN-(Ti,Cr,Al)N coating were investigated in our previous works [54,60].

Passage of a Crack through a Coating Layer
A coating layer on the 0xy plane, with thickness of l* (Figure 1) was considered. A crack develops along the x axis, depending on the time point t or the number of cycles, and it is a random function in general case. At the time point t, the crack tip is at the point x = l(t). At each time point, the crack either retains its state or develops along the x axis. Reverse motion and possible deflection of the crack is not considered. Thus, the state of the crack at the time point t + Δt is either l(t + Δt) or l + Δl, into which the crack tip passes with the probabilities: Let ( , | 0 , 0 ) be a conditional probability that the crack tip is within the interval of ( , + ) at the time point t, if at the time point t0, the crack tip is at the point l0. Let us assume that the process of the transformation of the crack tip from the given state into the next one only depends on this state and does not depend on the previous one (the history of the process), then the theory of Markov processes can be applied in order to describe the crack formation process, according to which all probabilistic information about the process (multivariate distribution) is expressed through the conditional probabilities ( , | 0 , 0 ), which satisfy the following equation: ( , | 0 , 0 ) = ( , − | 0 , 0 ) × (1 − ( )) + ( − , − | 0 , 0 ) × ( − ) (2) Let the mean expected value for the conditional increment of the crack tip coordinate (at the fixed l(t)) in the short time Δt be denoted through: Accordingly, the dispersion (spread) of the above increment is denoted by: The limiting values (3), (4) are denoted, respectively: Through the calculation of the mean value (3) of the conditional increment of the crack tip coordinate over the time t, we obtain: A crack develops along the x axis, depending on the time point t or the number of cycles, and it is a random function in general case. At the time point t, the crack tip is at the point x = l(t). At each time point, the crack either retains its state or develops along the x axis. Reverse motion and possible deflection of the crack is not considered. Thus, the state of the crack at the time point t + ∆t is either l(t + ∆t) or l + ∆l, into which the crack tip passes with the probabilities: Let Π(l, t|l 0 , t 0 ) be a conditional probability that the crack tip is within the interval of (l, l + ∆l) at the time point t, if at the time point t 0 , the crack tip is at the point l 0 . Let us assume that the process of the transformation of the crack tip from the given state into the next one only depends on this state and does not depend on the previous one (the history of the process), then the theory of Markov processes can be applied in order to describe the crack formation process, according to which all probabilistic information about the process (multivariate distribution) is expressed through the conditional probabilities Π(l, t|l 0 , t 0 ), which satisfy the following equation: Let the mean expected value for the conditional increment of the crack tip coordinate (at the fixed l(t)) in the short time ∆t be denoted through: Accordingly, the dispersion (spread) of the above increment is denoted by: The limiting values (3), (4) are denoted, respectively: Through the calculation of the mean value (3) of the conditional increment of the crack tip coordinate over the time t, we obtain: For dispersion, respectively: Subsequently Let us assume that (∆l) 2 = c∆t Afterwards For the random function of l(t) with the calculated mean A and the dispersion B, it is possible to get a differential equation for the probability Π(l, t|l 0 , t 0 ). This is a Fokker-Planck-Kolmogorov equation.

∂Π ∂t
The Equation (10) is a parabolic-type equation of the mathematical physics, for the solution of which the analytical and numerical methods are used. A discrete analogue of the crack propagation can be obtained if it is assumed that the time t and length l are quantized in time and states. The time takes on the values of t = 0, 1, 2, . . . n, l = j = 0, 1, 2, . . . n, l 0 = i = 0, 1, 2, . . . n.
Subsequently, a replacement can be made: Π(l, t|l 0 , t 0 ) = Π ij (m), m = t − t 0 . A discrete equation of Markov type, which Π ij satisfies, has a form, as follows: Here, Π ij (m) is the probability of transformation from the state i to the state j in m steps. It is clear that i, j are the discrete values of the crack length, m is the number of cycles, during which a crack develops from the length i to the length j.
In the theory of Markov processes, Equation (10) is called a direct Fokker-Planck-Kolmogorov (FPK) equation. An inverse FPK equation can be obtained in a similar way: The difference between the Equations (10) and (12) is that, in (10), the right-hand side of the equation contains the differentiation with respect to l, while, in (12), the differentiation with respect to l 0 .
The process of solving an applied problem involves equations of both types, depending on the formulation of the problem. When it is required to define the probability that a crack tip reaches the length l * at the time point T * , it is advisable to use an inverse equation, which takes the influence of the initial conditions into account. It is advisable to use the direct equation if at the initial time, the conditions are determined in a probabilistic manner, i.e., the crack initiation phase l 0 is a random variable. For the Equations (10) and (12), the initial condition is set for the deterministic initial crack length l 0 : The condition for non-negativeness of the probability is as follows: Materials 2021, 14, 260 6 of 20 The condition for normalization of the probability is: The boundary conditions can be different, depending on a specific problem to be solved. If the crack length l(t) at the time point t 0 is not constant and it is a random variable with the distribution density P 0 (l), then the condition is as follows: The probability densities should be consistent, i.e., they should satisfy the following equation: Π(l, t|l 0 , t 0 ) , makes it easier to solve the problem mathematically and, therefore, let us formulate the basic relations for finding P(l,t), taking into account (17).
There are different algorithms for finding P(l,t), depending on the problem to be solved. At first, a functional solution to the Equation (12) with the initial condition (13) can be found, and then the Equation (17) can be solved.
Alternatively, an FPK equation can be directly obtained for P(l,t): Additionally, the initial and boundary conditions can be formulated for it. The initial conditions are as follows: Let us consider a model of a rough upper layer and assume that the interface x = 0 is an average surface, while a real surface has microroughness, where the value is determined by the random variable l 0 , with 0 < l 0 l * (Figure 2).
The condition for non-negativeness of the probability is as follows: ( , | 0 , 0 ) ≥ 0 (14) The condition for normalization of the probability is: The boundary conditions can be different, depending on a specific problem to be solved.
If the crack length l(t) at the time point t0 is not constant and it is a random variable with the distribution density P0(l), then the condition is as follows: The probability densities should be consistent, i.e., they should satisfy the following equation: ( , | 0 , 0 ) , makes it easier to solve the problem mathematically and, therefore, let us formulate the basic relations for finding P(l,t), taking into account (17).
There are different algorithms for finding P(l,t), depending on the problem to be solved.
At first, a functional solution to the Equation (12) with the initial condition (13) can be found, and then the Equation (17) can be solved.
Alternatively, an FPK equation can be directly obtained for P(l,t): Additionally, the initial and boundary conditions can be formulated for it. The initial conditions are as follows: Let us consider a model of a rough upper layer and assume that the interface x=0 is an average surface, while a real surface has microroughness, where the value is determined by the random variable l0, with 0 < 0 ≪ * ( Figure 2).
The growth of a crack starts from the length l0 (the length of the microcrack), taking the micro roughness into account. The growth of a crack starts from the length l 0 (the length of the microcrack), taking the micro roughness into account.
The initial condition for the probability density has the form, as follows: Let F 0,l (t, l 0 ) be the probability that the microcrack with the initial value of l 0 reaches the layer interface during the time period t > 0: either l = 0 or l = l * . The case of reaching l = 0 corresponds to the motion of the crack tip towards the free surface, i.e., spallation, and the case of l = l * exhibits the through fracture by the crack of the upper layer with the thickness of l * .
A case when the crack does not reach the layer interface, which is, the crack is damping in a layer, is determined by the probability of non-fracture of the layer.
where T(0, l 0 , l * ) is a random moment when the crack reaches the interface l * , when the mean (expected value) and dispersions are constant: has the form, as follows: Let us consider a model for which a crack develops only rectilinearly into the depth of a layer, without spallation and deviation in the layer. In this case, the initial and boundary conditions are as follows: The following condition should be satisfied: Let us denote the stationary approximation of the solution of Equation (22) on reaching by the crack tip of the interface l = l* through F l * (l 0 ). The condition (24) means that the probability of the crack tip reaching the state of l* from the state of l 0 is, in practice, a persistent event.
The stationary approximation of F l * (l 0 ) satisfies the equation that follows from (22) at Under the boundary conditions of the solution of the Equation (25) under the condition (26) in the analytical form is as follows: It is clear that In order to further specify the problem, it is necessary to know the mean (expected value) A(x) and the B(x), a dispersion of the random function l(t), where the probability density satisfies the Fokker-Planck-Kolmogorov equations. As a rule, in this case, the l(t) satisfies the stochastic differential equation: Other forms of the Equation (30) are as follows: The presented formulas make it possible to use Markov models of random processes to describe the processes of fatigue fracture that occur due to crack development.
The Paris' law, based on the fact that all events on a crack tip, including the crack propagation velocity, depend on the stress intensity factor, can be written as: where n is the number of load cycles, ∆K = K max − K min , K is the stress intensity factor. While addressing the time dependence, let us denote n·τ = t, where τ is the cycle time, then the Equation (34) can be presented as: The Equation (35) describes a deterministic process of the crack growth. However, in fact, it is a stochastic process, in which only the average crack propagation velocity can be determined, with an error of n(t), where the average value is n(t) = 0. Subsequently, the Equation (35) can be randomized and, as a result, after adding a random error n(t) to the right-hand side of the Equation (35):

Calculation of Coating Reliability Based on a Model with Independently (Sequentially) Fracturing Layers
The state of each layer does not depend on the states of other layers. A fracture of the upper layer leads to the failure of the entire coating. Let Qς(t) be the probability of the operational condition of the entire coating at the time point t and T i -the period of time during which an ith layer stays not damaged (a random variable). Afterwards, The probability of simultaneous event written in brackets for independent events is: By denoting we obtain The probability of the coating fracture is: The probability density (failure rate) is introduced by the formula: In a case when where λ i is the fracture rate. Subsequently, the average time until the fracture of the three-layer structure is: Two layers, which were rigidly connected to each other, were considered, with the different physical and mechanical characteristics of G i , v i , ρ i , i = 1.2 (shear modulus, Poisson's ratio, density) ( Figure 3).
we obtain The probability of the coating fracture is: The probability density (failure rate) is introduced by the formula: In a case when where λi is the fracture rate. Subsequently, the average time until the fracture of the three-layer structure is: Two layers, which were rigidly connected to each other, were considered, with the different physical and mechanical characteristics of , , , = 1.2 (shear modulus, Poisson's ratio, density) ( Figure 3). The crack trajectory (Material 1) is curvilinear in a real medium, and it fluctuates near the x axis. The average time of the layer fracture depends on the crack length. The minimum time of the layer fracture will be obtained for a crack rectilinear in average, and this is the lowest estimate of the crack resistance of the layer. The maximum time of crack resistance will be obtained in the case when the crack trajectory significantly deflects from the straight line (for example, it is a fractal curve with a dimension over 1). As is known, the crack propagation is the process of transformation of the deformation energy into the The crack trajectory (Material 1) is curvilinear in a real medium, and it fluctuates near the x axis. The average time of the layer fracture depends on the crack length. The minimum time of the layer fracture will be obtained for a crack rectilinear in average, and this is the lowest estimate of the crack resistance of the layer. The maximum time of crack resistance will be obtained in the case when the crack trajectory significantly deflects from the straight line (for example, it is a fractal curve with a dimension over 1). As is known, the crack propagation is the process of transformation of the deformation energy into the bond-breaking power, i.e., the formation of a free surface [31]. Thus, the crack trajectory is a line (surface), with a density vector of the elastic energy flux (Umov-Poynting vector) that is directed along it. The optical-mechanical analogy can be applied here. Let us consider the transmission of energy through the layer interface, namely, a change in the direction of the energy density vector. Let us introduce the transmission coefficient T and the energy reflections at the layer interface. Let P 1 be the power that is supplied with the crack from Layer 1 to the layer interface, P 2 -the power, passing into Layer 2, and P 11 -the power, which did not pass into Layer 2 through the layer interface, and it is localized at the layer interface, then Let assume that the directing unit vector of power density e 1 forms the angle Ψ with the axis z in Layer 1, and the vector e 2 forms the angle θ with the axis z in Layer 2.
It is clear that |e 1 | = |e 1 | = 1. A Snel van Royen formula (used in optics and acoustics) can be written, depending on the ratio of the physical and mechanical characteristics of the layer.
where G i is the shear modulus of the ith layer material and ρ i is the density of the ith layer material. In general, the angle θ determines the initial direction of the crack in Layer 2. A special case corresponds to the scenario of the fracture of the layer interface (delamination). The situation takes place when Ψ = Ψ c , a critical angle.
In this case, the energy is localized around the interface between the layers.
In the case of a crack, close to a vertical one, Ψ = θ = π 2 .
T ≈ 4n 2 n 1 When the crack tip in Layer 1 approaches the interface layer at the angle of Ψ = Ψ c , close to the critical one, then In general case, the Umov-Poynting vector is described for the elastic layers through the following formula: where σ ij is the stress tensor, u i are the components of the displacement velocity vector, and P j are the components of the vector P.
The angle under which the crack from Medium 1 enters the interface with Medium 2 is not known, and it can be assumed that the angle of incidence is random. Accordingly, the angle under which the crack enters Medium 2 is also a random angle, due to the connection between these angles under the Snel van Royen law. Therefore, the process of the crack passing from Layer 1 into Layer 2 will be characterized by the conditional probability density f (t 2 , θ 2 |t 1 , θ 1 ), which satisfies the equation, as follows: The following expression should be taken as an initial condition: where δ(θ) is the Dirac function.
With the initial condition (53), the solution of the Equation (52) is presented as: where P n (cos θ) are Legendre polynomials. Using the probability densities (54), it is possible to solve the problem of determining the direction of crack propagation in Medium 2 near the interface with Medium 1 (Figure 4).
the angle under which the crack enters Medium 2 is also a random angle, due to the connection between these angles under the Snel van Royen law. Therefore, the process of the crack passing from Layer 1 into Layer 2 will be characterized by the conditional probability density ( 2 , 2 | 1 , 1 ), which satisfies the equation, as follows: ( 1 , | 0 , 0 ) = 2 sin (sin ) The following expression should be taken as an initial condition: (53) where ( ) is the Dirac function.
With the initial condition (53), the solution of the Equation (52) is presented as: where (cos ) are Legendre polynomials. Using the probability densities (54), it is possible to solve the problem of determining the direction of crack propagation in Medium 2 near the interface with Medium 1 (Figure  4).
In a general case, the problem is formulated, as follows: to find the probability that the angle θ is located between the values of θ + and θ − at the time point t 1 : The following formula calculates the probability that a crack will pass from Layer 1 into Layer 2 almost vertically: In a general case, the problem is formulated, as follows: to find the probability that the angle θ is located between the values of θ + and θ − at the time point t 1 : The following formula calculates the probability that a crack will pass from Layer 1 into Layer 2 almost vertically: where ξ is a small quantity. In case, when a crack from Layer 1 approaches the interface at a small angle and, in fact, propagates along the interface between the layers, the probability is calculated by the formula: Of interest is the most probable angle of the direction of crack propagation from Medium 1 into Medium 2 at the time point t = t 1 , to be found by the formula: The necessary condition for the minimum probability (58) has the form, from which the most probable θ* can be found: The sufficient condition for the minimum is as follows: In the case a crack stops at the boundary at the time point t 0 , and then at the random time point t* it passes from Layer 1 into Layer 2, it is possible to formulate the problem of finding the most probable time for the crack to pass from Layer 1 into Layer 2 and the most probable angle at this time moment.
In this case, a two-dimensional probability distribution function can be introduced through the formula, as follows: The necessary condition is The most probable θ* and t 1 are to be found from (62). The sufficient condition: the matrix of should be positively definite. Let us consider the expression for B in (54). B can be found in the form of: Thus, through B in (54), the probabilities depend on the physical and mechanical characteristics of Layer 1 (through n 1 ) and Layer 2 (through n 2 ) from the angle Ψ, under which the crack tip in Layer 1 approaches the interface with Layer 2.

Discussion
On the qualitative level, the proposed model can be confirmed by the observation of the crack formation process in multilayer and nanolayer coatings during the action of stochastic loading. In particular, the cutting process is one of such processes [59]. Under the influence of loads (cutting forces) that arise during the process of cutting, microcracks are formed in the coating structure. In [2,3,21,35,36,52,53], it is found that the multilayer, in particular, nanolayer structure of the coating significantly affects the crack propagation. While, within one nanolayer, a crack usually propagates without changing its direction (in most cases, perpendicular to the nanolayer boundaries), then, while passing an interface between two nanolayers, a crack can either completely stop (Figure 5a) or transform into delamination between nanolayers (Figure 5b-d). Such a transformation can stop the crack propagation (Figure 5b) or form a step (Figure 5c,d). At the same time, the delamination can re-transform into a transverse crack (Figure 5c) or form a transverse crack branching from the delamination (Figure 5d). Under other conditions, a transverse crack can cut through the nanolayer structure without forming delaminations and slightly changing its direction at the nanolayer interface (Figure 5e), or with forming a straight transverse crack, which cuts through the coating as a whole (Figure 5f).
On the qualitative level, the proposed model can be confirmed by the observation of the crack formation process in multilayer and nanolayer coatings during the action of stochastic loading. In particular, the cutting process is one of such processes [59]. Under the influence of loads (cutting forces) that arise during the process of cutting, microcracks are formed in the coating structure. In [2,3,21,35,36,52,53], it is found that the multilayer, in particular, nanolayer structure of the coating significantly affects the crack propagation. While, within one nanolayer, a crack usually propagates without changing its direction (in most cases, perpendicular to the nanolayer boundaries), then, while passing an interface between two nanolayers, a crack can either completely stop (Figure 5a) or transform into delamination between nanolayers (Figure 5b-d). Such a transformation can stop the crack propagation (Figure 5b) or form a step (Figure 5c,d). At the same time, the delamination can re-transform into a transverse crack (Figure 5c) or form a transverse crack branching from the delamination (Figure 5d). Under other conditions, a transverse crack can cut through the nanolayer structure without forming delaminations and slightly changing its direction at the nanolayer interface (Figure 5e), or with forming a straight transverse crack, which cuts through the coating as a whole (Figure 5f).  Let us consider the influence of the nanolayer coating structure on the pattern cracking using the example of the Zr,Nb-(Zr,Nb)N-(Zr,Nb,Al)N coating ( Figure 6). In this coating, the nanolayers are characterized by a balanced combination of considerable hardness and ductility [57][58][59]. Figure 6 exhibits the structure of the coating that was deposited on a carbide substrate after the turning of steel. Cracks and internanolayer delaminations were formed in the coating structure as a result of the stochastic action of the cutting forces. An extensive delamination transforming into a transverse crack (Areas A, B, and C) arose in the internal area of the coating.
Area A is of particular interest (Figure 7a), which is characterized by the transformation of the internanolayer delamination first into an inclined and then into a transverse crack. All three possible ways of crack propagation through the nanolayer interface are exhibited here: the transformation into delamination (Ψ = Ψ c ), the change of propagation direction (Ψ = θ), and the retardation (P2 = 0, P1 = P11, T = 0, condition (45)). Let us consider Area A1 with the multiple transformation of the transverse crack into the delamination and reverse transformation. This process continues until the deformation energy decreases to a value, when no new transformation occurs and the crack growth stops [31]. At the same time, a decrease in power can be noticed in accordance with Equation (45). If the conditions (47) and (50) are met, then the energy is localized near the layer interface, which leads to the transformation of the crack into delamination (see Area A2, Figure 7c). A reverse transformation of the delamination into a transverse crack is possible (see Area A3, Figure 7c). In this case, condition (49) is met. It should be noted that the layer material was assumed to be isomorphic during the modeling. In reality, the layer material is neither isomorphic nor isotropic, and the crystalline structure of the layer influences the crack propagation. In particular, Area A2 (Figure 7c) depicts crystalline formations, which influence the crack growth direction. A crack can grow along the crystalline boundaries or can break a crystal (see Area A4, Figure 7e). Let us consider the influence of the nanolayer coating structure on the pattern cracking using the example of the Zr,Nb-(Zr,Nb)N-(Zr,Nb,Al)N coating ( Figure 6). In this coating, the nanolayers are characterized by a balanced combination of considerable hardness and ductility [57][58][59]. Figure 6 exhibits the structure of the coating that was deposited on a carbide substrate after the turning of steel. Cracks and internanolayer delaminations were formed in the coating structure as a result of the stochastic action of the cutting forces. An extensive delamination transforming into a transverse crack (Areas A, B, and C) arose in the internal area of the coating. Area A is of particular interest (Figure 7a), which is characterized by the transformation of the internanolayer delamination first into an inclined and then into a transverse crack. All three possible ways of crack propagation through the nanolayer interface are exhibited here: the transformation into delamination ( = ), the change of propagation direction ( ≠ ), and the retardation (P2 = 0, P1 = P11, T = 0, condition (45)). Let us consider Area A1 with the multiple transformation of the transverse crack into the delamination and reverse transformation. This process continues until the deformation energy decreases to a value, when no new transformation occurs and the crack growth stops [31]. At the same time, a decrease in power can be noticed in accordance with Equation (45). If the conditions (47) and (50) are met, then the energy is localized near the layer interface, which leads to the transformation of the crack into delamination (see Area A2, Figure 7c). A reverse transformation of the delamination into a transverse crack is possible (see Area A3, Figure 7c). In this case, condition (49) is met. It should be noted that the layer material was assumed to be isomorphic during the modeling. In reality, the layer material is neither isomorphic nor isotropic, and the crystalline structure of the layer influences the crack In some cases, the pattern of crack propagation upon reaching the layer interface can be more complex and, in particular, a crack can branch. For example, in Area B (Figure 8a), the crack branches into an internanolayer delamination (which stops about 100 nm from the branch point) and a transverse crack, which, then, in turn, transforms into an interlayer delamination. In Area C (Figure 8b), the interlayer delamination transforms into a transverse crack, which then cuts through two interlayer interfaces without changing its direction (the condition (49)), and then stops at the third interlayer interface (P2 = 0, P1 = P11, T = 0, the condition (45)).
The area of crack retardation (Area D, Figure 9a) is also of certain interest. It can be noticed how the crack retards gradually, moves along the crystal boundaries, and then stops, resting against the crystal boundary. At this moment, the crack propagation energy is insufficient to break the crystal or change the crack growth direction, and the crack propagation stops. propagation. In particular, Area A2 (Figure 7c) depicts crystalline formations, which influence the crack growth direction. A crack can grow along the crystalline boundaries or can break a crystal (see Area A4, Figure 7e). In some cases, the pattern of crack propagation upon reaching the layer interface can be more complex and, in particular, a crack can branch. For example, in Area B (Figure  8a), the crack branches into an internanolayer delamination (which stops about 100 nm from the branch point) and a transverse crack, which, then, in turn, transforms into an interlayer delamination. In Area C (Figure 8b), the interlayer delamination transforms into a transverse crack, which then cuts through two interlayer interfaces without changing its direction (the condition (49)), and then stops at the third interlayer interface (P2 = 0, P1 = P11, T = 0, the condition (45)). The area of crack retardation (Area D, Figure 9a) is also of certain interest. It can be noticed how the crack retards gradually, moves along the crystal boundaries, and then stops, resting against the crystal boundary. At this moment, the crack propagation energy is insufficient to break the crystal or change the crack growth direction, and the crack propagation stops.

Conclusions
Thus, it can be noticed that, with a certain combination of properties of alternating nanolayers, it is possible to decelerate the crack development and, thus, increase the crack resistance and resistance to brittle fracture of the coating. The possibility of finding the most probable angle of the direction of the crack propagation from Medium 1 into Medium 2 while using the proposed technique makes it possible to predict the total resistance of a multilayer coating to brittle fracture, since the nanolayer coating usually consists of a repeating sequence of nanolayers that are identical in composition and thickness. It should Similar patterns of crack propagation can be noticed in another area of the coating under consideration (Area E, Figure 9b-d). In particular, there is a step-by-step transformation of the crack growth direction under the mechanism of "delamination-transverse crack cutting through a nanolayer-delamination" (Figure 9b). The closer examination reveals the influence of the crystalline structure of the layers on the crack growth direction. With sufficient growth energy, the crack continues to grow along the intercrystalline interfaces and then cuts through the internanolayer interfaces (Area E1, Figure 9c) and, when the growth energy decreases, the crack stops, resting against the crystal at the nanolayer interface (Area E2, Figure 9d).

Conclusions
Thus, it can be noticed that, with a certain combination of properties of alternating nanolayers, it is possible to decelerate the crack development and, thus, increase the crack resistance and resistance to brittle fracture of the coating. The possibility of finding the most probable angle of the direction of the crack propagation from Medium 1 into Medium 2 while using the proposed technique makes it possible to predict the total resistance of a multilayer coating to brittle fracture, since the nanolayer coating usually consists of a repeating sequence of nanolayers that are identical in composition and thickness. It should be noted that the proposed model is dedicated to specific material of coating-solid ceramic materials based on nitrates, carbides, oxides, and carbonitrides of metals.
Thus, the proposed model includes the following stages of predicting: • the determination of the crack propagation velocity in a layer by calculating the conditional probability for location of a crack tip within the interval of (l, l + ∆l) at the time point t, if at the time point t 0 , the crack tip stays at the point l 0 . In this case, the theory of Markov processes is used in order to describe the crack formation process. This model also takes the surface roughness of the coating layer under consideration into account.

•
Predicting the behavior of a crack when it passes through the layer interface based on predicting changes in the direction of the energy density vector. The presented modeling uses the optical-mechanical analogy, and the problem of determining the direction of crack propagation in Medium 2 near the interface with Medium 1 is to be solved by determining the probability densities. The probabilities that are to be determined depend on the physical and mechanical characteristics of Layers 1 and 2 and the angle Ψ, at which the crack tip in Layer 1 approaches the interface with Layer 2. • It can be noticed that, in addition to the nanolayer structure, the pattern of crack propagation can also be affected by the crystalline structure of the coating. With a decrease in the deformation energy, the intercrystalline interfaces have a greater influence on the crack growth direction, and the crack can stop, resting against a crystal boundary. Thus, during the further modeling, it is also important to take the influence of the crystalline structure of the nanolayers into account.
In particular, as a possible direction for the further development of this model, the factors with a significant influence on crack formation in a multilayer coating can be considered, as follows: • the layer interface is gradient, which is, the transition from one layer to another occurs smoothly with a gradual increase in the percentage composition of some elements and a decrease in others [21,35,36]; • a three-dimensional model more accurately reflects the specific features of the process, although it requires much more computational capabilities; • the temperature factor affects both the mechanical properties of the layer materials and the state of the layer interfaces (in particular, due to the difference in the coefficients of thermal expansion of conjugate layers); and, • internal and residual stresses (depending, in particular, on the total coating thickness and the thicknesses of nanolayers [59]) have a noticeable influence on the crack formation. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.

Conflicts of Interest:
The authors declare no conflict of interest.
Abbreviations ω oscillation frequency in a layer under the action of cyclic loads G shear modulus E Young's modulus ν Poisson's ratio σ 0 yield strength