Introduction to the level-set full field modeling of laths spheroidization phenomenon in α/β titanium alloys

The fragmentation of α lamellae and the subsequent spheroidization of α laths, in α/β titanium alloys, are well known phenomena, occurring during and after deformation. We will illustrate the development of a new finite element methodology to model these phenomena. This new methodology is based on a level set framework modeling the deformation and the ad hoc concurrent or subsequent interfaces kinetics. In the current paper, we will focus on the modeling of the surface diffusion at the α/β phase interfaces and the motion by mean curvature at the α/α grain interfaces.


Introduction
Two-phase α/β titanium alloys are materials with numerous applications in different industrial domains, mostly due to their attractive mechanical properties. They present high strength, fatigue and corrosion resistance, low density and good ductility [1].
Titanium alloys exhibit different microstructures depending on the applied thermomechanical path. Spheroidization is a phenomenon observed at α/β titanium alloys with initial stable microstructure of α lamellae inside β grains, during deformation and subsequent thermal treatment. More precisely, during this phenomenon, the fragmentation of α lamellae and the subsequent spheroidization and coarsening of α laths is observed [2].
Spheroidization has received considerable attention due to its importance in microstructural control. The new spheroidized microstructure shows enhanced strength and ductility; so evidently, this phenomenon raises high interest for the industrial applications [1,2].
In this paper, we will illustrate a new finite element (FE) methodology in order to model these microstructural evolutions. Our interest is focused on the predominant mechanisms occurring during the first stages of the spheroidization at the lamellae interfaces without considering the microstructure deformation modeling. The α/α grain interfaces are introduced arbitrarily leading to surface diffusion at the α/β phase interfaces and motion by mean curvature at the α/α grain interfaces. For the purpose of modeling efficiently these interfacial kinetics, a level set framework was introduced. If some studies concerning the full field modeling of coarsening mechanism exist [3], at the authors' knowledge, the 2D or 3D full field modeling of the first steps of laths spheroidization in α/β titanium alloys is a poorly researched subject on the state of the art. Some basic cases of surface diffusion will be detailed in order to introduce the numerical methodology. A first case of immersion of a real microstructure from experimental data will also be considered. Finally, a simple case of lamellae splitting due to the interaction of motion by surface diffusion at the α/β interface and motion by mean curvature at the α/α sub-boundaries will be described.

The physical mechanisms
According to Semiatin [2], several mechanisms are activated simultaneously during thermomechanical process leading to spheroidization. During deformation, it is observed that high angle misorientations are formed inside the α lamellae. Subsequently these misorientations lead to the formation of grooves which imply the penetration of the β phase and lead to the splitting of the α lamellae (see Fig. 1). After the splitting of the α lamellae in smaller parts, a microstructure with α pancake shape particles is obtained. In order to reduce the interfacial energy, the α particles tend towards a shape where the surface area is minimized for a given volume i.e. a spherical one. Then, coarsening (i.e. volume diffusion well known as Ostwald ripening [4]) can take place. To be more precise the spheroidization of the microstructure is a consequent event of the four combined following mechanisms: & crystal plasticity during deformation which introduces the α/α sub-boundaries inside the α lamellae, & motion by surface diffusion at the α/β interfaces, & motion by mean curvature at the α/α interfaces, & coarsening.
Except the mechanism of crystal plasticity that it obviously occurs during hot deformation, there is no clear distinction of when each one of these mechanisms occurs. At this point it is very important to specify that during spheroidization, we do not have phase fractions evolution.
In the following, we will focus on the mechanisms which are responsible for the splitting of the α lamellae. Grooving is usually initiated by atomic scale processes near the region of α/α grain boundaries and the α/β phases intersection. These kinetics can be described by motion by surface diffusion at α/β interfaces and motion by mean curvature at the α/α sub-boundaries (Fig. 2).
These two competitive mechanisms, occurring simultaneously, lead to the lamellae splitting. After the splitting of the α lamellae, surface diffusion at the α/β interfaces of the laths becomes the dominant mechanism before the coarsening by volume diffusion. Next section illustrates experimental data describing this sequence.

Experimental analysis
Some hot compression tests have been realized (TA-6 V alloy) and are detailed below. We have worked with double coned samples and three experiments have been performed (Fig. 3) at 950°C. After the described thermomechanical treatments, a middle cut is realized and the microstructure is observed by SEM-EBSD analysis exactly at the center of the samples corresponding to the maximum strain area.
All the experiments start with a thermal treatment of 30 min at 950°C in order to reach equilibrium of the volume fraction of the two phases. The first experiment involves only deformation at 950°C (ε=20%, ε =0.1s −1 ) in order to visualize the microstructure exactly after deformation. Additional 15 min (respectively 1 h) of annealing at 950°C is considered for the second (respectively the third) experiment. The purpose is then to quantify the effect of annealing time on the microstructure evolution. Figures 4,5,6 illustrate representative SEM images from each stage of the microstructure evolution after each experiment. Just after deformation, we can see a not so fragmented microstructure but from the variation of the grey scale color inside the lamellae we can guess a lot of sub-boundaries (Fig. 4). After 15 min of annealing we see changes on the shape but not so sever ones. The shape evolution of the lamellae is obvious but we cannot observe spheroidized microstructure (Fig. 5). After one hour of annealing we can see finally a more spheroidized and coarsening microstructure (Fig. 6).
In order to discuss quantitatively of the microstructure evolutions, an important number of such representative images in Fig. 1 Splitting of α lamellae into α laths terms of lath area and lath aspect ratio evolutions have been analyzed. Basically, by approximating the α laths as ellipses, the aspect ratio is then defined as the ratio between the major axis and the minor one.
The results obtained, by measuring approximately 500 α laths for each case, are summarized in Fig.7. They illustrate a fast evolution of the mean shape ratio during the first 15mins with a quasi-stable mean surface value. During 15 mins and 1 h, the mean shape ratio decreases slowly whereas the increase of the mean surface ratio is important. These experimental results are coherent with the involved mechanisms [1,2]: & After deformation, sub-boundaries inside the lamellae are present and sphereodization is not activated. & During the first minutes of the thermal treatment, surface diffusion and motion by mean curvature are the predominant mechanisms which lead to the laths appearance and sphereodization of the laths without a real growth of the mean lath surface (limited effect of the volume diffusion). & After this first step, coarsening becomes the predominant mechanism.
From the observations above, we believe that, modeling the first steps of spheroidization implies the necessity to simulate the boundary motion due to surface diffusion at α/β interfaces and due to the mean curvature at the α/α sub-boundaries. Subsequent modelling of coarsening is of course of prime importance to predict the laths volume time evolution but is only a natural perspective of the present works which are, at yet, dedicated to the two first mechanisms.
Next section illustrates the governing equations of surface diffusion and motion by mean curvature.

The physical equations
Motion by surface diffusion According to Mullins [5], in order to describe the atoms flow at the α/β interface we can consider a surface flux j ! : where ν denotes the number of drifting atoms per unit area and υ ! denotes the average velocity of these drifting atoms.
Assuming local equilibrium, we can express υ ! with the Nerst-Einstein formula as following:  where D αβ denotes the surface diffusivity of the α/β interface, μ the chemical potential, k the Boltzmann constant and T the absolute temperature. The ∇ s operator corresponds to the surface gradient operator defined as the tangential component of the gradient: with n the outward-pointing unit vector normal to the surface and P = I − n ⨂ n. By considering Eq.(1) and Eq.(2), the following equation is obtained: Assuming that there is mass conservation, the surface motion can then be described by: where v n ¼ υ * •n denotes the normal velocity of the surface and Ω the atomic volume. By combining Eq. (4) and Eq. (5), we obtain: with Δ s = ∇ s • ∇ s the surface Laplacian operator (or Laplace-Beltrami operator). Considering κ as the mean curvature (sum of the main curvatures in 3D) and γ αβ the α/β interface energy and by ignoring the possible effects of anisotropy, the following relationship is also verified: Thanks to Eq.(6) and Eq. (7), we obtain: , the kinetic coefficient. Eq. (8) describes the relation between the motion by surface diffusion and the surface Laplacian of the mean curvature [6][7][8][9].

Motion by mean curvature
In order to describe precisely the surface evolution of an α lamella, the influence of the mean curvature should be also considered [8].
The grain boundary energy is indeed important for the lamellae splitting. The grain boundary energy is given by the well-known Gibbs-Thompson relationship where the normal velocity v n of the grain boundary is described proportionally to the mean curvature κ: , where γ αa denotes the grain boundary energy, b is the burgers vector norm associated with the hoping event, νˇis the Debye frequency, R the gaz constant and Q the apparent activation energy.

Motion of the interfaces
Surface diffusion and mean curvature motions are taking part simultaneously during the phenomenon of spheroidization. A global velocity combining both of these motions can then be summarized as: with χ α/β the characteristic function of the α/β phase interfaces and χ α/α the characteristic function of the α/α grain interfaces (Fig. 2).

Level set formulation
A level-set (LS) model was formulated in order to deal with the topological changes at the α/β interfaces. The level set method was chosen due to its capability to immerse/describe/capture easily in a FE context the interfaces [10][11][12][13] and also due to the fact that geometrical quantities as the mean curvature κ and the outside normal n can be obtained as: and φ is then defined over the domain Ω as a signed distance function to the interface Γ of the subdomain of interest that we will denote Π.
The sign convention of Eq. (11) corresponds to a distance function negative inside Π and positive outside.
Thus the interface velocity can be rewritten in a level set form as: with: It can also be proved that in the considered level-set formulation [6,14], v n can be re-written as: The velocity is then defined in the entire domain and corresponds in the vicinity of the zero level-set function of φ, i.e. Γ, to the interfaces velocity [6,7]. A new numerical methodology of resolution is detailed in the following section concerning the surface diffusion mechanism.

A surface diffusion methodology
For the modelling of the induced flow from the surface diffusion mechanism at the α/β interface, a FE methodology is adopted. At any time t, the transport velocity v * is defined by: The B coefficient is defined as a constant and it is chosen to neglect any anisotropy concerning the interface energy and the diffusivity. Additionally, isothermal conditions are assumed. The time evolution of Γ(t) is then obtained by solving the following convective system: The interface can then be obtained at each time step as the 0-isovalue of the distance function and the velocity is updated by following Eq.(18) before the following time step. At the following subsections, more extensive details are given for the resolution algorithm.

Surface diffusion velocity identification and transport resolution
The methodology used to obtain the surface diffusion velocity is based on the FE based strategy introduced by Bruchon et al. in [6,15].
Indeed, as P1 description of the LS is considered in the proposed methodology, one of the basic issues in the problem of surface diffusion is that the velocity is defined by the Laplacian of the curvature, which means that the velocity is a function of the fourth order spatial derivative of φ. The numerical strategy proposed by Bruchon et al. consists to solve this problem by considering a Bregularized^formulation. More precisely Eq. (18) is solved in a weak form by using a FE formulation. Further information can be found in [6,15].

Convection-Reinitialization methodology
By assuming the appropriate calculation of the surface velocity, the traditional strategy of convection and subsequent reinitialization steps is used. The main idea is to solve the advection equation and to rebuild the metric properties of the LS function in order to keep a distance function (‖∇φ(x, t)‖ = 1) at least near the interface Γ(t). Classical approaches consist in solving, separately, the convective part and the reinitialization part thanks to the resolution of a classical Hamilton-Jacobi system [13] or to adopt an unified advection and renormalization methodology by solving one single equation based on a smooth description of the LS [16]. Here, a new approach is proposed.
This new approach is based on a separate resolution of the transport and of the reinitialization part. Convective equation is firstly solved thanks to a stabilized P1 solver (SUPG or RFB method). Then, a parallel and direct reinitialization algorithm detailed in [17], which has been proven to be extremely fast and accurate, is used. In this algorithm, the Γ(t) interface is firstly discretized into a collection of segments (respectively triangles in 3D) and the nodal values of the level-set function are then updated by finding the nearest element of the collection and calculating the distance between the considered node and this nearest element. This method takes advantage of a space-partitioning strategy using k-d tree and an efficient bounding box strategy enabling to maximize the numerical efficiency for parallel computations.
Moreover, this methodology presents two other benefits: & it enables to avoid the validation/calibration of unphysical parameters necessary to reinitialize the LS function [13,15], & it enables to obtain directly an exact P1 description of n [18] before to solve Eq. (18) rather than following the classical less precise methodology where the normal is computed by performing a P1 interpolation of the first derivative of the LS function.

First academic case
In this section we examine an ellipsoid shape under surface diffusion. We want to test the efficiency of the proposed formulation for a simple case with an analytical solution. Considered computational domain is a 1 mm × 1 mm square centered in (0,0). An initial ellipse (a = 0.3 mm and b = 0.2 mm) of equation x 2 /a 2 + y 2 /b 2 = 1 is considered. Of course, this shape is going to evolve towards a circle shape while conserving its area. Thus limit radius, i.e. limit value of a and b is given by the value ffiffiffiffiffiffiffi ffi πab p ≈0:43416mm. Initially, in order to test the efficiency of the approach without dealing with the effect of meshing/remeshing operations in the results, a fixed mesh is considered during simulations. An initial isotropic mesh adaptation is considered in a ring centered in (0,0) and defined as 0.19mm ≤ r ≤ 0.31mm in order to keep a very fine mesh size, defined as h in Table 1, in all the zone crossed by the zero isovalue of the LS function during the simulation. The Fig. 8 illustrates the FE mesh used (a), a zoom on the FE mesh (b), the initial distance function field (c), the curvature field near the interface (d) and the normal velocity field near the interface obtained initially thanks to the FE resolution of Eq. (18) (e). Red or white line corresponds to the initial ellipse interface (0 isovalue of the LS function).
In all simulations B is assumed to be homogeneous and equal to 1 mm 4 /s. Exact velocity of the point (a(t),0) is known [15] and given by: Hence a(t + dt) can be easily evaluated thanks to a forward Euler method: and b(t + dt) can be easily obtained by verifying the area conservation at anytime: This method is then used with a time step of 1 ms to evaluate the Bexact^evolution of a and b values during the shape evolution. Concerning FE simulations, at each time step, the positions of (a(t),0) and (0,b(t)) are determined on the zero-isovalue of the distance function and then compare to the Bexact^solution. Final time of all simulations is fixed to 1 s. Table 1 summarizes all the parameters (time step, mesh size in the fine mesh zone, number of elements of the used mesh, method used), the CPU time of the simulations and the corresponding precision of the results obtained concerning the positions of (a(t),0) for t ∈ [0, 1s] by using the unified convective-renormalized approach described in [15,16] and the new one proposed here with different time step and h values. These cases are representatives of an important number of other performed simulations. Errors are defined as: where i denotes the discretization in time. From different simulations (Cases 1 to 4 of Table 1 are representatives), we can summarize the following comments: -Considering the precision of a(t) predictions (see e 1 and e 2 errors on Table 1), both approaches are relevant to model surface diffusion. Figure 9 illustrates, at t = 0.2 s and t = 1 s, the difference of the exact interface and the results obtained with the Case1. -The time step of the Case4 (dt = 0.1 ms) is the maximal value usable for h = 1μm and the unified formulation to avoid numerical instabilities. Such instabilities were not identified for the new proposed approach even for important time step and mesh size (Case 3 for example). -As illustrated in Table 1, calculation time of the new proposed approach is then clearly very attractive comparatively to the unified approach. -Even for the new proposed approach, decreasing the dt and h values below, respectively, 10 ms and 2 μm seem not improve the results quality (e 2 ≈ 3 % ). It can be

Second academic case: volume loss and mesh adaptation
Next, we consider a more realistic shape of a long ellipse with a = 0.5 mm and b = 0.1 mm and mesh adaptation. Indeed, in order to obtain acceptable calculation time to model surface diffusion of complex microstructure, a fixed meshing strategy as the one of the previous section is not an option. A meshing and remeshing strategy must be used.
To begin with, two metric based meshing strategies associated with the MTC topological mesher were tested. MTC is a P1 automatic remesher based on elements topology improvement that was developed for Lagrangian simulations under large strains [19]. This tool was extended to anisotropic mesh adaptation [19], for which it was extensively used in context of FE microstructure description [6,12,13,[15][16][17].
The first metric considered, is the metric described in [16]. This metric enables to obtain isotropic or anisotropic (in the normal direction of the interface) fine mesh in the vicinity of the interface without considering its curvature. In this paper, isotropic adaptation was considered (Method 1 in Table 2).
The second metric considered, is based on an a priori error estimator linking the interpolation error on the LS function to its gradient vector and hessian matrix. As already described, these variables allow to obtain the normal vector to the interface and its main curvatures. Using these data, a metric field can be built in order to obtain a fine mesh size in the vicinity of the interface depending on local curvatures [20] (Method 2 in Table 2).
Simulation and data results are reported in Table 2 for both meshing/remeshing methods. The BConv + exact Reinit^strategy was used. Final time of both simulations is fixed to 1 s. As the mesh size near the interface is of the same order than the mesh   size used in Case1 and Case2, same precision could be expected. However, remeshing is also synonymous of diffusion concerning the FE fields carrying the interface at each remeshing operation necessary to follow interface motion. Then, volume conservation was tracked for these cases. Results described in Table 2 illustrate that the mesh adaptation techniques come with good precision and fast calculation times. Meshing adaptation based on the main curvatures enables to obtain a very good conservation of the volume. This aspect will be of course very important for real configurations where the initial thin shape of the α lamellae implies very high ratio between the minimal and the maximal values of the interface main curvatures. So, this meshing strategy in terms of metric seems particularly indicated.
If the results described previously in terms of calculation time and precision could be sufficient to study the surface diffusion of one thin ellipse, it seems clear that to further improve our numerical framework, it is important to deal with real 2D or 3D configurations.
In order to face this problem, we adopt a new topological mesher, Fitz, developed by Shakoor et al. [20]. With Fitz, a body fitted meshing and remeshing is possible while using complex metric as previously described. It was then proved in [21] that this new mesher associated with a volume conservation constraint which is compatible both with implicit and body-fitted interfaces, the modification of the interface due to remeshing can be delayed enough so that a Lagrangian LS method becomes more than interesting compared to an Eulerian LS method, even when large deformations or displacements occur.
First tests are very promising, allowing in the Case6 configuration to obtain a precision of 2% concerning the volume conservation with a calculation time of 30s on 12 CPU. Figure 10 illustrates a result obtained with this numerical framework for a Ba = 0.5 mm and b = 0.05 mm^configuration.
Moreover we anticipate, with this mixed implicit/explicit description of the interfaces, an easier coupling between surface diffusion at the α/β interfaces and motion by mean curvature at the α/α grain interfaces.

Surface diffusion in real α/β microstructures
Having tested the validity of the surface diffusion solver in simple configurations, we consider now a more complex one. The FE immersion of a real microstructure, from a SEM experimental image, is investigated.
The considered microstructure of size 1.1mm × 0.7mm, described by the Fig. 11, is one similar to the configuration of the  . This configuration is interesting as splitting can be considered as achieved (see section 3). Binarized version (see Fig. 12) of the microstructure described by Fig. 11 and the signed distance function of the laths have been obtained thanks to the ImageJ software and have been immersed in an initial FE mesh thanks to our FE C++ library CimLib [22]. Some laths, near of the domain boundaries, were deleted in order to avoid border effects. A first Eulerian LS framework BConv + exact Reinit^simulation has been considered with representative physical parameters [2]. We experienced some difficulties regarding the shape evolution of the interface during the simulation. Having numerous laths interacting very quickly under surface diffusion, lead to oscillations of the interfaces and severe volume loss.
From the other hand, performing the same simulation with the Lagrangian LS framework based on Fitz methodology, lead to much more reasonable results considering the shape evolution of the lamellae and volume loss as illustrated by the time evolution of Fig. 13. The size of the mesh elements was fixed to h min = 1μm close to the lath interfaces and to h max = 5μm otherwise (see Fig. 14). The time calculation was 2 h for 12 CPUs. The volume loss from the beginning of the process to the end was limited to 2%. More quantitative validations of these simulations and discussions of the A and B coefficient values will be detailed in a forthcoming publication.
Interestingly, laths coalescence events are observed during surface diffusion. It is logical as one LS function is considered for all the laths (as in [3] in context of phase field description of the interfaces). However, at the authors knowledge, this kind of evolutions was not reported in the state of the art of α/β titanium alloys during the first steps of the spheroidization. Once again, it seems that consider the role of α/α grain interfaces are of prime importance to predict a realistic evolution of the contacts between the laths during surface diffusion phenomenon.
First simple case of lamellae splitting: As already mentioned, two main competitive mechanisms are involved in lamellae splitting. Motion by surface diffusion at the α/β interfaces and motion by mean curvature at the α/α sub-boundaries (see Fig. 2). The numerical framework based on the Fitz meshing/remeshing tools and the enhanced lagrangian methodology is now definitively adopted. The simple configuration of a lath crossed by a α/α grain interface as described by Fig. 15 is considered.
Eq. (17) is taken into account in order to evaluate the global velocity which is used to convect the LS functions (L1 and L2) describing the two laths. Normalized velocity due to surface diffusion (respectively due to motion by capillarity) is described by the Fig. 16 (resp. the Fig. 17).
A realistic ratio between the A and B coefficients is used in order to obtain a representative normalized evolution of the two laths. Fig. 18 illustrates the mesh used for the simulation (conform mesh at the interfaces and a ratio of five between the size of the coarse mesh elements and the fine ones). This figure clearly illustrates the impact of the mean curvature at the α/α grain boundary concerning the fast appearance of dynamic triple junctions. Even in context of small misorientations and so in context of weak values of the γ αa parameter, the mean curvature remains extremely high at the initial Tjunction between the α/β interface and the α/α interface.
Finally Fig. 19 illustrates the splitting of the lamella in two laths due to surface diffusion and motion by capillarity.
This case demonstrates that the proposed formalism enables to deal with surface diffusion and motion by mean curvature. The volume loss during this simulation was around Fig. 19 Splitting modeling, from top to bottom: initial configuration, at t = 0.1t end , t = 0.5t end , t = t end

Conclusions and perspectives
The first steps of a new FE numerical framework dedicated to the modelling of the mechanisms of spheroidization in α/β titanium alloys have been detailed. This numerical framework has been illustrated with some basic numerical cases on surface diffusion of α laths for checking the efficiency of the methodology. Optimal algorithm of resolution and meshing strategies were proposed in terms of precision and calculation time. A new way was also opened by considering a mixed implicit/explicit description of the interfaces with the use of body fitted meshes.
A surface diffusion case on a real microstructure obtained from experimental images was validated in terms of numerical efficiency. We also demonstrate that the proposed formalism enables to simulate lamellae splitting due to the competitive mechanisms of motion by surface diffusion and motion by mean curvature.
The perspectives of this works are numerous. This article illustrates a first step to propose full field simulation of sphereodization and coarsening in α/β titanium alloys. First of all, more representative simulations in terms of domain size must be realized and the obtained results must be validated with in-situ experimental results. In this context, the coefficients A, B need to be finely calibrated in order to respect realistic kinetics. Concerning the involved mechanisms, volume diffusion must be added to the global velocity in order to predict realistic volume evolution of the laths during spheroidization. It seems also important to highlight that the proposed methodology is usable in 3D without additional developments. At the authors knowledge, 3D full field modeling of laths for α/β titanium alloys was never addressed. It is then an exciting perspective of these developments. Finally, we plan also to simulate the deformation step thanks to an existing crystal plasticity finite element approach developed in a LS context [23].