A new front-tracking Lagrangian model for the modeling of dynamic and post-dynamic recrystallization
AA new front-tracking Lagrangian model for themodeling of dynamic and post-dynamicrecrystallization
Sebastian Florez ∗ , Karen Alvarado , and Marc Bernacki Mines-ParisTech, PSL-Research University, CEMEF – Centre demise en forme des mat´eriaux, CNRS UMR 7635, CS 10207 rueClaude Daunesse, 06904 Sophia Antipolis Cedex, FranceSeptember 18, 2020 abstract
A new method for the simulation of evolving multi-domains problems has beenintroduced in previous works (RealIMotion), Florez et al. (2020) and furtherdeveloped in parallel in the context of isotropic Grain Growth (GG) with noconsideration for the effects of the Stored Energy (SE) due to dislocations. Themethodology consists in a new front-tracking approach where one of the origi-nality is that not only interfaces between grains are discretized but their bulksare also meshed and topological changes of the domains are driven by selectivelocal remeshing operations performed on the Finite Element (FE) mesh. In thisarticle, further developments and studies of the model will be presented, mainlyon the development of a model taking into account grain boundary migration by(GBM) SE. Further developments for the nucleation of new grains will be pre-sented, allowing to model Dynamic Recrystallization (DRX) and Post-DynamicRecrystallization (PDRX) phenomena. The accuracy and the performance ofthe numerical algorithms have been proven to be very promising in Florez etal. (2020). Here the results for multiple test cases will be given in order tovalidate the accuracy of the model taking into account GG and SE. The com-putational performance will be evaluated for the DRX and PDRX mechanismsand compared to a classical Finite Element (FE) framework using a Level-Set(LS) formulation. ∗ corresponding author a r X i v : . [ c s . C E ] S e p Introduction
The modeling, at the mesoscopic scale, of Grain Growth (GG) and recrys-tallization (ReX) in polycrystalline materials during thermal and mechanicaltreatments has been the focus of numerous studies in the last decades. Indeed,mechanical and functional properties of metals are strongly related to theirmicrostructures which are themselves inherited from thermal and mechanicalprocessing.When looking to the so-called full-field (FF) methods, based on a full de-scription of the microstructure topology and modeling of grain boundary mi-gration (GBM) at mesoscopic scale, main numerical frameworks involve: MonteCarlo (MC) [1, 2], Cellular Automata (CA) [3, 4, 5, 6], Multi Phase-Field(MPF) [7, 8, 9, 10], Vertex/Front-Tracking [11, 12, 13, 14, 15] or Level-Set(LS) [16, 17, 18, 19] models. These numerical methods are developed by manyresearchers [20]. All the mentioned methods have, of course, their own strengthsand weaknesses [20, 21].When large deformation have to be considered (common in metal formingcontext), LS or MPF approaches in context of unstructured FE mesh and FEremeshing strategies remain the main powerful and generic approaches but witha strong limitation in terms of computational cost.In this context vertex and front tracking approaches appear as interestingcandidates. An explicit description of the interfaces is considered and GBMis imposed at each increment by computing the velocity of the nodes describ-ing the interfaces. While having a deterministic resolution (solving of partialdifferential equation - PDE), this methodology is very efficient. However, theimplementation of the topological events is not straightforward and the factto not describe the bulk of the grains could be limiting for some metallurgicalmechanisms such as appearance of new grains (nucleation) or substructures in-side the grains.Previous works dedicated to the creation of an improved front-trackingmethod, solving these weaknesses, have been published in previous articles[22, 23]. The model denominated TOpological REmeshing in lAgrangian frame-work for Large interface MOTION (ToRealMotion, hereafter TRM) maintainsthe interior of grains meshed, handling with relative ease the topological changesof the grain microstructure and allowing the treatment of in-grain operationsand at a higher computational performance than classical FE-LS models for thesame accuracy. The objective of the present article is then to adapt the TRMmodel to handle DRX and PDRX phenomena.2
The TRM model : Isotropic Grain Growthcontext
The TRM model has been presented in a previous work in [22], then adaptedto a parallel computational environment in [23]. This model uses the logic be-hind front-tracking methods where the discretization of interfaces is the minimaltopological information allowing to model 2D-GBM. The TRM model goes astep further by implementing also a discretization of the interior of the grains inthe form of simplexes to allow the interaction of the grain boundaries with thebulk of the grains and preventing inconsistencies of the physical domain suchas the overlapping of regions. The data structure of the TRM model is thenbuilt on top of a mesh with element and nodes, enabling also the possibility tocompute FE problems on it. This data structure defines geometrical entitiessuch as points , lines and surfaces by grouping sets of nodes and elements: each point regroups a P-Node and a set of connections to other points and lines .Each line is defined by an ordered set of L-Nodes , an initial point and a final point . Finally, surfaces are defined as a set of S-Nodes , a set of elements anda set of delimiting lines and points .This data structure can be constructed by performing some preprocessingsteps, in [22], a preprocessor able to transform MPF or LS data to the datastructure of the TRM model has also been introduced. This preprocessor isbased on the works presented in [24] and in [25], where a remeshing proceduretransforms a typical LS configuration (a FE mesh with grain boundaries beingdefined by the interpolated zero iso-value of several LS fields as in [26, 27])into a body-fitted mesh (where the interpolated zero iso-value coincide withsome nodes of the FE mesh) with the help of a joining and fitting algorithm ,expliciting the nodes on the “front” to track (as in front-tracking methods).Subsequently on the preprocessing step, a classification is made for the nodesof the mesh based on their topological representation on the microstructuralspace. Finally, a geometric identification algorithm is performed to build thegeometrical entities mentioned before. The reader is referred to [22] for a com-plete description of the reconstruction process.After the preprocessor, the data in the form of LS fields in no longer neededas the microstructure is now defined by the identified geometric entities. Theclassification of geometric entities is also helpful when computing geometricproperties, the area of surfaces can be computed by adding the contribution ofeach element of the grain while the curvature κ and normal (cid:126)n of interfaces can beobtained by approximating the interface with a high order mathematical form Which defines a node of the mesh with a topology degree equal to 0 on the microstructuralframework, hence a multiple junction. Nodes with a topology degree equal to 1, or a node belonging to a simple grain boundaryin the microstructure. Nodes with a topology degree equal to 2, or nodes in the bulk of grains with no connectionto the S-Nodes on other grains. physical mechanism can be simulated.This physical mechanism represents how the different geometries are supposedto evolve based on their current state. The TRM model has been developed tomove the different nodes of the mesh based on a user defined velocity field (cid:126)v and a time step dt . Once a velocity is defined a new position for each node N i on the mesh can be obtained as: (cid:126)u i = (cid:126)u i + (cid:126)v i · dt, (1)where (cid:126)u i is the current position of the node N i .Here, each node displacement can potentially produce an overlap of someof the elements attached to the node. The TRM model hence ensures the localconformity of the mesh by means of a “local-iteratively movement-halving” thatfinds iteratively the approximated maximal displacement that a node is able tomake in the direction of the velocity v i when an overlap takes place. This pro-cedure ensures at all times that both, the mesh and the microstructural domainare valid.Once several steps of Lagrangian movement are performed, it is highly proba-ble that the quality of the mesh become too poor to continue with the GBM, thisis why the TRM model implements a particular remeshing procedure stronglyinfluenced by the works in [28, 29], that improves the mesh quality and allowstopological events such as grain disappearance. This remeshing procedure mustbe adapted to the data structure previously defined and it has to ensure thatthe final mesh after remeshing is valid to be used by the TRM model. Some ofthe geometrical entities must adapt their sets to take into account the changesmade on the mesh, to do this, several local selective remeshing operators havebeen developed that allow to make changes on the mesh by maintaining a validdata structure: selective vertex smoothing, selective node collapsing, selectiveedge splitting, selective edge swapping and selective vertex gliding, see [22] fora complete description of these operators and the global remeshing procedure. An overlap in a mesh is produced when an element is partially or completely superposedby another element hence disrupting the 1:1 mapping of the numerical domain to the physicaldomain, such a mesh can not be used in a Finite Element resolution. The word selective denotes a variation of the original remeshing operations when per-formed over the data structure of the TRM model, as each remeshing operation will be per-formed differently over nodes with different topology (P-Node, L-Node and S-Node). .1 Parallel Implementation As presented in [23], a distributed memory system, was chosen for our model,allowing the computation of simulation across single CPUs and over a compu-tational cluster where multiple CPUs are connected over the local network andnot over a single mother board. For this implementation, the standard com-munication protocol Message Passing Interface (MPI) [30] was used and thefree library METIS [31] was used in order to obtain the initial partitioning ofthe domain. Additionally, multiple new functions were added to the originalsequential approach in order to make it work in parallel. Among these func-tions, the implementation of a repartitioning algorithm was essential to obtain are-equilibrium of charges as well as a coherent remeshing approach. This repar-titioning algorithm is based in a user-defined ranking system, where a uniquerank is attributed to each processor (in [23] the number of elements was usedto define the ranks), then a Unidirectional Element Sending algorithm was de-veloped to exchange layers of elements between the partitions as illustrated inFig. 1.The remeshing strategy was developed based on the blocking of the bound-aries between partitions, meaning that any application of a remeshing operatorthat changes the boundaries between partition is discarded by default. Nor thecreation or the deletion of nodes (and edges) is allowed at the boundary, also vertex smoothing is not allowed for the Shared-Nodes . Finally, in order toobtain a complete remeshing over the hole domain, the Unidirectional ElementSending algorithm, ensures the motion of the boundaries between partitions,hence unblocking those edges before the element scattering. This is very conve-nient as not a lot of changes are necessary over the original implementation ofthe remeshing operations but inconvenient as the remeshing must be performedtwo times in a single increment.Some other algorithms are necessary for the complete parallel implementa-tion, these algorithms handle the identification and reconstruction of geometricentities across the segmented domain as well as the computation of geometricproperties and the Lagrangian movement of the Shared-Nodes, see [23] for acomplete description.
The simulation of microstructural evolutions are given by the addition of com-plex and different phenomena as GG [17, 18, 32, 19, 33], Recrystallization (ReX)[17, 34, 35, 26, 27, 36, 37] or Zener Pinning (ZP) [38, 39, 40, 41, 42]. In [22], Approach where each processor interacts with its own independent memory location. Nodes present in the memory of multiple processor at the same time, these nodes resideat the boundaries between partitions.
P art , elements 3, 4 and 5appear only to be sent to P art , and not to P art . Elements 1 and 2 appearinitially on the list to be sent to P art and P art but they are filtered in thelast part of the algorithm, as the higher rank of the nodes of these elementsbelongs to P art . c) Selected elements to be sent to P art , the intersection ofelements from P art to be sent to P art and P art is empty [23].6sotropic GG with no influence of SE was used to compare the TRM modelto other approaches (LS-FE [35, 18, 19]), the base model used to represent thisphenomenon is commonly known as migration by curvature flow. The velocity (cid:126)v at every point on the interfaces can be approximated by the following equation: (cid:126)v c = − M γκ(cid:126)n, (2)where M is the mobility of the interface, γ the grain boundary energy, κ the local magnitude of the curvature in 2D and (cid:126)n the unit normal to the graininterface pointing to its convex direction. In an isotropic context as consideredhere, the terms M and γ are supposed as invariant in space.Of course when post-dynamic phenomena such as Static ReX (SRX) or Meta-Dynamic ReX (MDRX) are considered, the SE will act as another driving pres-sure of the GBM. Note that the SE within a grain can be variant, as therecould be regions on the grain that have accumulated more or less dislocationsduring the considered thermomechanical treatment (TMT). At the mesoscopicscale, the SE can be discussed following different hypotheses. Crystal plasticitycalculations and EBSD experimental data can bring dislocation density fieldand so SE field with fine precision until intragranular heterogeneities. Whilethis information is directly usable in pixel/voxel based stochastic approachessuch as MS or CA methodologies, generally it is homogenized by consideringconstant value per grain in deterministic front-capturing (MPF, LS) and front-tracking approaches. If this choice seems quite natural for phenomena wherestored energy gradients and nucleation of new grains are mainly focused on GBlike for discontinuous DRX (DDRX), it could be a strong assumption for phe-nomena where the substructure evolution is important, like for continuous DRX(CDRX). This aspect was for example studied in [43] in context of SRX with aFE-LS numerical framework. It was concluded that intragranular gradients onthe stored energy could indeed have a big impact on the grain morphology andthat simulations taken into account such variations were more in accordancewith experimental observations than simulation using a constant value of storedenergy, but with an important numerical cost as the FE mesh must be thenadapted at the intragranular heterogeneities scale. In the following, a constanthomogenized energy per grain is assumed. Nonetheless, the approach presentedin this article to model GG with a stored energy field can be used in the contextof a heterogeneous intragranular energy field, this aspect will be investigated ina forthcoming publication.Thus, here SE can act on the displacement of the interface by consideringthe difference of SE at both sides of the interface. We will adopt a slightlymodified methodology to the one presented in [34] to quantify it: (cid:126)v e = − M δ (˙ (cid:15) ) [ E ] ij (cid:126)n, (3)where the term [ E ] ij defines the difference of stored energy E between thegrains i and j ( E i − E j ), the term δ (˙ (cid:15) ) is a mobility coupling factor whose nature7s explained in [36] appendix c and where the direction of the unit normal (cid:126)n sets, for a given node of the interface, the order of the indices as: first theindex i and then j . Note that this definition holds even if the direction of (cid:126)n isambiguous (in the case of a flat interface with no convex side) as the directionof the velocity (cid:126)v e will be then pointed, in all cases, from the lower to the highervalue of stored energy no matter what the direction of (cid:126)n is. Moreover the valueof stored energy can be computed using the equation: E = 12 µb ρ, (4)where b corresponds to the norm of the Burgers vector and µ correspond tothe elastic shear modulus of the material.Finally, the contribution of driving pressures due to SE and capillarity canbe accounted by linearly adding the two velocities as in [34, 27]: (cid:126)v = − M ( δ (˙ (cid:15) ) [ E ] ij (cid:126)n + γκ(cid:126)n ) , (5)where (cid:126)v denotes the final velocity of the interface during GBM when SEeffects are included. Equation 2 can only be used in a one-boundary problem, as in a more gen-eral context, the presence of multiple junctions (the intersection points of morethan 3 interfaces) makes it impossible to compute a curvature κ or a normal (cid:126)n at these points. As explained in previous works [22, 23], we have used analternative methodology to compute the velocity due to capillarity at multi-ple points: Model II of [11], where the product κ(cid:126)n is directly obtained from anapproximation of the free energy equation of the hole system in a vertex context.Similarly, Eq. 3 only holds in a one-boundary problem as neither the valueof [ E ] ij nor the value of (cid:126)n can be obtained at these points. To solve this, a differ-ent approach has been developed to compute a “resultant” velocity due to storeenergy (cid:126)v e at multiple junctions. This approach is illustrated in Fig. 2 where forthe sake of clarity, the value of M has been held constant and equal to 1. Fig. 2a) shows a typical configuration where the boundaries of three grains convergeto a single point, each grain i has its own stored energy E i where E > E > E .The values of the velocity for each normal boundary have been computed withEq. 3 and are shown as white arrows for each node in the boundary of Fig. 2 b),here the index on the normal (cid:126)n ij term are only representative of their directionand serve to set the indices of each [ E ] ij terms, these expressions do not followthe Einstein notation summation laws, all summations will be represented bythe conventional Σ operator. A mobility coupling factor function of the effective strain rate ˙ (cid:15) with δ (˙ (cid:15) ) = 1 when ˙ (cid:15) = 0. . . . S t o r ed E ne r gy Figure 2: Graphical demonstration of the obtention of Eq. 6, a) typical triplejunction configuration with values of SE homogenized on each grain and thenormal vectors n computed at the nodes of the interfaces pointing to theirconvex side, b) computation of the term − [ E ] ij · (cid:126)n ij for each node of the interfaceexcept for the node at the triple junction, c) definition of the same configurationas in a) but in a differential portion of radius dr , d) the resultant driving forcesare applied at the center of the segments on the differential portion, e) and f)the driving forces are distributed at the ends of each segment and an expressioncan be formulated at the triple junction for its resultant driving force.If a portion of differential size dr centered at the multiple point is evaluated(see Fig. 2 c)) the boundaries between grains will appear as flat, here the9ifference on the stored energy can be seen as a distributed difference of potential[ E ] ij applied on the length of the grain boundary of size dr (analog to a givenpressure acting as a resultant force on a given interface). A normal (cid:126)n (cid:48) ij canbe obtained and used to compute a velocity of each boundary (2 d)) appliedat its center. Note that the direction of (cid:126)n (cid:48) ij can be chosen ambiguously onthis linear segment, however, as mentioned before, an eventual ambiguity onthe direction of (cid:126)n do not represent an ambiguity on the term − [ E ] ij · (cid:126)n ij as[ E ] ij · (cid:126)n ij = [ E ] ji · (cid:126)n ji with [ E ] ji = − [ E ] ij and (cid:126)n ji = − (cid:126)n ij . These velocities canbe divided and applied at the ends of each boundary and finally added at thejunction point (2 e) and f) respectively) to obtain a valid velocity vector fieldat multiple junctions. The expression on Fig. 2 f) can be extended to the casewhere the values of M are neither constant nor equal to 1: (cid:126)v e = − Σ M δ (˙ (cid:15) ) [ E ] ij (cid:126)n , (6)of course, this expression can be also used in cases of multiple junctionsof any order, where more than three interfaces meet. Eq. 6 will be used tocompute the value of v e at multiple junctions as an approximation to the yetunknown behavior of such configurations under the influence of stored energyin a transient state. Multiple changes on the topology of the microstructure occur during GG andReX. In general, the topological changes during GG are given by the disappear-ance of grains: on a shrinking grain, each of their boundaries evolve until they collapse to multiple junctions. Eventually, all boundaries collapse to a singlemultiple junction and the domain occupied by the grain disappears. This behav-ior was implemented on the original TRM model presented in [22] by means ofthe application of the selective node collapse operator, where some restrictionswhere made regarding the order of collapsing.In [22] we had opted to use this methodology to produce coherent topologicalchanges on the microstructure, leading to a series of rules on the selective nodecollapse operator (see section 2.4.1 of [22]). These rules gave to the P-Nodesa higher influence over other kind of nodes and prevented the collapsing ofnon-consecutive nodes as illustrated in Figures 3 and 4 respectively.The implementation of such node collapsing strategy allows a high controlover the order on which the topological changes occur, unfortunately this kindof reasoning can only be used on isotropic GG and can not be used when SEor spatial heterogeneities of the mobility/interface energy must be taken intoaccount.When considering SE, the kinetics of the GB are not only led by the move-ment of multiple points; flat surfaces can evolve with a given velocity and it ispossible that the velocity of simple boundaries becomes much more importantthan the velocity of multiple junctions. Fig. 5 illustrates this behavior with six10 i NjCollapsing zone of Ni
Figure 3: Node Collapsing rule in GG by capillarity. Some nodes are within thecollapsing zone of N i : Two S-Nodes (yellow) will collapse, one L-Node (blue)cannot collapse and one P-Node (red) cannot collapse. N i cannot collapse P-Nodes (topological degree), N i can collapse L-Nodes but N j does not belong tothe same line (i.e. they do not belong to the same grain boundary). [22] N i N a Collapsing zone of Ni N b N c N d Figure 4: Node Collapsing rule in GG by capillarity. Some nodes are within thecollapsing zone of N i : four L-Nodes (blue) N a , N b , N c and N d . Only Nodes N b and N c can collapse as they are consecutive to N i within the same line. [22]grains with a specific SE state. The circular grain in the middle grows due toits low SE compared to the SE of its surrounding grains. The circular grainis indeed surrounded by an initially squared grain that starts shrinking by thecombined effects of capillarity at their external boundaries and the surface takenaway by the circular growing grain. Fig. 5 (right), shows the moment when theboundary of the circular grain and the external boundaries of the initially squaregrain collide, unchaining a series of topological changes on the microstructure.These topological changes are illustrated in Fig. 6, where new multiple junc-tions (Points) appear, grain boundaries (Lines) are split and grains (Surfaces)are divided.These several changes on the microstructure (contact of different grain bound-aries in non-convex grains) can not be accomplished by the TRM model if therules described above and illustrated by figures 3 and 4 are maintained. This iswhy these two rules need to be overridden and a new condition implemented:If two non-consecutive nodes (nodes not connected by the microstructural wire-11 . . . S t o r e d E n e r gy Figure 5: Six grains with a specific SE balance, the circular grain in the middlegrows due to its low SE compared to the SE of its surrounding grains. Theinitially squared grain shrinks by the combined effects of capillarity at theirexternal boundaries and the surface taken away by the circular growing grain.Left: initial state, center: the circular grain grows, right: the boundary of thecircular grain and the external boundaries of the initially square grain collide.Figure 6: Details of the final event of Fig. 5, highlighting the changes on themicrostructure. Here, Points describe multiple junctions, Lines grain boundariesand Surfaces grains.frame) collapse, the classification of the remaining node is a P-Node. Addi-tionally, the remaining node is moved to the barycenter of the initial nodesinvolved in the collapse and the surrounding geometrical entities (points, linesand surfaces) are checked and updated if necessary. Take for example the sameconfiguration shown in Fig. 3 now in Fig. 7: here the collapsing of nodes N i and N k is possible, the remaining node N i is placed in the middle of the edge N i N j and its classification is changed from L-Node (blue) to P-Node (red). Nowtheir surroundings need to be checked for possible changes on the topology: allthree Lines (grain boundaries in green) need to erase node N k as a final/initialpoint and put in its place P-Node N i , similarly, one of the lines has to add N k as a node in their sequence of L-Nodes hence N k changes also its classificationto L-Node.Similarly, the situation presented in Fig. 4 can be reproduced with the new12igure 7: Node Collapsing of Fig. 3 when allowing the collapse between nonconsecutive nodes, S-Nodes are displayed in yellow, L-Node in blue and P-Nodesin red. The collapsing of nodes N i and N j produces a new P-Node ( N i ) whilethe pre-existent P-Node N k needs to be reclassified as a L-Node. Left: initialstate, right: state after collapse.rules of collapsing: on the initial state (Fig. 8a)), the collapsing zone of N i putsL-Nodes N a , N b , N c and N d inside. The first node to be collapse is L-Node N a . After this initial collapse (Fig. 8b)) the collapse produces L-Node N i to bemoved to the center of the edge N i N a and to become a P-Node. Moreover, thecollapsing zone of N i changes its position and leave L-Node N c out hence it willnot be collapsed. Additionally, a new topological change is identified (Fig. 8c)), here the collapse is also responsible of the fact that a Surface is divided intwo new surfaces (cyan and orange Surfaces). Finally, the remaining collapsesare performed between P-Node N i and L-Nodes N b and N d (Fig. 8d)), thesecollapses are performed in the conventional way following the rules of collapsingbetween P-Nodes and L-Nodes presented in [22].The new node collapsing rules have been implemented in the TRM modeland will be used from this point forward in the cases where SE is present. Thisnode collapsing technique will be able to perform the majority of the topologicalchanges in the microstructure. However not all topological changes can be han-dled by this node collapsing mechanism, indeed one special topological changeneeds to be handled differently: The creation of boundaries by the decompo-sition of unstable multiple junctions with more than 3 intersected boundaries.This phenomenon has been addressed in [22] section 3.2.4 where isotropic grainboundaries were considered, by the implementation of a new operator used onthe mesh where such a configuration appears. This operator is responsible forthe successive dissociation of grain boundaries from multiple junctions until astable configuration is obtained. Here we will use the same algorithm denomi-nated “Split of multiple junctions”. 13igure 8: Node Collapsing of Fig. 4 when allowing the collapse between nonconsecutive nodes, S-Nodes are displayed in yellow, L-Node in blue and P-Nodesin red. a) Initial state; b) L-Node N i performs the first collapse with N a . Thisproduces L-Node N i to be moved to the center of the edge N i N a and to become aP-Node. Moreover, the collapsing zone of N i changes its position and leave Node N c out; c) The collapse also has produced a Surface to be divided in two (cyanand orange Surfaces); d) Additional collapses are performed between P-Node N i and L-Nodes N b and N d , these collapses are performed in the conventionalway without moving N i as it is a P-Node now. Finally the algorithm for a time step of the TRM model in the context ofisotropic grain growth under the influence of stored energy and capillarity ispresented in Algorithm 1, where the step “Perform Remeshing and ParallelSequence” corresponds to the parallel implementation of the TRM model pre-sented in [23].
In order to model Recrystallization (ReX) with the TRM model, two additionalcomponents are necessary: the first is a procedure allowing to change the topol-14 lgorithm 1
Isotropic Grain Growth TRM Algorithm for capillarity and storedenergy Perform Remeshing and Parallel Sequence for all Points: P i do while Number of Connections > do split multiple point P i . for all Lines : L i do Compute the natural spline approximation of L i . for all L-Nodes : LN i do Compute curvature and normal ( κ(cid:126)n ) over LN i then compute (cid:126)v c for LN i (Eq 2). Compute the (cid:126)v e for LN i (Eq 3) for all P-Nodes :
P N i do Compute the product κ(cid:126)n over
P N i using model II of [11] then compute (cid:126)v c for P N i (Eq 2). Compute (cid:126)v e for P N i (Eq. 6) Delete Temporal Nodes for all
L-Nodes and P-Nodes :
LP N i do Compute final velocity (cid:126)v of Node
LP N i (Eq. 5) Iterative movement with flipping check in parallel ogy of the microstructure and to introduce new grains (i.e. nuclei); the secondcomponent is a model of apparition of nucleus which depends thermomechani-cal conditions. Here discontinuous dynamic ReX (DDRX) context is considered.Of course post-dynamic ReX (PDRX) and subsequent GG phenomena can alsobe investigated by considering microstructure evolutions when the deformationis completed. The combination of these two mechanisms can describe multipleTMTs that are used today in the material forming industry.With the purpose of simplicity, in this article we will use the same method-ology presented in [36] for the laws governing the introduction of new nucleiduring the modeling of hot deformation:In [36] the evolution of the dislocation density is accounted by a Yoshie-Laasraoui-Jonas Law [44] as follows: ∂ρ∂(cid:15) peff = K − K ρ, (7)which can be evaluated in a discretized time space with an Euler explicitformulation for the next increment step as : ρ ( t +∆ t ) = K ∆ (cid:15) + (1 − K ∆ (cid:15) ) ρ ( t ) , (8)where ρ ( t ) is the value of the dislocation density at time t and where thevalue of ∆ (cid:15) can be computed as ˙ (cid:15) peff · ∆ t with ∆ t the time step.15s explained in [36], when a grain boundary migrates, the swept area isassumed almost free of dislocations. This aspect is modeled by attributing tothese areas a value of dislocation density equal to ρ , then, for the grains withpart of their domain presenting ρ , their dislocation density is homogenizedwithin the grain (as intragranular gradients on the stored energy are not takeninto account neither in [36] nor in the present work), the final value of ρ for thegrowing grains is computed as: ρ ( t +∆ t ) = ρ t S ( t ) + ∆ Sρ S ( t +∆ t ) , (9)where S t and ∆ S denote the surface at time t and the change of surface ofa given grain.Additionally to Eq. 9, in PDRX the annihilation of dislocations by recoverymust be taken into account. This is done thanks to the following evolution law: dρdt = − K s ρ, (10)where K s is a temperature dependent parameter representing the static re-covery term. This recovery law is only taken into account in PDRX as in DRX,Eq. 9 already takes into account the annihilation phenomenon. The procedure consists of introducing volume (surface in 2D) of nuclei at a rateof ˙ S , once the local value of dislocation density has reached a critical value: ρ c .In [27, 36, 45] this value was determinated by iterating until convergence thefollowing equation: ρ ( i +1) ∗ c = − bγ ˙ (cid:15) K M δ (˙ (cid:15) ) τ ln (1 − K K ρ ic )
12 with ρ ic = ρ i − c + c · ( ρ ( i ) ∗ c − ρ i − c ) , (11)where i represents the iteration number, c is a convergence factor ( c < c = 0 . K and K represents the strain hardeningand the material recovery terms in the Yoshie–Laasraoui–Jonas equation dis-cussed in [44], the term b = 1 in 2D and b = 2 in 3D, τ is the dislocation lineenergy and ˙ (cid:15) is the effective deformation rate used during the deformation ofthe material.When solving this equation, two special cases may produce an erroneouscomputation: the first is given when K /K ∗ ρ c > ρ c < K /K whenever this16ituation occurs. The second is when ˙ (cid:15) = 0 which correspond to the intervalswhere PDRX is considered. Two solutions may be considered for this situation:the first is to block the nucleation when it is not necessary (metadynamic evo-lution for example), and the second to supply value of ˙ (cid:15) > (cid:15) s is usedinstead ˙ (cid:15) in PDRX: ˙ (cid:15) s = (cid:82) t ˙ (cid:15) dt (cid:82) t ˙ (cid:15) dt , (12)which accounts for the instant mean value of the real effective strain rate.Once a value of ρ c is computed, the surface per unit of time ˙ S of nuclei tobe inserted can be computed with the following equation corresponding to avariant of the proportional nucleation model [46]:˙ S = K g P c , (13)where the term K g is a probability constant depending on the processingconditions and P c is the total perimeter of the grains whose dislocation densityis greater than ρ c .Another constraint is given by the minimal radius r ∗ of nucleation (the radiusat which the nuclei should be inserted in the domain so the capillarity forceswould not make it disappear) which can be computed thanks to the followingequation [47]: r ∗ = ω γ ( ρ c − ρ ) τ , (14)where ω > ω accounts for the non spherical shape ofa grain inserted in a discretized domain such as in the TRM model. In section5.1 a value for this factor will be obtained based on numerical tests. Till here we have defined the tools needed to obtain the kinetics of the grainboundaries, where the pressure behind such kinetics can be of different nature:capillarity, stored energy or both. However, in order to model ReX it is nec-essary to have a way to introduce new grains into the domain of the TRMmodel. Nucleation, similarly to boundary migration, is one of the ways of themicrostructure to relax the high gradients of the stored energy appearing duringor after a TMT. Nucleation has been addressed by several approaches for eachmethodology able to simulate such behavior: LS-FE methods, relay on the defi-nition of circular LS fields (different from the already defined LS fields occupyingthe same spatial domain) to form nuclei [48, 36], CA and MC methods changethe crystallographic orientation and SE value of some cells [49, 50, 51, 52] in17rder to nucleate and vertex models form new grains by redefining new vertexand interfaces in the shape of triangles around the pre-existent vertex points [15].In the present work, a remeshing-reidentification procedure will be performedaround a central node N i in order to introduce nuclei. A circular region withcenter N i will be drawn and all edges crossed by this circle will be split at theintersection in a similar manner as in [24, 25] by successively applying an edgesplitting operation, regardless of the classification of the nodes defining the edge(P-Node/L-Node/S-Node) (see Fig. 9). The classification of the new nodesbeing placed by the splitting algorithm are as L-Nodes unless the split edgerepresents a grain boundary, in which case the inserted node will be classifiedas a P-Node (see Fig. 10 b) middle and c) middle). Once all edges are split,a Surface Identification algorithm will be performed over node N i (see section2.2.5 of [22]), all identified elements and nodes will be inserted into a new emptySurface defining the nucleus, and extracted from their previous Surfaces (grains),new Lines (grain boundaries) will be built with their respective Points (multiplejunctions) if any were formed by the nucleation process and all remaining linesand points inside the new surface will be destroyed (remaining lines and pointscan appear if the nucleation took place near over one grain boundary) see Fig.10 a), b) and c) left. However in a parallel context, an additional constraintwas added: shared nodes can not be involved in the nucleation process, neitheras a central node nor one of the nodes of a split edge. This constraint wasadded because of the lack of information (position of the edges to cut) aroundshared nodes and the performance of the nucleation process (as a great amountof information would be necessary to be transferred to other processors). Thisconstraint should not have a great impact on the general behavior of the modelas the domain of each processor (and their shared nodes) is changed by the Unidirectional Element Sending algorithm presented in [23] every time step,hence constantly unblocking the restriction to nucleate over the same region.18igure 9: Remeshing steps for the nucleation process of the TRM model. Top:initial state with a selected node (cyan) and a circle drawn over the mesh,middle: successive edge splitting steps to form the interfaces of the nucleus,bottom: the elements inside the nucleus are identified and extracted from itsprevious Surface container. 19igure 10: Examples of the formation of nucleus over different types of nodes,a) S-Node at the center and no crossed lines, b) L-Node at the center and onecrossed line, two P-Nodes are created, the initial line is divided and two newlines are created, c) P-Node at the center, 3 P-Nodes are created, the P-Nodeat the center is detached from all its lines and converted to S-Node, 3 new linesare created. 20
Numerical Tests
In this section different academic tests will be performed to evaluate the per-formance of the TRM model when simulating GBM under the influence of cap-illarity and stored energy. For these academic tests, adimensional simulationswill be considered. Moreover, results of simulations using the DRX and PDRXframeworks described in section 4 will be given, these simulations will use thenucleation approach presented in section 4.2 specially developed for the TRMmodel. The different physical parameters will be take as representative of the304L stainless steel. Comparisons with LS-FE predictions will be discussed.
In this test case, it will be evaluated the accuracy of the model when the geo-metric configuration leads to a competition between the driving forces given bythe capillarity and the SE. Here we will adopt a value of boundary energy andmobility equal to γ = 1 and M = 1 respectively. A circular domain with a valueof stored energy E = α is immersed in a squared domain with an attributedvalue of stored energy E = β ( see Figure 11 a) and b) left) where β > α . Thedifference on the stored energy [ E ] = β − α at the boundary will try to makethe circle expand at a rate v e = M [ E ] = [ E ] while the capillarity effect willtry to make it shrink at a rate v c = M κ = κ where κ is the local curvature.The analytical model for this configuration can be put in terms of a non-linearordinary differential equation in terms of the radius r of the circle as follows: drdt = − r + [ E ] , (15)or in terms of the surface S of the circle: dSdt = 2( − π + √ πS [ E ]) , (16)We have used an Euler explicit approach to solve this equation and the re-sults are used to compare the response of the TRM model for different values of[ E ] for two cases: the first is given for an initial radius of r = 0 . r = 0 .
025 (see Figure 11 a).left and b).left). The initial mesh for each oneof the two cases is given in Figure 11 a).right and b).right respectively. Notehow in the first case, the initial circle boundary is discretized by an amount ofnodes sufficiently capable of capturing precisely the value of its curvature, henceit will serve to evaluate the accuracy on the kinetics of a typical curved bound-ary, while in the second case, the circle boundary is only defined by a few nodesallowing to evaluate the behavior of a nucleus when it is inserted on the domain.Results for this first case are given in Fig. 12 along with the solution of Eq. 16for different values of [ E ], Fig. 12.left illustrates how the references curves aresuperposed to the different simulated curves with a very low error (around 221igure 11: Circle Test, left: initial state and right: initial mesh a) r = 0 . . r = 0 .
025 radius (Surface=0 . E ] = 10 /
3) shows a very good behavior losing only 1 .
2% of its surface at t = 0 . E ] = 0 , , , E ] = 40 which corresponds analytically tothe metastable configuration. In TRM simulation, the grain disappears. Thisbehavior is due to the low number of nodes at the interface, producing anoverestimation on the computed value of its curvature, making it shrink fromthe very first increment. A value of [ E ] ≈ , .
2% accordingly to itsanalytical value). Moreover, for this value, the error on the prediction of theevolution of the surface was also the highest, going up to 92% after t = 0 . E ] = 50 , , ,
80 (for which the higher driving force is the stored energy)22igure 12: Evolution of the surface (left) and L2 error (right) for the circle testcase for an initial circle radius r = 0 . . h = 0 .
006 anda delta time dt = 3 e −
5, the analytical results (References) are shown superposedto the simulated curves in black dashed lines. The expected metastable curveis given for a [ E ] = 10 / E ] increases. This result can be usedon the determination of factor ω used in Equation 14 where the authors haveestimate that a value of ω = 1 . E ] for a metastable state, see the curve [ E ] = 60 infigure 13.) is sufficient in order to give the inserted nuclei a growing state andprevent its early disappearance.Figure 13: Evolution of the surface (left) and L2 error (right) for the circle testcase with an initial circle radius of r = 0 .
025 (Surface=0 . h = 0 .
025 and a delta time dt = 1 e −
5, the analytical results (References)are shown as dashed lines of the same color of their corresponding simulatedevolution. The expected metastable curve is given for a [ E ] = 40 (orange curve),but metastability was found for [ E ] = 48 .
491 (Red curve)23 .2 Triple junction : The capillarity effect on the quasi-stable shape of multiple junctions
In [53, 54] analytic solutions for the movement of multiple junctions in a quasisteady-state under the influence of stored energy were presented. In [53] the socalled ”Vanishing Surface Tension” (VST) test was introduced to demonstratethe non-uniqueness of the solution presented in [54] hereafter called the ”Sharp”solution, this test (the VST test) takes the form of the limit problem given by: (cid:126)v · (cid:126)n = − M ([ E ] ij + (cid:15)γκ ) , with (cid:15) → , (17)which has subjected to several 2D test cases and a perturbation analysisto demonstrate that the VST solution correspond to one of the solutions when (cid:15) = 0 and to the unique solution otherwise.These solutions were later studied in [17, 34] using a LS-FE model to obtainthe same behavior both in 2D and 3D. Here we have reproduced with the TRMmodel two tests that show the same behavior as in [53, 17, 34] for the 2D so-lutions. For all test the ”Sharp” solution was obtained when capillarity effectswhere taken into account (with (cid:15) = 1) and the VST solution when no capillaritywas introduced in the system (hence with a value of (cid:15) = 0). Furthermore wehave developed analytic equations for the evolution of the growing surface inour specific case (see Figure 14), these analytic evolutions are valid up to thepoint of contact of the multiple junction with the lower edge of the equilateraltriangle (the limits of out domain) and allow us to make a more quantitativecomparison in terms of error.Figure 14: Initial state for the triple junction test, three phases immersed in adomain in the shape of an equilateral triangle, this shape is intended to maintainan orthogonal position of the boundaries with respect of its limits while theconfiguration evolves. a) Initial configuration and b) initial mesh.For this test case, the initial conditions are those presented in figure 14.left,three phases immersed in a domain in the shape of an equilateral triangle, this24hape is intended to maintain an orthogonal position of the boundaries withrespect of its limits while the configuration evolves. Two of the phases (the twoin the lower part of the domain) will have a constant value of stored energy of α and the third phase a value of β < α , this configuration will produce a globalmovement of the triple junction downwards at a constant and normal velocityof the flat interfaces equals to α − β . Eventually the triple point will reach thebottom part of the domain making it to split and evolve towards a lower energystate; even though this portion of the simulation is showed in some of the resultsit is not relevant to our study, hence we will give quantitative results up to thepoint of splitting. The initial mesh for every test performed is shown in figure14.right corresponding to a mesh size parameter of h = 0 . γ = 1 and M = 1respectively.The analytic solution for the evolution of the surface of the upper phase (thegrowing phase) for the Sharp solution is given by: S Cap = (cid:18) a √ − y (cid:19) √
34 (18)where a is the length of one of the sides of the equilateral triangle (here a = 1) and y is the vertical position of the triple junction measured from thebase of the triangle and given by the following expression: y = a √ − | (cid:126)v · (cid:126)n | t √ t is the time and the expression | (cid:126)v · (cid:126)n | is the instant normal velocityof the flat phase boundaries, i.e. α − β .Similarly the analytic response of the VST solution in terms of surface forthe growing phase is given by S NoCap = S Cap + (cid:18) π − √ (cid:19) ( | (cid:126)v · (cid:126)n | t ) . (20)Two test were performed: one with β = 2 and α = 4, i.e. | (cid:126)v · (cid:126)n | = [ E ] = 2and one with β = 10 and α = 20, i.e. | (cid:126)v · (cid:126)n | = [ E ] = 10. The two tests wereperformed with a time step ∆ t = 1 · − . Results for the evolution of the meshand the surface are given in figures 15 and 16 for the first and the second caserespectively. It is clear that the accuracy on the scalability of the solution isvery good as figures 15 a), b) and c) are almost equal to the ones of figures 16 a),b) and c) respectively which were obtained for a velocity 5 times higher. Notethat the only different frame is given for Figures 15 d).right and 16 d).right ashere the capillarity effects over the limits of the domain are not negligible andin Figure 15 d).right the configuration have had 5 times more time to evolve toits given state. 25igure 15: States for the triple junction test case with a value of [ E ] = 2, left:with (cid:15) = 0 and right: with (cid:15) = 1 at a) t = 0 .
02 b) t = 0 .
04, c) t = 0 .
06, d) t = 0 . . E ] = 10, left:with (cid:15) = 0 and right: with (cid:15) = 1 at the instant a) t = 0 .
004 b) t = 0 . t = 0 . t = 0 . Here a simulation with a few initial grains will be performed using the recrystal-lization method mentioned in section 4: the initial tessellation will be realizedthanks to a Laguerre-Voronoi cells generation procedure [55, 56, 57] over a rect-27igure 17: Evolution of the surface of the growing phase of the triple junctiontest case, from top to bottom: (top) Evolution of the surface, (center) Zoom inthe red zone, (bottom) L2-error over the evolution of the surface. Left: resultsfor the test with [ E ] = 2 and right: with [ E ] = 10angular domain of initial dimensions 0 . × . mm (see figure 18) and thevalues for M , γ , τ and k s are chosen as representative of a 304L stainless steelat 1100 ◦ C (with M = M ∗ e − Q/RT where M is a constant M = 1 . · mm /Js , Q is the thermal activation energy Q = 2 . · J/mol , R is the idealgas constant, T is the absolute temperature T = 1353 K , γ = 6 · − J/mm , τ = 1 . · − J/mm and k s = 0 . s − [27, 36] ). Additionally, theparameters K , K , K g and δ are taken as dependent of the absolute value ofthe component xx of the strain rate tensor ˙ ε ( | ˙ ε xx | ) which is defined as corre-28igure 18: Initial State for the DRX/PDRX test case.sponding to a plane deformation case. These parameters will be obtained usinga linear interpolation of the values presented in Tab. 1.Table 1: Parameter data table for the DRX PDRX test case, when in range | ˙ ε xx | = [0 . , . s − the values are interpolated. If | ˙ ε xx | > . | ˙ ε xx | = 0 . s − , the samestrategy applies when | ˙ ε xx | < . | ˙ ε xx | s − K mm − K K g mm · s − δ · · − · · − | ˙ ε xx | = 0), the parameter δ will take the valueof 9 .
18 following the findings in [58]. Also, as explained in section 4.1, duringPDRX the parameter ρ c will be computed using the apparent effective strainrate ˙ (cid:15) s (see Eq. 12 and Fig. 20.right) instead of the effective strain rate ˙ (cid:15) (equalsto 0 in this regime). Finally, outside the range of interpolation, the values arecomputed as follows: if | ˙ ε xx | > . | ˙ ε xx | = 0 . s − , similarly, the same strategy applies when | ˙ ε xx | < .
01, using the values for | ˙ ε xx | = 0 . s − (see Fig. 19 for an illustrationof the values of K K Kg and δ in function of | ˙ ε xx | ).Four cycles of deformation/coarsening will be applied as illustrated in Fig.20, Fig. 20.right shows the computed values of the effective strain rates ˙ (cid:15) s and ˙ (cid:15) ,while Fig. 20.left shows the strain deformation component ˙ ε xx , where multiplemarkers have been displayed, these markers correspond to different states alongthe simulation that will be useful when analyzing the results.Statistical comparisons of the TRM model and the response obtained by a29igure 19: Evolution of the parameters of table 1 in function of | ˙ ε xx | .Figure 20: Deformation loading strategy for the DRX and PDRX case: right:the computed values of the effective strain rates ˙ (cid:15) s and ˙ (cid:15) , left: the strain defor-mation component ˙ ε xx , where multiple markers have been drawn, correspondingto different states during the simulations.FE-LS approach presented in [17, 18, 36, 59] will be given. This approach uses amore classic method of mesh adaptation during calculations where the interfacesare captured with an anisotropic non-conform local refined mesh. This method-30logy will be denoted in the following as the Anisotropic Meshing Adaptation(AMA) model.A well known behavior of full field simulations of microstructural evolutionsis that the reduced mobility ( γM product) is classically impacted by the choiceof the numerical method and is not only a universal physical parameter. In otherwords, the reduced mobility is a physical parameter that needs to be identifiedcomparatively to experimental data. This identification may lead to differentvalues depending generally of the numerical method used [52]. In [22], the re-duced mobility was adjusted in order to minimize the L2-difference between themean grain size evolution curves considering TRM or AMA numerical strate-gies. Same methodology was used here in the global thermomechanical pathsleading to an increase of 40% in the optimal identified reduced mobility.Here, we have chosen the AMA case as a reference even though there is noway to know which model gives the most accurate response to the given physicalproblem; this choice on the other hand is given as an example of how the TRMmodel can indeed obtain similar responses to well established models in the fieldof microstructural evolutions.Multiple microstructural states have been retrieved from the results givenby the TRM model. This states are marked with numbers corresponding to thestates of figure 20:States 1 to 3 are given in figure 21, here state 1 illustrates the apparitionof the firsts nuclei in the positions where the dislocation density field reachesits value ρ (cid:62) ρ c . Then state 2 gives the end of the first stage of deformationwhere more nuclei have appeared, note that the value of the stored energy insome of the small grains is different from others, these grains have been presentlonger in the domain and consequently have been subjected to strain hardening,contrary to the nucleus that have appeared later, during or at the end of thisdeformation stage. Finally state 3 shows the end of the first grain coarseningstage where nuclei have been given time to growth as a product of the highdifference in energy with their surroundings.States 4 to 6 are presented in Fig. 22. In stage 4 only a small percent of thedomain have a dislocation density of at least ρ c and nucleation is restricted tothese zones, contrary to stage 5, where a bigger part of the domain have reachedthe value of ρ c , consequently new grains appear everywhere. Finally, the end ofthe second deformation stage is given in state 6 where the first peak of numberof grains is reached (4250 grains).States 7 to 10 are given in figure 23, these steps are representative of theends of the third and fourth deformation/coarsening cycles, where during thedeformation the nucleation process increases the number of grains while in thegrain coarsening stages the high value of the parameter δ (2.245 to 9.18 its dy-31igure 21: States 1 to 3 (see figure 20) obtained with the TRM model. 1: thefirsts nuclei appear, 2: end of the first stage of deformation, 2: end of the firstgrain coarsening stage.namic vs its static value) makes the grain number decrease rapidly (see figure28.a) ) for the evolution of the number of grains).States 11 to 14 are provided in figure 24, these states correspond to a valueof time t of 30, 40, 50 and 60 seconds respectively. In this range of time no de-formation is considered. Note how the limits of the scale in figure 24 change asa product of the disappearance of high energetic grains and to the annihilationof dislocations simulated through equation 10.Statistical values for the states 4 to 6 and 11 to 14 are given in figures 26and 27 respectively. The grain size distributions for the TRM model without amobility increase and with a mobility increase of 40% have been plotted alongwith the response given by the AMA case. Similarly the evolution of the meangrain size are provided for all simulations in figure 25.left and the L2-differenceto the AMA case is given in figure 25.rightFinally, the evolution of some representative values are given in Fig. 28: theevolution of the number of grains, the recrystallized fraction, the mean value of32igure 22: States 4 to 6 (see figure 20) obtained with the TRM model. 4: thenuclei appear on the regions where ρ > ρ c , 5: all the domain is now above thevalue of ρ c hence the nucleation occurs everywhere, 6: end of the second stageof deformation, here the maximum number of grains is reached (4250 grains). ρ pondered in surface ( ρ ) and the total perimeter of the grains whose dislocationdensity is greater than ρ c ( P c ) are provided.These results show a good agreement between the general behavior of thethe TRM model and the behavior of the AMA simulation when an increase of40% is considered to the reduced mobility M γ value (following the findings in[22]). The computational cost for the different iterations of the TRM modelis given in figure 29, where for the slower simulation the time needed for itscompletion was of 25 min and for the fastest of 20 min, compared to the timeneeded for the AMA case (4 hours and 38 min).33igure 23: States 7 to 10 (see figure 20) obtained with the TRM model. 7: endof the second grain coarsening stage, the number of grains drops very quicklygiven by the increase of the value of δ from 2.245 to 9.18 (its dynamic vs itsstatic value) , 8: end of the third deformation stage, 9: end of the third graincoarsening stage, 10: end of the fourth deformation stage.34igure 24: States 11 to 14 (see figure 20) obtained with the TRM model. thesestate correspond to a value of time t of 30, 40, 50 and 60 seconds respectively.Figure 25: Evolution of the Mean grain size (left) for the TRM model and theL2-difference with the AMA simulation (right).35igure 26: Grain size distributions pondered in surface for the states 4 to 6(example of a deformation Stage). A peek on the nucleus size can be observed.36igure 27: Grain size distributions pondered in surface for the states 11, 13and 14 (example of a grain coarsening stage). The values are distributed moreevenly on the size range (x axis) as a product of the grain growth.37igure 28: different values as a function of time for the DRX and PDRX testcase, a) Number of grains, b) Recrystallized fraction, c) Mean value of ρ pon-dered by surface, and d) Critical perimeter for the computation of the nucleationrate in equation 13 38igure 29: CPU-time for the different simulations using the TRM model, thecomputational cost drops as the number of simulated grains decreases.39 Discussion, conclusion and perspectives
In this article the TRM model presented in previous works in the context ofisotropic grain growth by capillarity has been adapted in order to take intoaccount bulk terms due to the stored energy during plastic deformation. Thisadaptation has made possible the integration of a recrystallization model to theTRM approach, for which a nucleation procedure has also been presented.The algorithms presented in section 3.1 and represented by Eq. 6 for thecomputation of the velocity at multiple junctions, although intuitive have notbeen published before to the knowledge of the authors, only [37] shows a similar(more indirect) approach in the context of vertex simulations.Results for the circle test case and tripe junction case have demonstrated thehigh accuracy of the TRM model in the modeling of boundary migration dueto capillarity and stored energy, where in the normal context (for typical grainboundaries and multiple junctions), an error no greater than 2% was found.Also, the circle test case showed the typical behavior of a nucleus when sub-jected to a wide range of stored energy around its metastable point and helpeddefine the safety factor ω used in Eq. 14, defining the minimal radius to nucleatein the context of the TRM model.Finally, a DRX/PDRX test case was considered in order to test the recrys-tallization model provided in section 4 for 304L stainless steel at 1100 ◦ C . Areference test case using the same ReX model but with a FE-LS strategy wasalso considered (AMA case) [17, 18, 36, 59]. Following the findings in [22], anoptimal reduced mobility was calibrated to performed the tests (40% higherthan the mobility used in the AMA case context). Results shows a very goodagreement between the two models. Moreover the computational cost of theTRM model was lower being between 20 to 25 minutes against the 4 hours and38 minutes needed for the AMA case for its completion.Perspectives for the presented work and the TRM approach include the im-plementation of a model able to treat full anisotropic boundary properties, aswell as the study of the in-grain gradients of stored energy. The 3D implemen-tation of the TRM model will also be studied in future works. Acknowledgments
The authors thank the ArcelorMittal, ASCOMETAL, AUBERT & DUVAL,CEA, FRAMATOME, SAFRAN, TIMET, Constellium and TRANSVALORcompanies and the ANR for their financial support through the DIGIMU con-sortium and ANR industrial Chair (Grant No. ANR-16-CHIN-0001).40 ata availability
The raw data required to reproduce these findings cannot be shared at this timeas the data also forms part of an ongoing study. The processed data requiredto reproduce these findings cannot be shared at this time as the data also formspart of an ongoing study. 41 eferences [1] A. D. Rollett, D. J. Srolovitz, M. P. Anderson, Simulation and theory ofabnormal grain growth-anisotropic grain boundary energies and mobilities,Acta Metallurgica 37 (4) (1989) 1227–1240. doi:10.1016/0001-6160(89)90117-X .[2] A. D. Rollett, D. Raabe, A hybrid model for mesoscopic simulation ofrecrystallization, Computational Materials Science 21 (1) (2001) 69–78. doi:10.1016/S0927-0256(00)00216-0 .[3] D. Raabe, Cellular automata in materials science with particular refer-ence to recrystallization simulation, Annual Review of Materials Science32 (2002) 53–76. doi:10.1146/annurev.matsci.32.090601.152855 .[4] L. Barrales Mora, G. Gottstein, L. Shvindlerman, Three-dimensional graingrowth: Analytical approaches and computer simulations, Acta Materialia56 (20) (2008) 5915–5926. doi:10.1016/j.actamat.2008.08.006 .URL https://linkinghub.elsevier.com/retrieve/pii/S1359645408005661 [5] L. Rauch, L. Madej, P. Spytkowski, R. Golab, Development of the cellularautomata framework dedicated for metallic materials microstructure evo-lution models, Archives of Civil and Mechanical Engineering 15 (1) (2015)48–61. doi:10.1016/j.acme.2014.06.006 .URL http://dx.doi.org/10.1016/j.acme.2014.06.006 [6] L. Madej, M. Sitko, A. Legwand, K. Perzynski, K. Michalik, Developmentand evaluation of data transfer protocols in the fully coupled random cellu-lar automata finite element model of dynamic recrystallization, Journal ofComputational Science 26 (2018) 66–77. doi:10.1016/j.jocs.2018.03.007 .URL https://doi.org/10.1016/j.jocs.2018.03.007 [7] I. Steinbach, F. Pezzolla, B. Nestler, M. Seeßelberg, R. Prieler, G. J.Schmitz, J. L. Rezende, A phase field concept for multiphase systems,Physica D: Nonlinear Phenomena 94 (3) (1996) 135–147. doi:10.1016/0167-2789(95)00298-7 .[8] N. Moelans, B. Blanpain, P. Wollants, Quantitative analysis of grainboundary properties in a generalized phase field model for grain growthin anisotropic systems, Physical Review B - Condensed Matter and Mate-rials Physics 78 (2). doi:10.1103/PhysRevB.78.024113 .[9] C. E. Krill, L. Q. Chen, Computer simulation of 3-D grain growth usinga phase-field model, Acta Materialia 50 (12) (2002) 3057–3073. doi:10.1016/s1359-6454(02)00084-8 . 4210] H. K. Kim, S. G. Kim, W. Dong, I. Steinbach, B. J. Lee, Phase-field mod-eling for 3D grain growth based on a grain boundary energy database,Modelling and Simulation in Materials Science and Engineering 22 (3). doi:10.1088/0965-0393/22/3/034004 .[11] K. Kawasaki, T. Nagai, K. Nakashima, Vertex models for two-dimensional grain growth, Philosophical Magazine B 60 (3) (1989)399–421. doi:10.1080/13642818908205916 .URL [12] D. Weygand, Y. Br´echet, J. L´epinoux, A vertex dynamics simulation ofgrain growth in two dimensions, Philosophical Magazine B 78 (4) (1998)329–352. doi:10.1080/13642819808206731 .URL [13] J. L´epinoux, D. Weygand, M. Verdier, Modeling grain growth and relatedphenomena with vertex dynamics, Comptes Rendus Physique 11 (3-4)(2010) 265–273. doi:10.1016/j.crhy.2010.07.015 .URL https://linkinghub.elsevier.com/retrieve/pii/S1631070510000800 [14] L. A. Barrales Mora, 2D vertex modeling for the simulation of grain growthand related phenomena, Mathematics and Computers in Simulation 80 (7)(2010) 1411–1427. doi:10.1016/j.matcom.2009.08.005 .URL http://dx.doi.org/10.1016/j.matcom.2009.08.005 [15] Y. Mellbin, H. Hallberg, M. Ristinmaa, A combined crystal plasticityand graph-based vertex model of dynamic recrystallization at large defor-mations, Modelling and Simulation in Materials Science and Engineering23 (4). doi:10.1088/0965-0393/23/4/045011 .[16] B. Merriman, J. K. Bence, S. J. Osher, Motion of Multiple Junctions: ALevel Set Approach, Journal of Computational Physics 112 (2) (1994) 334–363. doi:10.1006/jcph.1994.1105 .[17] M. Bernacki, Y. Chastel, T. Coupez, R. E. Log´e, Level set frameworkfor the numerical modelling of primary recrystallization in polycrystallinematerials, Scripta Materialia 58 (12) (2008) 1129–1132. doi:10.1016/j.scriptamat.2008.02.016 .[18] A. L. Cruz-Fabiano, R. Log´e, M. Bernacki, Assessment of simplified 2Dgrain growth models from numerical experiments based on a level setframework, Computational Materials Science 92 (2014) 305–312. doi:10.1016/j.commatsci.2014.05.060 .[19] L. Maire, B. Scholtes, C. Moussa, D. Pino Mu˜noz, N. Bozzolo, M. Bernacki,Improvement of 3-D mean field models for pure grain growth based on fullfield simulations, Journal of Materials Science 51 (24) (2016) 10970–10981.4320] J. Humphreys, G. S. Rohrer, A. Rollett, Recrystallization and RelatedAnnealing Phenomena (Third Edition), Elsevier, 2017. doi:10.1016/B978-0-08-098235-9.00004-5 .[21] M. Bernacki, N. Bozzolo, P. De Micheli, B. Flipon, J. Fausty, L. Maire,S. Florez, Numerical Modeling of Recrystallization in a Level Set FiniteElement Framework for Application to Industrial Processes, in: Recrystal-lization: Types, Techniques and Applications, Nova, 2019.[22] S. Florez, K. Alvarado, D. P. Mu˜noz, M. Bernacki, A novel highly effi-cient Lagrangian model for massively multidomain simulation applied tomicrostructural evolutions, Computer Methods in Applied Mechanics andEngineering 367 (2020) 113107. doi:10.1016/j.cma.2020.113107 .URL https://doi.org/10.1016/j.cma.2020.113107 [23] S. Florez, J. Fausty, K. Alvarado, B. Murgas, M. Bernacki, A novel highlyefficient Lagrangian model for massively multidomain simulations: parallelcontext, arXiv preprint arXiv:2009.04424.[24] M. Shakoor, P.-O. Bouchard, M. Bernacki, An adaptive level-set methodwith enhanced volume conservation for simulations in multiphase domains,International Journal for Numerical Methods in Engineering 109 (4) (2017)555–576. doi:10.1002/nme.5297 .[25] S. Florez, M. Shakoor, T. Toulorge, M. Bernacki, A new finite elementstrategy to simulate microstructural evolutions, Computational MaterialsScience 172 (2020) 109335. doi:10.1016/J.COMMATSCI.2019.109335 .[26] B. Scholtes, M. Shakoor, A. Settefrati, P.-O. Bouchard, N. Bozzolo,M. Bernacki, New finite element developments for the full field modelingof microstructural evolutions using the level-set method, ComputationalMaterials Science 109 (2015) 388–398. doi:10.1016/j.commatsci.2015.07.042 .[27] B. Scholtes, R. Boulais-Sinou, A. Settefrati, D. Pino Mu˜noz, I. Poitrault,A. Montouchet, N. Bozzolo, M. Bernacki, 3D level set modeling of staticrecrystallization considering stored energy fields, Computational MaterialsScience 122 (2016) 57–71. doi:10.1016/j.commatsci.2016.04.045 .[28] G. Comp`ere, J. F. Remacle, E. Marchandise, Transient mesh adaptivitywith large rigid-body displacements, Proceedings of the 17th InternationalMeshing Roundtable, IMR 2008 (iMMC) (2008) 213–230. doi:10.1007/978-3-540-87921-3-13 .[29] G. Comp`ere, J.-F. Remacle, J. Jansson, J. Hoffman, A mesh adaptationframework for dealing with large deforming meshes, International Journalfor Numerical Methods in Engineering 82 (7) (2010) 843–867. arXiv:1010.1724 , doi:10.1002/nme.2788 . 4430] D. W. Walker, Standards for message-Passing in a Distributed Memory En-vironment, Tech. rep., Center for Research on Parallel Computing (CRPC),Oak Ridge National Lab., TN (United States) (1992).[31] G. Karypis, V. Kumar, A fast and high quality multilevel scheme for par-titioning irregular graphs, SIAM Journal of Scientific Computing 20 (1)(1998) 359–392. doi:10.1137/S1064827595287997 .[32] H. Hallberg, Influence of anisotropic grain boundary properties on the evo-lution of grain boundary character distribution during grain growth - A 2Dlevel set study, Modelling and Simulation in Materials Science and Engi-neering 22 (8). doi:10.1088/0965-0393/22/8/085005 .[33] J. Furstoss, M. Bernacki, C. Ganino, C. Petit, D. Pino-Mu˜noz, 2D and3D simulation of grain growth in olivine aggregates using a full field modelbased on the level set method, Physics of the Earth and Planetary Interiors283 (2018) 98–109.[34] M. Bernacki, H. Resk, T. Coupez, R. E. Log´e, Finite element model of pri-mary recrystallization in polycrystalline aggregates using a level set frame-work, Modelling and Simulation in Materials Science and Engineering 17 (6)(2009) 64006. doi:10.1088/0965-0393/17/6/064006 .[35] M. Bernacki, R. E. Log´e, T. Coupez, Level set framework for the finite-element modelling of recrystallization and grain growth in polycrystallinematerials, Scripta Materialia 64 (6) (2011) 525–528. doi:10.1016/j.scriptamat.2010.11.032 .[36] L. Maire, B. Scholtes, C. Moussa, N. Bozzolo, D. P. Mu˜noz, A. Settefrati,M. Bernacki, Modeling of dynamic and post-dynamic recrystallization bycoupling a full field approach to phenomenological laws, Materials and De-sign 133 (2017) 498–519. doi:10.1016/j.matdes.2017.08.015 .URL http://dx.doi.org/10.1016/j.matdes.2017.08.015 [37] H. Hallberg, A modified level set approach to 2D modeling of dynamicrecrystallization, Modelling and Simulation in Materials Science and Engi-neering 21 (8) (2013) 85012. doi:10.1088/0965-0393/21/8/085012 .[38] D. Weygand, Y. Br´echet, J. L´epinoux, Zener pinning and grain growth:a two-dimensional vertex computer simulation, Acta Materialia 47 (3)(1999) 961–970. doi:10.1016/S1359-6454(98)00383-8 .URL https://linkinghub.elsevier.com/retrieve/pii/S1359645498003838 [39] G. Couturier, C. Maurice, R. Fortunier, Three-dimensional finite-elementsimulation of Zener pinning dynamics, Philosophical Magazine 83 (30)(2003) 3387–3405. doi:10.1080/1478643031000152771 .4540] G. Couturier, C. Maurice, R. Fortunier, R. Doherty, J. H. Driver, Finiteelement simulations of 3D Zener pinning, in: Materials Science Forum,Vol. 467-470, 2004, pp. 1009–1018. .[41] G. Couturier, R. Doherty, C. Maurice, R. Fortunier, 3D finite element simu-lation of the inhibition of normal grain growth by particles, Acta Materialia53 (4) (2005) 977–989. doi:10.1016/j.actamat.2004.10.044 .[42] A. Agnoli, M. Bernacki, R. Log´e, J. M. Franchet, J. Laigo, N. Bozzolo, Se-lective Growth of Low Stored Energy Grains During δ Sub-solvus Anneal-ing in the Inconel 718 Nickel-Based Superalloy, Metallurgical and MaterialsTransactions A: Physical Metallurgy and Materials Science 46 (9) (2015)4405–4421. doi:10.1007/s11661-015-3035-9 .[43] D. N. Ilin, N. Bozzolo, T. Toulorge, M. Bernacki, Full field modeling ofrecrystallization: Effect of intragranular strain gradients on grain boundaryshape and kinetics, Computational Materials Science 150 (March) (2018)149–161. doi:10.1016/j.commatsci.2018.03.063 .[44] A. Laasraoui, J. J. Jonas, Prediction of steel flow stresses at high temper-atures and strain rates, Metallurgical Transactions A 22 (7) (1991) 1545–1558. doi:10.1007/BF02667368 .[45] A. L. Cruz-Fabiano, Modelling of crystal plasticity and grain boundary mi-gration of 304L steel at the mesoscopic scale, Ph.D. thesis, MINES Paris-Tech (2014).[46] P. Peczak, M. J. Luton, The effect of nucleation models on dynamic re-crystallization i. homogeneous stored energy distribution, PhilosophicalMagazine B 68 (1) (1993) 115–144. arXiv:https://doi.org/10.1080/13642819308215285 , doi:10.1080/13642819308215285 .URL https://doi.org/10.1080/13642819308215285 [47] J. Bailey, P. B. Hirsch, The recrystallization process in some polycrystallinemetals, Proceedings of the Royal Society of London. Series A. Mathematicaland Physical Sciences 267 (1328) (1962) 11–30.[48] B. Scholtes, Development of an efficient level set framework for the fullfield modeling of recrystallization in 3D., Ph.D. thesis, MINES ParisTech(2016). arXiv:1410.7235 , doi:10.1017/jfm.2015.566 .[49] H. W. Hesselbarth, I. R. G¨obel, Simulation of recrystallization by cellularautomata, Acta Metallurgica Et Materialia 39 (9) (1991) 2135–2143. doi:10.1016/0956-7151(91)90183-2 .[50] C. H. Davies, Growth of nuclei in a cellular automaton simulation ofrecrystallisation, Scripta Materialia 36 (1) (1997) 35–40. doi:10.1016/S1359-6462(96)00331-4 . 4651] L. Sieradzki, L. Madej, A perceptive comparison of the cellular automataand Monte Carlo techniques in application to static recrystallization model-ing in polycrystalline materials, Computational Materials Science 67 (2013)156–173. doi:10.1016/j.commatsci.2012.08.047 .URL http://dx.doi.org/10.1016/j.commatsci.2012.08.047 [52] F. Villaret, B. Hary, Y. de Carlan, T. Baudin, R. Log´e, L. Maire,M. Bernacki, Probabilistic and deterministic full field approaches to sim-ulate recrystallization in ODS steels, Computational Materials Science179 (March) (2020) 109646. doi:10.1016/j.commatsci.2020.109646 .URL https://doi.org/10.1016/j.commatsci.2020.109646 [53] F. Reitich, H. M. Soner, Three-phase boundary motions under constant ve-locities. I: The vanishing surface tension limit, Royal Society of Edinburgh -Proceedings A 126 (4) (1996) 837–865. doi:10.1017/S0308210500023106 .[54] J. E. Taylor, The motion of multiple-phase junctions under prescribedphase-boundary velocities, Journal of Differential Equations 119 (1) (1995)109–136. doi:10.1006/jdeq.1995.1085 .[55] H. Imai, M. Iri, K. Murota, Voronoi diagram in the Laguerre geometry andits applications, SIAM Journal on Computing 14 (1) (1985) 93–105.[56] K. Hitti, P. Laure, T. Coupez, L. Silva, M. Bernacki, Precise generationof complex statistical Representative Volume Elements (RVEs) in a finiteelement context, Computational Materials Science 61 (2012) 224–238. doi:10.1016/j.commatsci.2012.04.011 .[57] D. N. Ilin, M. Bernacki, Advancing layer algorithm of dense ellipse pack-ing for generating statistically equivalent polygonal structures, GranularMatter 18 (3) (2016) 43. doi:10.1007/s10035-016-0646-9 .[58] L. Maire, Full field and mean field modeling of dynamic and post-dynamicrecrystallization in 3D – Application to 304L steel, Ph.D. thesis, PSL,Mines-ParisTech (2019).[59] R. Log´e, M. Bernacki, H. Resk, L. Delannay, H. Digonnet, Y. Chastel,T. Coupez, Linking plastic deformation to recrystallization in metals usingdigital microstructures, Philosophical Magazine 88 (30-32) (2008) 3691–3712. doi:10.1080/14786430802502575doi:10.1080/14786430802502575