Recrystallization at Crack Surfaces as a Specific Fracture Mechanism at Elevated Temperatures—Cellular Automata Simulation

A new hybrid discrete-continuum cellular automata approach is proposed to simulate the process of new phase/grain nucleation and growth. The method couples classical thermomechanics and the logics of cellular automata switching. Within the framework of the hybrid discrete-continuum cellular automata method, the space occupied by the simulated specimen is represented as a cellular automaton—a set of ordered active elements. Every element imitates an immovable region of space related to a part of material being characterized by the certain numerical parameters. The proposed approach enables calculating the magnitude of the local force moments and simulating dissipation of torsion energy leading to the formation of new defect structures. To illustrate the capacity of the proposed hybrid discrete-continuum cellular automata approach, the numerical simulations of thermally activated recrystallization of pure titanium near crack faces were conducted. The 3D cellular automaton simulated the microstructure evolution of the V-notched specimen region that imitated the crack tip vicinity at high homologous temperatures. Calculation of heat expansion with simultaneous thermal stresses accumulation and microrotation initiation was incorporated in the simulations permitting thereby to evaluate the local entropy and to monitor the evolution of crystal defects from initiation to storage. Perspectives of the proposed algorithms for simulations of the mechanical behavior of materials experiencing thermally induced twining or phase transformations are discussed.


INTRODUCTION
Over the last 50 years, the rapid development of aircraft and space technology as well as ever-increasing demands for equipment in the steel industry emerge new problems for the design of materials and structures operating under extreme conditions. Among them are materials employed for manufacturing nozzles for jet engines, tuyeres and crystallizers of blast furnaces, etc. High temperature gradients and drops, powerful thermal impacts at engine start or filling melt into the crystallizer, etc., with the regard of cyclic repetition of the processes result in catastrophic degradation of the materials used. The material structural degradation is very complex in nature. It is associated with various interrelated concurrent processes which are accompanied by rearrangement of the internal structure of the PHYSICAL MESOMECHANICS Vol. 23 No. 1 2020 substrate material, cracking and spallation of the protective coating, nucleation of the brittle thermal grown oxide layer, etc.
Experimental studies of materials under extreme conditions and subsequent full-scale tests of semi-finished structures are very laborious and expensive. In this regard, theoretical predictions based on microstructurally informed computer simulations are very relevant. Increasing demands for new materials call for the necessity to design and develop new advanced methods of computer simulation of the behavior of heterogeneous materials under extreme loading conditions.
The above said is of particular importance at the design of steam turbine blades which are widely accepted by the industry. Titanium alloys are popular blade materials, particularly in low pressure steam turbines that are most influencing for the power output, size and efficiency of steam turbines. The most promising ways of further increasing the efficiency of steam turbines is by extending the lengths of low pressure steam turbine blades and by increasing the steam reheat temperature. However, this results in substantial centrifugal stresses and creep of blades and rotors. As a consequence, advanced titanium alloys and manufacturing processes are needed to design the blades with sufficient reserve strength. Moreover, one of the most significant tasks is the coupled problem of cracking and the simultaneous recrystallization process under heat loading.
Simulation of recrystallization processes, which occur in materials for steam turbines at high homologous temperatures, is still challenging. Successful examples of simulation of this kind are provided in recent studies [1,2]. In the current work, we propose a hybrid discrete-continuum cellular automata approach based on coupling the classical thermomechanics and logics of cellular automata switching to simulate new phase generation and grain growth [3]. On the basis of the hybrid discrete-continuum cellular automata approach the numerical simulation of thermally activated recrystallization were conducted for pure titanium in the vicinity of crack edges. In doing so, the 3D cellular automaton simulates the behavior of the V-notched specimen region that mimics the crack tip vicinity. Simulations calculate numerically heat expansion in the material with account for thermal stresses accumulation and microrotation initiation. The latter gives rise to the generation of new lattice defects and increasing the local entropy. Every newly-nucleated grain has zero dislocation density, so that the gradient of the dislocation density between the existing and newly forming grains gives rise to an additional driving force for the growth of new grains (along with the thermal gradient).
In the present paper we focus on developing a new advanced modeling algorithm. The latter ensures carrying out detailed studies of the processes related to rearrangement of the internal structure at the crack tip under intensive thermal loading. It should be noticed that, the next stage of this study will be related to finding out the model parameters corresponding to physical characteristics of the specific material and comparison of the simulation results with experimental ones.

SIMULATION TECHNIQUE AND ALGORITHMS
The new algorithm simulating the recrystallization process is supposed to account for the possibility of a twin-like growth of new phase crystallites. The discrete-continuum method of heat transfer was designed and implemented in finite element codes by the authors earlier [4,5]. The transfer of mechanical energy within the framework of the hybrid approach of cellular automata was reported in [610]. The ideas of the modeling were inspired by the results reported in the monograph [11] as well as in Refs. [12,13]. Within the framework of hybrid discrete-continuum cellular automata approach, the space occupied by the simulated specimen is represented as a cellular automaton, i.e. as a set of ordered active elements. Every element imitates an immovable region of space related to a part of material being characterized by the following numerical parameters: the heat capacity, thermal conductivity, density, etc. Additional material properties used in model will be described below when the algorithm details will be clarified. Note that the active element of the cellular automaton itself is stationary; however, the mass and energy flows change its properties.
In order to determine the temperature of the elements, a discrete model of heat transfer in a heterogeneous medium should be introduced [4,5]. The values of the temperature, strain, thermal conductivity, heat capacity and coefficient of linear thermal expansion of each element are specified as initial conditions. Then, for each nth time step, a new value of the element temperature is calculated taking into account the heat fluxes from each neighboring element Here ik λ stands for the coefficient of cross heat conduction, l is the distance between the centers of the elements under consideration, Ω is the area of the adjacent face, ∆τ denotes the time increment.
As a result of inhomogeneous thermal expansion, internal thermal stresses are accumulated in various structural elements, which in turn are transformed into flows of the mechanical energy being redistributed over the mesh of the cellular automaton. Grain boundaries are explicitly considered through including their curvature and angles of misorientation between the neighboring grains. The developed algorithms enable calculating the magnitude of the local force moments, the vorticity tensor and dissipation of torsion energy through the forming new defect structures. .ine details of the algorithm for the modeling of the elastic energy transport have been comprehensively reported earlier [610].
The angular velocity of the ith element under the action of the matter flux through the boundary of the kth and lth elements (each kth element is located at the 1st coordination sphere of the ith, while each lthat the intersection of the 1st coordination spheres of the ith and the corresponding kth elements) is calculated as (see .ig. 1 for vector notations) 2 .
The matter flux velocity at the boundary of the kth and lth elements is calculated using the Thornbull relation ( ) , p p denote the values of pressure (hydrostatic stress) in the kth and lth elements, respectively, kl n is the vector normal to the boundary of the kth and lth elements, and kl m is the mobility of the boundary between the kth and lth elements.
The total angular velocity of the ith element is defined as the following sum 1 1 .
Here K is the number of elements at the 1st coordination sphere of the ith element, and L is the number of elements at the intersection of the 1st coordination spheres of the ith element and each kth neighbor. The three-dimensional rotation angle of the ith element over the time ( ) i τ ∆ã is proportional to the total angular velocity .
i i ∆ = τ ã ù (6) The variation of the force moment of the ith element over the time ( ) i τ ∆M is calculated as follows: .ig. 1. A schematic for calculating the angular rate of material torsion for a cell automaton element (a); a schematic for the geometric model of torsion induced increment of the dislocation density (b) (color online).
Here G is the shear modulus of the material in the ith element, and c r is the radius of the element. The term being responsible for the accumulation of latent energy of defects is expressed as Here tors k denotes the defect accumulation rate (the quantity that can be measured experimentally).
Note that the expression (8) is quite general in nature and, depending on the constructed switching rules of the cellular automaton can bring the meaning of the energy necessary to run a reversible structural-phase transformation. Then, the coefficient tors k can be calculated on the basis of the atomic configurations of the initial and final states.
At the first stage of the recrystallization simulation a specimen is modeled by a cellular automaton network containing L active elements. The network is divided in M clusters representing individual grains. .or each ith element, the initial values of temperature 0 , is the temperature threshold below which the probability of nucleating a new grain is zero. If the element under consideration becomes a nucleation center of a new grain, then all adjacent elements lo-cated at the first coordination sphere automatically join it. Thus, the grain nucleus of the new phase actually represents a group of elements.
It should be noted that the algorithm offered was implemented as software code for calculations with the help of high-performance graphics cards. In this regard, the use of conditional if ... else statements is highly undesirable, and relation (9) is written as follows: where sgn x = 1 at x < 0, while sgn x = 1 at x ≥ 0. When the grain nucleus appears, it should be taken into account within each phase it occurs: either an old or a new one. If the nucleation takes place within the old phase, then the 1 2 , sw sw T T and sw p parameters from relationships (1) and (2) correspond to the material with the initial value of the defect density. Otherwise, it is assumed the grain appears inside the new phase and these parameters take the values corresponding to ones of the material with the minimum defect density.
Thus, the nucleation of a new grain entails a reduction of the defects density down to a minimum value for all the elements that have joined it. If this takes place, the variable L, which is responsible for the total number of grains in the specimen, increases by 1, while the index of the new grain becomes equal to L. When all the elements of the cellular automaton have been analyzed the growth of new grains is completed. .ollowing calculations are performed at each nth time step of the algorithm. The n i K parameter that is equal to the number of neighbors belonging to grains with orientation angles different from I θ (orientation angle of the Ith grain containing the ith element) by less than HAGB θ (maximal possible angle of boundary misorientation) is calculated for each ith element. The probability of transition of ith element from the current lth grain in each of the grains (with the index K) is computed: .irst of all, the specific energy g IK of the boundary between grains I and K is calculated as [11] HAGB HAGB HAGB 1 ln , where | | IK I K θ = ψ −ψ is the grain misorientation angle between Ith and Kth grains.
.urther, according to the Turnbull equation, the grain boundary velocity is calculated as Within the framework of this algorithm it was assumed that the increment of the dislocation density directly depends on the torsion angle of the material as d 2 , k − ∆γ π where ∆γ is the torsion angle determined by formula (6). Such an approximation is based on a simple assumption on the direct proportionality of the increment of the dislocation density from the sector area as a result of sweeping away the radius of the active element by the torsion angle (see .ig. 1b).
When calculating the probability of grain growth, the crystal anisotropy should be taken into account, which can be determined through the directions and relative lengths of the lattice vectors a, b, c. The maximum relative length of the lattice vector equals to unity.
Thus, the value of the grain boundary velocity n ik L should be corrected by taking into account the relative location of the expected grain growth vector u and vectors a, b, c: (i) Cosine moduli of the angles between the vector u and each of the a, b, c ones are computed as α = |cos(u, a)|, β = |cos(u, b)|, γ = |cos(u, c)|. (18) (ii) Among the vectors a, b, c, the one is chosen whose cosine with the vector u has the maximum absolute value δ = max(α, β, γ) = |cos (u, d)|.
(19) (iii) By introducing the anisotropy effect coefficient k, the value n ik L is redefined as ( 1) | | , n n k ik ik e δ− = d L L (20) where k ≥ 0 (at k = 0 the lattice is isotropic). It should be underlined that when the lattice anisotropy is absent, as well as when the vector u is codirected with one of the lattice vectors, the grain boundary movement rate attains its maximum value.
The switching probability n ik P of the active element is calculated as max (0, ) , where τ is the value of the time step. When a certain ith active element reaches the necessary level of switching, the index of state (grain index) should be found, which will acquire a new state of the element. Here the level of switching depends on the variation between the probability n ik P and the random number generated each time loop for each element.
Along with the grain that will catch the element, the index of neighboring active element belonging to the catching grain should be calculated. .or this purpose, one can construct a membership function, which accounts for all switching probabilities n ik P as This function equals to unity when the random number η belongs to the interval The relation (8) makes it possible to calculate the part of the microscopic energy flows of the material in the cell automaton element, which comes to conditions when the deformation defects start to generate. .rom the thermodynamics point of view this assumes a local variation of the entropy and temperature which occur according to the following relation: 1 1 [Ä ] Ä Ä .

NUMERICAL SIMULATIONS O. RECRYSTALLIZATION PROCESS
Hybrid cellular automata simulations were tuned for the thermally activated recrystallization process in pure titaniuma representative hexagonal material near the crack edges were conducted. The 3D cellular automaton simulates the V-notched specimen region resembling the vicinity of the crack. In the simulation procedure, the heat expansion in the material was taken into account along with the thermal stresses accumu-   lation and microrotation initiation. The latter gives rise to new defect generation and local entropy increase. As mentioned above, every newly-nucleated grain has zero dislocation density. Thus, the dislocation density gradient initiates an additional driving force for new grains growth, which is added to the thermal gradient.
The specimen was represented by the cellular automata consisting of elements with characteristic size of 1 µm. The specimen dimensions were 80 × 120 × 10 µm.
The initial temperature of each element was set at 300 K, and the initial values of strains and stresses were equal to zero. The time step was set at 1 ns. .or numerical calculations, the material constants were taken as those corresponding to pure Ti. In all simulations, the inner surface of the notch was heated to 1800 K. The scheme of the numerical simulation is illustrated in .ig. 2.
.igures 37 depict the evolution of the spatial distribution of the temperature, grain structure, elastic PHYSICAL MESOMECHANICS Vol. 23 No. 1 2020 energy, local force moments and defect accumulation rate.
As one can see from the grain structure distribution patterns, a high temperature gradient gives rise to a characteristic columnar-like growth of the grains. The growth is directed from the region with elevated temperatures towards the cooler bulk of the material. Along with the nucleation and the columnar-like growth of new grains within the bulk, the growth of old initial grains takes place. The thermal expansion of the material is simulated as discussed above. The inhomoge-neous distribution of thermal stresses generates the appearance of local force moments that serve as the sources of the crystal lattice curvature. .urther, the lattice curvature is relaxed at microrotations giving rise to the increasing defect density within the material. The detailed algorithm and simulation results of these processes have been provided in [4]. Since the present study is not focused on the analysis of the stress/strain distributions, in what follows we discuss the results showing the capacity of the proposed method in a qualitative illustrative manner and providing new insights for further developments of the algorithm. It should be particularly noticed that the value of thermal stresses is estimated through the ratio of the elastic energy over the volume of the local structural element. .igure 5 illustrates the thermal stress development in the surface layer.
In this regard, the pattern of spatial distribution of the local force moments in the specimen volume (.ig. 7) and its surface layer (.ig. 6) are of particular interest. These figures illustrate the distribution of the specific values of the local force moment. Despite the fact that these values are three orders of magnitude lower than the level of the effective stresses, it is quite enough to generate a substantial curvature of the crystalline lattice and distort its translational symmetry.
The spatial distribution of the ratio of the torsion energy increment over the elastic energy influx per the computational algorithm step is presented in .ig. 8. This parameter is directly related to the rate of new defective structures generation in the material.

SIMULATION O. TWO-STAGE RECRYSTALLIZATION PROCESS ALONG THE HOT CRACK .ACES, TAKING INTO ACCOUNT THE TWINNING PROCESS
The following modeling was inspired by results of laboratory tests [14], where the mechanism of structure rearrangement in the aged IN738 superalloy was revealed. In the experiment by .. Kuhn, .. Zeismann and A. Brückner-.oit, it was observed that There is a comparatively thin layer on the cracked face that is very porous and heavily damaged. It is possible that this is true for the direct contact of the hot gas with the unprotected material ... This thin layer is fine grained and has no twinning. In doing so, the algorithm was extended up to the model of twin formation in the form of new grains at the boundaries of newly formed and grown grains. The latter have the same structure but with a changed orientation of the crystalline lattice. .or sake of simplicity at the current stage two points limiting the temperature range of the twin structure formation were taken as the twinning criterion. The authors are fully aware on the fact that the temperature can be not the only factor determining this process but the intensity of plastic deformation in each local region of the specimen. The gradient of dislocation density is the driving force for the growth of the twin structure in this algorithm. The former comes from the appearance of rotational deformation modes (see formulas (7), (8), (16), and (17)).
At the boundary between the growing grain of the new phase with the one of the old phase, a twin nucleus can be formed, whose difference from the main grain lies only in the fact that the relative lengths of the basis vectors of the twins (a, b, c) differ from the ones of these vectors in the main grain. The probability of a twin nucleation tw ( ) n i p as well as in the case with the formation of a grain nucleus of a new phase, is determined by the current temperature of the material. If this temperature falls within a certain interval, then this probability becomes nonzero. This relationship can   p  p  T T   T T  T  T  T  T  T   T  T  T T   T  T  T T T where sgn x = 1 at x < 0, while sgn x = 1 at x ≥ 0.
By means of the stochastic exitable cellular automata method a numerical experiment was carried out, where two-layer recrystallization and twinning of grains in a titanium specimen are simulated. At the modeling of two-layer recrystallization, it is assumed that initially due to high-temperature impact grains of the new phase are generated within the grains of the old phase. The former differ from the old one by a reduced content of structural defects. Then the seeds of the new phase tend to grow because of the smaller concentration of defects. Moreover, in a certain temperature range the nucleation of one more generation of grains .ig. 9. Grain structure (a) and spatial distribution of specific elastic energy (b), moment of forces (c) and relative torsion energy (total work of torsion divided by total elastic energy influx) (d) (color online). On the other hand, in the range of lower temperatures, the twinning of new grains takes place. The twins possess the same type of crystalline lattice much like the grains of the new phase, but differ from them by the fact that the growth of the twin occurs along another direction of translational symmetry.

µs
The test specimen has dimensions of 120 × 120 × 6 µm and is simulated as a cell automaton with an .CC type of element packing whose size made 1 µm. The specimen consists of two grains: the average grain size is equal to 120 × 60 × 6 µm. The time step is chosen as 1 ns, the total time of the numerical experiment is equal to 50 µs.
The boundary conditions imitate the thermal load being applied to the V-shaped notch. In doing so, both faces of the notch during the entire numerical experiment are heated up to the temperature of 1800 K, while the end faces of the specimen are cooled down to 300 K. The energy exchange through the facets is absent. The initial values of stresses and strains are set equal to zero. Within the framework of the model, the thermal expansion of the material is taken into account. The latter gives rise to local rotations that are responsible for the accumulation of latent elastic energy (torsion energy). The most picturesque results attained at testing of the proposed algorithm are shown at .ig. 9.
When comparing these results with the ones of columnar structure modeling (.igs. 6, 7, without considering the twinning), one can see that the formation of a solid fine-crystalline layer along the crack faces gives rise to a partial relaxation of thermal stresses. In doing so, the local moments of forces in this layer are remained at high enough level. This will obviously lead to the flaking off the highly damaged material at the crack faces and its gradual growth (expansion). At the same time, unlike the case of the columnar structure, the local moments of forces in the twinning layer are very weakly pronounced. This indicates for the substantial dissipation of the elastic energy of torsion. The latter becomes the driving force for the nucleation and the growth of twin structures.

SUMMARY
The new proposed cellular automata approach for recrystallization simulation was developed with account for the possibility of a twin-like growth of new phase crystallites. The method combines cellular au-tomata switching mechanisms for structure transformation simulation, classical thermal transfer, thermal expansion and generation of lattice defects (dislocations).
It was shown that the material adjacent to the crack faces undergoes the changing of its crystal structure. New grains are generated at the crack faces and grow forming columnar patterns replacing gradually the initial grain structure. This scenario of structural changes leads finally to significant local variation of all mechanical properties of the material in the transformed regions. The present work represents an attempt to develop the adequate simulation framework enabling quantitative modeling of the deformation behavior of bulk parts containing fatigue cracks under large thermal gradients arising from the elevated temperatures at the crack faces to the cooler bulk.
It is shown that in the mostly heated regions of the specimen not only the nucleation of the new crystalline structure occurs, but the fraction of the elastic torsion energy increase is also very high. Under such conditions, the probability of twin-resembling structure nucleation is very high due to the deformation development by the torsion mechanism. In this paper some simple algorithm presented to show how the twins can occur ruled by temperature conditions only. Even in this case one can see that the generation of the twins significantly changes the spatial distribution of the torsion energy. However, the proposed simulation methodology opens up an exciting perspective of simulating stress-induced twinning in modern materials where this mode of deformation mediates (or dominates) the overall mechanical response. These materials include (but are not limited to) steels with the transformation induced plasticity or twinning induced plasticity effect, magnesium and its alloys, etc. The development twinning-specific models is on the way and the results will be reported elsewhere.
.UNDING This work was performed within the frame of direction of research III.23 of Basic Research Program of State Academies of Sciences for 20132020. The work was also supported by the R.BR grant No. 18-08-00516_a and the R. President Council Grant for the support of leading research schools NSh-2718. 2020.8.