Microstructural modelling of above -transus heat treatment of additively manufactured Ti-6Al-4V using cellular automata

A heat treatment is an essential part of the metal additive manufacturing process chain. If an additively manufactured part, made of Ti-6Al-4 V, is heated above its β transus temperature, the columnar prior-β grains will become equiaxed β grains. This work quantitatively models this transition and the subsequent cooling down to room temperature by using the well-established cellular automata (CA) technique. Using this microstructural model allows visualisation of the local variation in the microstructure. The final microstructure consists of both the equilibrium phase α and β, organised in laths. This paper shows that the developed CA is capable of modelling the microstructural evolution during the entire above-β transus heat treatment. In order to get an accurate simulation of the microstructural change during such a heat treatment, the nucleation and grain growth functions are dependent on temperature. Since there exists a thermal gradient throughout the simulated cube, the local values of these functions will vary, leading to spatial differences in the nucleation frequency and growth velocity of new β grains. The model is verified by comparing the transformed volume fraction with a typical Johnson-Mehl-Avrami-Kolmogorov (JMAK) equation for isothermal grain growth. However, the JMAK equation insufficiently describes the grain growth during the initial stage of the heat treatment, namely while heating up to above the β transus temperature. Finally, the simulations of the second half of the heat treatment show that there are underexplored mechanisms during the growth of α laths when cooling down to room temperature. The simulations show, that it is not a requirement to nucleate α in the centre of the former β grains to form basketweave α. Moreover, the basketweave morphology in the simulated microstructures is a result of the difference between the viewing plane of the microstructure and the plane in which the laths grow, with a pure Widmanstätten morphology only appearing when the planes are parallel.


Introduction
Additive manufacturing (AM) allows the production of near-net shape parts, with limited production overhead [1][2][3]. Laser-based powder bed fusion (LPBF) is one of the AM techniques. During LPBF, the part is built up layer by layer, from a digital file. The LPBF machine distributes a thin layer of powder, and melts a cross section of the desired part into the powder. A new layer is placed on top of the previous one, and the process repeats until the part is completed. Due to the fast cooling of the melt pool, and the layer-based nature of the process, the microstructure of a part produced by LPBF is typically columnar [4,5]. Specifically for Ti-6Al-4 V, the microstructure consists of α' laths organised in columnar prior β grains [4]. This microstructure is referred to as columnar is the remainder of this work. Heat treating at a sufficiently high temperature will transform this columnar microstructure to an equiaxed one [4,6].
One of the materials that has gained much traction for use in LPBF is Ti-6Al-4 V. Because of its low density, and two-phase nature, Ti-6Al-4 V has good specific mechanical properties. The formation of oxides on the surface leads to good corrosion resistance [7,8]. In general, heating a Ti-6Al-4 V part above the β transus temperature transforms it into β. For part made by metal AM, the initial microstructure is made up of α'. The equilibrium phases below the β transus temperature are a mixture of α and β phases. The transformation to equiaxed β, reportedly consists of a process of nucleation and growth [9]. There are two cooling paths of interest in this work: cooling down quickly leads to martensitic α', while cooling down slowly leads to the aforementioned two-phase microstructure. Because of the fast cooling of the laser-based process, the initial microstructure for a heat treatment consists of α' laths in columnar former β grains. Experimental investigations of the microstructural evolution during an above β-transus heat treatment can be found in the work by Vrancken et al. [4], where the effect of the https://doi.org/10.1016/j.mtcomm.2020.101031 Received 27 December 2019; Received in revised form 25 February 2020; Accepted 26 February 2020 different process parameters on the microstructure are analysed. Pilchak et al. [9] also performed experimental investigations, where it was found that there exists a relation between β-grain size, heat treatment time and temperature. Galarraga et al. [10] found correlations between the cooling rate and the lath thickness of the final α-β microstructure. A review of the different attempts to relate the different microstructures and mechanical properties to the heat treatment can be found in the work by Shipley et al. [11]. Although stress relief heat treatments occur at a temperature significantly below the β transus temperature, an increase in temperature might accelerate the stress relaxation process [12]. During a stress relaxation, heating above the β transus could occur accidentally, and therefore a numerical analysis of this phenomenon has application for stress relaxation heat treatments. Additionally, heating above the β transus temperature has been proposed as a technique to remove unwanted microstructural features [13].
These experimental investigations show the importance of a good control of the microstructure through the heat treatment. Currently it is not easy to analyse the microstructure during the heat treatment at high temperature [14]. Additionally, to optimise the heat treatment experimentally requires a relatively large number of expensive trials. An alternative for experiments is a numerical investigation, where the evolution of the microstructure is modelled during the heat treatment.
Modelling of heat treatment of Ti-6Al-4 V is currently limited to averaged simulations, where the different fractions of different phases are compared using Johnson-Mehl-Avrami-Kolmogorov (JMAK) kinetics. In the work by Charles Murgau et al. [15] an above β-transus temperature heat treatement was modelled this way. Similarly, Da Costa Teixeira et al. [16] and Appolaire et al. [17] used JMAK kinetics to model a two stage heat treatment: heating up above the β-transus temperature and cooling down to room temperature. Additionally, they differentiated between three different α morphologies: grain boundary α (α GB ), Widmanstätten α (α W ) and basketweave α (α BW ). First, the α GB grew until it was present along the entire grain boundary. Next, the α W grew in from the α GB in a needle-like morphology. At the same time, α BW nucleated and grew in the centre of the β grain. Their model was capable of simulating the phase fractions of these morphologies, but unable to provide information on the exact position or distribution of different morphologies inside of the grains.
There are different techniques to perform a microstructural simulation, but the currently most popular approaches are the Monte-Carlo-Potts model, a Cellular Automata (CA) method or Phase Field models. An overview of the different microstructural modelling techniques can be found in the review by Markl and Körner [18]. The present work uses CA, because it has a clear link with physical phenomena, is computationally relatively cheap (compared to phase field simulations), and can be massively parallelised for time efficient computations.
CA was first successfully applied for microstructural modelling in the work by Rappaz and Gandin [19]. Their work modelled the dendritic solidification of an Al-7 wt% Si alloy. They showed that the method is capable of capturing the different types of grain morphologies that are present in a casting: columnar near the cast wall, and equiaxed in the centre. They used a constant temperature, consistent with thin plates, and coupled the CA method with Kurz-Trivedi-Giovanola (KTG) kinetics [19]. The KTG kinetics are particularly suited for simulation of rapid solidification processes, and assume a marginally stable dendritic growth front [20]. The CA method developed by Rappaz and Gandin was applied to LPBF process by Zinoviev et al. [2] in 2D and by Zinovieva et al. [21] in 3D. In their works, the temperature was not constant, but was calculated using an implicit and explicit finite difference scheme respectively. They showed that the microstructure obtained from numerical simulation by LPBF is columnar, and that the grain diameter increases from the bottom of the computational domain to the top. A coupling of a CA with the finite element method for obtaining the temperature field was performed by Guillemot et al. [22]. They developed a one-dimensional model for unidirectional solidification, and showed good agreement between the CA and a front-tracking model.
The CA method has also been applied to problems other than solidification, with success. Examples are its application to crack growth [23], dynamic recrystallization [24,25] or welding [26]. Most of these works focus on two-dimensional simulations. Since microstructure is intrinsically three-dimensional [21], it is important to take into account that features can grow into and out of the viewing plane. However, since it is very computationally expensive to perform full 3D simulations, alternatives have to be developed. An example of this is the work by Wei et al. [27], where the grain growth is modelled purely based on the thermal calculations.
A heat treatment above the β-transus temperature leads to a solidstate transformation. Therefore, KGT kinetics, as applied by Rappaz and Gandin [19] might not be suitable, since they predict the movement speed of solidifying dendrite tips. In this work, a parabolic growth law, derived by Pilchak et al. [9] from experimental observations of the above β heat treatment, is used instead.
Due to a lack of similar works for heat treatment of parts produced by LPBF, the present work aims to apply the CA method developed by Rappaz and Gandin to simulate this process. For validation, real micrographs are compared to the simulated microstructures. The proposed heat treatment heats the part up above the β-transus temperature. The two phase transformations that are modelled in this work are the transition from α' in columnar prior-β grains to equiaxed β grains and the transition from β to an α-β lath structure. In this work, we couple a CA with a finite difference (FD) based thermal solver.

Model description
As previously stated, the method used in this work is a CA, as described by Rappaz and Gandin [19], modified to a frontal CA [28]. This leads to the following transition rule: where SI k is the status index in cell k and N k is the neighbourhood of cell k. All the cells are initiated with a status index of zero, indicating that the cell is untransformed. If any cell in the neighbourhood of cell k is transformed into a growing state, with status index 1, cell k will get a status index of 1. If all the cells in the neighbourhood are transformed, then cell k will get a status index of 2, which means that the cell is completely transformed. This indicates the importance of the neighbourhood, since it determines the direction grains can grow. Most CA's use a regular grid, either hexagonal [29,30], or rectangular [2,3,19,21,25,31]. In this work, a rectangular grid is used. Rappaz and Gandin already pointed out that the choice of neighbourhood has a profound effect on the grain shape [19]. For rectangular grains, they compared the Neumann neighbourhood (cells sharing a side), and the Moore neighbourhood (cells sharing a point). Expanding the neighbourhood, for example to a third order neighbourhood [2], reduces the effect. This work uses an alternative approach where the neighbourhood is based on a circle with radius r [24,32]. The different possible neighbourhoods are demonstrated in Fig. 1.
A mathematical description of the neighbourhood that is solved in this simulation is the following: where x and z is the difference between the x-and y-coordinates of centres of cell k and cell l respectively, and r is the radius of the neighbourhood. The coordinate system is shown in Fig. 2a, where the x-y plane is parallel to the bottom of the simulated cube, and the z axis is perpendicular to it. The radius of the neighbourhood of all the cells with a status index of 1 increase in every CA time step by a fixed amount of half the cell size. The choice of this fixed value is motivated by analysing the circular neighbourhood in Fig. 1: the distance between the current cell centre and the nearest ones is equal to the cell size itself. The distance between the centre cell and the second nearest cells is the cell size, multiplied by the square root of two. Choosing the increase of the neighbourhood radius as half the cell size ensures that the nearest cells are captured first, and second nearest subsequently.
When any cell in the computational domain is captured, the diameter of all grains in the microstructure is calculated, based on the position of the grain-matrix interface. The diameter of the grain is linked to the elapsed time through the empirical relation, derived by Semiatin et al. [33,34]: where d is the grain diameter, k is a velocity constant, Q is the transformation energy, R is the universal gas constant, T is the temperature, and t is time. The subscript () 0 denotes the value of the parameter at the beginning of the calculation. The values of all the fixed parameters are shown in Table 2 and the temperature calculation is outlined in the next paragraph. The grain with the largest total time is fixed in diameter, and all other grains are allowed to grow using the circular neighbourhood, until they have reached the same total elapsed real time as the largest one. In order to calculate the evolution of real time accurately, a time dependent temperature is necessary. Therefore, the evolution of the temperature field is included in this simulation. The transient heat conduction equation is the basis to obtain the transient temperature field: where k is the thermal conductivity, the density and C p the specific heat capacity. To obtain the temperature field, this equation is solved using an alternate direction implicit (ADI) finite difference (FD) scheme [35]. Fig. 3a illustrates how the different parts of the model interact with one another. Fig. 3b shows the prescribed temperature profile, and illustrates when the temperature dependent CA is active.

Modelling methodology
In this work, the previously mentioned, well-established CA is coupled to the simple FD-based ADI solver. Although none of these modelling techniques is novel in their own right, their application to investigate the microstructure of a heat treatment above the β transus temperature quantitatively and qualitatively allows new insights in the evolution of the microstructure during such a heat treatment.
The simulation can be divided in two steps: heat treatment at high temperature, and cooling down. The first step consists of a temperature analysis, and a microstructural simulation. In the second step, the temperature is assumed uniform throughout the domain, and decreases linearly with time. Therefore, a temperature calculation is not required in this step. This section explains how each of these steps is set up. For each model, the input, computational domain and output are discussed. In this section, the numerical building blocks for the thermal calculation are referred to as elements, the CA domain is comprised of cells.

Step 1a: temperature calculation during heating up to above-β-transus temperature
In this work, the coupling between the temperature and the microstructural change is considered to be one way, meaning that there is no effect of the microstructure on the temperature. This also means that effects like release of latent heat from solid state transformations are not taken into account, because they are significantly lower than the heat input from the furnace.
The domain in which the temperature is calculated is a cube with an edge length of four centimetres. The ADI-FD scheme used here requires a ghost node just outside of the domain to take into account the boundary conditions. Therefore, the computational domain is divided in 17 elements in the x-and z-direction, and 17 elements in the y-direction (including the ghost elements). To improve the accuracy of the interpolation for the CA, which is explained in step 1b, the centre row of elements in x-and z-directions is subdivided in ten. Therefore, there are a total of one hundred fine thermal elements. The final size of each element is summarised in Table 1. The domain is shown in Fig. 2a. At the bottom (z = 0), the part is insulated, so there is a zero heat flux: At the other surfaces, the temperature is fixed to the furnace temperature.
This work choses to simulate the microstructure only in the centre of the cube, to minimise the effect of the boundary conditions. In addition, the CA, which is presented in the previous section, is limited to two dimensions. Therefore, a two-dimensional cross section of the part is taken in the centre of the cube. This cross section is indicated in red in Fig. 2. The output of the temperature calculation is a temperature profile as functions of position in the cube and time.
The CA mesh is comprised of 1000 by 1000 cells, located in the centre of the computational domain. Each cell has a side length of 1 μm, covering a total area of 1 by 1 mm 2 . This domain is smaller than the area covered by the fine thermal elements, to avoid possible edge effects from the coarse thermal mesh near the outside of the simulated cube.

3.2.
Step 1b: microstructural simulation during heating up to above-βtransus temperature The second half of the first step in the heat treatment procedure is the microstructural model for the high temperature step. The input for this model is an initial microstructure from experiments [36], which is converted to an idealised columnar microstructure. On the micrographs in Fig. 4, a change in colour marks the edges of the columnar grains. After conversion of the images into grey-scale, this results in different intensity values. Using this information, the MATLAB image analysis toolbox is used to divide the microstructure into superpixels. These superpixels are regions which have an apparent structure at a scale larger than the pixel size itself [37]. Inside of each superpixel, the image analysis code averages the intensity. Since the columnar grains all have similar colours, the superpixels from the same grain also have similar colours. The group of all the contiguous superpixels with a colour difference below a certain threshold form one columnar grain. The threshold used for the analysis shown in Fig. 4 is 7.8 percent, or an absolute value of 20 on an 8-bit grey-scale image.
To build up the idealised microstructure, line diameters are measured in the micrographs after identification of the positions of the grain boundaries. Line diameters are the most suitable metric to evaluate the size of the grains for a columnar microstructure, since the width of the columnar grains influences the number of nucleation sites more than the length of the columnar grains. Additionally, the difference in intensity in the initial micrographs incidentally exceeds the threshold and therefore a single columnar prior β grain is identified as multiple shorter grains. However, a line diameter is largely unaffected by this vertical segmentation of single columnar grains, since it primarily measures the distance between grain boundaries perpendicular to the columnar prior-β grain boundaries. To evaluate the line diameter, the program draws a number of horizontal lines across the micrograph, and counts the number of intersections between the columnar grains and this line. The average grain diameter is then equal to the division of the line length by the number of intersections, and the width of each individual grain is the difference between the positions of two consecutive intersections. This image analysis reveals that the feature size is in the order of 0.1 mm, and therefore a micrometre resolution is required. The size of each CA cell is 1 μm, and therefore the temperature, which is calculated on a grid size in the order of 1 mm, is interpolated between the different temperature nodes.
In each time step, the previously computed temperature is linked with the grain growth as indicated by Eq. (3). To model this microstructural evolution, nucleation is required, as indicated by Rappaz and Gandin [19]. This work also uses a Gaussian distribution for the nucleation rate, as indicated in their article. The nucleation function in this work is: in which N is equal to the nucleation density as a function of temperature, n avg is the average nucleation density in the domain, T is the standard deviation of the nuclei distribution, T µ is the average undercooling at which nucleation happens, and T is equal to the undercooling, namely T T , the difference between the current temperature and the β transus temperature, and represents the driving force. The average nucleation density is fixed at 50e-18, based on the micrographs after the heat treatment. This is an estimated value for the average nucleation density for the above-β heat treatment. Rescaling the grain density with the average CA cell. The nucleation chance is  Fig. 3a shows the flow chart of the temperature dependent CA, and the interactions between the CA, the thermal ADI, and the real grain growth. Fig. 3b shows the prescribed temperature profile on the exterior of the cube. The thermal ADI is active during the heat treatment above the β transus temperature. Since the temperature decreases slowly during the second half of the heat treatment, no thermal gradient is present in the microstructural simulation domain.  Table 2 Simulation parameters. The average undercooling is based on the works in [40][41][42], while taking into account the difference in process and process conditions. calculated, based on the work by Rappaz and Gandin [19], is the nucleation density, divided by the area of the CA cells, which is 1e-18 m 3 . This work will use the terminology from the original research [19], and refer to T as the undercooling, even though, in this case temperature required to initiate nucleation is higher than the equilibrium temperature (and therefore undercooling is technically overheating). Every time step, this nucleation probability is integrated between the previous undercooling, and the current undercooling at each available nucleation site. It is assumed that nucleation of the β grains occurs preferentially on the columnar prior-β grains [4]. Other sites, such as the α' lath boundaries could be considered, but due to the small CA cell size require to resolve the individual α' laths, this nucleation site is omitted from the simulations. If the difference between a uniformly distributed, but randomly generated number and the nucleation probability is smaller than zero, nucleation occurs in the cell. Therefore, the status of the cell will change from 0 to 1. When the status changes to 1 after nucleation, the cell is also assigned a random label, which allows to distinguish the different cells, and allows identification of the grain boundary, when two cells, which share a side, have a different label. At the same time, the grain diameter of this nucleus is initiated at zero.
In each time step, the size of the neighbourhood is calculated, and all the CA cells inside of the neighbourhood are captured by the growing grain. The status variable of these cells is updated until all the cells have a status of 2.

Step 2: cooling down to room temperature
The second step is cooling down to room temperature. According to Appolaire et al. [17] and Da Costa Teixeira et al. [16], the transformation from β to α + β, can be considered as two processes, namely growth of α on the grain boundaries of the β grains, and the growth of the laths from the grain boundaries into the interior of the β grains. The grain boundary α growth is a simple athermal CA, with five randomly distributed nuclei per cell.
The lath growth module consists of two major parts: surface reconstruction and growth. In this work, the so-called efficient least-square VOF interface reconstruction algorithm [38] is used to reconstruct the grain boundary of the equiaxed β grains. The boundaries are approximated as linear segments. As pointed out by Comminal et al. [38], this method is second-order accurate, and therefore should reconstruct the line segments accurately, providing a good estimate for the local normal to the β grain boundary.
Due to the orientation relationship, that exists between the β matrix and the new α laths, new α laths will only grow in a limited number of directions in each β grain. More specifically, twelve variations for the α lath growth are predominant [39]. In this work, the α laths originate at the grain boundaries, and grow inwards. As a result, half of these variations are directed outwards, and only six variations for the direction of the α laths need to be considered. These six variations are referred to as the preferential growth directions. These growth directions are rotated by a random angle, to represent the different crystal orientations in the different β grains. This is illustrated in Fig. 5a, which shows the local, randomly rotated coordinate system in the β grain. On the different grain boundaries, all six possible preferential growth directions are indicated.
Each α lath is nucleated by changing the status index of a cell in the grain boundary from 0 to 1. The position of these changes in status matrix is governed similarly to Eq. (6). However, because the temperature is uniform throughout the domain, and decreases slowly overall, combined with the speed at which α laths grow, the temperature dependence of the nucleation function is negligible. As a result, nucleation of the α laths is almost instantaneous, and the nucleation density is given by: where N is once more the nuclei density, n the number of nucleation sites, and lath the average lath thickness, which is assumed to be six micrometre. Multiplying the nuclei density by the cell size results in the nucleation probability [19]. Identically to the nucleation when heating up, for every eligible nucleation site, a random uniformly distributed number between zero and one is generated, and compared to the   Fig. 4a shows one of the equiaxed β grains, together with its local coordinate system, and the six preferential growth directions for some of the grain boundaries. Fig. 4b illustrates how the difference of the growth direction of the laths, with respect to the observation plane are reflected by rescaling the real lath growth velocity to obtain the apparent growth direction. nucleation probability. If the random number is smaller than the calculated probability, a new nucleus is created in the β grain boundary. The CA used for growth of the α laths is the same as the one described in the one in Section 2 on the model description, with one modification. The radius, indicated in Eq. (2), is quadrupled in the direction parallel to one of the previously mentioned preferential growth directions. Since many laths will grow simultaneously, the ones growing perpendicular to the grain boundary will impinge the ones growing in other directions. Therefore, the quadrupled direction is the one closest to the normal to the grain boundary, as calculated by the surface reconstruction algorithm [38]. Each α lath inherits a colour similar to the one of the α gb, in order to be able to identify different α lath colonies.
Of course, even though the CA simulation is two-dimensional, laths can grow in three dimensions. This is accounted for by multiplying the preferential growth direction by the cosine of the angle between the lath and the observation plane. Since there is no physical relation between the arbitrary observation plane, and the crystal orientation of the grains, this angle is a random angle between zero and 90 degrees. Fig. 5b illustrates this rescaling schematically. It shows the β grain boundary, and a lath growing perpendicular to it. The apparent growth direction is indicated on the observation plane.
Each of these simulation steps requires different parameters. The values of these parameters are summarised in Table 2. For the thermal simulation, constant material properties are assumed, since the range where microstructural change could occur, namely between the β transus temperature and the heat treatment temperature, is limited. Fig. 6a shows the microstructure of the as-built part, which is used as input for the simulated initial microstructure, which can be seen in Fig. 6b. The input for this initial, simulated, microstructure is obtained by using the image analysis code, developed specifically for this work, and described in section 3.2. Two qualitative differences can be noted when comparing the experimental micrograph in Fig. 6a and the simulated initial microstructure. Firstly, there are no major defects present in the simulated microstructure. The simulation does not consider pores, which seem to be present throughout the experimental sample. Also, the simulated microstructure excludes nucleation, since it assumes perfect epitaxiallity, which has shown to represent some experiments [47]. However, the choice of the cutting plane for the micrograph can influence the appearance of the grains. Since the grains become wider near the top of the sample, some grains, which originally were located in front or behind the cutting plane, can grow into view. Conversely, grains, which started in the image plane, can grow out of it and seem to disappear. In the simulation, only the image plane is defined, and therefore, grains can only disappear due to impingement from neighbouring grains. The initial, simulated, microstructure is obtained by using the same CA algorithm outlined earlier in this work, but without considering nucleation, since grain growth is epitaxial, without the temperature dependence of Eq. (3) and with a Moor neighbourhood. Fig. 6c and d show distribution of the grain diameters. They illustrate, that although the specific appearance of the modelled microstructure is different, the distribution of the grains is not dissimilar. A more accurate model for obtaining this columnar microstructure can be found in the work by Zinovieva et al. [21], but since this work focusses on the heat treatment, the initial microstructure in Fig. 6b is used as input for the actual heat treatment.  Table 1. The microstructural domain, comprised of 1000 × 1000 CA cells, is completely enclosed by the region with fine elements. In Fig. 7a, the temperature is depicted at six different points in time, in the x-direction of the domain (and therefore along a line at z = 2 cm and y = 2 cm).

Fig. 6.
Initial microstructure for the simulation [36]. Fig. 6a shows the micrograph of the asbuilt sample, in Fig. 6b, the initial microstructure of the simulations is depicted. Figs. 6c and 6d show the distribution of the grain diameters of the experimental grains and the simulated initial microstructure respectively. Note that the minimum grain diameter in both cases is limited to approximately 46 μm. Fig. 7a shows that the temperature increases gradually, but stays approximately constant throughout the domain. Focussing on the blue line in Fig. 7a, one can see that there exists a thermal gradient from the centre of the microstructural domain to the side (in x-direction). This gradient is shown in Fig. 7b, when the temperature on this line in xdirection is displayed separately. Fig. 7c shows the temperature on a line in the z-direction, corresponding to x = 2 cm and y = 2 cm. This temperature profile also shows a clear gradient, from 1255.5 K near the bottom of the fine elements to 1258.7 K near the top.
The origin of these gradients is the following: The sides and top of the cube are exposed to the hot air in the furnace directly, and will thus heat up first. The bottom of the cube does not contact the furnace air, and will thus heat up more slowly. Since the thermal conductivity of Ti-6Al-4 V is relatively high, the temperature gradients inside of the cube will stay relatively small, although at t = 20 s, a clear gradient is present in the cube. Looking solely at the fine elements in the centre of the cube, Fig. 6c indicates a gradient of four kelvin between the top and the bottom. Since the β transus temperature is 1253 K, and an undercooling (which represents a deviation from the equilibrium temperature) of 3 K ± 0.1 K is required for nucleation, the probability of nucleation is initially considerably higher at the top of the domain compared to the bottom. As mentioned earlier, undercooling is used in literature [19] when describing microstructural modelling, since most investigations are concerned with solidification or cooling down from high temperature to room temperature. However, in this work, the driving force is higher than the equilibrium temperature. Therefore, when this work refers to undercooling, this ought to be understood as overheating: heating above the equilibrium temperature. When the topmost fine elements heat up, they exceed the range in which nucleation can be expected, while the nucleation probability is still significantly larger than zero near the cooler bottom of the domain. This leads to a staggering of the occurrence of nuclei. This temperature gradient remains observable until the outside of the domain reaches the heat treatment temperature, when the temperature levels out again, which can be seen clearly in Fig. 7d. This figure also shows that the temperature in the centre lags behind the outside of the cube. After approximately 65 s, the temperature in the centre of the cube has reached 99.9 % of the heat treatment temperature. Fig. 8a shows the microstructure obtained after a heat treatment at 1323 K. The former β grains are surrounded by a thin layer of α, which allows identification of these equiaxed grains. Inside of the grains, two different structures can be identified: Widmanstätten α, which grows in from the grain boundary, and basket weave α, which is typically located in the middle of the grains. The Widmanstätten structure grows in from the α gb , and is characterised by the parallel plates in colonies. Basketweave α looks like interlacing lathes or groups of laths, forming a weave-like pattern. Fig. 8b shows the microstructure after holding the cube at 1323 K. The grains have transformed from the columnar former β grains of Fig. 6b to equiaxed. Fig. 8c shows the microstructure after cooling down to room temperature, where the α phase on the grain boundary is coloured white, both indicating the positions of the former equiaxed β grain, and highlighting α gb . Fig. 8d shows the final microstructure after the heat treatment. The α gb phase is less clearly visible here, since the α laths inherit their colour from their parent grain boundary α phase.
In Fig. 9, the grain diameters from the β grains in Fig. 8a and b are depicted in histograms, to compare them quantitatively. For each of these histograms, 400 measurements of line diameters are performed. Fig. 9 visually shows the agreement between the measured line diameters. In Fig. 8, a lognormal distribution is fitted to the data, allowing a quantitative comparison of the grains. The average grain diameter in the experimental micrograph is 220.1 μm, with a standard deviation of Fig. 7. Temperature in the fine elements, indicated in Table 1. Fig. 7a shows the temperature profiles as function of time and position on a line along the x-direction in the middle of the domain (where z = 2 cm and y = 2 cm). Fig. 7b the temperature profile at t = 20 s along the same line. Fig. 7c show the temperature profile along a line in the z-direction (x = y = 2 cm) at t = 20 s. In Fig. 7d, the simulated temperature in the centre of the cube (x = y = z = 2 cm) is plotted together with the prescribed temperature on the outside of the cube, matching a preheated furnace. 115.6 μm, while the grain diameter for the simulated equiaxed β grains is 169.6 μm, and the standard deviation 99.1 μm.

Discussion
Because of the symmetrical thermal boundaries, the temperature profile, shown in Fig. 7, is symmetrical. It can also be observed that the temperature profile increases rapidly. Within 45 s, the minimum temperature inside of the computational domain has reached 95 % of the temperature inside of the furnace. Additionally, the temperature profile flattens out, demonstrating that there is indeed only a limited effect of the boundary conditions under isothermal heat treatment conditions. On average, nucleation requires an undercooling of 3 K, as can be seen in the earlier described nucleation function. This means that nucleation will first take place at the side of the part that is heated up first. However, since the temperature is linear until the furnace temperature is reached, nucleation will continue homogeneously until the maximum nucleation density is reached. Therefore, the distribution of the grain diameter is largely constant throughout the microstructure, and only dependent on the random nature of nucleation.
It has been observed that the transition temperature in AM parts is not necessarily the theoretical one, which is measured in the technical Fig. 8. Microstructure after the heat treatment has been applied. Fig. 8a shows the microstructure after the real heat treatment. In Fig. 8b, the simulated β grains can be seen. Fig. 8c contains the microstructure at room temperature after the entire heat treatment is over. In this image, the alpha on the grain boundaries is coloured white, to show the different former β grains more clearly. Fig. 8d is the final microstructure after the heat treatment. alloy [6]. Therefore, the T , which is used in this work, can vary and is typically lower than 1253 K. However, only the nucleation density depends on this β transus temperature. This nucleation temperature contains two other tuneable parameters (specifically the average undercooling and the standard deviation of the undercooling). In this work, these parameters are derived from literature [40][41][42], but little is known about the exact values for a heat treatment above the β transus temperature. In Fig. 9, the measured line diameters are shown, together with the fitted lognormal distribution function. It is commonly assumed that the distribution of grains after normal grain growth [48][49][50] or even the distributions of the laths [41] is lognormal. To support this assumption, the measured grain diameters are displayed in a lognormal probability plot (Fig. 10). It shows that most of the measured line diameters are approximately lognormally distributed, but that there are significant outliers at low and high grain diameters. The origin of these outliers could be possible defects, resolution issues or for the simulations edge effects of the computational domain. Alternatively, a different distribution can be proposed, but this is beyond the intention of this work. Fig. 9 also revealed that there is a slight difference between the line diameter in the real micrographs and the virtual, simulated ones. Assuming the sample mean and standard deviations represent the real ones (motivated by the large number of grains measured using the image analysis code), a new sample of ten line diameters is taken to investigate the difference in the mean of the distribution. , this results in a z-value of 1.4194, hence failing to reject the null hypothesis that the means are the same for both distributions with a standard significance level of 0.05. Therefore, it is not possible to disprove that both microstructures represent the same lognormal distribution.
To illustrate the importance of including an accurate thermal calculation and to verify the microstructural simulation, the fraction transformed during the first part of the heat treatment (heating up) is displayed as a function of time in Fig. 11a. Alongside with the nonisothermal heat treatment, which is simulated using the temperature dependent CA, and alternative simulation is performed, where the temperature is fixed throughout the domain at the heat treatment temperature (therefore there is no temperature gradient). The undercooling for the nucleation function remains at four Kelvin to ensure that nucleation takes place. Since the second simulation mimics an isothermal heat treatment, the shape of the transformed fraction versus time curve should follow the Johnson-Mehl-Avrami-Kolmogorov (JMAK) equation, which is used to model different aspects of transformations in Ti-6Al-4 V [15,51,52]. The JMAK equation itself is of the form [53]: where X is the transformed fraction, and K and n are constants. To find a value for these constants, the equation is rewritten as: = + X K n t ln( ln(1 )) ln ln (8) which is shown in Fig. 12. A linear least square fit provides the parameters for equation (8) using the model output. The slope and intercept from this line represent the estimate for the values of the parameter with the smallest possible squared error. The linear fit and the equivalent JMAK curves are represented by the dashed lines in both Figs. 11a and 12 . Fig. 12 shows clearly that the non-isothermal heat treatment does not give good agreement with a JMAK-type curve. Moreover, the parameters of the linear regression are -7.0 for the intercept and 1.15 for the exponent, but this only gives a coefficient of determination (R 2 ) of 0.49, indicating an extremely bad fit. On the other hand, the isothermal heat treatment, which uses the same model as the non-isothermal one, but fixes the temperature at 1323 K with a nucleation undercooling of 4 K has an R 2 -value of 99.3 for the regression parameters -15.5 and 1.83 for the intercept and slope of the regression curve respectively. The conclusions from this result are twofold. On the one hand, it serves as verification for the CA method as microstructural simulation tool for the proposed heat treatment: if the heat treatment is kept perfectly isothermal, the JMAK theory holds true, confirming that the included physics suffice to gain insight in the microstructural change.
On the other hand, this results indicates that for a real heat treatment, it is important to take into account the change of temperature during the process, not only in time, but also spatially. Analysing the non-isothermal curve in Fig. 8a reveals two distinct regions: a steep growth, which starts when nucleation can occur, and thus, when the highest temperature is approximately four degrees above the β-transus temperature, and a slower growth, which starts after the wave of nuclei has stopped, but before the temperature has levelled out at 1323 K.
This non-isothermal heat treatment can be compared to the curves obtained by Murgau et al. [15] and Salsi et al. [54], who both used the additivity rule to applied the JMAK kinetics to a non-isothermal heat treatment. They analyse the evolution of the phase fraction during continuous cooling. Heating up a cube in a furnace can be treated as a case of heating with a high rate: the lowest heating rate is located at the bottom in the centre of the cube, and has a magnitude of 15.8 K/s.
Since the investigated phases and phase transformation is different, it is impossible to compare the curves in the aforementioned studies and the present work quantitatively, but comparing the shapes of the curve reveals that the beginning of the non-isothermal heat treatment in Fig. 11a has the same shape as the transformation curves obtained by Salsi et al. [54]. A zoomed in version of this curve is provided in Fig. 11b.
The second half of the curve in Fig. 11a represents the evolution of the grain size when nucleation seizes to happen due to the high overheating. The grain diameter increases according to an inverse quadratic function, also illustrated by Eq. (3). This is in partial contradiction to the work by Pilchak et al. [9], which analysed the early stage transformation. When they analyse the transformation of their titanium alloy, they notice the transformed fraction follows JMAK kinetics perfectly. The difference is probably due to the difference in the size of the examined work piece. Since the analysed domain for Pilchak et al. has dimensions of 6.35 mm by 6.35 mm by 12.7 mm, the temperature gradient inside of the part will be a lot smaller, and the lowest heating rate will be higher than observed in this work. Unfortunately, only the heating profile is provided up to the β transus temperature, but no information is provided on the temperature evolution once the β transus temperature is surpassed. Due to the sample geometry, it is safe to assume that the above transus heat treatment was isothermal, leading to the observed JMAK curve. Fig. 8d reveals that, contrary to what is assumed in the work by Da Costa Teixeira et al. [16], that it is not necessary to include nucleation of α laths to obtain a basket-weave structure. As a consequence, this simulation reveals that basket-weave α is just the consequence of Widmanstätten α growing into and out of the observation plane. Comparing Fig. 8a and d, it can be seen that the interweaving of laths in the simulations mimics the real laths closely, indicating that this out-ofplane growth is a feasible explanation for the transition from α-Widmanstätten to α-basket-weave. This is an extension of the study by Appolaire et al. [17], who only distinguished two different allotropiomorphs of α after the heat treatment above the β transformation temperature, namely grain boundary α and Widmanstätten plates. Determining the correct sequence of transformations requires either insitu experiments, or a full three-dimensional simulation. However, extending the CA method to three dimensions is difficult due to the associated computational cost. Different microstructural modelling studies for 3D CA already exist [21,24,28,55,56], but their number of cells is typically limited to one million cells. To resolve to lath growth during cooling down, the number of cells in the current study already exceeds this -it is namely 2.25 million. Extending the current domain while keeping the required resolution would lead to more than three billion cells, which, even with an efficient implementation of the CA method will require days to compute, rather than a few hours [28].
One of the origins of this difference between the performed simulations and the model used by Da Costa Teixeira et al. could be due to the difference in the used titanium alloy, or the difference in the heat treatment. Bein et al. [57] experimentally investigated the transformation kinetics of different titanium alloys, and different isothermal heat treatments. They show that, when heating at a high temperature, α laths only nucleate on the β-grain boundaries, while when heat treating at lower temperatures, more specifically 873 K, nucleation of the α laths also occurs homogeneously throughout the former β grains. The experiments by Salib et al. [58] confirms this influence of the heat treatment by verifying the transformation of the former β-grains into Widmanstätten α, growing in from the grain boundary. The effect of the alloy compositions on the microstructure is illustrated by in-situ experiments from the study by Aeby-Gautier et al. [14]. Higher equilibrium β phase fractions require higher concentrations of stabilising elements, such as vanadium, molybdenum, iron or chromium. When the α Widmanstätten laths grow, these elements are pushed out of the α phase. For near-β alloys, the concentration reaches a critical nucleation concentration more rapidly, resulting in nucleation in the centre of the β grains.

Conclusions
In this paper, the microstructural evolution of Ti-6Al-4 V is modelled above the β transus temperature. A cellular automata model is used due to its relatively low computational cost and therefore ability to simulate the change in microstructure at a high resolution. Both the transformation from martensitic α' to β during heating up, and the reverse transformation to an equilibrium α-β lath structure upon slow cooling is modelled. Image analysis software is developed for this work to capture the initial grain morphology from experimental micrographs, allowing a direct link between the actual primary LPBF process and the simulation.
The primary findings of this work are: • A temperature dependent version of the CA method is devised, which includes nucleation and grain growth for the α' to β transition. Comparing the transformed fraction versus time reveals the importance of the temperature evolution for parts with sufficiently large sizes. Not all regions in the computational domain have the same temperature and therefore there exists a gradient in the distribution of the grains.
• The isothermal CA is capable of capturing the α-β lath growth in two dimensions according to the mechanisms suggested in literature: growth of α on the grain boundaries and subsequent growth of the laths towards the inside of the former β grains.
• The developed model is not only capable of simulating a heat treatment up to a temperature of 1323 K, and also the subsequent heat treatment, but the results indicate that the model predicts similar curves as expected from JMAK kinetics, when the temperature is kept constant throughout the domain.
• The model indicates that nucleation is not necessary for obtaining basket weave α. It is possible to obtain this microstructure purely from laths growing in from the grain boundaries, due to the varying Fig. 11. a -Fraction transformed versus time for both the model developed in this work (non-isothermal) and an isothermal heat treatment. b -The non-isothermal curve from Fig. 11a, focussed on the initial few time steps. growth speed in with respect to the viewing plane. If the viewing plane is perfectly parallel to the lath direction, the apparent speed is maximal, while a viewing plane perpendicular to the lath growth direction results in an apparent speed of zero.

Data availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.