A Novel Physics-Based and Data-Supported Microstructure Model for Part-Scale Simulation of Laser Powder Bed Fusion of Ti-6Al-4V

The elasto-plastic material behavior, material strength and failure modes of metals fabricated by additive manufacturing technologies are significantly determined by the underlying process-specific microstructure evolution. In this work a novel physics-based and data-supported phenomenological microstructure model for Ti-6Al-4V is proposed that is suitable for the part-scale simulation of selective laser melting processes. The model predicts spatially homogenized phase fractions of the most relevant microstructural species, namely the stable $\beta$-phase, the stable $\alpha_{\text{s}}$-phase as well as the metastable Martensite $\alpha_{\text{m}}$-phase, in a physically consistent manner. In particular, the modeled microstructure evolution, in form of diffusion-based and non-diffusional transformations, is a pure consequence of energy and mobility competitions among the different species, without the need for heuristic transformation criteria as often applied in existing models. The mathematically consistent formulation of the evolution equations in rate form renders the model suitable for the practically relevant scenario of temperature- or time-dependent diffusion coefficients, arbitrary temperature profiles, and multiple coexisting phases. Due to its physically motivated foundation, the proposed model requires only a minimal number of free parameters, which are determined in an inverse identification process considering a broad experimental data basis in form of time-temperature transformation diagrams. Subsequently, the predictive ability of the model is demonstrated by means of continuous cooling transformation diagrams, showing that experimentally observed characteristics such as critical cooling rates emerge naturally from the proposed microstructure model, instead of being enforced as heuristic transformation criteria.


Introduction
Additive manufacturing has become an enabler for next-generation mechanical designs with applications ranging from complex geometries for patient-specific implants to custom lightweight structures for the aerospace industry. Especially metal selective laser melting (SLM) has gained broad interest due to its high quality and flexibility in the manufacturing process of load-bearing structures. Still, the reliable certification of such parts is an open research field not least because of a multitude of complex phenomena requiring the modeling of interactions between several physical domains on macro-, meso-and micro-scale [1].
A significant impact on elastoplastic material behavior, failure modes and material strength is imposed by the evolving microstructural composition during the SLM process [2][3][4][5]. Modeling the microstructure evolution in selective laser melting is thus an important aspect for more reliable and accurate process simulations and a crucial step towards certifiable computer-based analysis for SLM parts.
In [3,6] microstructure models are divided into the categories statistical [7][8][9][10][11][12][13][14][15][16][17], phenomenological [2-4, 6, 18-27] and phase-field [28][29][30][31][32] models. Statistical models are either data-driven and infer statements of coarse-grained trends from experiments or apply local stochastic transformation rules and neighborhood dependencies, which might be based on physical principles. This category includes in the context of this work also data-based surrogates and machine learning approaches as well as Monte-Carlo simulations and (stochastic) cellular automaton approaches. Without physical foundation, the reliability and predictive ability of purely data-driven approaches is rather limited, especially in the case of very scarce and expensive experimental data (e.g., dynamic microstructure characteristics in the high-temperature regime) or if the available data does not contain certain physical phenomena at all, which might result in high generalization errors. In case the simulation is based on stochastic rules (e.g., Monte-Carlo simulations) one encounters often challenges in terms of computationally demands as a reliable response statistic requires a large number of simulation runs. Furthermore, a consistent conservation of global and local physical properties for the individual simulation runs remains an open challenge. Stochastic properties might additionally be space and time dependent or functionally dependent on further physical properties. The inference of suitable and generalizable parameterizations of these stochastic properties is especially problematic in the case of limited experimental data.
On the other end of the spectrum of available models, phase-field approaches offer the greatest insight into microstuctural evolution and provide a detailed resolution of the underlying physics-based phenomena of crystal formation and dissolution, such as crystal boundaries and lamellae orientation. However, the resolution of length scales below the size of single crystals comes at a considerable computational cost, which hampers their application for part-scale simulations.
A preferable cost-benefit ratio can be found in the category of phenomenological modeling approaches. Here, microstructure evolution is described in a spatially homogenized (macroscale) continuum sense by physically motivated, phenomenological phase fraction evolution laws that can be solved at negligible extra cost as compared to standard thermal (or thermo-mechanical) process simulations. In this work we will propose a novel physics-based and data-supported phenomenological microstructure model for Ti-6Al-4V that is suitable for the part-scale simulation of selective laser melting processes. We present several original contributions, compared to existing approaches of this type. Compared to existing approaches of this type several original contributions, both in terms of physical and mathematical consistency but also in terms of the underlying data basis, can be identified.
From a physical point of view, the phase fraction evolution equations proposed in this work are solely motivated by energy considerations, i.e. deviations from thermodynamic equilibrium configurations are considered as driving forces for diffusion-based and non-diffusional transformations. Thus, the evolution of the most relevant microstructural phases, namely the β-phase, the stable α s -phase as well as the metastable martensitic α m -phase, is purely driven by an energetic competition and the temperature-dependent mobility of these different specifies. This is in strong contrast to existing approaches, where e.g. the formation of meta-stable phases is triggered by heuristic rules for critical cooling rates, which are taken from experimental observations and explicitly prescribed in the model to match the former. In the present approach, however, there is no need to prescribe such critical cooling rates as criterion for phase formation.
Instead, phase formation is a pure consequence of the underlying energy and mobility competition. Critical cooling rates can be predicted as a result of the modeling approach, and show very good agreement to experimental observations. From a mathematical point of view, the diffusion-based transformations are described in a consistent manner by evolution equations in differential form, i.e. ordinary differential equations that are numerically integrated in time, which renders the model suitable for the practically relevant case of solid state transformations involving temperature or time-dependent diffusion coefficients, arbitrary temperature profiles, and multiple coexisting phases. Again, this is in contrast to existing approaches modeling the phase evolution with Johnson-Mehl-Avrami-Kolmogorov (JMAK) equations. In fact, JMAK equations can be identified as analytic solutions for differential equations of the aforementioned type, which are, however, not valid anymore in the considered case of non-constant (temperature-dependent) parameters.
From a data science point of view, unknown parameters in existing modeling approaches are typically calibrated on the basis of single experiment data. As a consequence, this single experiment can then be represented with very good agreement while an extrapolation of the calibration data, i.e. a truly predictive ability, is only possible within very narrow bounds. In the present approach, a broad basis of experimental data in form of time-temperature transformation (TTT) diagrams is considered for inverse parameter identification of the (small number of) unknown model parameters. Moreover, the predictive ability of the identified microstructure model is verified on an independent data set in form of continuous-cooling transformation (CCT) diagrams, showing very good agreement in the characteristics (e.g. critical cooling rates) of numerically predicted and experimentally measured data sets. Since experimental CCT data is very limited (to only a few discrete cooling curves), the prediction of these diagrams by numerical simulation is not only relevant for model verification. In fact, the proposed microstructure model allowed for the first time to predict CCT data of Ti-6Al-4V for such a broad and highly resolved range of cooling rates, thus providing an important data basis for other researchers in this field. Eventually, the proposed model is exploited to predict the microstructure evolution for a realistic SLM application scenario (employing a state-of-the-art macroscale SLM model) and for the cooling/quenching process of a Ti-6Al-4V cube with practically relevant size (side length 10 cm).
The structure of the paper is as follows: Section 2 briefly presents the relevant basics of Ti-6Al-4V crystallography, the basic assumptions and derivation of the proposed microstructure evolution laws and finally the temporal discretization and implementation of the model in form of a specific numerical algorithm. Section 3 depicts the data-supported inverse parameter identification on the basis of TTT diagrams and model verification in form of CCT diagrams. In Section 4, first the basics of a thermo-mechanical finite element model employed for the subsequent part-scale simulations are presented. Then, in Sections 4 and 5 applicability of the proposed microstructure model to part scale simulations is demonstrated by means of two practically relevant examples, a realistic SLM application scenario as well as the cooling/quenching process of a Ti-6Al-4V cube.
2 Derivation of a novel microstructure model for Ti-6Al-4V In the following, we will derive a model for the microstructure evolution in Ti-6Al-4V in terms of volume-averaged phase fractions. First, we introduce fundamental concepts and outline our basic assumptions in Section 2.1. Afterwards, equilibrium and pseudo-equilibrium compositions of the microstructure states are presented in Section 2.2 as a basis for the subsequent concepts for arbitrary microstructure changes in Section 2.3. The model will be presented in a continuous and discretized formulation. The latter is then used in the numerical demonstrations.

Crystallographic fundamentals and basic assumptions
The aspects of microstructure evolution in Ti-6Al-4V alloys as considered in this work are assumed to be determined by the current microstructural state, the current temperature T as well as the temperature rateṪ , i.e. its temporal derivative. An overview of characteristic temperatures are given in Table 1. When cooling down the alloy from the molten state, solidification takes place between liquidus temperature T liq and solidus temperature T sol . The co-existent liquid and solid phase fractions in this temperature interval shall be denoted as X liq and X sol = 1 − X liq . Below the solidus temperature T sol the microstructure of Ti-6Al-4V is characterized by body-centered-cubic (bcc) β-crystals and hexagonal-closed-packed (hcp) α-crystals. First, β-crystals will grow in direction of the maximum temperature gradient for T sol < T ≤ T liq [4]. Depending on the prevalent cooling conditions, the α-phase can be further subdivided into stable α s and metastable α m -phases (Martensite). For sufficiently slow cooling rates |Ṫ | |Ṫ αm,min | the microstructure evolution can follow the thermodynamic equilibrium, i.e. stable α s nucleates into prior β-grains. This diffusion-driven transformation between alpha-transus start temperature T αs,sta and alpha-transus end temperature T αs,end results in a temperature-dependent equilibrium composition X eq α (T ) (see Figure 1, left) characterized by 90% α s -and 10% β-phase, i.e. phase fractions X αs = 0.9 and X β = 0.1, for temperatures below T αs,end [4,20,28,35].
Under faster cooling conditions, the formation rates of the stable α s -phase, which are thermally activated and limited by the diffusion-driven nature of this transformation process, cannot follow the equilibrium composition X eq α (T ) anymore such that β-phase fractions higher than 10% remain below T αs,end . At temperatures below the Martensitestart-temperature T αm,sta , the metastable Martensite-phase α m becomes energetically more favorable than the excess (transformation-suppressed) β-phase fraction beyond 10%. Under such conditions the β-crystals collapse almost instantaneously into metastable Martensite following a temperature-dependent (pseudo-) equilibrium composition X eq αm (T ) [4,6,22,33,36]. The critical cooling rateṪ αm,min is defined as the cooling rate above which the formation of stable α s is completely suppressed (up to the precision of measurements). The resulting microstructure consists exclusively of βand α m -phase fractions. The temperature-dependent Martensite phase fraction is denoted as X eq αm,0 (T ) for this extreme case. It is typically reported that the Martensite-end-temperature T αm,end , i.e. the temperature when the Martensite formation is finished, is reached at room temperature T ∞ going along with a maximal Martensite phase fraction of X αm = 0.9 in this extreme case (see Figure 1, right). Figure 1: Left: Equilibrium compositions X αs , X β and X liq resulting from slow cooling rates |Ṫ | |Ṫ αm,min |, such that X αm = 0. The experimental data is based on Katzarov [28], Kelly [4], Malinov [20] and Pederson [35]. Right: Metastable (pseudo-)equilibrium composition due to Martensite transformation for |Ṫ | > |Ṫ αm,min |, such that X αs = 0.
In addition to these most essential microstructure features, which will be considered in the present work, some additional microstructural distinctions are often made in the literature [4,6,37]. For example, the stable alpha phase α s can be further classified into its morphologies grain-boundary-α gb and Widmanstätten-α w phase, with X αs = X αgb + X αw . During cooling, grain-boundary-α gb forms first between the individual β-crystals. Upon passing the intergranular nucleation temperature T ig , lamellae-shaped Widmanstätten-α w grow into the prior β-crystal starting from the grain boundary. The appearance of such α-morphologies can have different manifestations such as colony-and basketweaveα s or equiaxed-α s -grains [4].
The microstructure model proposed in this work will only consider the accumulated α s phase without further distinguishing α gb and α w morphologies. Given their similar mechanical properties [38,39], this approach seems to be justified, as the proposed microstructure model shall ultimately be employed to inform homogenized, macroscale constitutive laws for the part-scale simulation of metal powder bed fusion additive manufacturing (PBFAM) processes. In a similar fashion, the so-called massive transformation is often observed in the range of cooling rates that are sufficiently high but still belowṪ αm,min , which leads to microstructural properties laying between those of pure α s and pure α m . For similar reasons as argued above, this microstructural species is not explicitly resolved by a separate phase variable in the present work. Instead, the effective mechanical properties of this intermediate phases are captured implicitly by the co-existence of α s and α m phase resulting from the present model at these cooling conditions. In addition and in summary, the following basic assumptions are made for the proposed modeling approach: 1. The microstructure is described in terms of (volume-averaged) phase fractions, i.e. no explicit resolution of grains and grain boundaries. 2. Only the most important phase species β, α s and α m are considered. 3. The Martensite-start-and Martensite-end-temperature are considered as constant, i.e. independent of the current microstructure configuration. 4. Currently, no information about (volume-averaged) grain sizes, morphologies and orientations is provided by the model. 5. The influence of the mechanical stress state, microstructural imperfections (e.g. dislocations) as well as further morphologies is not considered.
Partly, these assumptions are motivated by a lack of corresponding experimental data. In our ongoing research work, we intend to address several of these limitations. Remark (Cooling rates during quenching). Note that the cooling rates during quenching experiments are not temporally constant in general. Thus, the critical cooling rateṪ αm,min measured in experiments is usually the cooling rate at one defined point in time, typically defined at a high temperature value such that the measured cooling rate is (close to) the maximal cooling rate reached during the quenching process. In such a manner, a value ofṪ αm,min = 410 K/s has been reported in the literature for Ti-6Al-4V [33] but was interpreted in several contributions [6,22,40] as a fixed constraint for Martensite transformation.

Equilibrium and pseudo-equilibrium compositions
In this Section the temperature-dependent, thermodynamic equilibrium and pseudo-equilibrium compositions of the β, α s and α m are described in a quantitative manner. These phase fractions X i ∈ [0; 1] have to fulfill the following continuity constraints: For simplicity, the solidification process between liquidus temperature T liq and solidus temperature T sol is modeled via a linear temperature-dependence of the solid phase fraction X sol : Moreover, we follow the standard approach to model the temperature-dependent stable equilibrium phase fraction X eq α (T ), towards which the α s -phase tends in the extreme case of very slow cooling rates |Ṫ | |Ṫ αm,min |, on the basis of an exponential Koistinen-Marburger law ( [41]; see black solid line in Figure 1): While the alpha-transus start temperature T αs,sta = 1273K has been taken from the literature [4,22], the parameters T αs,end = 935K and k eq α = 0.0068 K −1 in (3) have been determined via least-square fitting based on different experimental measurements as illustrated in Figure 1 on the left side. It has to be noted that the equilibrium composition X eq α = f (T ) in form of a temperature-dependent function f (T ) as given in (3) could alternatively be derived as the stationary point (∂Π(X αs , T )/∂X αs )| (Xα s =X eq α )= 0 ⇔ X eq α − f (T )= 0 of a generalized thermodynamic potential Π(X αs , T ) containing contributions, e.g., from the Gibb's free energies of the individual phases β and α s , from phase/grain boundary interface energies or from (transformation-induced) strain energies [2]. Here, the dependence of the potential Π(X αs , T ) on X β has been omitted since X β = 1 − X αs for solid material under equilibrium conditions, i.e. in the absence of Martensite. In the present work, for simplicity, the expression for the equilibrium composition X eq α = f (T ) has directly been postulated and calibrated on experimental data instead of formulating the individual contributions of a potential, which would involve additional unknown parameters. Still, from a mathematical point of view it shall be noted that the expression X eq α = f (T ) in (3) is integrable, i.e., a corresponding potential can be found in general, resulting in beneficial properties not only of the physical model but also of the numerical formulation.
Next, we consider the second extreme case of very fast cooling rates |Ṫ | ≥ |Ṫ αm,min | at which the diffusion-driven formation of X αs is completely suppressed. For this case, we model the metastable Martensite pseudo equilibrium fraction X eq αm,0 (T ), emerging in the absence of α s -phase, based on an exponential law [4,6,22,41,42]: While the value T αm,sta = 848K has been taken from the literature [4,6,33], we choose k eq αm = 0.00415 K −1 such that (4) yields a maximal Martensite fraction of X eq αm (T ∞ ) = 0.9 at room temperature, which is in agreement to corresponding experimental observations [22] (see Figure 1 on the right).
Finally, we want to consider the most general case of cooling rates that are too fast to complete the diffusion-driven formation of the stable α s phase before reaching the Martensite start temperature T αm,sta but still below the critical rate |Ṫ αm,min |, i.e. a certain amount of stable α s phase has still been formed and consequently a Martensite phase fraction below 90% is expected at room temperature. For this case, we postulate an effective pseudo equilibrium phase fraction X eq αm (T ) for the α m phase that accounts for the reduced amount of transformable β-phase at presence of a given phase fraction X αs of the stable α s phase according to It can easily be verified that (5) fulfills the important relation X eq αm (T )+X αs < 0.9 for arbitrary values of the current temperature T and α s -phase fraction X αs . This means, for any given value X αs , an instantaneous Martensite formation according to X eq αm (T ) will never result in a total α-phase fraction X α = X αm + X αs that exceeds the corresponding equilibrium composition X eq α (which takes on a value of X eq α = 0.9 in the relevant temperature range below T αs,end ). In the extreme case that the maximal α s -phase fraction of 90% has already been formed before reaching the Martensite start temperature T αm,sta , Equation (5) ensures that no additional Martensite is created during the ongoing cooling process. Again, the pseudo equilibrium composition X eq αm =f (T ) in form of a temperature-dependent functionf (T ) as given in (5) could alternatively be derived as the stationary point (∂Π(X αm , X αs , T )/∂X αm )| (Xα m =X eq αm )= 0 ⇔ X eq αm −f (T )=0 of a generalized thermodynamic potential Π(X αm , X αs , T ), in which the current α s phase fraction X αs = X eq α can be considered as a fixed parameter. Thus, X eq αm =f (T ) according to (5) is not a global minimum of this generalized potential but rather a local minimum with respect to X αm under the constraint of a given α s phase fraction X αs = X eq α . This model seems to be justified given the considerably slower formation rate of the α s phase as compared to the (almost) instantaneous Martensite formation (see also the next section).
For a given temperature T αm,sta ≥ T ≥ T ∞ and α s -phase fraction X αs ≤ 0.9 during a cooling experiment, Equation (5) will in general yield a Martensite phase fraction such that X α = X αs + X αm ≤ 0.9, i.e. the sum of stable and martensitic alpha phase fraction might be smaller than the equilibrium phase fraction X eq α according to (3). In this case of co-existing α s -and α m -phase, it is assumed that X eq α in (3) as well as its complement X eq β = 1 − X eq α represent the (pseudo-) equilibrium compositions for the total α phase fraction X α = X αs + X αm and for the β phase fraction X β = 1 − X α . In other words, the driving force for diffusion-based α s -formation, as discussed in the next section, is assumed to result in the following long-term behavior: While at low temperatures, Martensite is energetically more favorable than the β-phase, which is the driving force for the instantaneous Martensite formation, it is assumed to be less favorable than the α s -phase. Therefore, there exists a driving force for a diffusion-based dissolution of Martensite into α s -phase, resulting in the following long-term behavior: lim t→∞ X αm =X eq αm = 0.
However, as discussed in the next section, the diffusion rates for this thermally activated process drop to (almost) zero at low temperatures such that Martensite is retained as meta-stable phase at room temperature. Thus, (7) can only be considered as theoretical limiting case in this low temperature region.

Evolution equations
Since the melting and solidification process is completely described by Equation (2), this section focuses on solidstate phase transformations for temperatures T < T sol below the solidus temperature (i.e. X sol = 1). To model the formation and dissolution of the α s -, α m -and β-phase, we propose evolution equations in rate form with the following contributions to the total rates, i.e. to the total time derivativesẊ αs ,Ẋ αm andẊ β , of the three phases: Here, e.g.,Ẋ β→αs represents the formation rate of α s out of β whileẊ αs→β represents the dissolution rate of α s to β.
The meaning of the individual contributions in (8), the underlying transformation mechanisms (e.g. instantaneous vs. diffusion-based) as well as the proposed evolution laws will be discussed in the following. Since the formation and dissolution of phases might follow different physical mechanisms in general, we have intentionally distinguished the (positive) ratesẊ x→y ≥ 0 andẊ y→x ≥ 0 of two arbitrary phases x and y instead of describing both processes via positive and negative values of one shared variableẊ y↔x . It is obvious that (8) satisfies the continuity Equation (1) for temperatures T < T sol below the solidus temperature (with X sol = 1 andẊ sol = 0), which reads in differential form: Thus, the phase fraction X β = 1 − X αs − X αm ∀ T < T sol can be directly calculated from (1) and only the evolution equations for the phases α s and α m will be considered in the numerical algorithm presented in Section 2.3.2.

Time-continuous evolution equations in rate form
In the following, the individual contributions to the transformation rates in (8) will be discussed. One of the main assumptions for the following considerations is that α m ↔ β transformations take place on much shorter time scales than α s ↔ β transformations [4,5,21], which allows to consider the former as instantaneous processes while the latter are modeled as (time-delayed) diffusion processes. In a first step, the diffusion-based formation of the stable α s -phase out of the β-phase is considered [34,37]. Modified logistic differential Equations [43] can be considered as suitable model and powerful mathematical tool to model diffusion processes of this type. Based on this methodology, we propose the following model for the diffusion-based transformation β → α s : Diffusion equations of this type typically consist of three factors: i) The factor (X β −X eq β ) represents the driving force of the diffusion process in terms of transformable β phase (see (6)) and has a decelerating effect on the transformation during the ongoing diffusion process. With the continuity relations X β = 1−X α and X eq β = 1−X eq α this term could alternatively be written as (X eq α − X αs ). ii) The factor with X αs leads to a transformation rate that increases with increasing amount of created α s -phase, i.e., it has an accelerating effect on the transformation rate during the ongoing diffusion process. Physically, this term can be interpreted as representation of the diffusion interface between α s -and β-phase, which increases with increasing size of the α s -nuclei (and thus with increasing α s -phase fraction). iii) The factor k αs (T ) represents the temperature-dependent diffusion rate of this thermally activated process. From a physical point of view, this term considers the temperature-dependent mobility of the diffusing species. When plotting the phase fraction X αs over time (at constant temperature), the factors i) and ii) together result in the characteristic S-shape of such diffusion-based processes as exemplary depicted in Figure 2 for four combinations of k and c. Depending on the type of diffusion process, different values for the exponent c αs can be derived from the underlying physical mechanisms resulting in more process-specific types of diffusion equations. In its most general form, which is applied in this work, the exponent c αs of the modified logistic differential equation is kept as a free parameter, which allows an optimal inverse identification based on experimental data even for complex diffusion processes [43].
When cooling down (Ṫ < 0) the material at temperatures below the Martensite start temperature (T < T αm,sta ) and the equilibrium composition of the stable α s -phase has not been reached yet (X αs < X eq α ), an instantaneous Martensite formation out of the (excessive) β-phase is assumed following the Martensite pseudo equilibrium composition X eq αm according to (5). Mathematically, this modeling assumption can be expressed by an inequality constraint based on the following Karush-Kuhn-Tucker (KKT) conditions: The constraint X αm −X eq αm ≥ 0 (first inequality in (11)) states that X αm cannot fall below the equilibrium composition X eq αm since Martensite will be formed instantaneously out of the β-phase. Here, the formation rateẊ β→αm (second inequality in (11)) plays the role of a Lagrange multiplier enforcing the constraint X αm −X eq αm = 0 as long as Martensite is formed. According to the complementary condition (third equation in (11)), this formation rate vanishes in case of excessive X αm -phase (i.e.Ẋ β→αm = 0 if X αm − X eq αm > 0). This scenario can occur e.g. during heating of the Martensite material (i.e.Ẋ eq αm < 0) since the diffusion-based Martensite-dissolution process (see below) cannot follow the decreasing equilibrium composition X eq αm in an instantaneous manner. Since the α m -phase is energetically less favorable than the α s -phase, there is a driving force for the transformation α m → α s . In contrast to the instantaneous formation of Martensite, the α m -dissolution to α s is modeled as (time-delayed) diffusion process [6] according tȯ whereX eq αm = 0 represents the long-term equilibrium state (t → ∞) of the Martensite phase (see Equation (7)). Equations (10) and (12) have the same structure and share the same exponent c αs as well as the same diffusion rate k αs (T ) since both describe the diffusion-based formation of α s -phase, i.e., both result in the same daughter phase. The underlying modeling assumption is that the transformation α m → α s can be split according to α m → β → α s , i.e., it is assumed that Martensite first dissolves instantaneously into the intermediate phase β, which afterwards transforms to the α s -phase in a diffusion-based β → α s process similar to (10). The model for the temperature-dependence of k αs (T ) as described below will result in diffusion rates that drop to (almost) zero at low temperatures such that Martensite is retained as meta-stable phase at room temperature, which is in accordance to the corresponding experimental data.
Eventually, also the case of lacking β-phase (i.e., X β < X eq β ) shall be considered: Due to the nature of the β-dissolution processesẊ β→αs according to (10) andẊ β→αm according to (11), which only yield contributions as long as X β > X eq β , this scenario can only arise fromẊ eq β > 0, i.e., ifṪ > 0 and T ∈ [T αs,end ; T αs,sta ]. Experimental data by Elmer et al. [34] strongly supports a diffusional behavior of β-phase build up, respectively α s -dissolution during heating of Ti-6Al-4V. We thus chose a diffusional nucleation model in Equation (13) for the resulting α s → β transformation: Here,X β = X β − 0.1 defines a corrected β-phase fraction. This ensures that the second factor in (13), which can be interpreted as a measure for the increasing diffusion interface during the formation process, starts at a value of zero when heating material with initial equilibrium composition X β = X eq β above T αs,end . By this means, the temporal evolution ofX β (i.e., of the additional β-material beyond 10%) begins with a horizontal tangent when exceeding T αs,end as also observable in corresponding experiments [34]. From a physical point of view this experimental observation -as well as the corresponding diffusion model in (13) -rather correspond to a phase nucleation and subsequent growth process of a new phase fractionX β (starting atX β = 0) than a continued growth process of a pre-existing phase fraction X β (starting at X β = 0.1). While existing experimental data in terms of temporal phase fraction evolutions is very limited, there is still significant experimental evidence that the α s → β dissolution should be considered as a (time-delayed) diffusion-based process [4,34] rather than an instantaneous transformation, which for simplicity has been assumed in some existing modeling approaches, like e.g., [6].
Let us again consider the case of lacking β-phase (i.e. X β < X eq β ), which can only occur in the temperature interval T ∈ [T αs,end ; T αs,sta ] as discussed in the paragraph above. If in this scenario, a certain amount of remaining Martensite material (i.e. X αm > 0) exists, it is assumed that the α m -phase fraction is decreased in an instantaneous α m → β transformation such that X β can follow the energetically favorable equilibrium composition X eq β . Mathematically, this modeling assumption can be expressed by an inequality constraint with the following KKT conditions: The constraint X β −X eq β ≥ 0 (first inequality in (14)) states that X β cannot fall below the equilibrium composition X eq β as long as a remainder of Martensite, instantaneously transformable into X β , is present. Again, the formation ratė X αm→β (second inequality in (14)) plays the role of a Lagrange multiplier enforcing the constraint X β −X eq β = 0 as long as the α m → β transformation takes place. According to the complementary condition (third equation in (14)), this formation rate vanishes in case of excessive β-phase (i.e.Ẋ αm→β = 0 if X β −X eq β > 0). Since in the considered scenario, the remaining contributions toẊ β vanish, i.e.,Ẋ αs→β = 0andẊ β→αs = 0 due to X β = X eq β as well asẊ β→αm = 0 due to T > T αm,sta , (8c) together with the differential form of the constraintẊ β =Ẋ eq β allows to explicitly determine the corresponding Lagrange multiplier toẊ αm→β =Ẋ eq β . Once all Martensite is dissolved, i.e., X αm = 0, the β-phase fraction cannot follow the corresponding equilibrium composition X eq β in an instantaneous manner anymore, but rather in a time-delayed manner based on the diffusion process (13).
To close the system of model equations proposed in this section, a specific expression for the temperature-dependence of the diffusion rates k αs (T ) and k β (T ) required in (10), (12) and (13) has to be made. The mobility of the diffusing species, which is represented by these diffusion rates, is typically assumed to increase with temperature. However, it is assumed that the diffusion rates do not increase in a boundless manner but rather show a saturation at a high temperature level. Moreover, for the considered class of thermally-activated processes, these diffusion rates are assumed to drop to zero at room temperature. The following type of logistic functions represent a mathematical tool for the described system behavior: The free parameters k 1 , k 2 , k 3 and c αs governing the β → α s -diffusion processes (10) and (12) will be inversely determined in Section 3 based on numerical simulations and experimental data for time-temperature-transformations (TTT). The same temperature-dependent characteristic as in (15) is also assumed for the α s → β-diffusion process (13).
Since the dissolution of α s -into β-phase is reported to take place at higher rates [4,34] as compared to the α s -formation out of β-phase, we allow for an increased diffusion rate k β (T ) of the form: Thus, only the two free parameters f and c β are required for the α s → β-diffusion. These two parameters will be inversely determined in Section 3 based on heating experiments taken from [34].
Remark (Martensite cooling rate). In our model the critical cooling rateṪ αm,min = −410 K/s is not prescribed as an explicit condition for Martensite formation as done in existing microstructure modeling approaches [6,40]. Instead, the process of Martensite formation is a pure consequence of physically motivated energy balances and driving forces for diffusion processes. In Section 3.3 it will be demonstrated that a value very close tȯ T αm,min = −410 K/s results from the present modeling approach in a very natural manner when identifying the critical rate for pure Martensite formation from continuous-cooling-transformation (CCT) diagrams created numerically by means of this model. Remark (Johnson-Mehl-Avrami-Kolmogorov (JMAK) equations). In other publications [4,6,19,23,27,34,40,44], Johnson-Mehl-Avrami-Kolmogorov (JMAK) equations are used to predict the temporal evolution of the considered phase fractions. It has to be noted that JMAK equations are nothing else than analytic solutions of differential equations (for diffusion processes) very similar to (10), which are, however, only valid in case of constant parameters k αs , X αs and X eq β . Since these parameters (due to their temperature-dependence) are not constant for the considered class of melting problems, only a direct solution of the differential equations via numerical integration, as performed in this work, can be considered as mathematically consistent. We furthermore want to note that the mathematical form of JMAK equations, which involve logarithmic and exponential expressions, is prone to numerical instabilities in practical scenarios, especially for the extremely high temperature rates that appear in SLM processes. In contrast, the proposed algorithm in the following Section 2.3.2 has a simple and robust mathematical character.

Temporal discretization and numerical algorithm
For the numerical solution of the microstructural evolution laws from the previous section, we assume that a temporally discretized temperature field (T n ,Ṫ n ) based on a time step size ∆t is available at each discrete time step n (e.g. provided by a thermal finite element model as presented in Section 4.1). In principle, any time integration scheme can be employed for temporal discretization of the phase fraction evolution equations from the last section. Specifically, in the subsequent numerical examples either an implicit Crank-Nicolson scheme or an explicit forward Euler scheme have been applied. For simplicity, the general algorithmic realization of the time-discrete microstructure evolution model is demonstrated on the basis of a forward Euler scheme.
For the following time integration procedure of (8) it has to be noted that only the ratesẊ αm→αs ,Ẋ β→αs ,Ẋ αs→β corresponding to diffusion processes will be integrated in time. Instead of integrating the ratesẊ β→αm andẊ αm→β corresponding to instantaneous Martensite formation and dissolution processes, the associated constraints in (11) and (14) will be considered directly by means of algebraic constraint equations. In a first step, assume that the temperature data T n+1 ,Ṫ n+1 of the current time step n + 1 as well as the microstructure data X n αs , X n αm ,Ẋ n αm→αs ,Ẋ n β→αs ,Ẋ n αs→β of the last time step n is known. With this data, the microstructure update is performed: αs = X n αs + ∆t(Ẋ n β→αs +Ẋ n αm→αs −Ẋ n αs→β ) and X n+1 αm = X n αm − ∆tẊ n αm→αs .
The time integration error in (17) might lead to a violation of the scope X n+1 αs , X n+1 αm ∈ [0; 0.9] of the phase fraction variables. In this case the relevant phase fraction variable is simply limited to its corresponding minimal or maximal value, respectively. Similarly, if X n+1 α = X n+1 αs + X n+1 αm exceeds the maximum value of 0.9 the individual contributions X n+1 αs and X n+1 αm are reduced such that the maximum value X n+1 α = 0.9 is met and the ratio X n+1 αs /X n+1 αm is preserved. Subsequently, the β-phase fraction is calculated from the continuity equation X n+1 α . Afterwards, the updated equilibrium phase fractions X eq,n+1 αs and X eq,n+1 αm are calculated according to (3)-(5) with X n+1 αs and T n+1 before the corresponding equilibrium composition X eq,n+1 αs of the β-phase is updated. Next, a potential instantaneous Martensite formation out of the β-phase according to (11) is considered as follows: Similarly, a potential instantaneous Martensite dissolution into β-phase according to (14) is considered as follows: Again, if necessary the increment in (19a) is limited such that the updated phase fraction X n+1 αm does not become negative. As a last step, the diffusion-based transformation ratesẊ n+1 β→αs ,Ẋ n+1 αm→αs andẊ n+1 αs→β according to (10), (12) and (13), all evaluated at time step n + 1, are calculated. With these results, the next time step n + 2 can be calculated starting again with (17).
Remark (Initial conditions for explicit time integration). While the diffusion process according to Equation (12) will start at a configuration with X αm = 0 and X αs = 0, Equation (10) needs to be evaluated for X αs = 0 to initiate the diffusion process. However, the evolution of Equation (10) based on an explicit time integration scheme will remain identical to zero for all times for a starting value of X αs = 0. Therefore, during the first cooling period the phase fraction has to be initialized at the first time step t n where X eq α > 0. In the following, the initialization procedure considered in this work is briefly presented. In a first step, (10) is reformulated using the relations X β = 1 − X α and X eq β = 1 − X eq α as well as the approximate assumptions k αs (T ) = const., X eq α = const. and X αm = 0 (i.e., X α = X αs ) for the initial state [43]: Based on the analytic solution in [43], evaluated after one time step ∆t, the initialization for X n αs at t n reads: We compared this approach with an implicit Crank-Nicolson time integration (either used for the entire simulation or only for initialization of the first time step), where the initial α s -phase fraction X αs does not need to be set explicitly and found no differences in the resulting diffusion dynamics according to (10).

Inverse parameter identification and validation of microstructure model
The four parameters θ diff,αs = [c αs , k 1 , k 2 , k 3 ] T according to Equations (10), (12) and (15) as well as the two parameters θ diff,β = [c β , f ] T according to Equations (13) and (16), are so far still unknown and need to be inversely identified via experimental data sets. For the inverse identification of θ diff,αs , we use so-called time-temperature-transformation (TTT) experiments, a well-known experimental characterization procedure for microstructural evolutions [4,6,44,45].
As TTT-experiments only capture the dynamics of cooling processes, we will identify the parameters θ diff,β governing the heating dynamics of the microstructure via data from heating experiments taken from [34].

Inverse identification of α s -formation dynamics via TTT-data
Time-temperature transformation (TTT) experiments [46] are one of the most important and established procedures for (crystallographic) material characterization. The goal of the TTT investigations is to understand the isothermal transformation dynamics of an alloy by plotting the percentage volume transformation of its crystal phases over time. Thereto, the material is first equilibriated at high temperatures such that only the high-temperature phase is present. Afterwards, the material is rapidly cooled down to a target temperature at which it is then held constant over time so that the isothermal phase transformation at this temperature can be recorded. Rapid cooling refers here to a cooling rate that is fast enough so that diffusion-based transformations during the cooling itself can be neglected and can subsequently be studied under isothermal conditions at the chosen target temperature. The procedure is repeated for successively reduced target temperatures. The emerging diagram of phase contour-lines over the T × log(t) space is commonly referred to as TTT-diagram.
In the present work, the simulation of TTT-curves for Ti-6Al-4V was conducted as follows: We initialized the microstructure state at T = 1400 K > T αs,end with pure β-phase, such that X β = 1.0. Afterwards, the microstructure was quickly cooled down withṪ = −500 K/s 1 to a target temperature T target ∈ [350, 1300] and the evolution of the microstructure at this target temperature was recorded over time. The range of target temperatures was discretized in steps of 10 K such that 95 individual target temperatures and hence microstructure simulations were considered. Figure  3 depicts the isolines of simulated X αs -(left) and X αm -phase fractions (right), after identifying the model parameters via the experimental data [4,6,44,47] shown in the left figure. Please note that the phase-fractions in Figure 3 were normalized with X eq α (T ), which takes on a value of 0.9 for temperatures below T αs,end , as this was also the case in the underlying experimental investigations. We highlighted three temperatures T = {400, 800, 1000} K to discuss the microstucture evolution at these points.  Figure 3: Simulation of the TTT-diagram for the α s -and α m -phases using the maximum likelihood point estimate for the uncertain kinetic parameters θ * diff,αs of the microstructure evolution, along with experimental data by Malinov [47] and Kelly [37]. Left: Contour-lines for X αs ; Right: Contour-lines for X αm . Contour lines are shown for the 1%, 5%, 45%, 55%, 95% and 99% normalized phase fractions. Three temperatures are marked in red and discussed in the analysis.
At T = 1000 K (between T αs,sta and T αs,end , i.e. above T αm,sta ) we get a pure β → α s transformation. When looking at higher target temperatures T > 1000 K it can be observed that the isoline with 1% phase fraction for X αs is shifted to longer times, which results from a decreasing value of X eq α (i.e. a decreasing driving force) that slows down the initial dynamics of X αs -formation at elevated temperatures, even-though k αs is already saturated at its maximal value in this temperature range (see Figure 4). At T = 800 K we are now below T αm,sta , so that the initial cool-down results (instantaneously) in a Martensite phase fraction according to X eq αm . With ongoing waiting time, the remaining β-phase transforms in a diffusion-driven manner into stable α s -phase. In the TTT-diagram, we can already notice that the isoline with 1% phase fraction for X αs is shifted to longer times as compared to the higher temperature level T = 1000 K, which, this time, is caused by the decreasing value of k αs for lower temperatures, as depicted in Figure 4 and modeled in Equation (15). Moreover, the right-hand side of Figure 3 shows that the Martensite phase fraction decreases again for waiting times t > 100s, which represents the diffusion-based dissolution of Martensite into α s -phase according to (12).
Finally, the transformation at T = 400 K initially results in almost the maximal possible amount of Martensite (Remember: The 100%-isoline in Figure 3 corresponds to a phase fraction of 0.9 for T < T αs,end ). As the diffusion rate k αs is almost zero for such low temperatures, the Martensite phase cannot be dissolved to stable α s in finite times and the Martensite phase remains present as metastable phase at low temperatures, which agrees well with experimental observations. In the TTT-diagram this effect shifts the isolines asymptotically to t → ∞ when approaching the room temperature. All in all, the shift to longer times due to a low value of k αs for low temperatures and the delay effect due to a decreasing (driving force) value of X eq α at high temperatures leads to the typical C-shape of TTT-phase-isolines. Inverse parameter identification was conducted for the diffusion parameters θ diff,αs by maximizing the data's likelihood for θ diff,αs . We assumed a (conditionally independent) static Gaussian noise for the measurements on the log(t)-scale. This assumption is equivalent to a log-normal distributed noise in the data along the time-scale. The maximumlikelihood point estimate can then be determined by solving the following least-square optimization problem (see Remark below for more details): In Equation (22) the term X αs ,TTT (log(t i ), T i , θ diff,αs ) describes the simulated TTT-phase fractions (normalized by X eq α (T )) at time t i , Temperature T i and for diffusional parameters θ diff,αs . The index i marks here the specific temperatures and times for which the corresponding observed experimental data X αs ,TTT,exp.,i was recorded.
We utilized a Levenberg-Marquardt optimization routine [48], which is implemented in our in-house software framework QUEENS [49] to iteratively solve Equation (22). The result of this optimization procedure is given by the parameter set θ * diff,αs = [c αs , k 1 , k 2 , k 3 ] T = [2.51, 0.294, 850.0, 0.0337] T . Additionally, Figure 4 visualizes the temperaturedependent diffusion rate k αs (T, θ * diff,αs ) according to Equation (15) that results from these parameters.  (22) is the result of the following assumptions for the inverse identification task of θ * diff,αs : The simulation model is expressed by a mapping y(θ, c), with θ being model parameters that should be inferred from data via inverse analysis and c being coordinates on which the model output is recorded. Further model inputs (e.g., parameters that we want to keep fixed) are omitted to avoid cluttered notation. We assume that a vector of n obs scalar experimental observations y obs,C , recorded at coordinates C = {c i }, was disturbed by conditionally independent and static Gaussian noise with variance σ 2 n on the log(t)-space. This assumption is equivalent to log-normal distributed noise on the t-space. Statistically, the former assumption is expressed by the so called likelihood function l(θ) = p(y obs,C |y(θ, C)), which is the probability density for the observed data y obs,C , given a specific choice of the parameterized simulation model y(θ, C), evaluated at the same coordinates C. Gaussian conditional independent noise implies now that (theoretically) repeated observations of experimental data y obs,ci at each location c i will result in n obs Gaussian probability distributions with variance σ 2 n for y obs,ci . Independence in particular means that a disturbance by noise at c i does not influence the disturbance by noise at any other point. Ideally, the simulation output should reflect the mean of the observation noise at each c i . The latter can be expressed by the product (due to the conditional independence) of n obs Gaussian distributions with their mean values being the simulation outputs y(θ, c i ) at each c i and with variance σ 2 n : l(θ) = p(y obs,C |y(θ, C)) = nobs i=1 N y obs,ci |y(θ, c i ), σ 2 n ∝ exp − i (yobs,c i −y(θ,ci)) 2

Remark (Assumptions for inverse identification). The least-squares optimization problem in Equation
Despite l(θ) being a true probability density function for y obs,C , the latter is not the case w.r.t. the model parameters θ, such that l(θ) is mostly referred to as the likelihood function. The maximum likelihood (ML) point estimate θ * refers then to a value of θ that maximizes the likelihood function l(θ), respectively the probability density of the observations y obs,C for the specific choice θ * . The maximum of l(θ) can be found by minimizing the sum of the square terms in the argument of the exponential function, which leads ultimately to Equation 22.

Inverse identification and validation of diffusional heating dynamics
So far, the inverse identification of θ diff,αs via TTT-data only accounts for the cooling dynamics of the microstructure model. In the following, we inversely identify the parameter set θ diff,β via experimental data from [34], which consists of three temperature and β-phase data-sets for three positional measurements for an electron beam welding process with Ti-6Al-4V. Note that the specific measurement positions x = {4.5, 5.0, 5.5} mm defined in the original work [34] are not relevant here since we only aim at correlating temperature and phase fraction data. Please note also that this identification step reuses θ * diff,αs , respectively the temperature characteristics of k αs as identified above and only scales the latter by the factor f (see Equation (16)). The maximum likelihood estimate θ * diff,β was again calculated by solving a least-squares optimization problem with the Levenberg-Marquardt optimizer: Equation (23) accounts for all three phase-temperature measurements at locations = {4.5, 5.0, 5.5} mm simultaneously leading to a robust point estimate θ * diff,β with low generalization error.
The inverse identification yields the parameter set θ * diff,β = [c β , f ] T = [11.0, 3.8] T . The comparison of experimental data from [34] and the simulation-based prediction in Figure 5 show very good agreement. We summarize that the Elm er: Tem perat ure profile at x = 4.5 mm Figure 5: Experimental measurements for β-phase evolution in a welding process of Ti-6Al-4V by Elmer et al. [34] along with the corresponding simulation results based on the proposed microstructure model and the identified parameter set θ diff,β = θ * diff,β .
dissolution of the α s -phase in the case of X αs > X eq α can be successfully modeled by a diffusional approach for α s -nucleation by scaling k αs with a factor of f * = 3.8 and selecting c * β = 11.0. We will use this approach in the subsequent numerical demonstrations.

Validation of calibrated microstructure model via CCT-data
To validate the calibrated microstructure model, we will now compare the predicted microstructure evolutions with further experimental data. Another common experimental approach for the characterization of microstructural phase evolutions are the so-called continuous-cooling transformation (CCT) experiments [33,50]. Here, the microstructural probe is again equilibriated at a temperature above T αs,sta = 1273 K, such that X β = 1.0. Afterwards, the probe is cooled down to room temperature T ∞ = 293.15 K at different cooling ratesṪ CCT . Note that the true cooling rateṪ (t) is time-dependent andṪ CCT is only a representative descriptor. Eventually, the evolving microstructure of these cooling procedures is recorded on the t × T -space. In contrast to TTT-experiments, the dynamics of the cooling process itself and the thereof resulting microstuctural phase transformations are now the main aspect of these experimental procedures. Hence, the particular temperature profile of the cooling process is now of great importance and has a large impact on the emerging phases. The true cooling ratesṪ (t) in practical CCT-experiments are changing over time, being highest in the beginning of the cooling procedure and tending to zero for infinite long times. Following the procedure in [33], the characteristic cooling rateṪ CCT is in the following defined as the cooling rate at 900 • C, respectively 1173.15 K Ṫ CCT :=Ṫ T =1173.15 K . A schematic CCT-diagram for Ti-6Al-4V based on characteristic experimental cooling rates taken from [33] as well as the prediction by the proposed microstructure model are shown in Figure 6. In [33] two distinct CCT-cooling ratesṪ CCT are reported which divide the evolving microstructure for continuous Ti-6Al-4V cooling into three characteristic regimes. For cooling faster thanṪ CCT = −410 K/s, only martensitic transformation was observed. For cooling betweenṪ CCT = −410 K/s andṪ CCT = −20 K/s, coexisting Martensite and stable α s transformation were found. These coexisting α-morphologies are sometimes also described as massive-α [4]. Eventually, for cooling rates slower thanṪ CCT = −20 K/s only stable α s emerges. In the regime of fast cooling rates, the transformation of Martensite starts at temperatures below the Martensite-start-temperature T αm,sta . For slower cooling rates stable α s -nucleation is the dominating effect, which can already be observed at higher temperatures. From Figure 6 it can already be concluded that the characteristic kink in the isoline X α = 0.01 atṪ CCT = −410 K/s, which separates the pure Martensite from the massive-α regime, is predicted very well by the proposed model. A more detailed discussion and comparison will be presented in the following sections.

Inverse determination of cooling profiles
In the CCT-simulations presented in the subsequent section, we followed the experimental set-up and results presented in [33]. In their work, the authors realized different cooling rates through convectional air cooling for low and moderate cooling rates and through water quenching for high cooling rates. The material probe in [33] was a Ti-6Al-4V Jominy end quench bar [51] of 31.8 mm diameter cooled down by the mentioned media only at the front face side and thermally isolated at all remaining surfaces. We want to note that analogous CCT-investigations for steel, based on the numerical simulation of Jominy end quench tests, have been conducted in the past [52], nevertheless using a simpler microstructure model. To mimic the temperature evolution of the actual experiments, we approximated the temperature profile by the analytic solution for a semi-infinite solid body (sib) under surface convection, which is given in [53]: Here, erfc denotes the complimentary error function, T ∞ is the temperature of the cooling fluid and T 0 the initial temperature of the solid according to [33]. The variable x indicates the distance coordinate measured from the cooled front surface of the solid and t is the time since initiation of the cooling process. The parameter α describes the thermal diffusivity of the material, chosen according to α = 10 mm 2 s as an averaged value for Ti-6Al-4V based on [54]. The parameter g = h/k defines the ratio of the convective heat transfer coefficient h on the surface of the solid to its thermal conductivity k. Since the specific heat transfer coefficient h from the experiment is unknown, we determine the parameter g in the following in an inverse manner to optimally match four exemplary cooling curves provided in [33]. We relax the original assumption of g = const. underlying the analytic solution (24) slightly by introducing a quadratic temperature-dependence for g(T ). As discussed above, to make the numerically predicted microstructure evolution comparable to the experiments, the model requires cooling curves that match the cooling behavior from the experiments in good approximation. Considering a potentially temperature-dependent parameter g(T ) allows us to represent the experimental cooling curves with higher accuracy. Also from a physical point of view it seems reasonable that g(T ) is not necessarily constant across the large temperature spans relevant for these cooling/quenching experiments: Based on the experimental set-up as described in [33], we can directly set T 0 = 1323 K. Furthermore, as the continuous cooling of the Ti-6Al-4V probes was described to be conducted with air or water [33], we set T 0 = 293.15 K, assuming room temperature for the cooling fluid. The cooling rates in experiments can typically be adjusted by either different cooling air speeds or switching to water as a cooling medium for the extreme case of quenching. In order to keep the number of unknown parameters limited, we only considered air cooling in the numerical realization of the cooling curves. To vary the cooling rates as required for the CCT-experiments in the next section, we will adjust the cooling rate by means of an additional scaling parameter s g , which can be interpreted as a manipulation of the cooling air flow velocity. For the inverse parameter identification in this section, we fix the scaling parameter to the default value s g = 1. The remaining parameters θ thermo = [a g , b g , c g ] T in (25) are specific to the performed experiment and are inversely determined based on four temperature curves (a-, b-, c-and d-curve in Figure 7) experimentally measured at four different positions x along the bar axis at otherwise identical cooling conditions [33]. We use again the maximum   [33] along with the corresponding model-based curves resulting from the identified temperature-dependent parameter g(T ) from Equations (25) and (24). The a-, b-, c-and d-curve were measured at positions x a = 3.2 mm, x b = 9.5 mm, x c = 12 mm and x d = 15.2 mm at otherwise identical cooling conditions. likelihood (ML) point estimate θ * thermo as an optimal choice for the identified parameters required in Equations (24) and (25). Under the assumption of Gaussian measurement noise, the ML point estimate is found to be the minimizer of the squared loss function over all temperature curves: In Equation (26), T exp.,a (t), T exp.,b (t), T exp.,c (t), T exp.,d (t) refer to the experimental temperature measurements in [33] and T sib (x = 3.2, t, θ thermo ), T sib (x = 9.5, t, θ thermo ), T sib (x = 12, t, θ thermo ), T sib (x = 15.2, t, θ thermo ) are the corresponding simulated temperature profiles using Equations (24) and (25). It is emphasized that all four experimental and corresponding model-based temperature curves (differing only in the measurement/evaluation position x) are considered simultaneously in the inverse analysis. The broader data basis resulting from this combined approach is a pre-requisite for a robust inverse identification procedure and an accurate representation of the experimental cooling curves. We utilize again the Levenberg-Marquardt optimizer [48] to solve for θ * thermo in Equation (26), which results in the identified parameter set θ * thermo = [a g , b g , c g ] T = [73.8, −39.3, 6.3] T m −1 . Figure 7 shows a good agreement of the analytic model with the temperature measurements by Ahmed et al. [33]. Equation (24) based on the identified parameters can now be used to generate arbitrary cooling curves for the simulation of CCT-experiments by varying the scaling parameter s g in (25), as shown in the following section.

Simulation-based prediction of CCT-diagrams for Ti-6Al-4V
For the simulation-based creation of CCT-diagrams, we chose an equidistant range of scaling parameter values s g such that 150 cooling curves betweenṪ CCT = −1 K/s andṪ CCT = −600 K/s were realized. Afterwards, the evolving microstructure was again recorded over the t × T -coordinate space in form of contour-lines for the individual phase fractions. Figure 8 displays the results of the CCT-simulations, along with four characteristic cooling curves analytically reproduced by means of (24) as well as the aforementioned experimental a-curve from [33]. Note that in contrast to the   (24) and the experimental a-curve from [33]. Left: contour-lines for X αs . Middle: contour-lines for X αm . Right: contour lines for X β .
TTT-diagrams in Figure 3, the isolines in Figure 8 show the true/absolute phase fractions such that the phase fraction values of all three plots would sum up to one for any given data point. In [33] it was found that cooling rates faster thanṪ CCT = −410 K/s resulted in a fully martensitic transformation of the microstructure while cooling rates beloẇ T CCT = −20 K/s led to a fully diffusional α s -formation. According to Figure 8, the results predicted by our model are in very good agreement with these findings as the contour-line for X αs = 1% coincides almost exactly with the cooling line ofṪ CCT = −410 K/s (see dashed orange line in Figure 8, left) meaning that the proposed model predicts a fully martensitic transformation for cooling rates faster thanṪ CCT = −410 K/s. We want to emphasize that the inverse identification of the microstructure model parameters was solely based on the previously presented TTT-data and not on the CCT data. Thus, the accurate prediction of the CCT diagram confirms the validity of our microstructure model. In contrast to existing approaches, the proposed microstructure model does not explicitly enforce this critical cooling rate as transformation criteria, instead a lower formation bound at roughlyṪ = −410 K/s emerges naturally form the dynamics of the diffusion equations, which is a very strong argument for the generality and consistency of the proposed modeling approach. Furthermore, the predicted results show also a fully diffusional α s -formation (i.e. X αm is close to zero) for cooling rates slower thanṪ CCT = −20 K/s (see dashed green line in Figure 8), which is also in very good agreement with the experimental observations by Ahmed et al. [33]. Besides these characteristic quantitative values, also the overall qualitative appearance of the predicted CCT-diagrams matches theoretical presentations in [33,50] very well.
SLM model of this type [59,61] will be applied to provide the temperature field during SLM as required by the proposed microstructure model. We will first introduce the (thermal part of the) SLM model and the numerical set-up, then demonstrate the microstructure evolution during SLM for selected points of the domain over time and eventually show snapshots of micro-structural states for cross-sections and different base-plate temperatures.

Macroscale thermal model for SLM process
The microstructure model is implemented in the parallel, implicit finite element code Diablo [71]. The evolution of the microstructure is driven by a one way coupling of the emerging thermal field. The latter is the solution of the balance of thermal energy, and the associated boundary and initial conditions: where Γ T is the portion of the total boundary ∂Ω subject to essential boundary conditions, and Γ q is the portion of the total boundary subject to natural boundary conditions. The constitutive behavior is characterized by a temperaturedependent Fourier conduction of the form where k is a second-order tensor of thermal conductivities, which may also be a function of spatial coordinates. For the current problem, the heat source r in Equation (27a) plays an important role, that being to represent the deposition of laser energy into a powder. While there exist various models in the literature, the model in the current implementation was taken from [72] in an effort to provide the most natural description of the physical process.
The numerical implementation consists of an iterative nonlinear solver that uses consistent Fréchet derivatives, with the solution computed on first-order finite elements via an iterative method to solve the linear system of equations. Time integration is performed using the generalized trapezoidal rule. The code uses distributed memory parallelism to speed up the solution of the spatial problem. The solution of the thermal field is evaluated at the Gauss points of the finite element discretization. For a full description of the model and its implementation, see [59,61]. While Diablo contains additional physics, including solid mechanics and mass transfer, a detailed description is omitted at this point and we direct the interested reader to [71].
In the numerical demonstration, we are investigating the selective laser melting process of a onemillimeter-sided cube onto a base-plate subject to different values of constant Dirichlet boundary conditions T bp ∈ {303 K, 500 K, 700 K, 900 K, 1100 K, 1300 K} (which are applied solely on the "bottom" of the baseplate, Γ T = min (z)) and the associated effect on the emerging microstructure distribution. The domain is mostly insulated, with the exception of the free surface (Γ q = max (z)), for which there exists a Neumann boundary condition that accounts for the energy loss due to both radiation and evaporation. The initial condition is defined as T 0 = 303 K. The geometric set-up and one snapshot of the emerging temperature field and phase variable are depicted in Figure 9.
These problems were run using 128 cores, and typically took around 15 to 20 hours of wall clock time to complete. The significant resources needed to simulate laser powder bed fusion (LPBF) problems depend on several factors, including its multiscale nature, which can be relevant with respect to both the spatial and temporal domains. The spatial scales are resolved via the use of h-refinement, with refinement indicators tailored to the LPBF problem, as can be seen in Figure 9. With respect to the temporal scales, is it noted that the large 2 number of time steps (associated with the active simulation time and the high speed of the laser) necessitates more sophisticated methodological treatment in order to reduce the wall clock run times to the useful range without significantly degrading the solution accuracy. Both the geometric and temporal scale issues are currently being addressed by the Diablo development team, as addressed in the manuscripts [73,74]. Figure 9: Representative state of the thermal model. The top frame depicts the mesh, for which refinement is mediated by distance from the free surface, the bottom left frame depicts the temperature field, and the bottom right frame depicts the internal phase variable (with gray representing the baseplate, which is consolidated material throughout the simulation) as well as containing several temperature contours, with blue being T = 505 K, red being the melt temperature, and yellow being the vaporization temperature.

Microstructure evolution for the center node of the one millimeter-sided cube
We start the part-scale demonstrations by analyzing the microstructure evolution over time at the cube's center point P = [0.0, 0.0, 0.5] T mm. The layer that contains this point is processed at t ≈ 1.72 s. Figure 10a to 10f depict the respective temperature profiles along with the temperature rates and the corresponding microstructure evolution for all six investigated base-plate temperatures. The temperature profile is supported by characteristic temperatures such as the Martensite-start temperature T αm,sta , the α s -end and start temperatures, T αs,end and T αs,sta as well as the solidus temperature T sol and liquidus temperature T liq .
For T bp < 900 K, respectively in Figure 10a to 10c, we can observe the formation of Martensite directly after the first melting by the laser in the time period 1.72 s < t < 3 s and then again slightly before t = 4 s, followed by a subsequent stabilization of the Martensite phase The amount of formed Martensite phase for t > 4 s is dependent on the metastable Martensite equilibrium phase fractionX eq αm in Equation (5) and decreases with rising T bp . As stated in Equation (6), the metastable Martensite phase will transform to α s for t → ∞. Nevertheless, the diffusion rate for X αm transformation is at lower temperatures (T < 700 K) so small (see Figure 4) that Martensite dissolution, respectively the formation of stable α s remains unnoticeable. At a base-plate temperature of T bp = 900 K (see Figure 10d) no Martensite formation is recorded as the center node's temperature stays above the Martensite start temperature T αm,sta but below the α-transus end temperature such that T αm,sta < T bp < T αs,end . The latter condition will result in the maximum amount of stable α s -transformation as the α-equilibrium phase fraction X eq α = 0.9 at this temperature is still the same as at room temperature. Figure 10d only shows the first 13 seconds of the SLM process but the α s -nucleation will continue until X αs = X eq α = 0.9. Other points within the cube might nevertheless experience a martensitic transformation when the steady-state temperatures fall below the Martensite start temperature T αm,sta .
For a base-plate temperature of T bp = 1100 K (see Figure 10e), which fulfills T αm,sta < T αs,end < T bp < T αs,sta , we observe a similar characteristic but at slower α s -formation rates and a lower equilibrium phase fraction X eq α < 0.9. Eventually, for a base-plate temperature of T bp = 1300 K (Figure 10f), which is above T αs,sta , the microstructure at P stays in the β-regime for all times (as long as the base-plate temperature is held at T bp = 1300 K) and no α s -or α m -transformation is observed.
We can summarize the following points for the microstructure evolution at the cube's center node P : • Without preheating of the base-plate, respectively for moderate preheating (T bp < T αm,sta ), the Martensite phase is dominating the microstructure in the deposited material. • For T bp ≤ 900 K the center point experiences six liquid-solid transitions that are followed by martensitic transformations. For higher base-plate temperatures T bp > T αm,sta the material melts more often but the martensitic transformation is not triggered anymore. The visible temperature peaks in Figures 10a to 10f consist actually of several very closely packed peaks that are caused by the neighboring laser tracks that heat up the center node P again. This process happens so fast (the scanning speed in the simulation was 600 mm/s) that these peaks are not distinguishable on the presented time-scale. The 17 visible peaks corresponds to the number of layers that are processed above the center node P such that the cube consists in total of 34 layers. • The absence of martensitic transformation for T bp ≥ T αm,sta in combination with slower formation dynamics of the stable α s -phase is especially interesting from a mechanical perspective. The β-phase is known to be softer than the α s -or α m -phase [75], such that build-up of residual stresses is expected to be smaller if the microstructure is dominated by the β-phase for longer times. • Theoretically, Martensite will dissolve for t → ∞ but the transformation will only be noticeable for T > 700 K.

Microstructure distributions in center plane of the one millimeter-sided cube
In this section, we continue the investigations of the evolving microstructure in the SLM example of the one millimeter sided cube. In contrast to the former local analysis of the center node, we now want to investigate the resulting microstucture distributions within a complete cross-section in the vertical x-z-center plane that contains the center node P of the cube, for several selected snap-shots in time. We furthermore consider the exemplary base-plate temperatures T bp = 303 K and T bp = 900 K to compare the room-temperature case with the first preheated case that does not result in martensitic transformation. Snapshots are taken for the physical times t ∈ [3, 4, 15] s and phase fractions X αs , X αm and X β are depicted in Figures 11a to 11d.
For a base-plate at room temperature (T bp = 303 K) we observe spatially strongly heterogeneous microstructure distributions during the ongoing SLM process (see Figure 11a and 11b). The Martensite phase in these cases is propagating upwards (in positive z-direction) starting from the base-plate towards the last processed layer.
When looking at the results for the two investigated cases at t = 15 s (see Figure 11c and 11d), we find a rather homogeneous microstructure distribution in both cases with only a slight deviation for the case T bp = 900 K. For T bp = 300 K, we observe an almost complete Martensite transformation. On the contrary, for T bp = 900 K no Martensite is formed. Note, that for this case the phase transition from α s to β is still ongoing (until an equilibrium phase fraction of X αs ≈ 0.9 is reached) at the depicted snap shot at t = 15 s. Moreover, we notice a slightly higher α s -concentration at the bottom of the cube, i.e. a higher amount of β has already transformed into α s at t = 15 s, which is caused by the higher base-plate temperature that leads to higher nucleation rates k →αs (T ) (see Equation (15)).
Coming back to the case T bp = 300 K, it is emphasized that for larger geometries we would in general expect a more heterogeneous microstructure state due to a more heterogeneous temperature distribution with strong spatial gradients at the surfaces of the geometry. An extended dwelling time at elevated temperatures, due to the increased thermal mass of larger parts in combination with an increasing amount of absorbed laser energy, would result in longer times at a temperature band (T αm,sta < T < T αs,sta ) favoring the diffusional dissolution of α m into α s in the core region of such parts, while the higher cooling rates and lower temperature levels at the free surfaces are expected to foster remaining Martensite phase fractions. In the next example, a component of larger size will be analyzed to verify these considerations.
Figure 11a: Simulated microstructure distributions for a SLM process of a one millimeter-sided cube within its vertical x-z-center plane that contains the center node P . The microstructure state is shown for a base-plate temperature of Tbp = 303 K and the processing time t = 3 s. The gray area above marks the layers that are not yet processed by the laser at the time t = 3 s. In the right figure, which depicts the β-phase fraction X β , the melt-pool is indirectly visible through the decreased β-phase fraction in the right corner of the up-most layers due to the laser that has just previously scanned this plane in x-direction from left to right. In a similar fashion, the decreased amount X β in the upper left corner of the right Figure stems from the heat of the laser that is currently melting the subsequent track with a now 45 • rotated scanning direction (diagonal scanning), starting in a corner of the cube.
Figure 11b: Simulated microstructure distributions for a SLM process of a one millimeter-sided cube within its vertical x-z-center plane that contains the center node P . The microstructure state is shown for a base-plate temperature of Tbp = 303 K and the processing time t = 4 s. Martensitic transformation propagates from the bottom of the cube to its top, which can be seen by the higher Xα m -phase fractions towards the bottom of the cube in the left figure. x Figure 11c: Simulated microstructure distributions for a SLM process of a one millimeter-sided cube within its vertical x-z-center plane that contains the center node P . The microstructure state is shown for a base-plate temperature of Tbp = 303 K and the processing time t = 15 s. A homogeneous microstructure state with full martensitic transformation (Xα m = 0.9) can be observed throughout the cube.
Figure 11d: Simulated microstructure distributions for a SLM process of a one millimeter-sided cube within its vertical x-z-center plane that contains the center node P . The microstructure state is shown for a base-plate temperature of Tbp = 900 K and the processing time t = 15 s. Due to the preheated base-plate no Martensite is formed. The diffusional phase transition from αs to β is not finished at the depicted snap shot. Figure 11: Microstructure distributions within the vertical x-z-center plane that contains the center node P of the one millimeter-sided cube. The sub-figures show the resulting microstructure state of the phase-fractions X αs , X αm and X β at different times of the SLM process and or two different base-plate temperatures T bp = 303 K and T bp = 900 K.
5 Part-scale demonstration: Microstructure evolution during a quenching process As SLM simulations with consistently resolved laser heat source for larger geometries (one centimeter sided cubes and above) are still hampered by the associated computational costs, even when most modern HPC systems are considered, we now want to investigate geometry-scaling effects and the effect of the increased thermal mass on the microstructure through rapid convective cooling of a 100 mm -sided cube. The numerical examples do not intend to quantitatively mimic SLM processes but rather demonstrate qualitative microstructural characteristics that evolve for parts of practical relevant size under rapid cooling or quenching. These kind of heat treatments have broad applications in metal processing.
We model a preheated cube at T 0 = 1300 K which is in thermo-mechanical contact with a base-plate with a Dirichlet boundary condition at the bottom of first T bp = 300 K and then in a second investigation of T bp = 900 K (see Figure  12, red line). The cube is subject to a convective boundary condition with heat-transfer coefficient α c = 1000 W m 2 K on its free surfaces ( Figure 12, blue lines). The surrounding atmospheric temperature is set to room temperature at T ∞ = 300 K. The interface between base-plate and the bottom of the cube is modeled via thermo-mechanical contact interaction employing a numerical formulation recently developed in [76] and choosing a thermal contact resistance that is equivalent to an effective heat-transfer coefficient of α tc = 5 · 10 5 W m 2 K (Figure 12, bold black line). The remaining free surfaces of the base-plate are assumed to be adiabatic (see Figure 12, green lines). Note that the order of magnitude of the heat transfer coefficient α c = 1000 W m 2 K at the free surfaces of the cube was chosen to mimic convective cooling, e.g., due to forced air flow. As we only want to demonstrate qualitative results here, we omit a further detailed description of a specific cooling scenario. These quenching simulations were conducted with our in-house multi-physics framework BACI [77].
As a consequence of the convective cooling and the thermo-mechanical contact, the cube starts cooling down until thermodynamic equilibrium is reached. We simulated the microstructure evolution for this cooling process and investigated the resulting steady-state microstructure state at t = 5000 s. Figure 13 shows the resulting crystallographic distribution for a base-plate temperature of T bp = 300 K and Figure 14 demonstrates the result for T bp = 900 K.
For both investigated cases, we see a several millimeter strong Martensite coating, which is especially pronounced at the corners of the cube that are characterized by higher temperature rates and lower temperature levels. These qualitative findings are in good agreement with experimental results for quenching of Ti-6Al-4V as well as additively manufactured parts [4,78]. In case of T bp = 300 K ( Figure 13) we observe an even higher Martensite amount on the base-plate-cube interface due to the higher heat transfer at the base-plate. This effect is inverted when the base-plate is heated to T bp = 900 K ( Figure 14) and stable α s -phase can be found instead of the previous Martensite phase. Similar to the preceding SLM demonstration, a pre-heated base-plate resulted in a significant reduction of metastable Martensite phase for material close to the cube-base-plate interface. The core of the material in both demonstrations ( Figure 13 and  Figure 12: Schematic set-up for cooling simulation of a 100 mm -sided cube on a base-plate. The initial temperature of the cube is prescribed to T 0 = 1300 K. The atmospheric temperature is set to T ∞ = 300 K and the free surfaces of the cube have a convective boundary condition (blue) with heat transfer coefficient α c = 1000 W m 2 K . The base-plate has a Dirichlet boundary condition on the bottom (red) of first T bp = 300 K and, in a second simulation run, of T bp = 900 K and is assumed adiabatic on its free surfaces (green). The cube and the base-plate are modelled to be in thermo-mechanical contact (bold black line) using a thermal contact resistance that is equivalent to an effective heat-transfer coefficient of α tc = 5 · 10 5 W m 2 K .  Figure 13: Simulated microstructure distribution of α s -and α m -phases in the vertical center plane of the 10 cm -sided cube after 5000 s of cooling. The simulation used a thermal heat-transfer coefficient of α c = 1000 W m 2 K , a base-plate temperature of T bp = 300 K and a thermal contact resistance that is equivalent to an effective heat-transfer coefficient of α tc = 5 · 10 5 W m 2 K .  Figure 14: Simulated microstructure distribution of α s -and α m -phases in the vertical center plane of the 10 cm -sided cube after 5000 s of cooling. The simulation used a thermal heat transfer coefficient of α c = 1000 W m 2 K , a base-plate temperature of T bp = 900 K and a thermal contact resistance that is equivalent to an effective heat-transfer coefficient of α tc = 5 · 10 5 W m 2 K . Figure 14) is composed of stable α s -phase due to the increased thermal mass of the material when compared to our former SLM examples with a 1 mm sided cube in Section 4.3.
Also for larger SLM-manufactured geometries, Martensite can be expected mostly in proximity to surfaces in form of a Martensite coating, which is in compliance to experimental findings. During the SLM process, temperature rates in regions close to the melt-pool are extremely high such that we would initially expect a purely martensitic transformation as soon as the temperature decreases below the Martensite start-temperature T αm,sta . Nevertheless, the re-heating of the material during the processing of subsequent layers and tracks leads to lower cooling rates in subsequent thermal cycles on the other hand, and to an increasing overall temperature niveau within the deposited material. This effect might lead to longer dwelling times with conditions that are suitable for diffusion-based Martensite dissolution into α s -phase such that the core of larger SLM processed parts can also be dominated by the stable α s -phase, depending on the specific geometry and scanning strategy. With the development of more efficient part-scale simulation strategies for selective laser melting, we are planning on answering these questions in future investigations.

Conclusion and Outlook
In this work, we proposed a novel physics-based and data-supported phenomenological microstructure model for part-scale simulation of Ti-6Al-4V selective laser melting (SLM). The model predicts spatially homogenized phase fractions of the most relevant microstructural species, namely the stable β-phase, the stable α s -phase as well as the metastable Martensite α m -phase, in a physically consistent manner, i.e. on the basis of pure energy and mobility competitions among the different phases. The formulation of the underlying evolution equations in rate form allows to consider general heating/cooling scenarios with temperature or time-dependent diffusion coefficients, arbitrary temperature profiles, and multiple coexisting phases in a mathematically consistent manner, which is in contrast to existing approximations via JMAK-type closed-form solutions.
Altogether, the model contains six free (physically motivated) parameters determined in a robust inverse identification process on the basis of comprehensive experimental time-temperature transformation (TTT) and transient heating data sets. Subsequently, it has been demonstrated that the identified model predicts common experimental procedures such as continuous-cooling transformation (CCT) experiments [33] with high accuracy and reflects well-known dynamic characteristics such as long term equilibria or critical cooling rates naturally and without the need for heuristic transformation criteria as often applied in existing models.
The part-scale simulation of selective laser melting processes with consistently resolved laser heat source is still an open research questions. In contrast to existing microstructure models that resolve the length scale of individual crystals/grains, the proposed continuum approach has the potential for such part-scale application scenarios. To the best of the authors' knowledge, in the present contribution a homogenized microstructure model of this type has for the first time been applied to predict CCT-and TTT-diagrams, which are essential means of material/microstructure characterization, and to predict the microstructure evolution for a realistic SLM application scenario (employing a state-of-the-art macroscale SLM model) and for the cooling/quenching process of a Ti-6Al-4V cube with practically relevant dimensions.
In the investigated SLM process, Martensite could be identified as the dominating microstructure species due to the process-typical extreme cooling rates and the comparatively small part size considered in the present study. A preheating of the material (e.g. via a preheated base-plate) resulted in a decreased Martensite formation and higher phase fractions of the stable α s -phase. In a subsequent simulation, the rapid cooling/quenching process of a larger cube (side length 10 cm) resting on a cold metal plate and subject to free convection on the remaining surfaces was considered. It was demonstrated that the high cooling rates in near-surface domains can lead to a strong martensitic coating of several millimeters, a behavior that is well-known from practical quenching experiments. The slower average cooling rates and the higher thermal mass in this example resulted in a large core domain dominated by the stable α s -phase.
In future research work, the proposed microstructure model shall be employed to inform a nonlinear elasto-plastic constitutive model, thus contributing to the long-term vision of achieving accurate thermo-mechanical simulations of selective laser melting processes on part-scale. Furthermore, SLM process parameters shall be inversely adjusted to yield specific microstructural distributions and hence desired mechanical properties by deploying novel efficient multi-fidelity approaches for (inverse) uncertainty propagation [79].