Data workflow to incorporate thermodynamic energies from Calphad databases into grand-potential-based phase-field models

In order to approximate Gibbs energy functions, a semi-automated framework is introduced for binary and ternary material systems, using Calphad databases. To generate Gibbs energy formulations by means of second-order polynomials, the framework includes a precise approach. Furthermore, an optional extensional step enables the modeling of systems in which a direct generation leads to the unsatisfactory results in the representation of the thermodynamics. Furthermore, an optional extensional step enables the modeling of systems, in which a direct generation leads to the unsatisfactory results, when representing the thermodynamics. Within this extension, the commonly generated functions are modified to satisfy the equilibrium conditions in the observed material systems, leading to a better correlation with thermodynamic databases. The generated Gibbs energy formulations are verified by recalculating the equilibrium concentrations of the phases and rebuilding the phase diagrams in the considered concentration and temperature ranges, prior to the simulation studies. For all comparisons, a close match is achieved between the results and the Calphad databases. As practical examples of the method, phase-field simulation studies for the directional solidification of the binary Ni–35Mo and the ternary NiAl–10Mo eutectic systems are performed. Good agreements between the simulation results and the reported theoretical and experimental studies from literature are found, which indicates the applicability of the presented approaches.


Introduction
The phase-field method is widely used to investigate the microstructure evolution during the solidification processes in multicomponent material systems. To calibrate this method with physical data, thermodynamic information needs to be incorporated. The information required to describe the driving forces for the phase transition is mostly derived from CAL-PHAD databases [1]. In these databases, the Gibbs energy for each contained phase is stored as a function of pressure, temperature and concentrations.
To incorporate this thermodynamic information into the different phase-field models, based on the Allen-Cahn or Cahn-Hilliard model, several approaches have been established [2][3][4][5][6][7][8][9][10][11][12]. During the simulations, an obvious approach is the direct access to the information from the databases. This approach, for example, is used by Steinbach et al. [8], in their coupling study from 2007. To reduce the computational effort for the simulations in this work, the calculations of the thermodynamic information are only performed at certain time intervals. Based on these calculations, the quasi-equilibrium stages in between are extrapolated. A further approach, used by Qin and Wallach [9], is to precompute the thermodynamic information before the simulation, using MTDATA [13], and store the results in a data file. This allows a reduction in the computational effort, during the simulation. However, if not all required information is sufficiently precalculated, the data need to be approximated by means of interpolation, which can result in a loss of accuracy.
Another approach is to directly incorporate the Gibbs energy functions from the CALPHAD databases into the model, without calling external libraries. This approach was used by Böttger et al. [10] and Zhu et al. [11]. Their procedures allow the exact determination of thermodynamic information from the CAL-PHAD databases. As the functions in the CALPHAD databases are often stored in computationally intensive formulations, iterative methods are required to solve them. Therefore, further approaches, which use more simply approximated functions of the Gibbs energies, have been established to efficiently compute the thermodynamic conditions during the simulations. A commonly used procedure is the approximation of the Gibbs energies, using parabolic functions. The accuracy of this approach depends on the considered temperature and concentration ranges. With increasing ranges, the ability of the parabolic functions, to reproduce the thermodynamic information, becomes limited. Hence, it is challenging to represent a broad range of temperature-composition space with a high accuracy. In [3], Welland et al. compare their parabolically approximated functions with functions composed of an approximated minimizer approach. For both approximation approaches, the authors show a good accordance with the original data, in the vicinity of the equilibrium states. With a growing distance to the equilibrium, the parabolically approximated functions show a stronger deviation from the original data, whereas the approximated minimizer functions are still in good accordance. However, for most processes investigated with phase-field simulations, a limited range is sufficient to reproduce the phase transitions of the particular interest to a specific problem.
The advantages of the parabolically approximated functions are their lightweight computations and the proper representation of the thermodynamic information. In recent years, different binary [3,5,[14][15][16] and ternary [5,[17][18][19][20][21] alloys have therefore been modeled with parabolically approximated functions. A further advantage of this approach is that these functions are suitable for the application in free-energy-based [8,22] as well as in grand-potential-based phase-field models [23,24]. For the grand potential models, all required parameters, such as the concentrations, the chemical potentials and the grand potentials, can be calculated directly and uniquely from the approximated functions.
In most of the above-mentioned publications, the focuses of the works are on the microstructure evolution, which is highly affected by the system and the process parameters. In these works, the challenges included in the modeling of the Gibbs energy functions and the modeling process itself are rarely discussed in great detail. These challenges can be summarized as the accuracy of the modeled functions, in representing the thermodynamics and their suitability for the utilization in phase-field models. Therefore, an efficient and semi-automated framework is presented in this paper, so as to generate Gibbs energy formulations for the application in binary and ternary material systems, which is exploited in the introduced phase-field model of [17,23]. The functions in this framework are derived on the basis of a numerical least squares method, using near-equilibrium concentrations, so as to prevent their validity from being solely confined to the equilibrium conditions. To reduce the models' sensitivity to external influences, further thermodynamic information is therefore taken into account, in the vicinity of the equilibrium conditions. In order to generate formulations for the Gibbs energy, the presented framework primarily follows a common approach. For the case of a more complex material system, however, it additionally includes an optional extension of the general approach, for which the common approach leads to unsatisfactory results. To demonstrate the applicability of the framework, the derivation of parabolic Gibbs energy formulations is performed for the simulations of the directional solidifications in binary Ni-35Mo and ternary NiAl-10Mo eutectic systems, which are two practical examples of non-trivial material systems. In the literature, phase field simulations have not yet been carried out for these systems. The obtained microstructures are verified by the reported theoretical and experimental results of the individual systems.
In the upcoming sections, the utilized phase-field method is initially introduced in more detail. In the next section, both the framework for the generation of Gibbs energy formulations and its extension are presented in general. Afterward, the workflow of the framework is exploited in detail, in order to obtain the Gibbs energy information for the mentioned binary and ternary material systems. Finally, twodimensional simulation studies are conducted, in order to verify the procedure and the resulting Gibbs energy formulations.

Phase-field model
The used phase-field model is based on the grand potential functional [23,24]. A detailed description of its implementation is given in the literature [17,25,26], and further examples for its application to the directional solidification of eutectic systems are reported in [17,19,20,[27][28][29][30]. The local phase fractions of the phases are represented by the N order parameters /â, which are stored in the vector /. Similar to the works [20,31,32], the phases a; b; . . . differ from their indicesâ andb by a different labeling. The K chemical potentials l i are derived from the mass balance of the concentrations and are collected in the vector l. The time evolution equations o ot À Á of the coupled phase fields (based on an Allen-Cahn approach), the chemical potentials (based on Fick's law) and the temperature T read as: Equation (1) describes the phase transition, taking place during the simulations, with a diffuse interface. The shape of this interface is modeled by the gradient energy density að/; r/Þ and by the potential energy density xð/Þ. While the gradient energy density is a function of the phase fields and their gradients, the potential energy density is only a function of the phase fields. Both further depend on the interfacial energies câb between the evolved phases [17]. The interface thickness is related to the parameter e, and the movement of the interface is controlled by the kinetics of the interface and by the driving forces from the investigated system. The kinetics of the interface, which is a function of the diffusion coefficients, and the concentrations are described by the parameter s, based on [23], and the driving forces are described by the differences between the phase-dependent grand potentials wð/; l; TÞ. The term 1 N P N b¼1 ðrhs 1;b þ rhs 2;b Þ is the Lagrange multiplier K, which is introduced to fulfill the constraint P N a¼1 o/â ot ¼ 0. By calculating the evolution of the chemical potentials l, Eq. (2) describes the diffusion processes in the investigated system. The function Mð/; l; TÞ is the mobility of the interface and includes the information of the diffusion coefficients D of the involved phases [23]. In Eq. (2), the function hâ [33] interpolates between the different phases, and the anti-trapping current J at [23,34,35] balances the effects of the artificially enlarged interface. For the investigated phaseâ, the parameter câðl; TÞ represents the concentration vector of the K components. Equation (3) describes the evolution of the temperature T, in the growth direction x, using the initial base temperature T 0 , the applied gradient G and its velocity v G .
To ensure the thermodynamic consistency of the derived phase-field model, the thermodynamic properties, required to calculate the phase transitions, are derived from thermodynamic databases. These properties mainly consist of the grand potentials w a , the concentrations c a ðl; TÞ and the chemical potentials l. Following Kellner et al. [19], and by assuming a constant pressure and volume, the grand potentials can be defined as: w a ðl; TÞ ¼ g a ðc a ðl; TÞ; TÞ À l Á c a ðl; TÞ : As the concentrations c a ðl; TÞ and the chemical potentials l can be derived from the Gibbs energy g a ðc; TÞ, both the grand potentials and the driving force can be expressed as a function of the Gibbs energies and their derivatives, with respect to the element concentrations [19]. Formulations for the Gibbs energies of different material systems can be found in thermodynamic CALPHAD databases. In the databases, the Gibbs energies are stored in the form: The term g 0 of Eq. (5) describes the Gibbs energy, due to the mechanical mixing of the contained phases. The second term id g mix represents the Gibbs energies, due to the ideal mixing contribution, and describes a statistic distribution of the concentrations in the phase. During the calculations with CALPHAD software packages like FactSage [36], OpenCalphad [37], Pandat [38], Pycalphad [39] and Thermo-Calc [40], this part is automatically added to the calculated Gibbs energies. The last term of Eq. (5), xs g mix , is the excess Gibbs energy of mixing and describes all nonideal mixing contributions. The different terms can refer to either a single or multiple sublattice model [1,41]. As mentioned before, these formulations are often computationally intensive, especially J Mater Sci (2021) 56:11932-11952 the calculation with multiple sublattices. Furthermore, the required derivation of the chemical potentials l and the concentrations c a ðl; TÞ from the CALPHAD formulations can lead to solutions that are not unique. Therefore, a simpler approximation of these formulations, which leads to unique solutions of the derivations, is required for the used grand potential model. To generate these approximations, a framework has been developed, which transfers the computationally intensive CALPHAD formulations into simpler terms, with respect to the original data and the equilibrium conditions of the phases. A comparison between the required computational effort for Gibbs energy calculations, based on the CALPHAD formalism, and the later used parabolic approach is performed in A. Based on the result, and depending on the number of components in the system, the parabolic approach can reduce the number of necessary cycles in the simulation studies by up to 7 times, for a single sublattice. In case of multiple sublattices, this amount can increase several times.

Generation of Gibbs energy functions
For binary and ternary material systems, the semiautomated framework for the generation of the approximated Gibbs energy functions is schematically illustrated in Fig. 1 and is explained in this section. The tool is implemented to be utilized prior to the simulation studies. Before starting the calculations, the procedure requires a detailed definition of the material systems and the investigated phase transformation reaction. Apart from the CALPHAD database, the considered components and the reaction type have to be specified. Depending on the chosen numbers of components, an automated distinction between binary and ternary systems is realized. For a further definition of the phase transformations, the considered temperature range and the concentration of the melt, normally the eutectic composition of the melt, must be known. To ensure a solidification of the solid phases in the subsequently performed phase-field simulations, the temperatures should be chosen in such a way that they are below the melting temperatures of the individual phases. It has to be mentioned that the choice of the temperatures can have a critical impact on the accuracy of the generated Gibbs energy formulations and hence on the stability of the resulting phase-field simulations. In order to obtain a satisfactory formulation of the Gibbs energies for all involved phases, multiple calculation runs with different temperature ranges can be required.
The generation process of the approximated temperature-dependent Gibbs energies can be divided into five steps, with an optional modification loop between the third and the fifth step. STEP 1 The equilibrium concentrations for each phase c eq a , evolved in the examined reaction, are calculated based on the considered temperatures and concentrations of the melt, by a software that can handle Calphad databases. In this work, the Thermo-Calc software package [40] is used. STEP  The presented workflow within the Pace3D framework [42] consists of several C-based tools, used to generate macrofiles for the employed CAL-PHAD software package Thermo-Calc and to subsequently verify the calculation results from the database. By using a numerical least squares method, the exact equilibrium conditions are usually not fully reproduced. However, by adjusting both the concentration range and the number of calculation points around the equilibrium condition in STEP 2, the accuracy of the calculations can be increased. Thus, a close match between the original and the derived thermodynamic data can be achieved in a certain range around the equilibrium conditions. Furthermore, different binary (B) and ternary (T) approaches can be chosen for the approximated formulations. Apart from the mentioned parabolic approaches, further approximation function approaches are implemented into the Gibbs energy framework. For binary and ternary systems, a summary of the implemented approaches is given in Table 1. For the required derivations, all implemented functions lead to unique solutions.
In STEP 5, the approximated Gibbs energy functions are checked to verify the accordance with the CAL-PHAD data and the satisfaction of the thermodynamic equilibrium, at the desired temperatures and concentrations of the elements. Therefore, the original and the derived data of the Gibbs energies are compared visually, as exemplified later in Figs. 2 and 4. For a quantitative comparison, the maximum and average deviations between both are calculated, and the chemical potentials at the equilibrium concentrations are checked for equality, so as to ensure a correct representation of the equilibrium conditions. For nonvariant reactions, a phase diagram is additionally derived from the fitted Gibbs energy formulations and compared with the phase diagram from CALPHAD, as shown later in Figs. 3 and 5. If the functions do not satisfactorily reproduce the original data, the mathematically obtained Gibbs energy formulations of STEP 3 are modified by including more thermodynamic criteria. These criteria can be stated as: 1. Equality of the resultant chemical potentials, at the equilibrium concentrations, for all involving phases. 2. Existence of a common tangent between the equilibrium concentrations of the solid phases, at the regarded temperature.
Considering a total number of N phases (including liquid as the last phase) and a chemical potential of phase i as l i ¼ og i ox for binary and l ij ¼ og i ox j , j ¼ 1; 2 for ternary systems, the following equations are obtained to ensure equal chemical potentials: and a common tangent: J Mater Sci (2021) 56:11932-11952 binary systems After the modifications, steps 4 and 5 are repeated to validate the newly generated Gibbs energy formulations. As the conditions for the validation are already considered in the modeling of the new Gibbs energy formulations, an accordance between the modeled and the original CALPHAD material system can be ensured.
The modification of the fitted Gibbs energy functions should not be considered as the adjustment of the CALPHAD data. In fact, the approximation process enables a simplified and computationally efficient reproduction of the original data, which can be used in the simulation studies. The additionally performed modifications in this process increase the accordance with the expected physical properties of the systems.
It has to be mentioned that with the presented approach, a thermodynamic modeling of the Gibbs energies cannot be realized for all material systems. The tool chain is limited to binary and ternary material systems and, in addition, cannot handle stoichiometric phases. However, as discussed in the following, the approach for systems with well-defined equilibrium conditions, without stoichiometric phases, directly leads to a satisfying result, for the approximated Gibbs energy formulations, even without a modification of the originally derived functions. In order to generate Gibbs energy functions for the eutectic systems Ni-35Mo and NiAl-10Mo, respectively, exemplary calculation chains are presented in the following sections. Ni-35Mo Table 1 Implemented approaches for approximated functions   Label Type Equation Comparison of the CALPHAD data points [43] and the fitted functions in the representation of the Gibbs energies, at 1590 K. Points: CALPHAD data, solid lines: fitted curves, dashed line: common tangent of the solid phases.

Figure 3
Comparison of the phase diagram, based on the CALPHAD database [43], and the fitted Gibbs energy functions used in this work, for the binary Ni-35Mo system.
describes a eutectic reaction in the binary Ni--Mo system and NiAl-10Mo describes a eutectic reaction in an isopleth section of the ternary system Ni--Al--Mo.
In addition to the validations of the generated Gibbs energy formulations, performed in STEP 5 of the tool chain, the suitability of the generated functions is investigated for the conduction of directional solidification processes in phase-field simulations. Therefore, the results of the phase-field simulations are subsequently presented in the mentioned systems.
Generation of Gibbs energy functions for the binary Ni-35Mo system In this section, all relevant steps of the framework, required to derive the Gibbs energy functions for the binary eutectic Ni-35Mo system, are explained to demonstrate the applicability of the presented method for a eutectic reaction. The thermodynamic database, utilized for the binary Ni-Mo system, is developed by Yaqoob et al. [43]. Starting from STEP 1 of the tool chain (see Fig. 1), Ni and Mo are defined, while the mole fraction of the molybdenum is set to the eutectic point, i.e., x Mo ð Þ eut ¼ 35:3 mol-%. The pressure of the system is expressed by the atmospheric pressure (101325 Pa), while for the existing components, a total number of 1 mol is defined in the computations. In order to reproduce the phase diagram of the upcoming calculations, the temperatures 1580, 1585 and 1590 are chosen to define a temperature range beneath the eutectic point. Using a macrofile, this information is transferred to the software package Thermo-Calc [40], so as to calculate the thermodynamic properties. In Table 2, the results of the equilibrium calculation from STEP 1 are listed for the middle temperature T ¼ 1585 K.
In order to calculate the Gibbs energies at the temperatures 1580 and 1590, in STEP 2, a user-specified number of points, surrounding the equilibrium concentrations, are selected for the involving phases, so as to obtain the amounts of the Gibbs energies. These points are specified in regular steps, in the vicinity of the equilibrium points. In this case, two hundred data points are selected for each phase, in a radius of 20 mol-% Mo and a step width of 0:1 mol-% Mo. The Gibbs energy values calculated at these points are the basis of the approximation of the Gibbs energy functions in STEP 3. In this case, binary second-order polynomials (B p ), as given in Table 1, are used to approximate the data points. The approximated functions are numerically calculated, by using a least squares method for the individual temperatures. Through a linear interpolation, carried out in STEP 4, the obtained coefficients of the approximated functions have become temperaturedependent. The resultant temperature-dependent functions are listed in Table 5, which can be found in the appendix.
In order to test the accuracy of the fitted functions, with regard to the reproduction of the thermodynamics, different procedures are carried out in STEP 5. In the vicinity of the equilibrium concentrations, the original amounts of the Gibbs energies from the CALPHAD database are first visually compared with the fitted amounts. This comparison is illustrated in Fig. 2. The average and maximum deviations between the original CALPHAD data g cal and the approximated Gibbs energies g app are reported in Table 5, showing a deviation of less than \0:04%, for all three phases. For each phase, the exact amounts of the Gibbs energies of the modeled functions and of the CALPHAD databases, at the equilibrium concentrations, are also summarized in Table 5, for a temperature of 1590 K. Next, the resultant chemical potentials for the involving phases are compared and checked for their equality (Eq. (6)). It can be noted that a maximum deviation does exist between the chemical potentials of 6%, which is negligible in the authors' opinion. The checking procedure is repeated for different temperatures below the melting point, which leads to similar correlations between the CALPHAD data and the approximated functions. Therefore, it can be concluded that the fitted functions behave well at the desired temperature, in terms of the representation of the Gibbs energies and the corresponding chemical potentials. As a final accuracy test, the phase diagram, based on the fitted functions, is rebuilt and compared with the CALPHAD data, by means of the common tangent rule.
In the system, the common tangents between the Gibbs energy curves are constructed for any possible two-phase combinations, i.e., d-fcc, fcc-liquid and dliquid. The points where the tangent lines and the Gibbs energy curves meet describe the equilibrium concentrations of the corresponding phases. This procedure is exemplarily illustrated in Fig. 1 of [3] and in Fig. 8 of [44]. By repeating the procedure for different temperatures in the considered temperature range, the touch points referring to the lowest energy levels are transferred to the temperature-concentration plot of the original phase diagram. By connecting these points, the transition lines of the rebuilt phase diagram are constructed and can subsequently be compared with the original phase diagram. For Ni-35Mo, the result of this reconstruction is given in Fig. 3, together with the original phase diagram of [43].
In the top diagram of Fig. 3, a good accordance is observed between the original and the reconstructed phase diagram. For the single areas around the equilibrium points, a closer look at the enlargements I to III, taken below, shows that the reproduced phase transition lines of the solids deviate slightly from the original lines. With this, the concentration of the eutectic reaction is slightly changed from 35.3 to 35:4 mol-% Mo. However, for all reproduced transition lines in the diagram, a maximum deviation of 0:4 mol-% Mo is found, which results in a maximum percentage deviation of 1%. This good correlation indicates the accuracy of the fitted Gibbs energy functions, concerning the representation of the thermodynamics for the binary Ni-35Mo system. This can also be seen in an enlarged temperature range , compared with the predefined temperature range.
As all validations show a good accordance between the modeled and the original system, a modification of the Gibbs energy functions, as described in Fig. 1, can be omitted. The subsequently performed study of the directional solidification process of the modeled system Ni-35Mo, taking place within phase-field simulations, is presented in section 5.

Generation of Gibbs energy functions, for the ternary NiAl-10Mo system
The thermodynamic properties for the NiAl-10Mo system are calculated on the basis of the database of Peng et al. [45], which describes the quaternary Ni-Al-Cr-Mo system and contains the ternary Ni-Al-Mo subsystem. NiAl-10Mo describes a nontrivial eutectic reaction in an isopleth section of the ternary system Ni-Al-Mo. To obtain different solidification velocities for this system, experimental investigations of the directional solidification process are performed by Zhang et al. [46]. In this experimental work, Mo-rich fibers are observed, embedded in NiAl-rich matrices.
Starting from STEP 1 of the tool chain, the system components Al, Mo and Ni are defined and the mole fractions of the components at the eutectic point are set to 43:9 mol-% Al, 9:4 mol-% Mo and 46:7 mol-% Ni. Based on the database, the eutectic temperature is 1875:47 K. For the calculations, the atmospheric pressure and a total component number of 1 mol are used again. In a procedure similar to the binary Ni-35Mo system, several temperatures below the eutectic point are selected for the calculation of the equilibria, which later results in a reproduction of the considered phase diagram section. In Table 3, the results of the equilibrium calculation from STEP 1 are exemplarily listed for the temperature T ¼ 1871 K. To calculate the Gibbs energies in STEP 2, a user-specified number of points, surrounding the equilibrium concentrations, is selected. Despite the binary case, the composition of a ternary system is defined by the concentrations of two components. Therefore, a regular mesh of data points, with a step width of Dx ¼ 0:1 mol-% Mo, Dy ¼ 0:1 mol-% Al, in a radius of 0:5 mol-%, around the equilibrium concentration, is used to determine the required concentration values for the involving phases. It has to be mentioned that only points in the simplex are selected. In STEP 3, ternary second-order polynomials (T p ), as given in Table 1, are used for the approximation with the least squares method. In STEP 4, the temperature fitting is done by considering the approximation results of the different temperatures. The obtained functions are depicted in Table 7 and are referred to as the initial fitting functions for the Gibbs energies. Based on the results, a maximum deviation of 0:01% is found between the approximated Gibbs energies g app and the original data g cal from CALPHAD. This deviation is found for the fiber phase at the exemplarily listed temperature of T ¼ 1871 K. At the equilibrium concentrations, however, comparisons of the chemical potentials lead to deviations of up to 18% and 13%, for l 0 and l 1 , respectively. To analyze the impact of these deviations, the equilibrium concentrations, which are based on the approximated Gibbs energy functions g app , are recalculated, leading to an unphysical negative amount for the concentration of Mo, in the matrix phase. Hence, the Gibbs energy functions g app do not reproduce the accurate material system. Therefore, the optional modification procedure is used to adjust the fitting coefficients (a 0 ; :::; a 5 ) for all phases, as described in Fig. 1. In this procedure, as mentioned before, the coefficients are modified as variables, in order to satisfy the Eqs. (6) and (7). Consequently, the chemical potentials, obtained by the modified functions, are similar for all phases, with minimum deviations of the resultant Gibbs energies from the original CALPHAD data. To establish the equilibrium conditions, the existence of a common tangent plane is furthermore ensured for the solid phases. The outcome of this procedure is compiled in Table 8 and is labeled as the modified fitting functions of the Gibbs energies g mod . Based on these functions, the maximum deviation between the modified Gibbs energies and the original CALPHAD data at the equilibrium concentrations has increased to 0:02%, for the previously shown exemplary temperature of 1871 K. However, for the mentioned temperature, the chemical potentials are almost equal, showing a maximum deviation of 0.18%. This indicates a good representation of the thermodynamics, by the approximated functions.
For a visual comparison of the modified fitted curves with the CALPHAD database [45], the obtained results are illustrated in Fig. 4. In this figure, isoconcentration lines (x Mo =cte), based on the fitted (a) (b) (c) Figure 4 Comparison between the modified Gibbs energy fitting curves, for the matrix phase in (a), the fiber phase in (b), the liquid phase in (c) and the CALPHAD database [45] of the NiAl-10Mo system. The temperature is 1871 K, and the concentration of Mo is shown next to each curve. functions, are plotted and the amounts of Mo concentrations are shown beside each line. According to the CALPHAD data, the Gibbs energies are sketched with single squares. For the involved phases, the comparison shows a good correlation between the Gibbs energies, in the vicinity of the equilibrium points. The maximum deviation with 0.02% is found in the comparison of the Gibbs energies from the NiAl matrix phase, shown in Fig. 4a.
By using the common tangent rule, the phase diagram around the equilibrium concentrations is reconstructed in a similar procedure as the binary Ni-35Mo system. For the isopleth section NiAl-10Mo, the reconstructed phase diagram in Fig. 5 is compared with the originally derived phase diagram from the CALPHAD database. The enlargements I to III, depicted below the diagram, indicate a good correlation between the phase diagram of the CALPHAD data (black lines) and the reconstructed diagram of the fitted Gibbs energy functions (red marks).
Based on the generated Gibbs energy functions, two-dimensional phase-field simulation studies are performed in the next section, for the previously presented binary Ni-35Mo system and the ternary NiAl-10Mo system, to show that the generated material systems are suitable for the application in phase-field simulations.

Simulation results
To proof the suitability of the generated thermodynamic models, in order to investigate the directional solidification process with the phase-field model, simulation studies for the material systems Ni-35Mo and NiAl-10Mo are conducted. As the simulation results have a high dependency on the used material properties (like interfacial energies, diffusion coefficients, etc.), simple simulation settings are used, in terms of domain and phase fillings, so as to reduce the impact of these values on the microstructure adjustments. To demonstrate the capability of the Gibbs energy fitting method in the current work, all phase-field simulations are performed in two-dimensional rectangular domains, with defined settings of the two solid phases next to each other, beneath the liquid phases. During the solidification, the solid phases grow into the liquid domain. To ensure an independent growth of the solid phases, an infinite liquid domain is ensured, by using a moving window technique [47]. In the domains, periodic boundary conditions are applied in the directions perpendicular to the growth direction. On the solidified side of the domain, a Neumann zero boundary condition is imposed, while on the opposite side, a Dirichlet boundary condition is set. The utilized parameters of the simulations are listed in Tables 6 and 9. In Fig. 6a, b, the exemplary results of two-  dimensional phase-field simulations are, respectively, presented for the systems Ni-35Mo and NiAl-10Mo. Both simulations show the directional solidification process for a growth velocity of 30 lm s À1 . Considering Fig. 6a, a stable growth of the d phase (red color ffi 47 mol-% Mo) is observed next to the fcc phase (dark blue ffi 27 mol-% Mo). In Fig. 6b, which indicates the ternary system, the same color style represents ffi 90 mol-% Mo for the Mo-rich fiber phase and ffi 0.08 mol-% Mo for the NiAl-rich matrix phase. For both simulated material systems, similar stable growth patterns are also obtained at the velocities 15 lm s À1 and 25 lm s À1 . In both cases, the simulation results are in good agreement with the CALPHAD databases, regarding the reproduction of the concentration of the elements in the involving phases.
For both material systems, simulation studies with different growth velocities v growth are performed to study the relationship between the adjusting undercooling DT and the spacing k. For this purpose, the width of the simulation domain for Ni-35Mo is varied between 36 and 76 cells (1:8 lm to 3:8 lm), whereas the domain width for the simulations of NiAl-10Mo is fixed to 150 cells. This configuration is required, due to the adjusting phase fractions in NiAl-10Mo. To resolve the fine Mo-rich fibers in 2D, a minimum domain width of 100 cells is required. By choosing a domain width of 150 cells, it can be ensured that the initial oscillations of the domain width, taking place during the establishment of a stable growth from the start settings, are also sufficiently resolved. For each simulation, however, a constant length scale would require domain widths between 150 and 200 cells, for NiAl-10Mo. To reduce the computational effort of this study, the domain width is set constant and the length scales of a simulation cell are varied instead, in order to reproduce the physical lengths from 1.106 to 1:483 lm.
In Fig. 7a , for NiAl-10Mo. In these plots, the profiles of the resulting undercoolings show the typical Jackson-Hunt-type shape for all six studies. The theory of Jackson and Hunt explains the relationships between the undercoolings, the rod spacings and the growth velocities of the solidifying materials in the eutectic reactions [48]. This theory is well known for its convenient achievements in the explanation of the growth dynamics in the mentioned processes [49,50]. The minima of the curves represent the spacing k ext: , for which the smallest undercooling has been measured. Based on the experimental work of [46], the corresponding experimental ranges for k ext: , shown in Fig. 7b, are additionally indicated by dashed lines, for the observed velocities. While the left and right lines indicate the minimum and maximum value of k ext: , the middle line represents its average value. In Fig. 7b, the undercooling minima between the simulations and the experiments show a good accordance for all three velocities. The maximum deviation of 13% is found for the study Sim NiAlÀ10Mo v15 , with the smallest growth velocity. For the system Ni-35Mo, the authors could not find any corresponding experimental results in the literature. Therefore, no comparison with experimental results is shown for the binary system.
In Fig. 7c, d, the spacings k ext: are plotted against the corresponding growth velocities v, in a double logarithmic graph. In [19], a similar comparison is performed for simulations and experiments of the system NiAl-34Cr. The black curves correspond to the analytical description k ¼ ffiffiffiffiffi ffi , taken from [19]. The parameters K 1 and K 2 are material-specific values and can be calculated by following the descriptions in [23]. For both material systems, good accordances are found between the simulated and the analytically calculated relationships of v and k ext: . The largest deviation with 3.3% can be observed for the simulation series Sim NiAlÀ10Mo

v25
. The additional comparison with the average experimental values from [46] shows a larger scattering of the data points around the analytic curve. Here, the largest deviation with 12.5% is found for the experimental studies with the growth velocity 15 lm s À1 .
In Table 4, a concluding analysis of the phase fractions, adjusting in the simulations, is shown. For the system NiAl-10Mo, the evolving phase fractions of the simulations are compared with the experimental results from Zhang et al. [46]. As there are no published experimental results for Ni-35Mo, the  phase fractions from the simulations are compared with the theoretically calculated phase fractions from the phase diagram of Yaqoob et al. [43]. The theoretical values are calculated by the lever rule. Depending on the used simulation domains, the phase fractions within one study, with a constant growth velocity, differ from each other by less than 2%. In Table 4, the phase fractions of the simulation with the smallest undercoolings thus are chosen for each study, respectively. The spacings k ext: , at which these minima occur in the undercooling-spacing curves, are additionally labeled in Table 4. For all studies, it can be seen that the phase fractions in the simulations are smaller, when compared to their theoretical or experimental reference. For the comparison of the binary system Ni-35Mo, a decrease of the phase-fraction deviations is found with increasing growth velocities. For the system NiAl-10Mo, a contrary trend is observed. While for NiAl-10Mo, a good accordance is found between the simulated and the experimentally observed phase fractions, larger deviations are observed for the comparison between the simulated and the theoretical phase fractions of Ni-35Mo. For the binary system, the maximum deviation is found in the simulations Sim NiÀ35Mo v15 , with 5.1%, and for the ternary system, it is found in the simulations Sim NiAlÀ10Mo v30 , with 6.0%. As can be observed, the spacing-undercooling relationships of the simulations are found to be in good accordance with the analytically expected results, indicating that the generated Gibbs energy functions are applicable for the phase-field simulation of the directional solidification in both systems.

Discussion and conclusion
For the generation of Gibbs energy formulations, which are to be used in grand-potential-based phasefield models, an efficient and semi-automated framework is introduced in this work. Thus, the focus of this paper is rather on the generation of the Gibbs energy formulations and on the challenges occurring during this process than on the microstructure evolution of the subsequently performed phase-field simulations. First, using a least squares method and second-order polynomials, a general approach to modeling binary and ternary systems is presented. The usage of a fitting approach, instead of an approximation, solely based on the thermodynamic values at the equilibrium conditions, leads to a stable and computationally efficient formulation of the Gibbs energies. Furthermore, as the functions are valid for a certain range around the equilibrium conditions, they can be calculated prior to the simulation studies, and no recalculations of the Gibbs energy functions are required during the simulations. This approach is known in the phase-field community, but most publications lack a detailed description of the process currently being used. In addition, these general approaches often require a subsequent adjustment and optimization, in order to meet the conditions of the observed systems. For this purpose, the framework also includes an optional extension for material systems, in which the results of the common approach turn out to be unsatisfactory. Within this extension, the generated functions are modified to satisfy the equilibrium conditions in the observed material systems. By taking the equilibrium conditions into account, the deviations resulting from the fitting process of the general approach are reduced, which leads to the generation of stable and robust Gibbs energy formulations. In order to adjust the Gibbs energy functions during this extended modeling process, the description of the criteria allows the optimization and fitting procedure to be transferred to other approximation methods.
By investigating the approximation process for the two different material systems Ni-35Mo and NiAl-10Mo, it is shown that the framework is suitable for the reproduction of the commonly known approach, which generates computationally efficient Gibbs energy functions from thermodynamic CALPHAD databases. For the ternary system NiAl-10Mo, the generated Gibbs energy functions from the common approach lead to unphysical results, whereas the generated Gibbs energy formulations from the extended approach are in good agreement with the thermodynamic information from the CALPHAD database. For the binary system Ni-35Mo, the common approach already produces satisfying results; therefore, the extended approach has not been applied.
For both material systems, the subsequently performed phase-field simulations of the directional solidification process demonstrate the suitability of the generated formulations. Despite the fact that the approximated Gibbs energy formulations of Ni-35Mo are in good agreement with the thermodynamic data, the phase fractions in the simulations differ from the expected values. In contrast, it can be said that the system NiAl-10Mo generally shows smaller deviations in the phase fractions between the simulations and the experiments. In the simulations of the system Ni-35Mo, the larger deviations of the phase fractions indicate that the calculated equilibria between the Gibbs energy functions differ from the expected equilibrium states of the phase diagram. While the expected phase fractions from the phase diagram of ), the undercooling temperature differs significantly from the considered middle temperature. Due to the data used for the modeling, the equilibrium conditions for these temperatures are not well represented by the generated functions, which results in larger deviations of the phase fractions. By using the extended approximation approach for the modeling of NiAl-10Mo, the equilibrium conditions of other undercooling temperatures are also taken into account, when generating the Gibbs energy functions. This leads to a maximum difference of 0.7 %, between the values of the phase fractions from the simulations of NiAl-10Mo, and also to a better agreement with the measured values from the experiments.
On the one hand, this shows the importance of the phase-field simulations, when validating the approximation process, and on the other hand, this also shows the benefits of the presented extension. However, the deviations of the phase fractions of Ni-35Mo are still in an acceptable range. Hence, the generated Gibbs energy functions can be used for further investigations, performed within phase-field simulations.
Despite the fact that the framework is limited to binary and ternary material systems without stoichiometric phases, the validation results show that the presented approaches can be used to generate computationally efficient and accurate Gibbs energy formulations, which then can be applied in directional solidification simulations, using the phase-field method. By using these approximated functions, the thermodynamic conditions, occurring during phasefield simulations, can be reproduced in a certain range around their equilibrium states. Thus, the simulation with different off-eutectic melt compositions or with other disturbances of the equilibrium states are viable. Furthermore, by using these approximated functions, the runtime of the simulations can be reduced significantly, as shown in ''Appendix.''

Acknowledgements
Early concepts have been funded within the BMBF project SKAMPY (Ultra-scalable multiphysics simulations for solidification processes in metals), the cooperative graduate school ''Gefü geanalyse und Prozessbewertung'' of the Ministry of Baden-Wü rttemberg and the Helmholtz graduate school ''Integrated Materials Development for Novel High Temperature Alloys''. The support is gratefully acknowledged. In continuation, the completion of the data workflow has been achieved through the DFG project PAK 983/1, with the grant number MU 959/48 and the title ''Phase-field simulation and experimental microstructure research''. Further, we are grateful for the provided computational resources on the bwUniCluster at KIT. The general methods to compute computationally efficient expressions for Gibbs energy representations are funded through the Science Data Center ''MoMaF'' of the state of Baden-Wü rttemberg.

Funding
Open Access funding enabled and organized by Projekt DEAL.

Declarations
Conflict of interest The authors declare that there are no financial interests or personal relationships to affect the reported scientific results of this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licen ses/by/4.0/.

Appendix 1: Computational effort
In order to determine the computational effort for a free energy calculation, Eq. (5) from section 2 needs to be expanded. For the sake of simplicity, it is assumed that only a single sublattice is present. If there were multiple sublattices, the estimated computational effort for the CALPHAD method would rise substantially. The excess Gibbs energy of mixing is expanded with a Redlich-Kister polynomial [51]. This yields the following expression for the Gibbs free energy of a single phase: The grand potential formulation of the phase-field method does not inherently require the Gibbs energy of the entire system, but only the Gibbs energy of each order parameter, corresponding to a phase. This is due to the assumption of a local equilibrium between the present phases. Hence, the analysis can be restricted to the Gibbs energy of a single phase. The mole fractions of the K individual components are given by the variables x i . To evaluate a Gibbs energy value at a given temperature and mole fraction, the calculation effort P can be written as: with P i as the effort of one specific operation of N operations and n i as the number of times this operation occurs in the Gibbs energy formulation. Equation (8) contains additions, subtractions, multiplications and logarithms, which form the set of considered operations with their efforts P i . Counting these operations for a regular solution (k ¼ 0) yields 2K þ ðK 2 þ KÞ=2 þ 1 multiplications, 2K þ ðK 2 þ KÞ=2 þ 2 additions and K logarithms. Subregular solutions (k ¼ 1) additionally require ðK 2 þ KÞ multiplications, ðK 2 þ KÞ=2 additions and ðK 2 þ KÞ=2 subtractions. The term g 0 i is a temperature-dependent expression, containing polynomial and logarithmic terms. Its formulation is highly dependent on the considered phase. As this paper studies the alloy systems Ni-Al-Mo and Ni-35Mo, the contribution of Mo to the Mo-rich phase with a BCC crystal structure, taken from [45], is used as an example for a g 0 : In Eq. (10), six additions, seven multiplications, one division, three powers and one logarithm have to be calculated. The corresponding phase description from [43] includes the same number of operations. As the powers are known in advance, they can be translated into multiplications, instead of using the more expensive power function. Hence, the three powers can be replaced by seven additional multiplications. Each of the K components has a contribution to the mechanical mixing energy, and thus, its evaluation requires 6K additions, 14K multiplications, K divisions and K logarithms, assuming that the calculation is similarly expensive in other components. Compared to this, a parabolic approach of the form requires ðK þ 1Þ 2 þ ðK þ 1Þ þ K multiplications and ðK þ 1Þ 2 þ ðK þ 1Þ =2 þ K þ 2 additions. The coefficients A i;j ðTÞ; B i ðTÞ and C(T) are assumed to be linear, as a compromise between accuracy and efficiency. Thus, their calculation involves ðK 2 þ KÞ=2 þ K þ 1 additions and multiplications. Determining the computational effort P i requires knowledge of the specific computer architecture. As an estimate for the effort required for the calculation time, we can use the reported latency as a specific instruction, as this represents the time required for doing the operation in series. Here, Agner Fog's [52] instruction table for the Haswell architecture is used as a representation of a modern CPU. For additions as well as subtractions, the Haswell architecture has a latency of 3 cycles and 5 cycles for multiplications. Division takes between 10 and 24 cycles. For a lower limit, which is in favor of the CALPHAD approach, 10 cycles are used. The natural logarithm is calculated indirectly and has a latency range of 58 to 630 cycles. Again, the lowest number of cycles is used.
Based on these considerations, Fig. 8 shows the number of cycles, required for the considered approaches, as a function of the components in the solution. For a ternary system, the parabolic approximation takes about 10 times fewer cycles, implying a possible speedup of 10. For the generation of approximated Gibbs energy formulations, this shows the benefit of the introduced framework. By using parabolic approaches, the required computational effort of the presented phase-field simulations could be reduced dramatically.

Appendix 2: Fitted Gibbs energy functions and simulation parameters for Ni-35Mo
See Tables 5 and 6.  As an exemplary temperature below the melting point, the temperature at 1590 K is used to calculate and compare the amounts of the Gibbs energies (g) and the chemical potentials (l) for different phases. app: approximated, cal: CALPHAD, eq: equilibrium and simulation parameters for NiAl-10Mo See Tables 7, 8 and 9.  As an exemplary temperature below the melting point, the temperature at 1871 K is used to calculate and compare the amounts of the Gibbs energies (g) and the chemical potentials (l) for different phases . app: approximated, cal: CALPHAD, eq: equilibrium The exemplary temperature of 1871 K, below the melting point, is utilized to calculate and compare the amounts of the Gibbs energies (g) and the chemical potentials (l) for different phases. mod: modified, cal: CALPHAD, eq: equilibrium