Recrystallization at crack surfaces as a specific fracture mechanism at elevated temperatures—Cellular automata simulation Текст научной статьи по специальности «Физика»

УДК 539.3, 539.5

Recrystallization at crack surfaces as a specific fracture mechanism at elevated temperatures—Cellular automata simulation

D.D. Moiseenko1, P.V. Maksimov1, S.V. Panin12, S. Schmauder3, V.E. Panin1, D.S. Babich4, F. Berto5, A.Yu. Vinogradov56, and A. Brückner-Foit7

1 Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634055, Russia 2 National Research Tomsk Polytechnic University, Tomsk, 634050, Russia 3 Institute for Materials Testing, Materials Science and Strength of Materials (IMWF), University of Stuttgart, Stuttgart, 70569, Germany

4 National Research Tomsk State University, Tomsk, 634050, Russia 5 Norwegian University of Science and Technology, Trondheim, 7491, Norway

6 Togliatti State University, Togliatti, 445020, Russia

7 Institut für Werkstofftechnik Qualität und Zuverlässigkeit, University ofKassel, Kassel, 34125, Germany

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.

Keywords: hot crack, recrystallization, fracture, computational solid mechanics, multiscale modeling, cellular automata

DOI 10.24411/1683-805X-2018-15003

Рекристаллизация материала вблизи берегов трещины как особый механизм разрушения при высоких температурах. Моделирование методом клеточных автоматов

Д.Д. Моисеенко1, П.В. Максимов1, С.В. Панин1,2, S. Schmauder3, В.Е. Панин1, Д.С. Бабич4, F. Berto5, А.Ю. Виноградов5,6, A. Brückner-Foit7

1 Институт физики прочности и материаловедения СО РАН, Томск, 634055, Россия 2 Национальный исследовательский Томский политехнический университет, Томск, 634050, Россия

3 Университет Штутгарта, Штутгарт, 70569, Германия 4 Национальный исследовательский Томский государственный университет, Томск, 634050, Россия 5 Норвежский университет естественных и технических наук, Тронхейм, 7491, Норвегия 6 Тольятинский государственный университет, Тольяти, 445020, Россия

7 Университет Касселя, Кассель, 34125, Германия

Предложен новый дискретно-континуальный подход клеточных автоматов для моделирования процессов зародышеобразования и последующего роста новых фаз (зерен). Метод объединяет подходы классической термоупругости и логику переключений клеточных автоматов. В рамках дискретно-континуального подхода клеточных автоматов пространство, занимаемое моделируемым образцом, представляется в виде клеточного автомата — множества упорядоченных активных элементов. Каждый элемент имитирует неподвижную область пространства, содержащую часть материала, которая характеризуется определенным набором числовых параметров. Предложенный алгоритм позволяет вычислять величину локальных моментов сил и моделировать диссипацию энергии кручения, определяющую образование новых дефектных структур. Для иллюстрации возможностей предлагаемого метода были проведены численные эксперименты по термически активируемой рекристаллизации вблизи берегов трещины в чистом титане. С использованием трехмерных клеточных автоматов выполнено моделирование эволюци микроструктуры образца с V-образным надрезом, имитирующим вершину трещины, под действием интенсивных термических нагрузок. Расчет термических напряжений в результате теплового расширения под действием ротационных мод деформации позволил оценить величину локальной энтропии и исследовать эволюцию процесса дефектообразования — от зарождения исходных дефектных структур до стадии их интенсивного накопления. В работе обсуждаются пути дальнейшего развития метода с позиции моделирования термически активируемого двойникования и иных структурно-фазовых превращений.

Ключевые слова: горячая трещина, рекристаллизация, разрушение, вычислительная механика твердых тел, многоуровневое моделирование, клеточные автоматы

© Moiseenko D.D., Maksimov P.V., Panin S.V., Schmauder S., Panin V.E., Babich D.S., Berto F., Vinogradov A.Yu., Brückner-Foit A., 2018

1. 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 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.

2. Simulation technique and algorithms

The new algorithm simulating the recrystallization process is supposed to account for the possibility of a twinlike 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

t" = Tn-1+-1- iQ,, (1)

CiPiV k =1

where T^1, T" are the temperatures of ith element at

(n - 1)th and nth time steps, c is the heat capacity of ith element, p. is its density, V is the element volume, Q^ denotes the flux of thermal energy from a neighboring element with an index k into the element under consideration with index i at nth time step, and Nis the number of neighbors. Hereinafter, all upper indices mean the step numbering incrementing in time, unless otherwise is specified. The change in the thermal energy Qik is calculated from the Fourier law

Qk= ^ (TT1 - Tn -1) AT.

(2)

Here Xik stands for the coefficient of cross heat conduction, l is the distance between the centers of the elements under consideration, Q is the area of the adjacent face, At 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. Fine details of the algorithm for the modeling of the elastic energy transport have been comprehensively reported earlier [6-10].

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 1th—at the intersection of the 1st coordination spheres of the ith and the corresponding kth elements) is calculated as (see Fig. 1 for vector notations)

=

rikl'

(3)

The matter flux velocity at the boundary of the kth and lth elements is calculated using the Thornbull relation

vkl = mkl(Pl- Pk )nkl, (4)

where pj, pk denote the values of pressure (hydrostatic stress) in the kth and Ith elements, respectively, nkl is the vector normal to the boundary of the kth and lth elements, and mkl 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

K L

®i = XX®ikl. (5)

k =1l =1

Here Kis 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 T(Ayi) is proportional to the total angular velocity

Ay, = ©,. T. (6)

The variation of the force moment of the ith element over the time t(AM i) is calculated as follows:

AM: =

Gn rc.Ayi _ Gn rc ©^t

(7)

2 2

Here G is the shear modulus of the material in the ith element, and rc is the radius of the element. The term being responsible for the accumulation of latent energy of defects is expressed as

[AA,

. ktorsG/ | Ay" | nrc

Aji.

(8)

rikl

Here ktors 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

Fig. 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)

necessary to run a reversible structural-phase transformation. Then, the coefficient ktors 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. For each ith element, the initial values of temperature T°> thermal energy Qf, elastic energy E0, and dislocation energy O0 are set. Each ith grain is also characterized by Euler orientation angles y t, rip The following algorithm is proposed for the simulation of the nucleation and growth of a new phase under thermal loading. For each element with index i the probability of a new grain nucleation inside the ith grain with the initial phase (psw )n is based on the temperature of the element T" which is calculated at every time loop of the algorithm. When the temperature of a certain element reaches the critical value Tsw1, (psw)" takes a nonzero value according to the following equation:

(Psw )" =

o, t" < ts

2(Psw )B

sw1,

sw1

T - T

sw2 sw1

Tsw1 - Ti" < (Tsw1 + Tswl)/2, 2(Psw )B

T T- Tn

sw2 i

T - T

sw2 sw1

(9)

(Tsw1 + Tsw2 V2 - Ti" < Tsw2 ,

[0, Tsw2 - Tn .

Here Tsw2 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 located 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:

(p )" = (p )

\rswn \rswf n

(1 + sign(T" - Tw1)) X

X (1 + sign (Tsw1 + Tsw2 - 2Tn)) T" Tw1 ■

2(Tsw2 T sw1)

+ (1 + sign (TSw2 - Tn ))(1 + sign (2Tn -

— T - T )) Tsw2 Ti Tsw1 Tsw2)) ~(rp _rp \ 2(Tsw2 Tsw1)

where signx = -1 at x <0, while signx =1 at x > 0.

(10)

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 Tsw1, Tsw2 and psw 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. Following calculations are performed at each nth time step of the algorithm. The K" parameter that is equal to the number of neighbors belonging to grains with orientation angles different from 67 (orientation angle of the Ith grain containing the ith element) by less than 6HAGB (maximal possible angle of boundary misorientation) is calculated for each ith element. The probability of transition of ith element from the current Zth grain in each of the grains (with the index K) is computed:

p" = p" (EK-1 - e"-1), 1 - k - k" . (11)

First of all, the specific energy gK of the boundary between grains I and K is calculated as [11]

( 6 >

Y IK = Y HAGB "

IK

HAGB

1 - ln-

IK

HAGB

(12)

where 6IK = |yI -yK | is the grain misorientation angle between Ith and Kth grains.

Further, according to the Turnbull equation, the grain boundary velocity is calculated as

2^yn

v"k = (mkp"k + mkPk + m" pk )e kB(T"+T">,

where

pik =

Pik =

Pik =

Qi-1 - Qk"-1

V

E" 4 - E"-1

V

O" - o"-1

V

(13)

(14)

(15)

(16)

are the driving forces initiated by the gradients of thermal energy, elastic energy and defect energy, correspondingly. In Eqs. (14)-(16) m*k is the grain boundary mobility in the thermally activated recrystallization process, m" is the grain boundary mobility in the mechanically driven ("cold") dynamic recrystallization process, m" is the grain boundary mobility in recrystallization process initiated by the dislocation density gradient, Q is the area of adjacent active elements, and kB is the Boltzmann constant.

The energy of defects in the ith element ©n at the nth time loop is calculated by the following relation

©; = -dipniGb, (17)

where is the coefficient depending on the spatial distribution of the dislocations, p" is their density, and bi is the Burgers vector.

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 - kd Ay/2n, where Ay 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 Fig. 1, b).

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 v;k 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

a = |cos(u, a)|, P = |cos(u, b)|, y = |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

8 = max(a, P, y) = |cos (u, d)|. (19)

(iii) By introducing the anisotropy effect coefficient k, the value vink is redefined as

Vk =| d| vikek(8-1),

(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 P" of the active element is calculated as:

max (0, -v") • t

n Pik

l

(21)

where t 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 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. For this purpose, one can

construct a membership function, which accounts for all switching probabilities Pikn as

f k-1 Y k \

f" (n, k)=2 -

n-E p;

j=1

n-E p;

;=1

k -1

n-E p;

;=1

n-E p;

;=1

(22)

This function equals to unity when the random number n belongs to the interval

f k -1 k \

\

E p; ; E Pjn j=0 j=0

(23)

and it is equal to zero when n is beyond the interval (23). Moreover, when the element stays in the old grain,

k = K; +1, (24)

and the probability becomes equal to

P K n+1 = 1 - E P. (25)

j =1

Summing up the aforesaid reasoning, the new state (grain index) of ith element at the next time loop will be as follows:

k; +1

7n+l

= E k • fk (n, k).

(26)

k=1

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. From the thermodynamics point of view this assumes a local variation of the entropy and temperature which occur according to the following relation:

[ AEd ] n = Tin-1AS; + S;-1AT". (27)

These relationships will be further used for the criteria of twin structures nucleation as well as subsequent failure of the material.

3. Numerical simulations of recrystallization process

Hybrid cellular automata simulations were tuned for the thermally activated recrystallization process in pure tita-nium—a 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 accumulation 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 p.m. The

Fig. 2. Initial structure of the specimen simulated (a) and the loading scheme (b) (color online)

1800 1425

Î1050

675

300

Fig. 3. Evolution of the temperature field during heat loading in the vicinity of the notch tip: 0.5 (a), 1.0 (b), 3.0 (c), 5.0 |is (d) (color online)

Fig. 4. The patterns of evolving grain structure in the V-notched specimen at different time steps: 0.5 (a), 1.0 (b), 3.0 (c), 5.0 |is (d) (color online)

Fig. 5. The evolution of the specific elastic energy at the surface of the V-notched specimen at various time steps: 0.5 (a), 1.0 (b), 3.0 (c), 5.0 |is (d) (color online)

specimen dimensions were 80x120x10 p.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. For 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 Fig. 2.

Figures 3-7 depict the evolution of the spatial distribution of the temperature, grain structure, elastic 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 inhomogeneous distribution of thermal stresses generates the appearance of local force moments that serve as the sources of the crystal lattice curvature. Further, the lat-

Fig. 6. The spatial distribution of the specific force moments at the specimen surface: 0.5 (a), 1.0 (b), 3.0 (c), 5.0 |is (d) (color online)

|M|, kPa

40

30

20

- 10

Fig. 7. The spatial distribution of the specific force moments in the V-notched specimen: 0.5 (a), 1.0 (b), 3.0 (c), 5.0 |is (d) (color online)

tice 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. Fig-

ure 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 (Fig. 7) and its surface layer (Fig. 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 trans-lational symmetry.

Fig. 8. The evolution of the defect accumulation rate in the specimen surface layer in the vicinity to the crack face: 0.5 (a), 1.0 (b), 3.0 (c), 5.0 |is (d) (color online)

Fig. 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)

the experiment by F. Kuhn, F. Zeismann and A. Bruckner-Foit, 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. For sake of simplicity at the current stage two points limiting the temperature range of the twin structure forma-

The spatial distribution of the ratio of the torsion energy increment over the elastic energy influx per the computational algorithm step is presented in Fig. 8. This parameter is directly related to the rate of new defective structures generation in the material.

4. Simulation of two-stage recrystallization process along the hot crack faces, 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

tion 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 (pw )" 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 be written as follows:

(PtwX" = (PtwX,

(1 + sign(T" - Twi)) X

T" - T

tw1

X (1 + sign(Ttw1 + Ttw2 2T"))- „ x

2(Ttw2 Ttw1)

+ (1 + sign(7^2 - Tn ))(1 + sign(2T" - Tw! -

- Tw2 ))

Ttw2 Tf

2(Ttw2 Ttw1)

(28)

where signx = -1 at x <0, while signx =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 recrys-tallization, 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 becomes possible between the grains of new phase. 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.

The test specimen has dimensions of 120x 120x 6 ^m and is simulated as a cell automaton with an FCC type of element packing whose size made 1 ^m. The specimen consists of two grains: the average grain size is equal to 120 x 60 x 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 1800K, 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 Fig. 9.

When comparing these results with the ones of columnar structure modeling (Figs. 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.

5. Summary

The new proposed cellular automata approach for re-crystallization simulation was developed with account for the possibility of a twin-like growth of "new phase" crystallites. The method combines cellular automata 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.

One of the authors (A. Vinogradov) acknowledges the partial support from the Ministry of Education and Science of RF through the State Assignment according to the contract No. 11.5281.2017/8.9. All co-authors have contributed substantially to the paper.

References

1. Ren F., Chen F., Chen J. Investigation on dynamic recrystallization behavior of martensitic stainless steel // Adv. Mater. Sci. Eng. - 2014. -doi 10.1155/2014/986928.

2. Wang Y., Shao W.Z., Zhen L., Lin L., Cui Y.X. Investigation on dynamic recrystallization behavior in hot deformed superalloy Inconel 718 // Mater. Sci. Forum. - 2007. - P. 546-549, 1297-1300. - doi 10.4028/www. scientific.net/MSF.546-549.1297.

3. Panin V.E., Egorushkin V.E., Moiseenko D.D., Maksimov P.V., Kul-kov S.N., Panin S.V Functional role of polycrystal grain boundaries and interfaces in micromechanics of metal ceramic composites under loading // Comp. Mater. Sci. - 2016. - V. 116. - P. 74-81. - doi 10.1016/j.commatsci.2015.10.045.

4. Moiseenko D.D., Maksimov P. V., Panin S.V, Panin V.E. Defect Accumulation in Nanoporous Wear-Resistant Coatings under Collective Recrystallization. Simulation by Hybrid Cellular Automaton Method // Proc. VII Eur. Cong. Comput. Meth. Appl. Sci. Eng. / Ed. by

M. Papadrakakis, V. Papadopoulos, G. Stefanou, V. Plevris. - Crete Island, Greece, 5-10 June 2016, V. 1. - P. 2080-2098.

5. Moiseenko D.D., Panin S.V., Maksimov P. V., Panin V.E., Berto F. Behavior of nanoporous thermal barrier coatings under cyclic thermal loading. Computer-aided simulation // AIP Conf. Proc. - 2015. -V. 1683. - P. 020152-1-020152-5. - doi 10.1063/1.4932842.

6. Moiseenko D.D., Panin V.E., Elsukova T.F. Role of local curvature in grain boundary sliding in a deformed polycrystal // Phys. Mesomech. -2013. - V. 16. - No. 4. - P. 335-347. - doi 10.1134/S10299599-13040073.

7. Moiseenko D.D., Pochivalov Yu.l, Maksimov P. V., Panin V.E. Rotational deformation modes in near-boundary regions of grain structure in a loaded polycrystal // Phys. Mesomech. - 2013. - V. 16. - No. 3. -P. 248-258. - doi 10.1134/S1029959913030077.

8. Panin S.V., Vinogradov A., Moiseenko D.D., Maksimov P. V., Berto F., Byakov A.V., Eremin A.V., Narkevich N.A. Numerical and experimental study of strain localization in notched specimens of a ductile steel on meso- and macroscales // Adv. Eng. Mater. - 2016. - V. 18. -P. 2095-2106. - doi 10.1002/adem.201600206.

9. Moiseenko D.D., Panin V.E., Maksimov P.V., Panin S.V., Berto F. Material fragmentation as dissipative process of micro rotation sequence formation: hybrid model of excitable cellular automata // AIP Conf. Proc. - 2014. - V. 1623. - P. 427-430. - doi 10.1063/1.4898973.

10. Moiseenko D.D., Panin S.V., Maksimov P.V., Panin V.E., Babich D.S., Berto F. Computer simulation of material behavior at the notch tip: Effect of microrotations on elastic energy release // AIP Conf. Proc. -2016. - V. 1783. - P. 020157-1-020157-5. - doi 10.1063/1.4966450.

11. Humphreys F.J., Hatherly M. Recrystallization and Related Annealing Phenomena. - Oxford: Elsevier, 2004. - P. 95-97.

12. Kroc J. Application of cellular automata simulations to modelling of dynamic recrystallization // Lect. Notes Comput. Sci. - 2002. -V. 2329. - P. 773-782. - doi 10.1007/3-540-46043-8_78.

13. Godara A., Raabe D. Mesoscale simulation of the kinetics and topology of spherulite growth during crystallization of isotactic polypropylene (iPP) by using a cellular automaton // Model. Simul. Mater. Sci. Eng. - 2005. - V. 13. - P. 733-751. - doi 10.1088/0965-0393/ 13/5/007.

14. Kuhn F., Zeismann F., Brueckner-Foit A. Crack growth mechanisms in an aged superalloy at high temperature // Int. J. Fatig. - 2014. -V. 65. - P. 86-92.

nocTynaaa b peaaKUHro 07.08.2018 r.

CeedeHua 06 aemopax

Dmitry D. Moiseenko, Cand. Sci. (Phys.-Math.), Senior Researcher, ISPMS SB RAS, [email protected]

Pavel V. Maksimov, Cand. Sci. (Phys.-Math.), Researcher, ISPMS SB RAS, [email protected]

Sergey V. Panin, Dr. Sci. (Eng.), Prof. of RAS, Deputy Director, ISPMS SB RAS, Prof., TPU, [email protected]

Siegfried Schmauder, Prof. Dr. Rer. Nat., Prof., University of Stuttgart, Germany, [email protected]

Victor E. Panin, Dr. Sci. (Phys.-Math.), Academician of RAS, Head of Laboratory, ISPMS SB RAS, [email protected]

Dmitriy S. Babich, Graduate Student, Tomsk State University, [email protected]

Filippo Berto, Prof., Norwegian University of Science and Technology, Norway, [email protected]

Alexey Yu. Vinogradov, Ph.D., Dr. Eng., Prof., Norwegian University of Science and Technology, Norway, Head of Laboratory, Deputy Director,

Institute of Advanced Technologies, Togliatti State University, [email protected]

Angelika Bruckner-Foit, Prof. Dr. Rer. Nat., UniversityofKassel,Germany, [email protected]