Modeling subgrain structure evolution during heat treatment

The fabrication of metal and alloy products often involves multi-stage production technologies [1]. In most cases, thermomechanical processing of structural alloy workpieces utilizes pressure processing methods over a wide temperature range [1, 2]. One of the key technological stages is the final annealing process [1, 3]. Annealing involves heating the material to $0.4-0.6$$0.4-0.6$0.4-0.60.4-0.6$0.4-0.6$ homologous temperatures, holding it for a certain time and then cooling slowly [3, 4]. This process results in the relaxation of residual stresses, leads to an increase in the average size of grains and subgrains and a decrease in the density of defects (mainly, dislocations), which significantly influences the mechanical properties of a final product [3].
Two key physical processes, recovery and recrystallization, drive the reorganization of a microstructure during the annealing of deformed metals and alloys [3, 5, 6]. Because of the recovery, the density of defects reduces gradually, the defects are reorganized with formation and evolution of a low-energy subgrain structure, and finally the properties of the material are partially restored [3-5]. A further change in the material structure is realized through recrystallization. In this process, new low-defect grains are formed, grow and absorb old defective grains, resulting in the formation of a novel grain structure with low density of defects $[3,5,6]$$[3,5,6]$[3,5,6][3,5,6]$[3,5,6]$. The driving force behind these processes is the reduction of the energy stored in defects, therefore they are competing processes [3-5]. Recovery starts at lower temperatures and/or smaller pre-strain values compared to recrystallization, because less activation energy is required [3-5]. During recovery, subgrains grow and their mutual misorientation increases due to subgrain migration and subgrain coalescence [5, 7, 8]. Migration occurs as a result of the displacement of a low-angle boundary and is mainly implemented through diffusion [ 5,8 ]. Coalescence involves the merging of neighboring subgrains undergoing rotation and the disappearance of a low-angle boundary between them [5,8]. Since larger potential nuclei (or subgrains) have a growth advantage during recrystallization, subgrain size is a dominating factor for subsequent grain structure evolution [3, 5]. Thus, the recovery process prepares subgrain structure for recrystallization [3, 5].
Heat-resistant nickel alloys are widely used in diverse industrial applications such as the aerospace, petrochemical, nuclear power industries [9, 10]. Nickel alloys are utilized in the aerospace industry to manufacture gas turbine engine components [9, 11]: blades, discs, casings and liners, in the petrochemical industry to produce valve and pipeline parts, and heat exchangers, and in the nuclear power industry to make fuel elements, parts of steam generators and cooling systems. Nickel alloys are known for their unique combination of properties such as high heat-resistance, toughness and resistance to oxidation and creep, alongside with high temperature stability of the phase composition [9, 10, 12]. The enhanced strength properties of these alloys are linked to their multiphase structure, which is formed through a sequence of technological processes, including the annealing process [9, 10, 12].
In a nickel-based superalloy, the two main phases are the matrix of the $\gamma$$\gamma$gamma\gamma$\gamma$-phase with a face-centered cubic Ni - Cr lattice and the ${\gamma }^{\prime }$${\gamma }^{\prime }$gamma^(')\gamma^{\prime}${\gamma }^{\prime }$ strengthening phase based on ${\mathrm{Ni}}_3\mathrm{Al}$${\mathrm{Ni}}_3\mathrm{Al}$Ni_(3)Al\mathrm{Ni}_{3} \mathrm{Al}${\mathrm{Ni}}_3\mathrm{Al}$ [12]. Depending on the type of thermomechanical treatment, several secondary phases ${\gamma }^{\prime \prime },\delta ,\sigma ,\mu$${\gamma }^{\prime \prime },\delta ,\sigma ,\mu$gamma^(''),delta,sigma,mu\gamma^{\prime \prime}, \delta, \sigma, \mu${\gamma }^{\prime \prime },\delta ,\sigma ,\mu$, Laves phases and some others can form [12, 13]. The main strengthening effect is provided by the ${\gamma }^{\prime }$${\gamma }^{\prime }$gamma^(')\gamma^{\prime}${\gamma }^{\prime }$ and ${\gamma }^{\prime \prime }$${\gamma }^{\prime \prime }$gamma^('')\gamma^{\prime \prime}${\gamma }^{\prime \prime }$ phases [12, 13]. In addition to secondary phases, the grain and subgrain structures, in particular their sizes, strongly
affect the strength of materials [12, 14]. It is worth noting that materials with larger grains exhibit better crack and creep resistance, while the fine-grained structure of materials enhances resistance to low-cycle fatigue and yield strength [14, 15]. During thermomechanical treatment, including annealing, material structure and phase composition can be purposefully changed, which provides ample opportunities to create functional material-products [14, 15].
In order to provide effective management of recovery and recrystallization processes needed to design the material structure suitable for specific applications, it is crucial to develop mathematical models that can explicitly account for material's structure at different scale levels [11, 16]. An effective approach to conducting numerical studies is a multilevel approach with internal variables, which makes it possible to explicitly describe deformation mechanisms and structure evolution at different scales of a material [11, 16]. In the previous study, the authors offered an advanced multilevel statistical model that takes into account the subgrain structure topology for describing the coalescence process [17, 18]. The purpose of this work is to develop further and modify the available statistical model for describing the recovery process driven by the subgrain boundary migration and to apply it to a representative volume of Inconel 718 nickel-based alloy subgrains to describe the annealing process and the interaction between migration and coalescence processes.

Mathematical modeling of the evolution of a subgrain structure in polycrystalline materials makes it possible to study in detail the deformation and strain hardening mechanisms and the recrystallization processes [16]. Taking into account the hierarchy of material structure in mathematical modeling enhances the accuracy and capabilities of digital twins of materials and products and expands their applications. The subgrain structure is generally described by approaches and models that vary in their level of detail. In the most common approach, the average subgrain size is described by the macrophenomenological evolutionary relation [7, 16]. The model developed by Sandström is one of the most popular models in this class [7]. It is based on the consideration of the rearrangement of dislocations within subgrains and their boundaries during static annealing. Another way to describe grain and subgrain structures implies using mean-field models, where the change in the size of structural elements is modeled on the assumption that effective grain or subgrain are "immersed" into the matrix of the surrounding material [11]. Mean-field models cannot be used to describe such processes as coalescence because of the need to consider local interactions between structural elements [17, 18]. The most accurate, but also computationally resource-intensive models are direct ones, which explicitly take into account the topology of grains or subgrains [11].
In direct models, the geometry (close to the real geometry) of the representative volume of subgrains or grains is described using polyhedral; this is a labor-intensive and complex task [5, 16, 19]. The main mathematical and computational difficulty arises due to the non-affine reorganization of the structure driven by the occurrence of new elements and the disappearance of old ones. This leads to an increase in computational costs, which increase many times with an increase in the number of structural elements of the model, which is typical for the recrystallization process. The evolution of polyhedral as a result of recrystallization and recovery processes can be described by applying different approaches [16, 20, 21]: the Monte Carlo model technique, the cellular automata method, the phase field method, the graph theory method, etc. Because of the computational complexity, explicit structure reorganization is carried out locally in some areas; complete reorganization with preservation of the main moments of statistical sampling of a representative volume is also possible [22].
An effective tool for studying the evolution of material structure, including a subgrain one, is a multilevel approach with internal variables [11, 16]. In [17, 19, 22], the authors proposed a novel way to model the dynamic recrystallization and coalescence processes using the advanced multilevel statistical model. This method permits describing the material structure evolution in the framework of a statistical approach involving the geometric characteristics of Laguerre polyhedra. The advantage of the above-mentioned models is their computational efficiency and ability to describe in detail the material structure evolution taking into account the topology (shape, size, mutual arrangement of elements) of grains or subgrains, as well as the interaction of adjacent structural elements. Introducing internal variables that correspond to the structural elements of the model provides a physically-based description of the processes of its reconstruction. The proposed model was applied to describe the process of coalescence of subgrains during dynamic recrystallization [17, 18]. In this paper, we perform generalization to model the recovery process with regard to the migration of subgrain boundaries. The developed model is used to describe the subgrain structure evolution caused by the combined action of the recovery mechanisms (subgrain boundary migration and subgrain coalescence) during annealing.
The material structure of metals and alloys is hierarchic [11, 16], therefore, in the study, three scale levels are distinguished: macrolevel, mesolevel-I, and mesolevel-II [17, 18]. The macrolevel is associated with a representative volume of the polycrystal where the effective macroproperties of material are determined. At mesolevel-I, individual grains having a
characteristic size of $10-100\mu \mathrm{m}$$10-100\mu \mathrm{m}$10-100 mum10-100 \mu \mathrm{~m}$10-100\mu \mathrm{m}$ and separated by high-angle boundaries are considered. Each grain consists of subgrains, which have a characteristic size ranging from 0.1 to $1\mu \mathrm{m}$$1\mu \mathrm{m}$1mum1 \mu \mathrm{~m}$1\mu \mathrm{m}$ and are separated by low-angle boundaries. At mesolevel-II, a homogeneous subgrain is investigated. In the model case studied, we consider a representative volume of subgrains with polyhedral structure.
In order to model the evolution of subgrain structure during the recovery process, the structure topology should be determined before a material undergoes annealing [17]. The initial state of a polyhedral structure was modeled using the algorithm of Laguerre polyhedra implemented in the free software package Neper [23]. In the initial configuration, it was assumed that all the energy stored in defects after preliminary plastic deformation is concentrated at subgrain boundaries, forming dislocation walls [5, 24, 25]. The computational difficulties listed above can be overcome using a computationally efficient method for describing subboundary migration that shares conceptual similarities with a coalescence model [17, 18]. The migration of boundaries is modeled using a discrete approach, as a result of an elementary act, the absorbed subgrain merges completely with the growing subgrain. The description of the computational scheme of this process is given below for an arbitrarily chosen boundary (facet) between two subgrains and can be naturally generalized.
The subgrain coalescence model being applied is based on the authors' previous works [17, 18]. According to the model [17], the coalescence of a pair of adjacent subgrains is realized through the rotation $\mathbf{r}$$\mathbf{r}$r\mathbf{r}$\mathbf{r}$ of the crystal lattice of one subgrain to align with its neighbor, followed by the dissociation of their boundaries. This leads to the formation of a new subgrain that inherits the total volume and topology of the two parent subgrains. The coalescence of subgrains undergoing rotation is described by applying two criteria, energy and time. According to the energy criterion, under the rotation $\mathbf{r}$$\mathbf{r}$r\mathbf{r}$\mathbf{r}$ of the crystal lattice, the total energy of the boundaries of two subgrains $E_{sb}$$E_{sb}$E_(sb)E_{s b}$E_{sb}$ should reduce, that is, the inequality $E_{sb}\ge E_{sb}^{\mathbf{r}}$$E_{sb}\ge E_{sb}^{\mathbf{r}}$E_(sb) >= E_(sb)^(r)E_{s b} \geq E_{s b}^{\mathbf{r}}$E_{sb}\ge E_{sb}^{\mathbf{r}}$ (where $E_{sb}^{\mathbf{r}}$$E_{sb}^{\mathbf{r}}$E_(sb)^(r)E_{s b}^{\mathbf{r}}$E_{sb}^{\mathbf{r}}$ is the total energy of boundaries after rotation) holds true. According to the time criterion, the elapsed calculation time $t$$t$tt$t$ should exceed the time required for the subgrain boundary to dissociate $t_c$$t_c$t_(c)t_{c}$t_c$, that is, the inequality $t\ge t_c$$t\ge t_c$t >= t_(c)t \geq t_{c}$t\ge t_c$ should be satisfied. The physical basis of a migration model stems from the works of Humphreys F.J. [5]. The migration rate $v_{ij}$$v_{ij}$v_(ij)v_{i j}$v_{ij}$ of the boundary between the subgrains $i$$i$ii$i$ and $j$$j$jj$j$ is determined by the following relation [5]:

where $m_{\text{lag}}$$m_{\text{lag }}$m_("lag ")m_{\text {lag }}$m_{\text{lag }}$ is the mobility of a low-angle subgrain boundary, and $p_{ij}$$p_{ij}$p_(ij)p_{i j}$p_{ij}$ is the driving pressure due to the difference in the energies stored in defects. The mobility of a subgrain boundary $m_{\text{lag}}$$m_{\text{lag }}$m_("lag ")m_{\text {lag }}$m_{\text{lag }}$ is determined by the mutual misorientation angle $\theta$$\theta$theta\theta$\theta$ and manifests itself as the mobility of a high-angle boundary $m_{\text{bag}}[5]$$m_{\text{bag }}[5]$m_("bag ")[5]m_{\text {bag }}[5]$m_{\text{bag }}[5]$ :

where ${\theta }_m$${\theta }_m$theta_(m)\theta_{m}${\theta }_m$ is the maximum misorientation angle of subgrains, which corresponds to the transition from low-angle to highangle boundary. The value of $m_{\text{bag}}$$m_{\text{bag }}$m_("bag ")m_{\text {bag }}$m_{\text{bag }}$ depends on temperature and is calculated by using the following relation [5, 26]:

where $m_{hag,0}$$m_{hag,0}$m_(hag,0)m_{h a g, 0}$m_{hag,0}$ is the pre-exponential term, $Q_b$$Q_b$Q_(b)Q_{b}$Q_b$ is the activation energy of boundary diffusion, $k_B$$k_B$k_(B)k_{B}$k_B$ is the Boltzmann's constant, and $T$$T$TT$T$ is the absolute temperature.
The subgrain boundary migration is determined by a decrease in the surface energy; in multicomponent alloys, the secondphase dispersed particles have a strong impact on the process [5,26]. The pressure $p_{ij}$$p_{ij}$p_(ij)p_{i j}$p_{ij}$ on the subgrain boundary on the side of the $i$$i$ii$i$-th subgrain toward the $j$$j$jj$j$-th subgrain is calculated by applying the relation [5]:

where $\alpha$$\alpha$alpha\alpha$\alpha$ is the subgrain geometry (form) factor, $e_{sb}$$e_{sb}$e_(sb)e_{s b}$e_{sb}$ is the subgrain boundary energy, $r_i$$r_i$r_(i)r_{i}$r_i$ and $r_j$$r_j$r_(j)r_{j}$r_j$ are the radii of spheres having a volume equal to the volume of subgrains, $p_z$$p_z$p_(z)p_{z}$p_z$ is the braking pressure due to the dispersed particles of the material. If the value of $p_{ij}$$p_{ij}$p_(ij)p_{i j}$p_{ij}$ is positive, then the boundary migrates from the $i$$i$ii$i$-th subgrain into the depth of the $j$$j$jj$j$-th subgrain. To describe the braking force $p_z$$p_z$p_(z)p_{z}$p_z$ of multicomponent alloys [26], the classical Zener relationship (formulated assuming spherical particles) is used [5]:

where $f_v$$f_v$f_(v)f_{v}$f_v$ is the volume fraction of particles, and ${\varrho }_{av}$${\varrho }_{av}$ϱ_(av)\varrho_{a v}${\varrho }_{av}$ is the average radius of particles. The low-angle subgrain boundary energy $e_{sb}$$e_{sb}$e_(sb)e_{s b}$e_{sb}$ is described by the Read-Shockley relationship [5, 27]:

where $e_{sb,m}$$e_{sb,m}$e_(sb,m)e_{s b, m}$e_{sb,m}$ is the high-angle boundary energy.
Below we present a description of the concept of modeling the low-angle subgrain boundary migration for the arbitrarily chosen facet between two subgrains. For this facet, the pressure $p_{ij}$$p_{ij}$p_(ij)p_{i j}$p_{ij}$ acting on the boundary is calculated by using formula (4) at each moment in time. When $p_{ij}>0$$p_{ij}>0$p_(ij) > 0p_{i j}>0$p_{ij}>0$, the volume $\mathrm{\Delta }{v}_{ij}$$\mathrm{\Delta }{v}_{ij}$Deltav_(ij)\Delta v_{i j}$\mathrm{\Delta }{v}_{ij}$ absorbed by the growing subgrain is calculated. The calculations rely on the known geometric characteristics like the facet area $s_{ij}$$s_{ij}$s_(ij)s_{i j}$s_{ij}$ and the size of adjacent subgrains $r_i$$r_i$r_(i)r_{i}$r_i$ and $r_j$$r_j$r_(j)r_{j}$r_j$. It is assumed that the facet moves exclusively along its normal, and its area $s_{ij}$$s_{ij}$s_(ij)s_{i j}$s_{ij}$ remains unchanged (Fig. 1a). The geometric inaccuracy is corrected through the introduction of a multiplier $\beta$$\beta$beta\beta$\beta$. Thus, the volume $\mathrm{\Delta }{v}_{ij}$$\mathrm{\Delta }{v}_{ij}$Deltav_(ij)\Delta v_{i j}$\mathrm{\Delta }{v}_{ij}$ of the $j$$j$jj$j$-th subgrain absorbed by the $i$$i$ii$i$-th in a time $\mathrm{\Delta }t$$\mathrm{\Delta }t$Delta t\Delta t$\mathrm{\Delta }t$ is written as:

The migration process of the arbitrary $i$$i$ii$i$-th facet is represented schematically in Fig. 1a.
The migration variables listed above are calculated at each step of the computational procedure. The geometry of subgrains changes discretely: the $j$$j$jj$j$-th subgrain is considered absorbed when $\mathrm{\Delta }{v}_{ij}$$\mathrm{\Delta }{v}_{ij}$Deltav_(ij)\Delta v_{i j}$\mathrm{\Delta }{v}_{ij}$ exceeds its initial volume $v_j$$v_j$v_(j)v_{j}$v_j$, that is, the inequality $\mathrm{\Delta }{v}_{ij}\ge v_j$$\mathrm{\Delta }{v}_{ij}\ge v_j$Deltav_(ij) >= v_(j)\Delta v_{i j} \geq v_{j}$\mathrm{\Delta }{v}_{ij}\ge v_j$ holds true. In this context, by absorption is meant the moving boundary dissociation with appropriate structural reorganization (Fig. 1a). In the following, the described procedure is repeated for the renewed subgrain structure.

Figure 1: Schematic diagram illustrating the subgrain structure reorganization during the migration of (a) arbitrary facet and (b) facets of neighboring subgrains.

The process of migration typically occurs simultaneously at many facets of the considered area. The continuous movement of these facets inevitably gives rise to the situations where the volume of a subgrain is completely absorbed by the adjacent subgrains. A schematic representation of one of these situations is given in Fig. 1b. In the developed model, the $j$$j$jj$j$-th subgrain is assumed to be absorbed if the total volume of its absorbed part from all adjacent subgrains, $\sum _{i=1}^N\mathrm{\Delta }{v}_{ij}$$\sum _{i=1}^N \mathrm{\Delta }{v}_{ij}$sum_(i=1)^(N)Deltav_(ij)\sum_{i=1}^{N} \Delta v_{i j}$\sum _{i=1}^N\mathrm{\Delta }{v}_{ij}$ ( $N$$N$NN$N$ denotes the number of neighboring moving subboundaries), exceeds the volume of the absorbed subgrain $v_j$$v_j$v_(j)v_{j}$v_j$, that is, the inequality $\sum _{i=1}^N\mathrm{\Delta }{v}_{ij}\ge v_j$$\sum _{i=1}^N \mathrm{\Delta }{v}_{ij}\ge v_j$sum_(i=1)^(N)Deltav_(ij) >= v_(j)\sum_{i=1}^{N} \Delta v_{i j} \geq v_{j}$\sum _{i=1}^N\mathrm{\Delta }{v}_{ij}\ge v_j$ holds true. When this criterion is satisfied, the reorganization of a polyhedral structure is realized through the merging of the absorbed $j$$j$jj$j$-th subgrain and the adjacent neighbor, which contributes the greatest amount of $\mathrm{\Delta }{v}_{ij}$$\mathrm{\Delta }{v}_{ij}$Deltav_(ij)\Delta v_{i j}$\mathrm{\Delta }{v}_{ij}$. The volumes of other subgrains, $\mathrm{\Delta }{v}_{ij}$$\mathrm{\Delta }{v}_{ij}$Deltav_(ij)\Delta v_{i j}$\mathrm{\Delta }{v}_{ij}$, except for the absorbed subgrain, are dropped from further calculations. The physical justification for this assumption is its consistency with the competitive growth of subgrains: the subgrain with the highest driving force and migration rate dominates the absorption process [5]. Thus, the subgrain migration on the initial polyhedral structure is realized discretely; the act of migration is the complete absorption of the neighboring subgrain. This description of migration is consistent with the description of coalescence from [17]; this ensures simple integration of the calculation moduli for these processes.
In the case of simultaneous local realization of both subgrain growth mechanisms - coalescence and migration - preference is given to the more intensive process, as a result of which the subgrain volume is completely absorbed more rapidly.
The discrete approach applied to modeling recovery processes is conceptually analogous to the cellular automata method [16, 20], in which microstructure evolution is realized through discrete switching events. The main advantage of the proposed approach lies in its computational efficiency and explicit description of the topology of the polyhedral structure. This makes it possible to describe the interaction between neighboring subgrains. The accuracy of the description can be further improved through additional partitioning of the polyhedra into smaller elements (clusters, voxels) [23] or through optimization-based methods for regenerating the subgrain structure [22]. An explicit description of boundary motion would require continuous reconstruction of the polyhedral geometry at each time step, which entails the use of computationally intensive methods similar to the phase-field method [21]. A detailed description of the polyhedral structure evolution lies beyond the scope of the present study.

At present, there is a noticeable lack of experimental studies focusing on the evolution of subgrain structure, especially for multicomponent alloys including nickel-based superalloys [3, 5]. This refers to difficulties in studying the influence of recrystallization and recovery processes at elevated temperatures, particularly at the subgrain level [5, 28]. Most experimental studies of subgrain structure have been performed on simple metals and alloys such as pure metals, binary and ternary alloy systems, and alloys involving a dominant matrix phase. The lack of sufficient natural data on the behavior of grain and subgrain structures in multicomponent alloys and the effect of second phases on boundary migration,
coalescence, new grain formation, dislocation movement, etc. [28] causes the complexity of the identification procedure. For this reason, the model parameters were decomposed into two sets: one where the secondary phase particles affect calculations, and another where they do not. Since the matrix of the Inconel 718 alloy is the $\gamma$$\gamma$gamma\gamma$\gamma$-phase [13,29,30], then the identification of some of the parameters from the second set was carried out using the experimental results obtained for pure nickel; these parameters are given below.
At the subgrain scale, the matrix contains particles of second phases ${\gamma }^{\prime }$${\gamma }^{\prime }$gamma^(')\gamma^{\prime}${\gamma }^{\prime }$ and ${\gamma }^{\prime \prime }$${\gamma }^{\prime \prime }$gamma^('')\gamma^{\prime \prime}${\gamma }^{\prime \prime }$, which are taken into account through the Zener force (5). The ${\gamma }^{\prime \prime }$${\gamma }^{\prime \prime }$gamma^('')\gamma^{\prime \prime}${\gamma }^{\prime \prime }$ phase particles are elongated discs with dimensions of $20-30\mathrm{nm}$$20-30\mathrm{nm}$20-30nm20-30 \mathrm{~nm}$20-30\mathrm{nm}$, and the ${\gamma }^{\prime }$${\gamma }^{\prime }$gamma^(')\gamma^{\prime}${\gamma }^{\prime }$ phase particles are spherical particles with dimensions of $10-15\mathrm{nm}$$10-15\mathrm{nm}$10-15nm10-15 \mathrm{~nm}$10-15\mathrm{nm}$ [29-31]. The formation of these particles occurs during alloy manufacturing; in the temperature range of $550-{660}^{\circ }\mathrm{C}$$550-{660}^{\circ }\mathrm{C}$550-660^(@)C550-660^{\circ} \mathrm{C}$550-{660}^{\circ }\mathrm{C}$, the alloy typically undergoes long-term aging and, in the temperature range of $700-{900}^{\circ }\mathrm{C}$$700-{900}^{\circ }\mathrm{C}$700-900^(@)C700-900^{\circ} \mathrm{C}$700-{900}^{\circ }\mathrm{C}$, short-term aging [29, 30]. At temperatures above $650-750{}^{\circ }\mathrm{C}$$650-750{}^{\circ }\mathrm{C}$650-750^(@)C650-750{ }^{\circ} \mathrm{C}$650-750{}^{\circ }\mathrm{C}$, the ${\gamma }^{\prime \prime }$${\gamma }^{\prime \prime }$gamma^('')\gamma^{\prime \prime}${\gamma }^{\prime \prime }$ phase is transformed into the $\delta$$\delta$delta\delta$\delta$ phase, which is placed at grain boundaries [28-30]. In this paper, it was assumed that the investigated alloy was subjected to standard heat treatment before annealing [30, 31]. Being in this state, the volume fractions of particles of the ${\gamma }^{\prime }$${\gamma }^{\prime }$gamma^(')\gamma^{\prime}${\gamma }^{\prime }$ and ${\gamma }^{\prime \prime }$${\gamma }^{\prime \prime }$gamma^('')\gamma^{\prime \prime}${\gamma }^{\prime \prime }$ phases were approximately equal, 0.04 and 0.13 , respectively, while coarse $\delta$$\delta$delta\delta$\delta$-phase particles are not observed [30-32]. The characteristic sizes of the ${\gamma }^{\prime }$${\gamma }^{\prime }$gamma^(')\gamma^{\prime}${\gamma }^{\prime }$ and ${\gamma }^{\prime \prime }$${\gamma }^{\prime \prime }$gamma^('')\gamma^{\prime \prime}${\gamma }^{\prime \prime }$ phase particles are significantly smaller than the average subgrain size. This justifies applying the classical Zener relationship (5) for an averaged description of their effect on migration [5]. An explicit description of individual fine particles of the ${\gamma }^{\prime }$${\gamma }^{\prime }$gamma^(')\gamma^{\prime}${\gamma }^{\prime }$ and ${\gamma }^{\prime \prime }$${\gamma }^{\prime \prime }$gamma^('')\gamma^{\prime \prime}${\gamma }^{\prime \prime }$ phases within the subgrain volume is impractical due to the significant computational costs. Thus, the average value of $p_z/e_{sb}$$p_z/e_{sb}$p_(z)//e_(sb)p_{z} / e_{s b}$p_z/e_{sb}$ in relation (5) from particles of the ${\gamma }^{\prime }$${\gamma }^{\prime }$gamma^(')\gamma^{\prime}${\gamma }^{\prime }$ and ${\gamma }^{\prime \prime }$${\gamma }^{\prime \prime }$gamma^('')\gamma^{\prime \prime}${\gamma }^{\prime \prime }$ phases was estimated as $2.1\cdot {10}^7{\mathrm{m}}^{-1}$$2.1\cdot {10}^7{\mathrm{m}}^{-1}$2.1*10^(7)m^(-1)2.1 \cdot 10^{7} \mathrm{~m}^{-1}$2.1\cdot {10}^7{\mathrm{m}}^{-1}$.
A computational annealing experiment was conducted on a representative volume of subgrains ( $N_0=10000$$N_0=10000$N_(0)=10000N_{0}=10000$N_0=10000$ ) within a single grain. The state of the material prior to annealing was assumed to be the same as that after rolling at a strain value of $\epsilon =2.25$$\epsilon =2.25$epsi=2.25\varepsilon=2.25$\epsilon =2.25$. The initial average size (diameter) of subgrains $d_{av0}=0.25\mu \mathrm{m}$$d_{av0}=0.25\mu \mathrm{m}$d_(av0)=0.25 mumd_{a v 0}=0.25 \mu \mathrm{~m}$d_{av0}=0.25\mu \mathrm{m}$ after preliminary deformation was determined based on the experimental results for pure nickel [33]. The initial sizes of subgrains are specified by a one-parameter Rayleigh distribution [34]. Analysis of the obtained data [5, 24] indicated that subgrains exhibit a predominantly equiaxed shape, therefore the computational experiment started by assuming that the range of sphericity $\psi$$\psi$psi\psi$\psi$, when uniformly distributed, was between 0.8 and 0.9 . The available experimental data on subgrain size and sphericity distribution made it possible to create a polyhedral structure using Neper [23]. The geometric correction factor $\beta$$\beta$beta\beta$\beta$, which accounts for the difference between the idealized absorbed volume and the actual shape of the absorbed region, was set to 1.0 . The simulation results and the comprehensive sensitivity analysis performed (see the results of modeling) show that the adopted estimate is an acceptable approximation and remains within the scope of the formulated problem of an approximate description of the subgrain structure.
Later, the problem of determining the mutual misorientation angle $\theta$$\theta$theta\theta$\theta$ between adjacent subgrains was solved. Using the Read-Shockley relation (6), the angle $\theta$$\theta$theta\theta$\theta$ uniquely determines the surface energy $e_{sb}$$e_{sb}$e_(sb)e_{s b}$e_{sb}$. It was assumed that the energy $e_{st}$$e_{st}$e_(st)e_{s t}$e_{st}$ stored during plastic deformation is primarily concentrated at subgrain boundaries. For rolling deformation $\epsilon =2.25$$\epsilon =2.25$epsi=2.25\varepsilon=2.25$\epsilon =2.25$, the value of this energy, experimentally found in [25] as $e_{st}^{\text{exp}}=3.3\mathrm{MJ}/{\mathrm{m}}^3$$e_{st}^{\text{exp }}=3.3\mathrm{MJ}/{\mathrm{m}}^3$e_(st)^("exp ")=3.3MJ//m^(3)e_{s t}^{\text {exp }}=3.3 \mathrm{MJ} / \mathrm{m}^{3}$e_{st}^{\text{exp }}=3.3\mathrm{MJ}/{\mathrm{m}}^3$, was compared with our calculated values. The calculations were performed on the polyhedral structure generated in this study, and the angle $\mathit{\theta }$$\mathit{\theta }$theta\boldsymbol{\theta}$\mathit{\theta }$ was assumed to be distributed under a Rayleigh law [35, 36]; the rotation axis was set in a random way through the uniform distribution over the sphere. The angle $\theta$$\theta$theta\theta$\theta$ was assigned to the sample of subgrains in a random way, and the value of $e_{st}$$e_{st}$e_(st)e_{s t}$e_{st}$ was found by integrating over all subboundaries under a Read-Shockley law (6) and attributed to the entire representative volume of subgrains $v_{sb}$$v_{sb}$v_(sb)v_{s b}$v_{sb}$ :

where $N_{sb}$$N_{sb}$N_(sb)N_{s b}$N_{sb}$ is the number of subgrain boundaries in a representative volume, and $s^{(i)}$$s^{(i)}$s^((i))s^{(i)}$s^{(i)}$ is the boundary area. Thus, the functional $f_1(\theta )$$f_1(\theta )$f_(1)(theta)f_{1}(\theta)$f_1(\theta )$ to be minimized had the following form:

Analysis of the numerical experiments carried out in this study yields the Rayleigh distribution with scale parameter ${\theta }_{av}$${\theta }_{av}$theta_(av)\theta_{a v}${\theta }_{av}$ corresponding to the average misorientation angle of ${2.29}^{\circ }$${2.29}^{\circ }$2.29^(@)2.29^{\circ}${2.29}^{\circ }$, and the value of the minimized functional was $f_1=0.001$$f_1=0.001$f_(1)=0.001f_{1}=0.001$f_1=0.001$. The final stage of identification involves determining the mobility of the high-angle boundary $m_{\text{bag}}$$m_{\text{bag }}$m_("bag ")m_{\text {bag }}$m_{\text{bag }}$ depending on temperature, relation (3). To this end, the annealing of pure nickel was modeled with a consideration of migration and coalescence processes. To determine the parameters $m_{hag,0}$$m_{hag,0}$m_(hag,0)m_{h a g, 0}$m_{hag,0}$ and $Q_b$$Q_b$Q_(b)Q_{b}$Q_b$, the minimization problem was solved so that the numerical data on the average size of subgrains best correspond to the experimental results obtained at 200 and ${340}^{\circ }\mathrm{C}$${340}^{\circ }\mathrm{C}$340^(@)C340^{\circ} \mathrm{C}${340}^{\circ }\mathrm{C}$ (Fig. 2) [33]. The corresponding objective functional was formulated as:

where $d_{av}^{\text{exp}}$$d_{av}^{\text{exp }}$d_(av)^("exp ")d_{a v}^{\text {exp }}$d_{av}^{\text{exp }}$ is the experimental average subgrain size at the annealing time $t_i^{\text{exp}},d_{av}^{\text{num}}$$t_i^{\text{exp }},d_{av}^{\text{num }}$t_(i)^("exp "),d_(av)^("num ")t_{i}^{\text {exp }}, d_{a v}^{\text {num }}$t_i^{\text{exp }},d_{av}^{\text{num }}$ is the corresponding value predicted by the model at $t_i^{\text{exp}}$$t_i^{\text{exp }}$t_(i)^("exp ")t_{i}^{\text {exp }}$t_i^{\text{exp }}$, and $N^{\text{exp}}$$N^{\text{exp }}$N^("exp ")N^{\text {exp }}$N^{\text{exp }}$ is the number of experimental points. The verification of the subgrain size in the model at 260 and ${300}^{\circ }\mathrm{C}$${300}^{\circ }\mathrm{C}$300^(@)C300^{\circ} \mathrm{C}${300}^{\circ }\mathrm{C}$ is shown in Fig. 2. The maximum mean deviation of the simulation results from the experimental data (determined according to relation (10)) amounted to $0.42\mathrm{\%}$$0.42\mathrm{\%}$0.42%0.42 \%$0.42\mathrm{\%}$ for identification and $1.15\mathrm{\%}$$1.15\mathrm{\%}$1.15%1.15 \%$1.15\mathrm{\%}$ for verification. Further computational experiments were carried out in the temperature range $200-340{}^{\circ }\mathrm{C}$$200-340{}^{\circ }\mathrm{C}$200-340^(@)C200-340{ }^{\circ} \mathrm{C}$200-340{}^{\circ }\mathrm{C}$. All identified parameters are given in Tab. 1.

(-) Experiment $T={200}^{\circ }\mathrm{C}$$T={200}^{\circ }\mathrm{C}$T=200^(@)CT=200^{\circ} \mathrm{C}$T={200}^{\circ }\mathrm{C}$

◇ Experiment $T={300}^{\circ }\mathrm{C}$$T={300}^{\circ }\mathrm{C}$T=300^(@)CT=300^{\circ} \mathrm{C}$T={300}^{\circ }\mathrm{C}$
△ Experiment $T={340}^{\circ }\mathrm{C}$$T={340}^{\circ }\mathrm{C}$T=340^(@)CT=340^{\circ} \mathrm{C}$T={340}^{\circ }\mathrm{C}$

To justify the assumption of assigning the absorbed volume to the most active subgrain, a series of test calculations was performed. To obtain an upper-bound estimate, the case of pure nickel annealing at $340{}^{\circ }\mathrm{C}$$340{}^{\circ }\mathrm{C}$340^(@)C340{ }^{\circ} \mathrm{C}$340{}^{\circ }\mathrm{C}$ was considered (Fig. 1). The total cumulative fraction of the absorbed volume attributed to non-dominant neighbors (i.e., the fraction discarded when the absorbed subgrain is assigned to the dominant neighbor) amounted to $0.44\mathrm{\%}$$0.44\mathrm{\%}$0.44%0.44 \%$0.44\mathrm{\%}$ of the entire volume considered over 40000 s . This estimate consists of two contributions: (a) $0.41\mathrm{\%}-$$0.41\mathrm{\%}-$0.41%-0.41 \%-$0.41\mathrm{\%}-$ the cumulative volume fraction of subgrains discarded in cases where subgrain merging is realized through the migration mechanism, and (b) $0.03\mathrm{\%}$$0.03\mathrm{\%}$0.03%0.03 \%$0.03\mathrm{\%}$ - the volume of subgrains accumulated by migration but discarded in cases where merging is realized through the competing coalescence mechanism. The obtained estimates confirm the adopted assumption: for both mechanisms considered, the dominant subgrain contributes the overwhelming majority to the absorbed volume. At the same time, the accuracy of describing the structural reorganization can be further improved through additional discretization of the polyhedral structure into smaller elements. In addition, the representativeness of the volume considered was verified by doubling the number of subgrains in the sample: the overall deviation in the average subgrain size between the two calculations (10000 and 20000 subgrains) amounted to $0.44\mathrm{\%}$$0.44\mathrm{\%}$0.44%0.44 \%$0.44\mathrm{\%}$. The average error of a single discrete subgrain merging act is about $4.3\mathrm{\%}$$4.3\mathrm{\%}$4.3%4.3 \%$4.3\mathrm{\%}$, with the largest contribution coming from acts at the initial stage of annealing; this value decays thereafter because the recovery process is directed toward subgrain coarsening [5, 40]. For this reason, the influence of this error on the final result is negligible.

This section includes the calculated results obtained using the developed model for the Inconel 718 alloy during annealing in the temperature range $200-{340}^{\circ }\mathrm{C}$$200-{340}^{\circ }\mathrm{C}$200-340^(@)C200-340^{\circ} \mathrm{C}$200-{340}^{\circ }\mathrm{C}$. Fig. 3 demonstrates how the average subgrain size (diameter) increases with annealing time. Higher temperatures accelerate recovery processes, and, consequently, increase the rate of subgrain growth. Analysis of the results obtained for the Inconel 718 alloy revealed a significant slowdown in the growth of subgrains (Fig. 3) in comparison with pure nickel (Fig. 2). This can be attributed to the retarding effect of the second-phase particles on the migration of subgrain boundaries. This effect is consistent with the previously described Zener's mechanism [5, 26]. It is easy to see that the obtained results (Fig. 3) are well described by the following classical relation [5, 33]:

where $d$$d$dd$d$ is the size of grains or subgrains after annealing over time $t,d_0$$t,d_0$t,d_(0)t, d_{0}$t,d_0$ is the size of grains or subgrains before annealing, $n$$n$nn$n$ is the power exponent, and $k$$k$kk$k$ is the temperature-dependent material constant. It is shown that the developed model has great potential for detailed study of subgrain structure; the obtained results are given below.

Fig. 4 illustrates a comparison between the contributions of coalescence and boundary migration processes to the subgrain structure evolution. The contributions were assessed by determining the fraction of subgrains $f_N=N/N_0$$f_N=N/N_0$f_(N)=N//N_(0)f_{N}=N / N_{0}$f_N=N/N_0$, where $N$$N$NN$N$ is the total number of absorbed subgrains as a result of migration and coalescence processes. The activation energy of subgrain coalescence is higher than the activation energy of low-angle boundary migration; therefore, migration occurs at lower temperatures. The results obtained in the temperature range $200-{300}^{\circ }\mathrm{C}$$200-{300}^{\circ }\mathrm{C}$200-300^(@)C200-300^{\circ} \mathrm{C}$200-{300}^{\circ }\mathrm{C}$ show that the primary mechanism of recovery is the migration of low-angle subgrain boundaries. Subgrain coalescence is more active at higher temperatures ( $300-340{}^{\circ }\mathrm{C})$$300-340\left.{}^{\circ }\mathrm{C}\right)$300-340{:^(@)C)300-340 \left.{ }^{\circ} \mathrm{C}\right)$300-340{}^{\circ }\mathrm{C})$. The results highlight the competition between realized recovery mechanisms and show their impact on subgrain size changes.

Figure 4: Evaluation of the contributions of coalescence and migration processes to the growth of subgrain sizes; the ordinate axis shows the quantitative fraction $f_N$$f_N$f_(N)f_{N}$f_N$ of absorbed subgrains by the mechanisms of migration and coalescence.

Fig. 5a shows the dependence of the released specific energy in the considered volume of subgrains vs time. The stored energy decreases over time, which aligns with previously described physical mechanisms of migration and coalescence of subgrains, which reduce it. When the process temperature increases, the energy release rate (a tangent line to this graph) increases as well. Fig. 5b shows the dependence of the average angle of mutual misorientation of subgrains vs time. The mobility of low-angle boundaries increases with increasing the angle of misorientation of adjacent subgrains (formula (2)); therefore, migration is more active at boundaries with a high angle of mutual misorientation of subgrains. Thus, the average angle of mutual misorientation in the volume of subgrains decreases. In turn, coalescence manifests itself more actively at low angles of misorientation between neighboring subgrains, since the common boundary dissociation can be realized over shorter periods of time. As can be seen from Fig. 5b, at low temperatures ( $200-{260}^{\circ }\mathrm{C}$$200-{260}^{\circ }\mathrm{C}$200-260^(@)C200-260^{\circ} \mathrm{C}$200-{260}^{\circ }\mathrm{C}$ ), the angle of mutual misorientation continuously decreases, which is coherent with high intensity of the migration process. At high temperatures ( $300-{340}^{\circ }\mathrm{C}$$300-{340}^{\circ }\mathrm{C}$300-340^(@)C300-340^{\circ} \mathrm{C}$300-{340}^{\circ }\mathrm{C}$ ), it decreases at the initial stage of annealing and then gradually increases, which is coherent with high intensity of the coalescence process. Note that the migration of a subgrain boundary enlarges its defectiveness and, as a result, its misorientation increases [41, 42]. This phenomenon is not addressed within the structure of the model proposed.

The histograms presented in Fig. 6 show the evolution of subgrain size distribution for the Inconel 718 alloy for the moments in time $t=0,2\cdot {10}^4$$t=0,2\cdot {10}^4$t=0,2*10^(4)t=0,2 \cdot 10^{4}$t=0,2\cdot {10}^4$ and $4\cdot {10}^4\mathrm{s}$$4\cdot {10}^4\mathrm{s}$4*10^(4)s4 \cdot 10^{4} \mathrm{~s}$4\cdot {10}^4\mathrm{s}$ at temperature ${340}^{\circ }\mathrm{C}$${340}^{\circ }\mathrm{C}$340^(@)C340^{\circ} \mathrm{C}${340}^{\circ }\mathrm{C}$ (theoretical distribution - solid lines). Based on the results, it can be concluded that the average subgrain size increases, followed by an increase in distribution dispersion. However, before annealing, the initial subgrain size distribution resembles the given Rayleigh distribution, and, during annealing, it tends to evolve into a lognormal distribution, which is consistent with the experimental results from [5, 34, 43].

Fig. 7 displays the evolution of the largest subgrain conglomerate formed at the final stage of annealing Inconel 718 at 340 ${}^{\circ }\mathrm{C}$${}^{\circ }\mathrm{C}$^(@)C{ }^{\circ} \mathrm{C}${}^{\circ }\mathrm{C}$ for different moments in time $t=0,2\cdot {10}^4$$t=0,2\cdot {10}^4$t=0,2*10^(4)t=0,2 \cdot 10^{4}$t=0,2\cdot {10}^4$ and $4\cdot {10}^4\mathrm{s}$$4\cdot {10}^4\mathrm{s}$4*10^(4)s4 \cdot 10^{4} \mathrm{~s}$4\cdot {10}^4\mathrm{s}$. The obtained results indicate the presence of abnormally large subgrains in the alloy under certain conditions, which is consistent with the experimental observations described before [41, 44]. The abnormal growth of subgrains occurs in materials when the structure contains special (favorably oriented) subgrains that have an increased value of the driving force and, correspondingly, mobility and size relative to their neighbors. For pure nickel, the number of such subgrains is great regardless of the generated sample. For Inconel 718, the braking force from the second-phase particles has a strong impact on the growth of subgrains. Therefore, in Inconel 718, the number of subgrains prone to anomalous growth is significantly reduced. The anomalous growth of subgrains may or may not occur in the material structure depending on the samples of subgrains generated for Inconel 718.
According to the Large Angle Boundary Migration (LABM) mechanism, the anomalous subgrains, characterized by a large size and increased misorientation, are considered to be the potential recrystallization nuclei [ $15,41,43$$15,41,43$15,41,4315,41,43$15,41,43$ ]. The model presented does not consider the formation of recrystallized grains; therefore, the subgrain size stops growing when it reaches a limiting value.
An effective method for assessing the adequacy of a mathematical model is a comprehensive analysis of its sensitivity to perturbations of parameters and input data [45]. The methodology used in this procedure is described in detail in [45]. In the present study, a sensitivity analysis of the model was carried out with respect to the parameters $f_v$$f_v$f_(v)f_{v}$f_v$ and ${\varrho }_{av}$${\varrho }_{av}$ϱ_(av)\varrho_{a v}${\varrho }_{av}$ of the Zener relationship (5), the initial distribution of subgrain misorientation angles $\theta$$\theta$theta\theta$\theta$, the parameters $m_{hag,0}$$m_{hag,0}$m_(hag,0)m_{h a g, 0}$m_{hag,0}$ and $Q_b$$Q_b$Q_(b)Q_{b}$Q_b$ of the migration relationship (3), and the geometric coefficient $\beta$$\beta$beta\beta$\beta$ of relationship (7). This set of parameters most fully characterizes the behavior of the Inconel 718 alloy system under consideration. Following the methodology of [45], at the first stage of the analysis, a baseline (unperturbed) solution was determined, relative to which the robustness of the developed model was assessed. This solution was taken as the dependence $d_{av}(t)$$d_{av}(t)$d_(av)(t)d_{a v}(t)$d_{av}(t)$ obtained at a temperature of ${300}^{\circ }\mathrm{C}$${300}^{\circ }\mathrm{C}$300^(@)C300^{\circ} \mathrm{C}${300}^{\circ }\mathrm{C}$ for the Inconel 718 alloy with a sample of 10000 subgrains and the baseline set of parameters (Tab. 1). At the second stage, the parameter under investigation was perturbed. In the general case, the perturbed parameter $A^{\ast }$$A^{\ast }$A^(**)A^{*}$A^{\ast }$ is defined as $A^{\ast }(0)=A(0)(1+\varphi )$$A^{\ast }(0)=A(0)(1+\varphi )$A^(**)(0)=A(0)(1+phi)A^{*}(0)=A(0)(1+\phi)$A^{\ast }(0)=A(0)(1+\varphi )$, where $A$$A$AA$A$ is the baseline value and $\varphi$$\varphi$phi\phi$\varphi$ is a random variable uniformly distributed within the interval $[-\delta ,\delta ]$$[-\delta ,\delta ]$[-delta,delta][-\delta, \delta]$[-\delta ,\delta ]$. At the third stage, a plan of computational experiments was developed to estimate the deviations from the baseline solution, and the ranges of relative perturbations $\delta$$\delta$delta\delta$\delta$ were chosen. At the fourth stage, a series of computational experiments was performed according to the developed plan. At the fifth stage, the relative norms of the deviations of the perturbed parameters ${\mathrm{\Delta }}_A=\frac{\Vert A^{\ast }(0)-A(0)\Vert }{\Vert A(0)\Vert }$${\mathrm{\Delta }}_A=\frac{‖A^{\ast }(0)-A(0)‖}{\Vert A(0)\Vert }$Delta_(A)=(||A^(**)(0)-A(0)||)/(||A(0)||)\Delta_{A}=\frac{\left\|A^{*}(0)-A(0)\right\|}{\|A(0)\|}${\mathrm{\Delta }}_A=\frac{\Vert A^{\ast }(0)-A(0)\Vert }{\Vert A(0)\Vert }$ and of the response ${\mathrm{\Delta }}_{d_{av}}=\frac{\Vert d_{av}^{\ast }(t)-d_{av}(t)\Vert }{\Vert d_{av}(t)\Vert }$${\mathrm{\Delta }}_{d_{av}}=\frac{‖d_{av}^{\ast }(t)-d_{av}(t)‖}{‖d_{av}(t)‖}$Delta_(d_(av))=(||d_(av)^(**)(t)-d_(av)(t)||)/(||d_(av)(t)||)\Delta_{d_{a v}}=\frac{\left\|d_{a v}^{*}(t)-d_{a v}(t)\right\|}{\left\|d_{a v}(t)\right\|}${\mathrm{\Delta }}_{d_{av}}=\frac{\Vert d_{av}^{\ast }(t)-d_{av}(t)\Vert }{\Vert d_{av}(t)\Vert }$ were determined [45], and a conclusion regarding the model's robustness was drawn.