Phase field simulation of α/β microstructure in titanium alloy welds

The microstructure formation during β → α + β transformation in heat affected zone of titanium alloy welds is studied using an approach combining a phase field solid-solid transformation model with a heat transfer finite element method (FEM). FEM is used to model macroscopic heat transfer during welding cycle and to compute the thermal history at several points across the weld. The thermal history is subsequently used as input to the phase field model describing microstructure evolution. The chemical component of Gibbs free energy, atomic mobility and elastic tensors are parameterized for Ti–6Al–4V using available literature data. A classical nucleation theory based model is parameterized using recent continuous cooling experiments and is used to account for nucleation events. The study explains the graded microstructures observed in titanium alloy welds and provides insights into the underlying processes that occur during welding.


Introduction
Combination of light weight, high strength, good corrosion resistance and ability to sustain high temperatures makes titanium alloys increasingly attractive for a multitude of applications-especially so in the aerospace industry [1][2][3][4][5][6]. Many of the technological applications of titanium alloys require joining of parts by welding. Welding not only plays an important role during fabrication, but also in maintenance and repair of components. During the welding process, a region of the base material close to the weld is substantially heated such that in a fusion zone (FZ) the material melts and solidifies, and in an adjacent heat affected zone (HAZ) no melting occurs but the microstructure can be affected [7][8][9][10]. The rapid heating and cooling in the HAZ may alter the grain structures and also change the relative proportion of present phases. Moreover, different regions of welded parts, depending on the distance from weld line, experience different thermal history, which leads to a graded microstructure in the vicinity of weld. These changes may be responsible for degrading the mechanical properties of the weld [11], hence detrimental to performance of the entire component. In this work we focus on investigating the processes that occur in the HAZ of titanium alloys using a framework that combines a heat transfer finite element model and a phase-field solid-solid transformation model.
Depending on heating/cooling rates, α/β titanium alloys such as Ti-6Al-4V may undergo a transition from β phase (body centered cubic -BCC) to the α phase (hexagonal closed pack -HCP). As a result, a two phase α/β lamellar microstructure forms. It is known that the lamellar α/β microstructure contributes to the strength of the material. Furthermore, when the cooling rates are fast, diffusionless transformation leads to the formation of a¢ martensite phase [12][13][14], which may be detrimental for mechanical properties of the weld. Therefore, understanding the processes related to solid state transformations in the HAZ is necessary to predict the mechanical properties in the weld. The heating and cooling associated with welding process are fast and therefore challenging to study in situ microstructural evolution experimentally, leaving simulations as an attractive alternative option. While typical size of the HAZ is a few millimetres, the typical length scale associated with the α/β microstructure is of the order of microns. Therefore, a framework that combines a macroscale heat transfer model [15] with a phase field model that resolves the length scale relevant for the α/β microstructure [16][17][18][19][20][21][22] is used in the present work.
The calculations of the thermal history performed in the macroscale apply established numerical framework of Computational Welding Mechanics (CWM) [23][24][25]. The mesoscale phase field model describes evolution of the α phase precipitates in β phase matrix with all possible variants and includes elastic effects arising from transformation strains. The evolution of compositions of titanium, aluminium and vanadium elements is studied using multi component diffusion equations. Although martensite microstructure is not explicitly modelled, regions where martensite can form are identified based on the temperature and local value of the alloy composition. The nucleation phenomena is incorporated using the classical nucleation theory, parameterized to experimental critical cooling rates required to avoid martensite formation [13,14]. The microstructural evolution is simulated in representative volume elements (RVE) placed at preselected regions in HAZ. Each of these simulations uses the temperature versus time data (thermal history) obtained from the heat transfer FEM calculations.
The paper is organized as follows. Section 2 presents details of the FEM heat transfer model and the calculation of the thermal history for a butt weld. The phase field model along with a nucleation model is described in section 3. An approach to incorporate martensite formation in the simulation framework is also discussed. In section 4, results of our phase field simulations of microstructural evolution in the HAZ are presented. The limitations of the approach and significance of the results are discussed in section 5. Section 6 concludes the paper.

Heat transfer model
In order to predict the microstructures in the weld, it is important to compute the thermal history in the weld for a given set of processing parameters. The method for computing thermal history at some representative points across the weld follows methodology used by Deng and Murakawa [26] and Yaghi et al [27] amongst many others. Figure 1(a) depicts the geometry of butt-weld. The part in grey represents the base metal of the titanium alloy to be welded together, while that in blue represents the weld material added during welding. The direction of weld is highlighted by the red arrow and the welding speed is taken to be 1.0 mm s −1 . During welding, the governing partial differential equation for the three-dimensional heat dissipation is described by where K is the thermal conductivity, ρ the density, C p the specific heat, T the temperature, t the time and  Q is the rate of moving heat generation per unit volume due to the weld torch. The heat generation per unit volume, or heat flux Q, is related to the applied voltage V (volt) and current I (amps) of the weld torch through h = Q VI where η is the efficiency of the welding process. In this thermal analysis, V, I and η are taken to be respectively, 25 volts, 30 amps and 0.9. The temperature-dependent physical properties required for the thermal analysis, specifically thermal conductivity (W/(m-K)) and specific heat (J/(kg-K)), are obtained from MPDB v7.37 [28]. The solidus and liquidus and latent heat per unit mass, are taken to be, respectively,1604°C, 1668°C and 392 kJ kg −1 . It is assumed that the physical properties of the weld material are the same as the base material. Both convection and radiation heat transfer on all free surfaces during the welding process have been taken into account, employing a lumped heat transfer coefficient [29,30]. Figure 2 shows the thermal history obtained at different locations in the weld along the line on the mid plane shown in figure 1(b). The model described in section 3 uses temperature dependent thermodynamic and kinetic data obtained using CALPHAD method to simulate the microstructure. At each location in the weld, the time evolution of the temperature shown in figure 2 is used as an input to the phase field model.

Microstructure model
Titanium alloys can exhibit a two-phase coexistence of the HCP a and BCC b phases. The a phase becomes unstable above a critical temperature b T (β transus). There are twelve equivalent ways to form HCP product phase from BCC parent structure (variants). It is also known that under certain conditions a¢ martensite phase may appear. This section describes the ingredients of the phase field model, the scheme to identify regions with martensite formation, and the nucleation model of α precipitates in β matrix.  The β phase corresponds to the case when all components of the order parameter vanish. The free energy describing the transition is given as where f int represents the interfacial contribution to the free energy, f chem is the chemical free energy and f el is the elastic energy associated with the transformation. The interfacial contribution is expressed in terms of the order parameter components as Here, the first term is a double well potential that represents an activation barrier between α and β. The second term is the usual gradient term that describes the energy penalty of creating interfaces. The barrier height W and the gradient coefficient h K are related to the α-β interfacial energy and interfacial width Δ with The interfacial width is assumed to be d D = 5 , where d is the grid spacing, taken to be d m = 0.05 m. The third term in equation (3) ensures that for a large enough value of x, only one of the twelve order parameter components is non-zero.
The chemical free energy takes standard form: Al, V, is an interpolating function. V m is the molar volume, taken to be = -V 10 m . m 5 3 G α and G β are temperature dependent Gibbs free energies obtained from thermodynamic databases (The reader is referred to supplementary information is available online at stacks.iop.org/MRX/7/ 046517/mmedia for more details). X Al,α/β and X V,α/β are the aluminium and vanadium composition variables that serve as auxiliary variables following the method introduced by Kim, Kim and Suzuki [32]. Note that we have utilized the condition X Ti +X Al +X V =1 and expressed the Gibbs energies only in terms of the aluminum and vanadium compositions. The auxiliary compositions are computing by using the conditions ⎛ where μ Al,Ti =μ Al -μ Ti and μ V,Ti =μ V -μ Ti are the diffusion potentials of aluminium and vanadium respectively. The precipitation of the HCP α phase in the BCC β matrix leads to generation of transformation strains due to the change in lattice parameters. The elastic energy is expressed assuming homogenous elasticity:  The kinetics of the transformation is described by set of equations: Equation (11) describes the evolution of variants and (12)(13)(14) are the multi-component diffusion equations describing the evolution of the compositions [33]. G h is the interfacial mobility and M j represents the atomic mobility of each component (see supplementary information for further information). The equations of motion (11)(12)(13)(14) must be solved subject to the mechanical equilibrium constraint given in equation (15). The elastic displacements and strains are also obtained by solving equation (15). For computational convenience equations (11)- (15) are casted in a dimensionless form. Rescaled variables are introduced as = ¢ t t t , The constant x is rescaled as m 0 A large value x¢ = 1000 is set for this constant to ensure that only one of the order parameters is non-zero. By setting ( The G α and G β are represented in the Redlich-Kister (RK) polynomial form following CALPHAD approach. In the present work, we use the published thermodynamic database for Ti-6Al-4V [34]. In order to obtain auxiliary compositions, a set of non-linear equations (5)-(8) has to be solved. In order to gain computational efficiency, we linearize the numerical problem by adopting a parabolic approximation to the Gibbs free energies. The details of this method are provided in the supplementary materials.
For the diffusion equation, atomic mobilities M j are specified for both α and β phases. For the β phase mobilities, the published database for the Ti-Al-V system [35] is used. To the best of our knowledge, there is no complete mobility data for this alloy for α phase. In reference [19], it was assumed that the diffusivities of the alpha stabilizer aluminium are the same for both α and β phases while the diffusion of the vanadium was not modelled. A similar approach is adopted here but all atomic mobilities in α and β phases are taken to be equal.
A typical welding cycle involves rapid temperature changes, it is therefore important to allow interfacial mobility G h to vary as a function of temperature. To the best of our knowledge, there is no experimental data for this quantity. Therefore, it is chosen such that during slow continuous cooling, the fraction of a single precipitate is close to the thermodynamically predicted fractions over a range of temperatures. Following references [36,37], Arrhenius form is adopted for the rescaled interfacial mobility ⁎ G The reader is referred to the supplementary for further details.

Martensite
It has been shown that fast cooling leads to formation of a¢ martensite by a diffusionless martensitic transformation [12][13][14]. In principle, there is a range of cooling rates which lead to mixed microstructures where a¢ coexists with the / a b microstructure. The a¢ martensite also has 12 crystallographic equivalent variants and forms a complex microstructure. In the present work, we do not model the martensitic microstructure explicitly. The main concern of the present work is to predict whether a given region will form martensite or not. Here it is assumed that martensite will form if a given location has not undergone any change in composition by the time the temperature reaches the martensite start temperature during continuous cooling. Therefore, it is possible to identify regions where the martensite can form, purely based on the local value of temperature and composition. In order to do this, a simple condition is applied: where, is a local fraction of α phase, The variable f m denotes whether or not martensite can form at a given location. T ms is the martensite start temperature equal to 810°C [14], X i is the local composition of a component (i=Al, V), and X i 0 is the corresponding initial (nominal composition). This condition simply assumes that during a given heat treatment, in case the composition at a given spatial point has not changed by the time temperature reaches T , ms that point instantaneously transforms to martensite. If the condition f = 1 m is satisfied at a given time step and grid point (the point is assumed to transform to martensite), the order parameters and diffusion equations are not evolved any further.

Nucleation
According to classical nucleation theory [38], the rate at which α phase is seeded in β matrix can be computed by following expression N V is the density of nuclei and the quantities Z and b* are given by Here D l =M l β RT, where M l β is the atomic mobility, calculated using [35]. X l β is the local composition and X l α,eq is the equilibrium composition of α. ΔG nuc is the driving force for nucleation that is approximated by   is the corresponding chemical potential of the a phase, evaluated at the equilibrium composition of a. The chemical potentials are evaluated by using the full Gibbs free energy with RK form (see supplementary information). The nucleation barrier DG* is given as For homogeneous nucleation, the factor f=16π/3. In the present work, we consider the nucleation to be heterogeneous and consider f as an adjustable parameter. The nucleation density N V is also regarded as adjustable.
The probability of nucleation event at a given site is related to rate J by where t 0 is a characteristic time scale, Δt is numerical discretization of the time step and δ is the smallest length scale. Spherical nuclei of radius 0.2 μm are added at each time with probability given by equation (23). Each of the nuclei is assigned with an orientation chosen randomly from one of the α phase variants. It is ensured that there is no overlap of nuclei with other nuclei as well as pre-existing α and martensite regions are excluded from nucleation procedure. The parameters N V and f are obtained such that the critical cooling rate required to avoid martensite formation is close to that observed in recent continuous cooling experiments [13,14]. The experiments report a critical cooling rate in the range of 2°C-2.5°C s −1 . In the present work, the nucleation parameters are set such that the critical cooling rate is ∼2°C s −1 . The phase field equations are solved using finite difference method for a stress-free system with periodic boundary conditions. A 3-dimensional grid of size´128 128 128 corresponding to a m m ḿ6.4 m 6.4 m 6.4 m volume element is considered. The RVE represents a region deep inside a grain of a polycrystal. Note that due to small size of the RVE, grain boundary nucleation is not taken into account. The initial condition corresponds to the nominal composition of Ti-6Al-4V ). The temperature is continuously varied at different cooling rates, starting from the beta transus temperature =  b T 940 C and finishing at T F =500°C. Several simulations with different choices of the nucleation parameters N V and f were tried. The phase fraction of α is computed as f α =〈f α 〉, the martensite phase fraction f m is calculated as fraction of grid points where the martensite formation condition f m =1 is satisfied and the β phase fraction is calculated as f β =1-f m -f α . A critical cooling rate of 2°C s −1 is obtained using f=0.0001 and N v ∼7.83×10 9 m −3 and so these parameters are used to compute the nucleation rate for all simulations in this paper. The reader is referred to the supplementary materials for a discussion on the cooling rate dependence of phase fractions and microstructures obtained by the continuous cooling simulations.

Simulation of initial (base metal) microstructure
The first step of the simulations is to generate an initial α/β morphology that represents the base metal microstructure. In the present work, the initial microstructure is generated by simulating continuous cooling from the β phase with a very slow cooling rate (0.01°C s −1 ). The simulation conditions are same as those used in section 3 and the nucleation rate is calculated using the parameters f=0.0001 and N v ∼7.83×10 9 m −3 . Figure 3 shows the evolution of the microstructure during the continuous cooling. The corresponding composition images in figure 3 show that at early times, although α domains have started to form, the compositions are still far from equilibrium. As cooling proceeds, the domains grow and the compositions approach the equilibrium values. The microstructure stops evolving at low temperatures. The final temperature microstructure has coarse α domains with thin aluminium depleted and vanadium rich β regions between them. The microstructure at T=500°C and t=44000 s in figure 3 is taken as the base metal microstructure and is used as an initial condition for the welding simulations at different locations.

Microstructure evolution during welding cycle
As discussed in section 2 and shown in figure 2, we have selected six points across the weld, each with its own thermal history. The temperature in each RVE associated with those points is considered homogeneous. The initial microstructure (prior to welding) in each of the RVE represents statistically similar region of the material. Therefore, it is reasonable to select the same initial state for each RVE. The phase field equations for a stress-free system with periodic boundary conditions are solved at each of these points, starting from the initial microstructure and using the thermal history for each RVE.
Composition of aluminium X Al at the end of the welding cycle (Time=100 s) at all selected points is presented in figure 4. RVE (6), the closest to the weld line (the edge of the HAZ) shows a large fraction of martensite (dark shaded regions), with only a few spherical α domains, surrounded by β regions observed. Such spherical α domains in a β matrix are also observed at location (5) but the fraction of martensite is significantly decreased. A completely different microstructure is observed at location (4), where martensite is not present, instead we see an α/β microstructure with non-equilibrium compositions (regions marked by arrows). Those regions are partially dissolved α domains. Naturally, no martensite is present at locations (1)-(3). At locations (1) and (2), the microstructure is essentially the same as the initial one. In RVE (3), β regions have grown slightly compared to the base metal microstructure with some partially dissolved domains. Such graded microstructures are indeed observed in a large class of titanium alloy welds. For example, in a recent study on inertial friction welding of Ti-6Al-4V, a similar sequence of microstructures has been reported [40]. Far away from the weld interface, the morphology is same as the initial base metal microstructure. For HAZ regions close to the interface, the martensite forms where as in regions further away, partially dissolved alpha regions referred to as 'ghosts' are formed. Figure 4 also shows that a mixture of martensite, alpha and beta can be observed in HAZ regions close to the weld center. A similar coexistence of phases has also been reported in Ti-6Al-4V fusion welds [41].
The simplest way to understand the mechanism behind gradation of microstructure in HAZ is to analyse the evolution of microstructure in time. Figure 5 captures RVE (6) at selected time points. At this location, the peak temperature during the welding is much higher than the β transus temperature b T and the system is kept above T β a considerable amount of time. Although α is thermodynamically unstable above T β , at time point A, some amount of α domains are still observed due to the fast heating rate. By the time the temperature reaches its maximum value, all α domains disappear and a homogeneous β phase with the nominal composition = = X X 0.102, 0.036 . During the cooling phase, below b T , some α nuclei appear (point C). Since the cooling rate is relatively fast, the precipitates are small and equiaxed. As soon as the temperature reaches the martensite start temperature T ms , regions that meet the martensitic transformation condition f m =1 are assumed to transform to martensite (dark shaded regions in time point D). Below this temperature no further nucleation takes place and the microstructure does not evolve significantly (compare points D and E). These results are fully consistent with experimental observation of martensite close to the weld line in titanium alloys [8]. Furthermore, [8] reports that the locations in the HAZ that transformed completely to β during heating, eventually formed martensite on cooling. Note that, the evolution at location (5) (not shown here) is quite similar to that for location (6). A complete dissolution occurs during heating. However, since cooling rate is slower, number of nucleation events is higher and there is more time available for diffusion and a relatively smaller region meets the martensite formation criterion. As a result the amount of martensite is expected to be much smaller (see figure 4).
The maximum temperature for thermal history at location (4) is only slightly above T β , as shown in figure 6. Since the microstructure starts to evolve only close to T β (point A) and the sample temperature is above T β only for a short time, the complete transformation to β for this case does not happen. At the peak (point B), domains of α phase are well preserved, although some of these domains have started to dissolve. As soon as temperature falls below T β , α domains have a tendency to grow. However, due to the rapid cooling and slow kinetics at low temperatures, there is insufficient time for any significant change to happen. As a result, non-equilibrium microstructure with partially dissolved α domains will be retained (points C, D and E). Considering above, it is evident that for regions in the weld where the peak temperature only slightly exceeds the transformation temperature, partially dissolved α regions may be observed. As discussed earlier, a similar phenomenon is indeed observed in welding experiments on titanium alloys, and reported as 'ghost domains' [10,11]. At location (3) the peak temperature is about 875°C (see figure 2). Although T β is not reached, the temperature is high enough to cause some evolution of the microstructure during heating. As seen in figure 4, the β domains at this location have grown slightly compared to the initial microstructure (for example, see the white arrow in the image corresponding to location 3). At locations (1) and (2), the peak temperatures are much lower than T β . Due to the low temperatures, the kinetics is too slow to allow for significant motion of the domains within the short time interval associated with the heating/cooling cycle.
The phase fractions plotted as function of the distance from the weld line at the end of the welding cycle are shown in figure 7. Essentially all phases (α, β and martensite) coexist up to 7 mm from the weld line. The α phase volume fraction starts decreasing from its nominal value of the base metal at around 11 mm, and falls to zero at approximately 6 mm from weld line. Martensite forms in a region between the weld line to a distance of 7 mm from the weld line. Close to the weld line (4 mm) about 90 percent of the system transforms to martensite. As a result, β phase volume fraction shows a maximum at around 7 mm from weld line. Figure 7 shows how microstructure changes dramatically over a distance of a few millimetres.

Discussion
The interfaces of α/β microstructures in Ti-6Al-4V are barriers to dislocation motion and they contribute to the yield strength of the alloy. At the same time, it has been observed that crack propagates by the formation of interfacial microcracks and the joining up of these microcracks through the a plates [42]. The resistance to crack propagation through these alloys and hence fracture toughness are found to depend on the thickness of these a plates. In addition, fatigue crack initiation is affected by the crystallography of the / a b lamellar microstructure and the dislocation slip transmission through the interface [43]. Therefore, in the spirit of Integrated The dark domains correspond to martensite regions where the condition f m =1 is satisfied.
Computational Materials Engineering (ICME), virtual simulations to relate microstructure and mechanical properties would require the morphology of the microstructure of the welds. There have been some studies of α/β microstructures evolution in titanium alloy welds [44][45][46][47] which successfully describe the volume fraction of different phases in the welds and mechanical properties can be predicted using pre-postulated material models. However, mesoscale microstructure and the complex morphologies that arise due to the transformation are not resolved. For example, the geometry of the a plates as well as compositional inhomogeneities that may occur due to welding process are not accounted for. The phase field approach simulates non-equilibrium phenomena that occur inside the weld and resolves morphological changes. with relative ease since explicit interface tracking is not required. The resulting microstructure morphology can serve as microstructure input for subsequent models to study microstructure-property relationship. The combined FEM/phase-field framework implemented in this paper can be applied to a large class of manufacturing processes as well as other alloy systems. However, there are some limitations that must be pointed out and should be addressed in future work. Since the size of RVE is only a few microns, the model does not take into account grain boundary nucleation, which can have consequences for the simulated morphologies. Here, we should also remark that although we have used available published data for Ti-6Al-4V as an input to our calculations, some assumptions had to be applied. For instance, there are limited data available for important kinetic parameters, atomic mobilities are not available for all the phases and there are no experimental measurements of the interfacial mobilities. Nevertheless, our results explain many aspects of the weld microstructure evolution which are in agreement with experimental observations.

Summary
We have developed and applied a combined FEM/phase field framework to simulate the evolution of α/β microstructure in heat affected zones of titanium alloy welds. The framework combines finite element heat transfer calculations with a phase-field microstructure model. The finite element calculations provide the thermal history at different locations in the simulated weld at the mm scale. The phase-field model is designed to simulate the microstructure evolution at selected representative points taking obtained thermal history as input. Regions where martensite can form are inferred from the local temperature and compositions during the simulations. The model is parameterized from available literature data for Ti-6Al-4V alloy. The simulations explain the graded microstructures that are commonly observed in HAZ of titanium alloy welds. It is shown that at locations closest to the weld line, due to rapid heating and cooling associated with welding process, mixed microstructures of a, b and martensite form. Further away from the weld line, non-equilibrium microstructures with partially dissolved a regions are observed. The model points to conditions which lead to such non-equilibrium configurations, which can be useful to design process conditions to achieve desired microstructures. The computational methodology presented in this paper may be used to study phase transformations and microstructures in other materials and manufacturing processes.