Modeling of the transition from transgranular to intergranular fracture at elevated temperatures in EI698 nickel alloy

Nickel alloy EI698 is utilized in the fabrication of gas turbine blades, disk components, deflectors, and other engine parts that operate at temperatures reaching up to ${800}^{\circ }\mathrm{C}$${800}^{\circ }\mathrm{C}$800^(@)C800^{\circ} \mathrm{C}${800}^{\circ }\mathrm{C}$. This alloy belongs to the family of heat- and corrosion-resistant nickel-based alloys. While it lacks a direct analogue in the international market, its composition and mechanical properties are comparable to UNS07750 (USA) and NCF750 (Japan). Components manufactured from such alloys typically function under conditions of significant variable mechanical and thermal loads. Understanding the mechanisms of damage accumulation in nickel-based alloys is essential for predicting the durability of highly loaded structural elements, where damage tolerance design prevails. The failure resistance of a material and its capacity for plastic deformation are fundamentally based on dislocation transfer mechanisms. Modern theories for fracture prediction in nickel-based alloys comprise complex compilations of models that account for microstructure, temperature, loading conditions, fatigue effects, creep, and aggressive environments $[1-3]$$[1-3]$[1-3][1-3]$[1-3]$. Studying the deformation behavior of nickel-based alloys, researchers found that as the temperature increases, the failure mechanisms changes from predominantly transgranular to intergranular cracking [4-6]. Researchers attribute this to changes in the mechanism of dislocation transfer across the grain boundary. The molecular dynamics methods can be considered as most detailed and complete [7]. However, at this stage of development, the problems solved by molecular dynamics methods relate to problems described by units of angstroms and picoseconds. Multi-level material models show good results in damage accumulation processes modeling [8,9]. These models allow to explicitly describing fracture at different scale levels, but they are too resource-intensive for solving problems of real structures. The reason is that there is a great variety of scenarios describing changes in material response at the micro level, which leads to a wide range of parameters at the macro scale.
Due to the different orientations of crystal lattices in the nucleation sites at grain boundaries, dislocations and discontinuities inevitably arise. The width of these boundaries is comparable to the size of atoms. Therefore, direct modeling of the grain boundary is not feasible using the finite element method. The phase-field theory belongs to nonlocal damage models and allows describing crack behavior at various scales. The considered scale in such models is defined by a characteristic length within which all processes are approximated by continuous and differentiable functions. Thus, phase-field models can be classified as multi-scale models, where the characteristic length determines the energy levels we operate with. On the other side, any modeling of discretized problems is always limited by computational resources. For scale depended damage accumulation modeling, the phase field fracture method has recently gained wide popularity. This theory is based on the Griffiths balance with the addition of a non-local region of damaged material [10]. It is a powerful tool for predicting intricate fracture behaviors and significantly expands the applicability of the finite element method for durability prediction of structural elements. In this paper, an efficient computational method for modeling the transition from transgranular to intergranular fracture mechanisms based on phase field fracture theory is discussed.

Operating conditions of heat-resistant nickel alloys determine the list of effects that must be taken into account for numerical modeling. For an acceptable result from the point of view of engineering calculations, this list must take into account such effects as isotropic and kinematic hardening, creep and thermal expansion of the material. Well-known continuum mechanical formalisms are used to describe the generalized model. In this study the analytical representation of a nonlinear material involves a yield surface. Inside the yield surface the material obeys Hooke's law. The surface can expand (isotropic hardening) and shift relative to the center (kinematic hardening) in the space of principal stresses. To describe the yield surface radius, a three-parameter exponential equation proposed by Voice is used. The shift of the yield surface center description is based on the Chaboche approaches. The calculations take into account the thermal expansion and viscous response of the material in the nonlinear region. The viscous response is described by the exponential creep law with strain hardening.
The model parameters characterizing the mechanical properties of the material were determined based on the stress-strain curve obtained during monotonic static loading and the analysis of hysteresis loops obtained during harmonic cyclic loading over a wide range of temperatures ( $20,370,400,450,550,650,700$$20,370,400,450,550,650,700$20,370,400,450,550,650,70020,370,400,450,550,650,700$20,370,400,450,550,650,700$ and ${750}^{\circ }\mathrm{C}$${750}^{\circ }\mathrm{C}$750^(@)C750^{\circ} \mathrm{C}${750}^{\circ }\mathrm{C}$ for static loading). The uniaxial tension tests were carried out in accordance with the requirements and limitations of ASTM standards E8 and E21. It should be noted that the characteristic parameters of the model are not directly related to the classical mechanical properties, except for the material's elastic modulus, which is determined from the linear portion of the true stress-strain curve in accordance with ASTM E111. The yield stress here is the elastic limit defined as $10\mathrm{\%}$$10\mathrm{\%}$10%10 \%$10\mathrm{\%}$ deviation of the tangent. The evolution of back stresses and kinematic hardening parameters were determined from analysis of hysteresis loops obtained from low-cycle fatigue tests conducted in accordance with the requirements of ASTM E606 at temperatures 23, 450 and ${650}^{\circ }\mathrm{C}$${650}^{\circ }\mathrm{C}$650^(@)C650^{\circ} \mathrm{C}${650}^{\circ }\mathrm{C}$ and strain range $0.4,0.5,0.6,0.7,0.8$$0.4,0.5,0.6,0.7,0.8$0.4,0.5,0.6,0.7,0.80.4,0.5,0.6,0.7,0.8$0.4,0.5,0.6,0.7,0.8$ and $1\mathrm{\%}$$1\mathrm{\%}$1%1 \%$1\mathrm{\%}$ for fully symmetric tension-compression cycle. At the next step, the isotropic hardening parameters were determined by maximizing the correspondence of the true stress-strain diagrams after subtracting the back stresses. The creep model parameters were determined based on the results of tests conducted in accordance with ASTM E139 at temperatures of $400,550,650,700$$400,550,650,700$400,550,650,700400,550,650,700$400,550,650,700$, and ${750}^{\circ }\mathrm{C}$${750}^{\circ }\mathrm{C}$750^(@)C750^{\circ} \mathrm{C}${750}^{\circ }\mathrm{C}$ and applied stress levels ranging from 200 to 830 MPa . The values of the coefficients of thermal expansion were provided by the material supplier. The ranges of parameter variations in the conducted experimental studies define the applicability domain of the model presented below.
As a result of testing smooth cylindrical specimens, all required material parameters were obtained. The dependences of Young's modulus and elastic limit on temperature were specified in finite element calculations by a piecewise linear function (Fig.1). For such effects as isotropic and kinematic hardening, creep, and thermal expansion, the most common models were used. The equations and obtained parameters of the used models are entered into the Tab. 1.

This paper does not focus on the yield function detailed description, since it can be changed by choose another set of equations presented in the final element software. For example, the exponential law of isotropic hardening can be replaced by a power law with the appropriate correction of the parameters. For numerical calculations, only models available in ANSYS for describing nonlinear material effects of continuum mechanics were used. A detailed description of all possible combinations of nonlinear effects is presented in the nonlinear materials section of the ANSYS help documentation [11]. The method presented below based on the new finite element will be valid for all allowed.

Voronoi tessellation is the most commonly used tool for modeling the structural heterogeneity of crystalline materials (Fig. 2a). The average grain size, determined from the analysis of fractographic studies of the fracture surfaces of specimens made from the material under consideration, is taken as 100 microns. The width of the grain boundary region was determined by the characteristic length of the phase field model. This parameter also determines the mesh size of the finite element model. Convergence of the results is achieved provided that there are 8 elements per characteristic length. Relatively high values of the ratio of grain size to grain boundary size (0.05-0.1) are used to reduce the number of elements of the finite element model (Fig. 2b). Parameters characterizing the mechanical properties of the material for the intergranular and transgranular space are the same for models of continuum mechanics. The difference in the material behavior will be regulated only by the energy required to crack propagation in phase field fracture model.

Thus, in finite element modeling, the material and element type for both the grain and the grain boundary are chosen to be the same, with differences specified by different sets of element real constants.

The phase field fracture approach is implemented using ANSYS user programmable features. A new custom element was created that implements phase field fracture theory for 2D and 3D problems. According to the original phase field fracture theory [12,13] the balance of a potential energy of a cracked body $\mathrm{\Psi }$$\mathrm{\Psi }$Psi\Psi$\mathrm{\Psi }$ is presented as the sum of the stored bulk energy ${\mathrm{\Psi }}_b$${\mathrm{\Psi }}_b$Psi_(b)\Psi_{b}${\mathrm{\Psi }}_b$ and the fracture energy dissipated by the formation of a new crack surface ${\mathrm{\Psi }}_f$${\mathrm{\Psi }}_f$Psi_(f)\Psi_{f}${\mathrm{\Psi }}_f$ :

where $G_c$$G_c$G_(c)G_{c}$G_c$ - critical value of the energy release rate, ${\psi }_s$${\psi }_s$psi_(s)\psi_{s}${\psi }_s$ - is the stored strain energy density, the sum of the elastic, plastic and creep energy that is not dissipated into heat, $\phi$$\phi$varphi\varphi$\phi$ - the phase field variable. The crack surface density function in perpendicular to the crack surface direction $\gamma (\phi ,\mathrm{\nabla }\phi )$$\gamma (\phi ,\mathrm{\nabla }\phi )$gamma(varphi,grad varphi)\gamma(\varphi, \nabla \varphi)$\gamma (\phi ,\mathrm{\nabla }\phi )$ is presented in the following form [14]:

where $l$$l$ll$l$ - is the length scale parameter.
The minimization problem of $\mathrm{\Psi }$$\mathrm{\Psi }$Psi\Psi$\mathrm{\Psi }$ governs crack initiation and propagation. Coupled balances of the external ( $\partial W_{\text{ext}}$$\partial W_{\text{ext }}$delW_("ext ")\partial W_{\text {ext }}$\partial W_{\text{ext }}$ ) and internal ( $\partial W_{\text{int}}$$\partial W_{\text{int }}$delW_("int ")\partial W_{\text {int }}$\partial W_{\text{int }}$ ) virtual works must be satisfied:

The external mechanical loading is defined by the variation of the external work increment as:

where $\mathbf{u}$$\mathbf{u}$u\mathbf{u}$\mathbf{u}$-displacement field, $\mathbf{b}$$\mathbf{b}$b\mathbf{b}$\mathbf{b}$-is a prescribed body force field per unit volume, while $\mathbf{h}$$\mathbf{h}$h\mathbf{h}$\mathbf{h}$ - is a boundary traction field per unit area.
The general energy functional (1) is quadratic and convex for $u$$u$uu$u$ and $\phi$$\phi$varphi\varphi$\phi$ separately. For a fixed value of $\phi$$\phi$varphi\varphi$\phi$, the functional $\mathrm{\Psi }(u$$\mathrm{\Psi }(u$Psi(u\Psi(u$\mathrm{\Psi }(u$,$)canbeefficientlyminimizedbysolvingalinearsystemofequations.Forfixedu$$)canbeefficientlyminimizedbysolvingalinearsystemofequations.Forfixedu$)canbeefficientlyminimizedbysolvingalinearsystemofequations.Forfixedu) can be efficiently minimized by solving a linear system of equations. For fixed u$)canbeefficientlyminimizedbysolvingalinearsystemofequations.Forfixedu$, the functional $\mathrm{\Psi }(,\phi )$$\mathrm{\Psi }(,\phi )$Psi(,varphi)\Psi(, \varphi)$\mathrm{\Psi }(,\phi )$ is minimized in a similar manner. Thus, the numerical implementation of the phase field model is relatively simple and can be done using a robust algorithm that minimizes each field alternately. The fracture energy balance equation for finite element implementation in simplest elastic case is [13,14]:

where ${\psi }_0(\epsilon )$${\psi }_0(\epsilon )$psi_(0)(epsi)\psi_{0}(\varepsilon)${\psi }_0(\epsilon )$ - the elastic strain energy density, $g(\phi )=(1-\phi {)}^2+k$$g(\phi )=(1-\phi {)}^2+k$g(varphi)=(1-varphi)^(2)+kg(\varphi)=(1-\varphi)^{2}+k$g(\phi )=(1-\phi {)}^2+k$ - material stiffness degradation function, $k$$k$kk$k$ - is a parameter to prevent numerical instabilities due to zero stiffness. For numerical calculations presented below $k={10}^{-7}$$k={10}^{-7}$k=10^(-7)k=10^{-7}$k={10}^{-7}$. More detailed information about the numerical implementation methods and their verification in ANSYS finite element software can be found in [14].
In this study, the both sides of the balance equation are modified. The right side is expanded to account for elastic ( ${\mathit{\psi }}_{el}$${\mathit{\psi }}_{el}$psi_(el)\boldsymbol{\psi}_{e l}${\mathit{\psi }}_{el}$ ), plastic ( ${\psi }_{pl}$${\psi }_{pl}$psi_(pl)\psi_{p l}${\psi }_{pl}$ ) and creep ( ${\psi }_{cr}$${\psi }_{cr}$psi_(cr)\psi_{c r}${\psi }_{cr}$ ) parts of the total strain energy density. On the left side it is assumed that the critical energy release rate $G_c$$G_c$G_(c)G_{c}$G_c$ depends on temperature $T_K$$T_K$T_(K)T_{K}$T_K$. Finally, the fracture energy balance equation is written in the following form:

where $h(\phi )$$h(\phi )$h(varphi)h(\varphi)$h(\phi )$ - plastic degradation function. This function can take in to account the energy dissipation as heat during plastic deformation. The necessity of introducing this function is due to the fact that the phase-field model in the presented formulation does not distinguish between plastic deformation along slip lines and the breaking of interatomic bonds, as it operates solely on the values of the local strain energy density. This, in turn, leads to incorrect crack propagation trajectories when most of the energy is consumed in the formation of new slip bands. The plastic degradation function follows the same quadratic form as the stiffness degradation function when dissipation is not considered and it is represented for consistency in the following form [15]:

where ${\beta }_p$${\beta }_p$beta_(p)\beta_{p}${\beta }_p$ - weight factor for the plastic contribution to damage. For ${\beta }_p=1$${\beta }_p=1$beta_(p)=1\beta_{p}=1${\beta }_p=1$ all plastic work going to new crack formation and $h(\phi )=g(\phi )$$h(\phi )=g(\phi )$h(varphi)=g(varphi)h(\varphi)=g(\varphi)$h(\phi )=g(\phi )$, in case ${\beta }_p=0$${\beta }_p=0$beta_(p)=0\beta_{p}=0${\beta }_p=0$ all plastic work dissipated as heat.

The interaction between the developed element implementing phase-field fracture and the nonlinear material can be schematically represented by the simplified algorithm shown in Fig. 3. Proposed realization of user defined element allows for the use of any nonlinear material available in ANSYS finite element software, since the interaction between the element and the material is based on the function calls of the respective libraries.

It is well known that crack growth resistance characteristics are temperature dependent. The methods to introduce the dependence of critical energy release rate $G_c$$G_c$G_(c)G_{c}$G_c$ from temperature in phase field fracture models is discussed in [16-18]. For the nickel alloy under consideration, these dependencies were experimentally obtained in two ways. In the first case, the method proposed in the work [6] was used. Moreover, the same material was considered in the mentioned work, therefore, in this study there will be many references to it. In the second case, a correlation between the critical value of the strain energy density and critical energy release rate is assumed. Critical strain energy density here is the area under the stress-strain curve in true stress-strain coordinates.
According to [6] the change in the work $\mathrm{\Delta }W$$\mathrm{\Delta }W$Delta W\Delta W$\mathrm{\Delta }W$ in each cycle of deformation can be calculated as the area of the hysteresis loop created in load ( $P$$P$PP$P$ ) and crack opening displacement ( $u_{COD}$$u_{COD}$u_(COD)u_{C O D}$u_{COD}$ ) coordinates (Fig.4):

For experimental values of $P$$P$PP$P$ and $u_{COD}$$u_{COD}$u_(COD)u_{C O D}$u_{COD}$ the Gauss equation was used to obtain the area of the each hysteresis loop:

where $P_i$$P_i$P_(i)P_{i}$P_i$ and $u_i$$u_i$u_(i)u_{i}$u_i$ experimentally obtained from tests of the precracked specimens under cycle loading coordinates of hysteresis loop points, $n$$n$nn$n$ - points count in the cycle $N$$N$NN$N$. Finally, the critical energy release rate for compact tension specimen calculated from:

where $\mathrm{\Delta }{a}_N=a_N-a_{N-1}$$\mathrm{\Delta }{a}_N=a_N-a_{N-1}$Deltaa_(N)=a_(N)-a_(N-1)\Delta a_{N}=a_{N}-a_{N-1}$\mathrm{\Delta }{a}_N=a_N-a_{N-1}$ is the crack length increment, $B$$B$BB$B$ - specimen thickness.
In the second case there is a direct correlation between critical fracture energy release rate $G_c$$G_c$G_(c)G_{c}$G_c$ and critical strain energy density ${\psi }_c$${\psi }_c$psi_(c)\psi_{c}${\psi }_c$ for initial conditions $\phi =0$$\phi =0$varphi=0\varphi=0$\phi =0$. For true stresses and true strains on the strain energy determination, the accumulation of the phase field before the initialization of the macrocrack can be neglected:

Considering the dimensions in (11) determines the form of the function $f\left({\mathit{\psi }}_c\right)$$f\left({\mathit{\psi }}_c\right)$f(psi_(c))f\left(\boldsymbol{\psi}_{c}\right)$f\left({\mathit{\psi }}_c\right)$. It is necessary to bring the strain energy density to a dimensional strain. Then, taking into account the extensometer base in tensile tests of smooth specimens $L_0$$L_0$L_(0)L_{0}$L_0$ and specimen thickness $B(10)$$B(10)$B(10)B(10)$B(10)$ :

The energy release rate during crack growth depends on the presence of a stress concentrator. However, in the context of this study, only the critical energy release rate is considered, which is treated as a material property dependent on temperature. At the same time, the critical value of the strain energy density is also a temperature-dependent material property. In Fig. 5, the red circles show the values of the critical energy release rate determined by (10) for the first cycles of harmonic loading at different temperatures [6]. The black rhombs correspond to the critical release rate determined from tests of smooth cylindrical specimens through the strain energy density (12) at the same temperatures. The obtained results showed that the qualitative distributions of these parameters coincide and, when normalized to a common scale, can be approximated by a single curve. A simple power law was used to approximate these results:

where $T_K$$T_K$T_(K)T_{K}$T_K$ is a temperature in Kelvin scale.

The crystallization process leads to the fact that the dislocation density is maximum at the grain boundary. Moreover, with increasing temperature, the mechanism of dislocation transfer across the grain boundary slows down, leading to an even greater accumulation of dislocations at the grain boundary during deformation. Based on this, it is assumed that the value of the critical energy release rate for the grain boundary will be significantly lower in contrast to grain body. Parametric studies conducted by the author showed, when a crack going from one grain to other thru grain boundary the additional energy is needed to transfer the crack growth process into the grain boundary field. Thus, in order for the crack to always bypass the grains and grow only in the intergranular space, it is necessary for the critical value of energy release rate to be at least 2 times lower than that for the grain body for monotonic static loading. This is well aligned with the fact that even with a concentration of dislocations at the grain boundary, the crack does not always propagate into the intergranular
space. According to fractographic studies of heat-resistant nickel alloys [ $5,6,19$$5,6,19$5,6,195,6,19$5,6,19$ ], a transgranular fracture mechanism is observed at room temperatures for all types of cyclic loading. For these conditions, it is reasonable to neglect the difference in critical energy release rate and consider the material as isotropic. At high temperatures, the material predominantly exhibits an intergranular fracture mechanism. In this regard, the following dependence for the critical energy release rate is proposed by analogy with (13):

The coefficients in Eqn. (14) were chosen so that the transition temperature from intragranular to intergranular fracture mechanism would be around ${400}^{\circ }\mathrm{C}$${400}^{\circ }\mathrm{C}$400^(@)C400^{\circ} \mathrm{C}${400}^{\circ }\mathrm{C}$. The relative positions of the approximating functions for intergranular and transgranular space are shown in Fig. 6. Finally, to minimize the deviations of the calculated values from the experimental ones, a dependence of the characteristic length of the phase field on temperature has also been introduced:

As previously noted, all the presented equations types are included in the custom element for the ANSYS software package. The source code and documentation of the created finite element can be found here [20] as an open source project.

When definition the parameters of the generalized model of continuum mechanics, the conformity between the experimental uniaxial stress-strain diagrams and the numerical results was checked first. Fig. 7 shows the curves corresponding to the two extreme points of the temperature range under study. The analysis of the modeling accuracy for tests under monotonic static loading showed a root mean square deviation of 32.79 MPa . The general trends in the change of the material behavior with increasing temperature are in good agreement with the changes of the modulus of elasticity and the elastic limit shown in Fig. 1. In this regard, intermediate results are not presented here and below in order to focus on the main effect.

The finite element modeling of the smooth cylindrical specimen in a 2D axisymmetric configuration made from the alloy under study is presented on Fig. 8. The results demonstrate that by specifying the overall mechanical properties of the grain and the intergranular boundary, the desired difference in material behavior can be achieved solely by manipulating the critical value of the energy release rate. For the above parameters at room temperature (Fig. 8a) the location of crack initiation is the grain body, at the same time, for high temperature (Fig. 8b), location of crack initiation is the grain border. The presented results do not account for the expenditure of plastic energy on heating the material, and crack growth occurs in a plane that coincides with the plastic slip planes. Parametric studies have shown that, in order for the macro crack growth plane to coincide with the experimental one, it is usually necessary to dissipate about 90 percent of the
plastic energy as heat. The value of the coefficient ${\beta }_p=0.1$${\beta }_p=0.1$beta_(p)=0.1\beta_{p}=0.1${\beta }_p=0.1$ was chosen according to the recommendations presented in [15].

Fig. 9 shows the calculated crack trajectories in precracked compact tension specimen under conditions where all plastic work contributes to the formation of new interfaces (Fig. 9a) and under conditions where 90 percent of the energy is dissipated as heat (Fig. 9b). One of the important aspects of the obtained results is that the process of transitioning from intragranular failure to intergranular failure can be absent with a significant difference in the values of the critical energy release rate. For transition temperatures (Fig. 9c), the difference in absolute values of the critical energy release rate was about $35\mathrm{\%}$$35\mathrm{\%}$35%35 \%$35\mathrm{\%}$, yet a purely intragranular failure mechanism can still be observed. This allows us to conclude that the investigated method takes into account the fact that additional energy is required to change the crack trajectory from its initial path. At the same time, the initiation of damage sites in the grain boundary region, which are not directly connected to the main crack, can be found. This also aligns well with material structure studies, where a field of microcracks can be observed around the tip of the main crack.

An efficient computational method for modeling the transition from transgranular to intergranular fracture mechanisms based on phase field fracture theory is proposed. The method is based on the implementation of a new finite element for phase fields fracture model. The proposed finite element implementation allows for accounting the decrease in the critical energy release rate with increasing temperature, as well as the dissipation of energy as heat during plastic deformation. The obtained results correlate well with previously conducted fractographic studies.

[1] Jiang, R., Song, Y.D., Reed, P.A. (2020). Fatigue crack growth mechanisms in powder metallurgy Ni-based superalloys—A review, Int. J. Fatigue, 141, 105887, DOI: https://doi.org/10.1016/j.ijfatigue.2020.105887.
[2] Pineau, A., Benzerga, A.A., Pardoen, T. (2016). Failure of metals I: Brittle and ductile fracture, Acta Mater., 107, pp. 424-483, DOI: https://doi.org/10.1016/j.actamat.2015.12.034.
[3] Pineau, A., McDowell, D.L., Busso, E.P., Antolovich, S.D. (2016). Failure of metals II: Fatigue, Acta Mater., 107, pp. 484-507, DOI: https://doi.org/10.1016/j.actamat.2015.05.050.
[4] Wang, K., Jing, H., Xu, L., Zhao, L., Han, Y., Song, K., Qi, X., Su, M. (2022). Fracture mechanism of a Ni-base alloy under high-temperature cyclic deformation: Experiments and microstructure characterization, Mater. Charact., 189, 111944, DOI: https://doi.org/10.1016/j.matchar.2022.111944.
[5] Shlyannikov, V., Sulamanidze, A., Kosov, D. (2024). Isothermal and thermo-mechanical fatigue-crack-growth analysis of XH73M nickel alloy, Theor. Appl. Fract. Mech., 129, pp. 104182, DOI: https://doi.org/10.1016/j.tafmec.2023.104182.
[6] Shlyannikov, V., Sulamanidze, A., Kosov, D. (2024). Generalization of crack growth mechanisms under isothermal and thermomechanical fatigue by COD and ERR parameters, Theor. Appl. Fract. Mech., 131, pp. 104392, DOI: https://doi.org/10.1016/j.tafmec.2024.104392.
[7] Wang, M., Incecik, A., Yang, C., Gupta, M.K., Królczyk, G., Andriukaitis, D., Li, Z. (2023). A critical review on molecular dynamics applied to structure fracture and failure analysis, Eng. Anal. Bound. Elem., 150, pp. 413-422, DOI: https://doi.org/10.1016/j.enganabound.2023.02.028.
[8] Trusov, P., Shveykin, A., Kondratev, N., Makarevich, E. (2020).Thermomechanical Processing of Steels and Alloys: Multilevel Modeling. Encyclopedia of Continuum Mechanics, Berlin, Heidelberg, Springer Berlin Heidelberg, pp. 2496-2511.
[9] Kurmoiartseva, K., Kotelnikova, N., Trusov, P. (2020).Modeling of Polycrystalline Materials Deformation with Dislocation Structure Evolution and Transition to Fracture., pp. 80-94.
[10] Wu, J.-Y., Nguyen, V.P., Nguyen, C.T., Sutula, D., Sinaie, S., Bordas, S.P.A. (2020).Phase-field modeling of fracture., pp. 1-183.
[11] ANSYS Inc. (2021). ANSYS 2021 R1 User Documentation.
[12] Simoes, M., Martínez-Pañeda, E. (2021). Phase field modelling of fracture and fatigue in Shape Memory Alloys, Comput. Methods Appl. Mech. Eng., 373, pp. 113504, DOI: https://doi.org/10.1016/j.cma.2020.113504.
[13] Martínez-Pañeda, E., Golahmar, A., Niordson, C.F. (2018). A phase field formulation for hydrogen assisted cracking, Comput. Methods Appl. Mech. Eng., 342, pp. 742-761, DOI: https://doi.org/10.1016/j.cma.2018.07.021.
[14] Kosov, D., Tumanov, A., Shlyannikov, V. (2024). ANSYS implementation of the phase field fracture approach, Frat. Ed Integrità Strutt., 18(70), pp. 133-156, DOI: https://doi.org/10.3221/IGF-ESIS.70.08.
[15] Díaz, A., Alegre, J.M., Cuesta, I., Martínez-Pañeda, E. (2025). A COMSOL framework for predicting hydrogen embrittlement-Part II: phase field fracture, Eng. Fract. Mech., 319, pp. 111008, DOI: https://doi.org/10.1016/j.engfracmech.2025.111008.
[16] Zhang, Z., Zhang, R., Sun, S., Yi, M. (2024). A thermo-mechanically coupled phase-field fatigue fracture model, Acta Mech., DOI: https://doi.org/10.1007/s00707-024-04163-y.
[17] Ruan, H., Peng, X.-L., Yang, Y., Gross, D., Xu, B.-X. (2024). Phase-field ductile fracture simulations of thermal cracking in additive manufacturing, J. Mech. Phys. Solids, 191, pp. 105756, DOI: https://doi.org/10.1016/j.jmps.2024.105756.
[18] Grossman-Ponemon, B.E., Mesgarnejad, A., Karma, A. (2022). Phase-field modeling of continuous fatigue via toughness degradation, Eng. Fract. Mech., 264, 108255. DOI: https://doi.org/10.1016/j.engfracmech.2022.108255.
[19] Shanyavskiy, A., Shlyannikov, V., Soldatenkov, A., Rubtsov, V. (2023). Micromechanics of fatigue, creep-fatigue
interaction and thermo-mechanical crack growth of XH73M nickel alloy, Procedia Struct. Integr., 43, pp. 215-220, DOI: https://doi.org/10.1016/j.prostr.2022.12.261.
[20] Tumanov, A.V. Tumanov, A.V.(2025). Github profile. Available at: https://github.com/Andrey-Fog.