- Email: [email protected]
Level set framework for the numerical modelling of primary recrystallization in polycrystalline materials
Available online at www.sciencedirect.com
Scripta Materialia 58 (2008) 1129–1132 www.elsevier.com/locate/scriptamat
Level set framework for the numerical modelling of primary recrystallization in polycrystalline materials M. Bernacki,* Y. Chastel, T. Coupez and R.E. Loge´ Centre de mise en forme des mate´riaux (CEMEF), Ecole Nationale Supe´rieure des Mines de Paris, UMR CNRS 7635, BP 207, 06904 Sophia Antipolis Cedex, France Received 7 December 2007; revised 9 February 2008; accepted 11 February 2008 Available online 26 February 2008
This letter describes a level set framework for the numerical modelling of primary static recrystallization in a polycrystalline material. The topological evolution of the grain structure is simulated in two and three dimensions, based on a kinetic law relating the velocity of the boundary to the thermodynamic driving force. The adopted finite element approach is described, discussed and tested from simple to more complex configurations. The possibility to accurately describe nucleation and growth is illustrated in three dimensions. Ó 2008 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Recrystallization; Finite element analysis; Level set framework; Primary recrystallization; Nucleation
Recrystallization phenomena inevitably occur during thermal and mechanical processes and have a major impact on the final in-use properties of the material. Their accurate control is therefore of prime importance. Theories for recrystallization that provide quantitatively correct predictions of crystallographic orientation and grain size distributions have long been sought to fill a critical missing link to model material processing from start to finish. To date, no such theory exists. The phenomena of grain growth or recrystallization seem simple, but their mechanisms are not very well understood. Multiscale models are in principle needed to fully describe recrystallization phenomena in a generic way [1–4]. Over the last decade, considerable progress has been made in the numerical simulation of primary recrystallization [5]. Common approaches include the Monte Carlo (MC) method [6], the cellular automaton (CA) methods [6], the phase field method [7] and the level set method [8]. The first two are probabilistic techniques which deliver grain structures with kinetics. They are associated with a two-dimensional (2-D) or 3-D geometric representation of the microstructure, discretized on a regular grid made of ‘‘cells” which are allocated to the grains. The standard MC method as derived from the Potts model (multistate Ising model) applies probabilistic rules to each cell at each time step of the simulation. * Corresponding author. E-mail: [email protected]
The use of this model in three dimensions is relatively easy and efficient [9]. However, the comparison between MC results and experiments is not straightforward [1]. Furthermore, the standard form of the model does not result in a linear relationship between stored energy and migration rate. The CA method uses physically based rules to determine the propagation rate of a transformation from one cell to a neighbouring cell [10], and can therefore be related to the microstructure and kinetics of a real system. In the case of primary recrystallization the switch rule is simple: an unrecrystallized cell will switch to being recrystallized if one of its neighbours is recrystallized. A major problem of the CA method to model recrystallization is the absence of effective methods for the treatment of nucleation phenomena [1]. The two others methods, i.e. the phase field and level set methods, have many common points. They both have the advantage of avoiding the difficult problem of tracking interfaces. The principle of the phase-field model consists in describing the location of phases by introducing an order parameter (the phase field) which varies smoothly from one to zero (or from minus one to one) through a diffuse interface [11]. The concept has been extended to deal with more complex problems involving more than two phases and for modelling microstructure evolution [12,13]. Each grain orientation is then described by a non-conserved order parameter field and the free energy density of a grain is formulated as a Landau expansion in terms of the structural order parameters. As for the
1359-6462/$ - see front matter Ó 2008 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.scriptamat.2008.02.016
1130
M. Bernacki et al. / Scripta Materialia 58 (2008) 1129–1132
MC or CA methods, the topological events are treated in a natural way as a result of energy minimization. In the case of 2-D ideal normal grain growth, published results illustrate the potential of the approach [12]. However, the main difficulty of this method remains the construction of the free energy density function. Furthermore, the energy minimization of each order parameter can involve very expensive and intensive calculations, particularly for 3-D systems. In contrast, the level set method [14,15] is now commonly used to follow propagating fronts in various numerical models [16]. In Ref. [17], the authors have extended the level set method to model the motion of multiple junctions when more than two regions or grains intersect. They developed an appropriate method for the grain growth problem, i.e. when the driving force for microstructure evolution is grain boundary energy (zero bulk stored energy). Each region has its own private level set function, and this function moves each level set with a normal velocity defined by the conditions at the nearest interface, and a reassignment step is used to avoid kinematic incompatibilities, i.e. the development of vacuum and overlapping regions. In this work, a new level set framework is proposed, in order to take into account bulk stored energies and to avoid the development of vacuum and overlapping regions in two or three dimensions. Nucleation and growth phenomena can also be described accurately. It is generally assumed, for relatively pure metals [2,18], that the kinetic law for grain boundary motion is described well by: ~ v ¼ MDf~ n;
ð1Þ
where M corresponds to the grain boundary mobility, Df is the driving force per unit area and ~ n is the outside unit normal to the grain boundary. The primary recrystallization driving force Df is approximated by: Df ¼ sDq;
ð2Þ
where s corresponds to the dislocation line energy and D q is the total dislocation density difference across the interface. If each interface Cij separates grains Gi and Gj, the oriented normal velocity from Gi to Gj can be defined as: ~ nij ¼ M ij ðej ei Þ; vij ~
ð3Þ
with ei = sqi the average stored energy of grain Gi. These equations are simple, but the associated prediction of topological evolution is complex, in particular for the treatment of multiple junctions. It has been shown that if no further conditions on the motion are imposed the solution of the problem is not unique [17,19]. A function /, defined over a domain X, is called level set of an interface C if, at any points of X, it corresponds to the distance from C. In turn, the interface C is then given by the level 0 of the function /: /ðxÞ ¼ dðx; CÞ; x 2 X; ð4Þ C ¼ fx 2 X; /ðxÞ ¼ 0g: Assuming that the domain contains NG grains, one level set function per grain {/i,1 6 i 6 NG} is considered. A sign convention is used: /i P 0 inside grain Gi and v is /i 6 0 outside grain Gi. If a common velocity field ~
defined over X, the motion of the interfaces (or any isovalue of the level set functions) is given by a convection equation: ( o/ i þ~ v r/i ¼ 0; ot ð5Þ 8i 2 f1; . . . ; N G g: /i ðt ¼ 0; xÞ ¼ /0i ðxÞ; At any time t the interface of grain Gi is given implicitly v in X will by the equation /i(t,x) = 0. The expression of ~ be detailed further. If the initial value of /i in Eq. (5) is a level set function, it is a priori not the case of the solution, i.e. k$/i(t,x)k 6¼ 1. Even if the interface remains the zero level of the solution, the gradient of /i can take any value and can generate instabilities for the finite element formulations adapted to the problem of convection. The problem is treated with a re-initialization technique [14]. The usual approach consists in re-initializing the level set functions by solving the following Hamilton–Jacobi equation: ( obi þ si ðkrbi k 1Þ ¼ 0; os with si ¼ signðbi Þ: ð6Þ bi ðs ¼ 0; xÞ ¼ /i ðt; xÞ; In practice, Eq. (6) is solved periodically, typically every few time steps, when the function /i becomes too irregular. Eq. (6) can be seen as a pure convection equation. ~ i ¼ si rbi , Eq. (6) becomes: Indeed, if we write U krbi k ( obi ~ i :rbi ¼ si ; þU os ð7Þ bi ðs ¼ 0; xÞ ¼ /i ðt; xÞ: In practice, to define the fictitious time step s, a stability ~ i kDs ¼ Ds ¼ h, where h correcondition is used: kU sponds to the mesh size. So, in fact, the fictitious time step is chosen as h. Instead of using this classical re-initialization method, the formulation developed in Ref. [20] is chosen. The idea proposed by this author is to modify Eq. (5) in order to keep the property of a level set function. The parameter k = ds/dt is introduced (in practice k h/Dt) to write the following equality o ot o ot o ¼ os þ ox r ¼ os ðot þ~ v rÞ and to transform os ot os Eq. (6): obi þ~ v:rbi þ ksi ðkrbi k 1Þ ¼ 0: ð8Þ ot The new formulation incorporating both convection and re-initialization becomes: ( o/ i þ~ v r/i þ ksi ðkr/i k 1Þ ¼ 0; ot 8i 2 f1; . . . ; N G g; /i ðt ¼ 0; xÞ ¼ /0i ðxÞ; ð9Þ which can be written as a pure convection problem: ( o/ i ~ i Þ r/i ¼ ksi ; þ ð~ v þ kU ot 8i 2 f1; . . . ; N G g: /i ðt ¼ 0; xÞ ¼ /0i ðxÞ; ð10Þ The problem is now defined by the system (10), but the common velocity field in X must be defined. In order to avoid kinematic incompatibilities, it is necessary to work with one same velocity field for all the level set functions, whereas the kinetic law is built according to
M. Bernacki et al. / Scripta Materialia 58 (2008) 1129–1132
parameters which are specific to each grain (bulk stored energy and unit normals). Geometry of the interfaces is known to be well described by the level set method [15], so it is well adapted to the calculation of the unit normal, or the curvature if it were needed. Considering grain Gi, the outward unitary normal of any isovalue of /i is defined by: ~ ni ¼
r/i kr/i k
¼
kr/i k¼1
r/i :
ð11Þ
At any time t, the level set approach allows, very easily, to describe the position of grains. Indeed, since the set {Gi nCi,1 6 i 6 NG} corresponds to a partition of the domain X, and, keeping in mind the sign convention described previously, the following results are verified for a point M with coordinates x: M 2 Ci () /i ðt; xÞ ¼ 0;
ð12Þ
M 2 Gi Ci () /i ðt; xÞ > 0 () 8j 6¼ i/j ðt; xÞ 6 0:
ð13Þ
The description of an interface by a level set function remains a ‘‘fuzzy” description of the interface, which means that the level 0 of the level set function (the interface) does not necessarily correspond to nodes of the mesh. The point M is located using the following relationship: M 2 Gi () /i ðt; xÞ ¼ max ð/k ðt; xÞÞ: 16k6N G
ð14Þ
A classical way to build a velocity field from the stored energy field for any node of the mesh is as follows: – find i with /i ðt; xÞ ¼ max16k6N G ð/k ðt; xÞÞ – compare /i(t,x) to a positive fixed parameter ‘ which defines a length scale related to the proximity of the node to the interface:
1131
where a is a positive fixed parameter calibrated to obtain a negligible exponential term outside the anisotropic part of the mesh (see Fig. 2). Expression (16) leads to a smoother velocity field and avoids the need to identify the neighbouring grains of each grain. Figure 1 illustrates the difference between Eqs. (15) and (16) for a 2-D triple junction using a common unstructured mesh. The chosen stored energy field is described in Figure 1. ~ v is shown at the boundaries; the smoothing effect of Eq. (16) appears clearly. To solve the problem defined by Eqs. (10) and (16), a finite element method based on the C++ library ‘CimLib’ [20] is used, with unstructured meshes and a stabilized P1 solver (SUPG or RFB method). The quality of the described approach is strongly related to the accuracy of level set function calculations around the interfaces. Indeed, a little disturbance of the levels around an interface leads to an error in the velocity estimation and consequently in their own evolution. For a given computational cost, optimal accuracy is obtained by using anisotropic meshes, with refinement close to the boundaries. The technique used to generate anisotropic meshes adapted to the grain structures is presented in Ref. [3]. Figure 2 illustrates the evolution of the triple junction presented in Figure 1 after 0, 90, 145 and 240 time steps, with automatic adaptive anisotropic remeshing. Results are very accurate and correspond to one of the exact solutions of this grain configuration (see Ref. [19]). More precisely, the L2 error between the exact and approached function defined by max06t6T end ðmax16k6N G ð/k ðt; xÞÞÞ is less than 3%. A second test case corresponds to a five grain microstructure, as shown in Figure 3. The energy of the circular grain is zero, whereas the energy of the four surrounding grains is unity. As for the previous test case, the interest of this configuration is that it is possible to
*
If/i ðt; xÞ > ‘ ) v ðt; xÞ ¼~ 0 else find j with /j ðt; xÞ ¼
*
ð/k Þ ) v ðt; xÞ max 1 6 k 6 NG k 6¼ i
¼ M ij f ð/i ðt; xÞ; ‘Þðej ei Þ~ ni ðt; xÞ
ð15Þ
with f a decreasing continuous function varying from 1 to 0 when /i(t,x) varies from 0 to ‘, and Mij the mobility at the interface between grain i and grain j. The main disadvantage of this approach lies in the management of multiple junctions: discontinuous velocity fields are generated, which then lead to convergence problems in the resolution of the convection equation. We propose another algorithm which considers, at any point x, all NG level set functions: Find i with
Figure 1. Representation of ~ v (left) for the formulation defined by Eq. (15) and (right) for the formulation defined by Eq. (16).
vðt; xÞ /i ðt; xÞ ¼ max ð/k ðt; xÞÞ ) ~ 16k6N G
¼
NG X
nj ðt; xÞ; M ij expðaj/j ðt; xÞjÞðei ej Þ~
j ¼ 1j 6¼ i ð16Þ
Figure 2. Evolution of the ±5 103 isovalues of the three distance functions with adapted anisotropic mesh.
1132
M. Bernacki et al. / Scripta Materialia 58 (2008) 1129–1132
Figure 5. Evolution of the ±5 103 isovalues of 10 distance functions. Figure 3. Evolution of the ±5 103 isovalues of the five distance functions with adapted anisotropic mesh.
compare the numerical results with the exact solution (defined by the growth of the circular grain with a unit normal velocity). Figure 3 shows the topology after 0, 90, 145 and 240 time steps using adaptive anisotropic remeshing. A very good agreement is obtained between the ±5 103 isovalues of the five distance functions and the exact solution (in black). More precisely, the L2 error defined as for the previous test case is less than 2%. One of the prominent advantages in using front capturing methods for describing interface motions is that there is no need for additional features when some regions (grains) disappear. Complex topological evolutions are handled automatically. In a similar way, it is possible to introduce new regions (grains), based on given criteria. A very simple method to create a nucleation site is to build a new level set function at the desired spatial location x and at a specific time increment. For example, a new distance function can be built such that the boundary of the nucleus is circular (in two dimensions) or spherical (in three dimensions). The circle or sphere can be centered around one node of the mesh with zero stored energy inside such that spontaneous growth occurs. A new level set function can be associated to a new nucleus, which evolves subsequently according to the principles described by Eq. (16). Figure 4 describes the creation and the evolution in three dimensions of one recrystallization nucleus, using the same triple junction configuration as in Figure 2. The chosen location is the node which is the closest to the triple junction. The criterion for nucleation is geometric here, but could be based on the distribution of any mechanical/physical quantity that has been stored in the finite element mesh, e.g. from previous modelling of the polycrystal deformation (crystal plasticity, etc.).
Figure 4. Example of nucleation in a 3-D three-grain case.
Another initial configuration of ten grains with an adapted anisotropic mesh is presented in Figure 5. A different stored energy is assumed in each grain. Three triple junctions are considered as possible nucleation sites. The results illustrate the capability of the level set method to model primary recrystallization in polycrystals. It has been shown that a level set framework, associated with adaptive anisotropic automatic remeshing, is a promising tool to describe primary recrystallization in a polycrystalline material. Full kinematic compatibilities can be enforced. Nucleation and growth phenomena can be implemented in this framework. Various studies are in process, including simulations with a statistical number of grains, the prediction of recrystallization kinetics and the comparison with other models. Future work will account for anisotropic mobility of grain boundaries (as related to the crystallography of the interface). Stored energy due to prior deformation of the polycrystalline aggregate will be used from crystal plasticity finite element simulations, and physical nucleation criteria will be studied. [1] A.D. Rollett, Prog. Mater. Sci. 42 (1997) 77. [2] G. Kugler, R. Turk, Comput. Mater. Sci. 37 (3) (2006) 284. [3] M. Bernacki, Y. Chastel, H. Digonnet, H. Resk, T. Coupez, R.E. Loge´, Comput. Methods Mater. Sci. 7 (1) (2007) 141. [4] R.E. Loge´, M. Bernacki, H. Resk, H. Digonnet, T. Coupez, in: Recrystallization and Grain Growth III Proceedings, Korea, 2007, pp. 1133–1138. [5] M.A. Miodownik, J. Light Met. 2 (2002) 125. [6] A.D. Rollett, D. Raabe, Comput. Mater. Sci. 21 (2001) 69. [7] L.Q. Chen, Scripta Metall. Mater. 32 (1995) 115. [8] H.-K. Zhao et al., J. Comput. Phys. 127 (1996) 179. [9] G.N. Hassold, E.A. Holm, Comput. Phys. 7 (1993) 97. [10] D. Raabe, Philos. Mag. A 79 (1999) 2339. [11] J.B. Collins, H. Levine, Phys. Rev. B31 (1985) 6119. [12] L.Q. Chen, Annu. Rev. Mater. Res. 32 (2002) 113. [13] A. Karma, Phys. Rev. Lett. 8711 (11) (2001) 115701. [14] J.A. Sethian, Level Set Methods, Cambridge University Press, 1996. [15] S. Osher, J.A. Sethian, J. Comput. Phys. 79 (1988) 12. [16] M. Sussman et al., J. Comput. Phys. 114 (1994) 146. [17] B. Merriman, J.K. Bence, S.J. Osher, J. Comput. Phys. 112 (2) (1994) 334. [18] G. Kugler, R. Turk, Acta Mater. 52 (2004) 4659. [19] F. Reitich, H.M. Soner, Proc. Roy. Soc. Ed. 126A (1996) 837–865. [20] T. Coupez, in: Numiform’07 Proceedings, Plenary Lecture, Porto, 2007, pp. 61–66.
Scripta Materialia 58 (2008) 1129–1132 www.elsevier.com/locate/scriptamat
Level set framework for the numerical modelling of primary recrystallization in polycrystalline materials M. Bernacki,* Y. Chastel, T. Coupez and R.E. Loge´ Centre de mise en forme des mate´riaux (CEMEF), Ecole Nationale Supe´rieure des Mines de Paris, UMR CNRS 7635, BP 207, 06904 Sophia Antipolis Cedex, France Received 7 December 2007; revised 9 February 2008; accepted 11 February 2008 Available online 26 February 2008
This letter describes a level set framework for the numerical modelling of primary static recrystallization in a polycrystalline material. The topological evolution of the grain structure is simulated in two and three dimensions, based on a kinetic law relating the velocity of the boundary to the thermodynamic driving force. The adopted finite element approach is described, discussed and tested from simple to more complex configurations. The possibility to accurately describe nucleation and growth is illustrated in three dimensions. Ó 2008 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Recrystallization; Finite element analysis; Level set framework; Primary recrystallization; Nucleation
Recrystallization phenomena inevitably occur during thermal and mechanical processes and have a major impact on the final in-use properties of the material. Their accurate control is therefore of prime importance. Theories for recrystallization that provide quantitatively correct predictions of crystallographic orientation and grain size distributions have long been sought to fill a critical missing link to model material processing from start to finish. To date, no such theory exists. The phenomena of grain growth or recrystallization seem simple, but their mechanisms are not very well understood. Multiscale models are in principle needed to fully describe recrystallization phenomena in a generic way [1–4]. Over the last decade, considerable progress has been made in the numerical simulation of primary recrystallization [5]. Common approaches include the Monte Carlo (MC) method [6], the cellular automaton (CA) methods [6], the phase field method [7] and the level set method [8]. The first two are probabilistic techniques which deliver grain structures with kinetics. They are associated with a two-dimensional (2-D) or 3-D geometric representation of the microstructure, discretized on a regular grid made of ‘‘cells” which are allocated to the grains. The standard MC method as derived from the Potts model (multistate Ising model) applies probabilistic rules to each cell at each time step of the simulation. * Corresponding author. E-mail: [email protected]
The use of this model in three dimensions is relatively easy and efficient [9]. However, the comparison between MC results and experiments is not straightforward [1]. Furthermore, the standard form of the model does not result in a linear relationship between stored energy and migration rate. The CA method uses physically based rules to determine the propagation rate of a transformation from one cell to a neighbouring cell [10], and can therefore be related to the microstructure and kinetics of a real system. In the case of primary recrystallization the switch rule is simple: an unrecrystallized cell will switch to being recrystallized if one of its neighbours is recrystallized. A major problem of the CA method to model recrystallization is the absence of effective methods for the treatment of nucleation phenomena [1]. The two others methods, i.e. the phase field and level set methods, have many common points. They both have the advantage of avoiding the difficult problem of tracking interfaces. The principle of the phase-field model consists in describing the location of phases by introducing an order parameter (the phase field) which varies smoothly from one to zero (or from minus one to one) through a diffuse interface [11]. The concept has been extended to deal with more complex problems involving more than two phases and for modelling microstructure evolution [12,13]. Each grain orientation is then described by a non-conserved order parameter field and the free energy density of a grain is formulated as a Landau expansion in terms of the structural order parameters. As for the
1359-6462/$ - see front matter Ó 2008 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.scriptamat.2008.02.016
1130
M. Bernacki et al. / Scripta Materialia 58 (2008) 1129–1132
MC or CA methods, the topological events are treated in a natural way as a result of energy minimization. In the case of 2-D ideal normal grain growth, published results illustrate the potential of the approach [12]. However, the main difficulty of this method remains the construction of the free energy density function. Furthermore, the energy minimization of each order parameter can involve very expensive and intensive calculations, particularly for 3-D systems. In contrast, the level set method [14,15] is now commonly used to follow propagating fronts in various numerical models [16]. In Ref. [17], the authors have extended the level set method to model the motion of multiple junctions when more than two regions or grains intersect. They developed an appropriate method for the grain growth problem, i.e. when the driving force for microstructure evolution is grain boundary energy (zero bulk stored energy). Each region has its own private level set function, and this function moves each level set with a normal velocity defined by the conditions at the nearest interface, and a reassignment step is used to avoid kinematic incompatibilities, i.e. the development of vacuum and overlapping regions. In this work, a new level set framework is proposed, in order to take into account bulk stored energies and to avoid the development of vacuum and overlapping regions in two or three dimensions. Nucleation and growth phenomena can also be described accurately. It is generally assumed, for relatively pure metals [2,18], that the kinetic law for grain boundary motion is described well by: ~ v ¼ MDf~ n;
ð1Þ
where M corresponds to the grain boundary mobility, Df is the driving force per unit area and ~ n is the outside unit normal to the grain boundary. The primary recrystallization driving force Df is approximated by: Df ¼ sDq;
ð2Þ
where s corresponds to the dislocation line energy and D q is the total dislocation density difference across the interface. If each interface Cij separates grains Gi and Gj, the oriented normal velocity from Gi to Gj can be defined as: ~ nij ¼ M ij ðej ei Þ; vij ~
ð3Þ
with ei = sqi the average stored energy of grain Gi. These equations are simple, but the associated prediction of topological evolution is complex, in particular for the treatment of multiple junctions. It has been shown that if no further conditions on the motion are imposed the solution of the problem is not unique [17,19]. A function /, defined over a domain X, is called level set of an interface C if, at any points of X, it corresponds to the distance from C. In turn, the interface C is then given by the level 0 of the function /: /ðxÞ ¼ dðx; CÞ; x 2 X; ð4Þ C ¼ fx 2 X; /ðxÞ ¼ 0g: Assuming that the domain contains NG grains, one level set function per grain {/i,1 6 i 6 NG} is considered. A sign convention is used: /i P 0 inside grain Gi and v is /i 6 0 outside grain Gi. If a common velocity field ~
defined over X, the motion of the interfaces (or any isovalue of the level set functions) is given by a convection equation: ( o/ i þ~ v r/i ¼ 0; ot ð5Þ 8i 2 f1; . . . ; N G g: /i ðt ¼ 0; xÞ ¼ /0i ðxÞ; At any time t the interface of grain Gi is given implicitly v in X will by the equation /i(t,x) = 0. The expression of ~ be detailed further. If the initial value of /i in Eq. (5) is a level set function, it is a priori not the case of the solution, i.e. k$/i(t,x)k 6¼ 1. Even if the interface remains the zero level of the solution, the gradient of /i can take any value and can generate instabilities for the finite element formulations adapted to the problem of convection. The problem is treated with a re-initialization technique [14]. The usual approach consists in re-initializing the level set functions by solving the following Hamilton–Jacobi equation: ( obi þ si ðkrbi k 1Þ ¼ 0; os with si ¼ signðbi Þ: ð6Þ bi ðs ¼ 0; xÞ ¼ /i ðt; xÞ; In practice, Eq. (6) is solved periodically, typically every few time steps, when the function /i becomes too irregular. Eq. (6) can be seen as a pure convection equation. ~ i ¼ si rbi , Eq. (6) becomes: Indeed, if we write U krbi k ( obi ~ i :rbi ¼ si ; þU os ð7Þ bi ðs ¼ 0; xÞ ¼ /i ðt; xÞ: In practice, to define the fictitious time step s, a stability ~ i kDs ¼ Ds ¼ h, where h correcondition is used: kU sponds to the mesh size. So, in fact, the fictitious time step is chosen as h. Instead of using this classical re-initialization method, the formulation developed in Ref. [20] is chosen. The idea proposed by this author is to modify Eq. (5) in order to keep the property of a level set function. The parameter k = ds/dt is introduced (in practice k h/Dt) to write the following equality o ot o ot o ¼ os þ ox r ¼ os ðot þ~ v rÞ and to transform os ot os Eq. (6): obi þ~ v:rbi þ ksi ðkrbi k 1Þ ¼ 0: ð8Þ ot The new formulation incorporating both convection and re-initialization becomes: ( o/ i þ~ v r/i þ ksi ðkr/i k 1Þ ¼ 0; ot 8i 2 f1; . . . ; N G g; /i ðt ¼ 0; xÞ ¼ /0i ðxÞ; ð9Þ which can be written as a pure convection problem: ( o/ i ~ i Þ r/i ¼ ksi ; þ ð~ v þ kU ot 8i 2 f1; . . . ; N G g: /i ðt ¼ 0; xÞ ¼ /0i ðxÞ; ð10Þ The problem is now defined by the system (10), but the common velocity field in X must be defined. In order to avoid kinematic incompatibilities, it is necessary to work with one same velocity field for all the level set functions, whereas the kinetic law is built according to
M. Bernacki et al. / Scripta Materialia 58 (2008) 1129–1132
parameters which are specific to each grain (bulk stored energy and unit normals). Geometry of the interfaces is known to be well described by the level set method [15], so it is well adapted to the calculation of the unit normal, or the curvature if it were needed. Considering grain Gi, the outward unitary normal of any isovalue of /i is defined by: ~ ni ¼
r/i kr/i k
¼
kr/i k¼1
r/i :
ð11Þ
At any time t, the level set approach allows, very easily, to describe the position of grains. Indeed, since the set {Gi nCi,1 6 i 6 NG} corresponds to a partition of the domain X, and, keeping in mind the sign convention described previously, the following results are verified for a point M with coordinates x: M 2 Ci () /i ðt; xÞ ¼ 0;
ð12Þ
M 2 Gi Ci () /i ðt; xÞ > 0 () 8j 6¼ i/j ðt; xÞ 6 0:
ð13Þ
The description of an interface by a level set function remains a ‘‘fuzzy” description of the interface, which means that the level 0 of the level set function (the interface) does not necessarily correspond to nodes of the mesh. The point M is located using the following relationship: M 2 Gi () /i ðt; xÞ ¼ max ð/k ðt; xÞÞ: 16k6N G
ð14Þ
A classical way to build a velocity field from the stored energy field for any node of the mesh is as follows: – find i with /i ðt; xÞ ¼ max16k6N G ð/k ðt; xÞÞ – compare /i(t,x) to a positive fixed parameter ‘ which defines a length scale related to the proximity of the node to the interface:
1131
where a is a positive fixed parameter calibrated to obtain a negligible exponential term outside the anisotropic part of the mesh (see Fig. 2). Expression (16) leads to a smoother velocity field and avoids the need to identify the neighbouring grains of each grain. Figure 1 illustrates the difference between Eqs. (15) and (16) for a 2-D triple junction using a common unstructured mesh. The chosen stored energy field is described in Figure 1. ~ v is shown at the boundaries; the smoothing effect of Eq. (16) appears clearly. To solve the problem defined by Eqs. (10) and (16), a finite element method based on the C++ library ‘CimLib’ [20] is used, with unstructured meshes and a stabilized P1 solver (SUPG or RFB method). The quality of the described approach is strongly related to the accuracy of level set function calculations around the interfaces. Indeed, a little disturbance of the levels around an interface leads to an error in the velocity estimation and consequently in their own evolution. For a given computational cost, optimal accuracy is obtained by using anisotropic meshes, with refinement close to the boundaries. The technique used to generate anisotropic meshes adapted to the grain structures is presented in Ref. [3]. Figure 2 illustrates the evolution of the triple junction presented in Figure 1 after 0, 90, 145 and 240 time steps, with automatic adaptive anisotropic remeshing. Results are very accurate and correspond to one of the exact solutions of this grain configuration (see Ref. [19]). More precisely, the L2 error between the exact and approached function defined by max06t6T end ðmax16k6N G ð/k ðt; xÞÞÞ is less than 3%. A second test case corresponds to a five grain microstructure, as shown in Figure 3. The energy of the circular grain is zero, whereas the energy of the four surrounding grains is unity. As for the previous test case, the interest of this configuration is that it is possible to
*
If/i ðt; xÞ > ‘ ) v ðt; xÞ ¼~ 0 else find j with /j ðt; xÞ ¼
*
ð/k Þ ) v ðt; xÞ max 1 6 k 6 NG k 6¼ i
¼ M ij f ð/i ðt; xÞ; ‘Þðej ei Þ~ ni ðt; xÞ
ð15Þ
with f a decreasing continuous function varying from 1 to 0 when /i(t,x) varies from 0 to ‘, and Mij the mobility at the interface between grain i and grain j. The main disadvantage of this approach lies in the management of multiple junctions: discontinuous velocity fields are generated, which then lead to convergence problems in the resolution of the convection equation. We propose another algorithm which considers, at any point x, all NG level set functions: Find i with
Figure 1. Representation of ~ v (left) for the formulation defined by Eq. (15) and (right) for the formulation defined by Eq. (16).
vðt; xÞ /i ðt; xÞ ¼ max ð/k ðt; xÞÞ ) ~ 16k6N G
¼
NG X
nj ðt; xÞ; M ij expðaj/j ðt; xÞjÞðei ej Þ~
j ¼ 1j 6¼ i ð16Þ
Figure 2. Evolution of the ±5 103 isovalues of the three distance functions with adapted anisotropic mesh.
1132
M. Bernacki et al. / Scripta Materialia 58 (2008) 1129–1132
Figure 5. Evolution of the ±5 103 isovalues of 10 distance functions. Figure 3. Evolution of the ±5 103 isovalues of the five distance functions with adapted anisotropic mesh.
compare the numerical results with the exact solution (defined by the growth of the circular grain with a unit normal velocity). Figure 3 shows the topology after 0, 90, 145 and 240 time steps using adaptive anisotropic remeshing. A very good agreement is obtained between the ±5 103 isovalues of the five distance functions and the exact solution (in black). More precisely, the L2 error defined as for the previous test case is less than 2%. One of the prominent advantages in using front capturing methods for describing interface motions is that there is no need for additional features when some regions (grains) disappear. Complex topological evolutions are handled automatically. In a similar way, it is possible to introduce new regions (grains), based on given criteria. A very simple method to create a nucleation site is to build a new level set function at the desired spatial location x and at a specific time increment. For example, a new distance function can be built such that the boundary of the nucleus is circular (in two dimensions) or spherical (in three dimensions). The circle or sphere can be centered around one node of the mesh with zero stored energy inside such that spontaneous growth occurs. A new level set function can be associated to a new nucleus, which evolves subsequently according to the principles described by Eq. (16). Figure 4 describes the creation and the evolution in three dimensions of one recrystallization nucleus, using the same triple junction configuration as in Figure 2. The chosen location is the node which is the closest to the triple junction. The criterion for nucleation is geometric here, but could be based on the distribution of any mechanical/physical quantity that has been stored in the finite element mesh, e.g. from previous modelling of the polycrystal deformation (crystal plasticity, etc.).
Figure 4. Example of nucleation in a 3-D three-grain case.
Another initial configuration of ten grains with an adapted anisotropic mesh is presented in Figure 5. A different stored energy is assumed in each grain. Three triple junctions are considered as possible nucleation sites. The results illustrate the capability of the level set method to model primary recrystallization in polycrystals. It has been shown that a level set framework, associated with adaptive anisotropic automatic remeshing, is a promising tool to describe primary recrystallization in a polycrystalline material. Full kinematic compatibilities can be enforced. Nucleation and growth phenomena can be implemented in this framework. Various studies are in process, including simulations with a statistical number of grains, the prediction of recrystallization kinetics and the comparison with other models. Future work will account for anisotropic mobility of grain boundaries (as related to the crystallography of the interface). Stored energy due to prior deformation of the polycrystalline aggregate will be used from crystal plasticity finite element simulations, and physical nucleation criteria will be studied. [1] A.D. Rollett, Prog. Mater. Sci. 42 (1997) 77. [2] G. Kugler, R. Turk, Comput. Mater. Sci. 37 (3) (2006) 284. [3] M. Bernacki, Y. Chastel, H. Digonnet, H. Resk, T. Coupez, R.E. Loge´, Comput. Methods Mater. Sci. 7 (1) (2007) 141. [4] R.E. Loge´, M. Bernacki, H. Resk, H. Digonnet, T. Coupez, in: Recrystallization and Grain Growth III Proceedings, Korea, 2007, pp. 1133–1138. [5] M.A. Miodownik, J. Light Met. 2 (2002) 125. [6] A.D. Rollett, D. Raabe, Comput. Mater. Sci. 21 (2001) 69. [7] L.Q. Chen, Scripta Metall. Mater. 32 (1995) 115. [8] H.-K. Zhao et al., J. Comput. Phys. 127 (1996) 179. [9] G.N. Hassold, E.A. Holm, Comput. Phys. 7 (1993) 97. [10] D. Raabe, Philos. Mag. A 79 (1999) 2339. [11] J.B. Collins, H. Levine, Phys. Rev. B31 (1985) 6119. [12] L.Q. Chen, Annu. Rev. Mater. Res. 32 (2002) 113. [13] A. Karma, Phys. Rev. Lett. 8711 (11) (2001) 115701. [14] J.A. Sethian, Level Set Methods, Cambridge University Press, 1996. [15] S. Osher, J.A. Sethian, J. Comput. Phys. 79 (1988) 12. [16] M. Sussman et al., J. Comput. Phys. 114 (1994) 146. [17] B. Merriman, J.K. Bence, S.J. Osher, J. Comput. Phys. 112 (2) (1994) 334. [18] G. Kugler, R. Turk, Acta Mater. 52 (2004) 4659. [19] F. Reitich, H.M. Soner, Proc. Roy. Soc. Ed. 126A (1996) 837–865. [20] T. Coupez, in: Numiform’07 Proceedings, Plenary Lecture, Porto, 2007, pp. 61–66.