Phase-Field Simulation of the Microstructure Evolution in the Eutectic Alloy NiAl-31Cr-3Mo

: The directionally solidiﬁed eutectic alloy NiAl-(Cr,Mo) is a promising candidate for structural applications at high temperatures, due to its increased creep resistance compared to its single phase B2 ordered NiAl counterpart. This system yields an eutectic trough connecting the invariant reactions of the ternary alloys NiAl-Cr and NiAl-Mo . During directional solidiﬁcation (DS) along this trough the evolved microstructures of the two-phase eutectic is changing from ﬁbrous to lamellar and back to ﬁbrous morphology while increasing and decreasing the amounts of Mo and Cr , respectively. To investigate these effects in the morphology, the phase-ﬁeld method has proven to be predestined in the last decades. However, as the modeling of quaternary systems is challenging for the simulation with a grand potential based phase-ﬁeld model, the focus of this work is on the generation of a material model for one deﬁned compound namely NiAl-31Cr-3Mo . The modeling is validated by investigating the microstructure evolution in two-and three-dimensional simulations of the DS process for two different growth velocities and by investigating their undercooling spacing relationships. The evolving microstructures obtained from three-dimensional large-scale simulations are presented and validated with corresponding micrographs from scanning electron microscopy (SEM) of directionally solidiﬁed samples with the same growth velocities. The simulation results show the theoretically expected behaviors and are in qualitative and quantitative accordance with DS experiments. The study of NiAl-31Cr-3Mo serves as the basis for a comprehensive data-driven analysis of microstructure properties and system quantities of the entire quaternary material NiAl-(Cr,Mo) . With this, an accelerated design of advanced materials is promoted.


Introduction
Over the last decades, simulations have established as a third essential tool in material science investigations, next to experiments and theoretical works.Especially when examining the solidification process, simulation methods have proven to be helpful tools.A big advantages of this technique is that intermediate states of the solidification can be recorded and examined just like the end state of the material.With this the evolution of the microstucture within the material can be investigated during the entire solidification process.Further, simulations enable a specific adjustment of the material and the process conditions.Through the advances in the computer technologies in the recent years [1], large-scale simulations of three-dimensional structures allow to study the evolution of microstructure formation phenomena in statistical volume elements on high performance clusters [2].
However, modern materials for industrial applications mainly consist of more than two or three components.Unfortunately, phase-field simulation studies of higher ordered materials are missing in literature.Therefore, this work is intended to be a further contribution to the investigation of multi-component systems using the phase field method.In a first extension to previous works, a new approach to model a quaternary system based on thermodynamic CALPHAD databases is introduced.To show the validity of the approach, it is executed for a eutectic composition of the real material system Ni-Al-Cr-Mo, which is a potential candidate for high temperature applications.Through the directional solidification at this composition, namely NiAl-(Cr,Mo), a two-phase eutectic pattern is evolving in the microstructure.In a first study, two-dimensional simulations of the directional solidification process are performed to investigate according to the theory of Jackson and Hunt [16] the relationship between the solidification velocity, the evolving lamellar spacings of the microstructure and the adjusting undercoolings at the solidification front.Next, three-dimensional large-scale simulations are performed to compare the adjusting microstructure patterns with respective experimental results, which have been also performed for this work.Finally, the conclusions of the simulation results and the comparison with theory and experiment are drawn.Based on these conclusions the presented modelling approach for eutectic quaternary systems is validated.

The Quaternary System NiAl-(Cr,Mo)
The intermetallic compound NiAl is of particular interest for engineering applications at elevated temperatures, such as turbine blades in jet engines or stationary gas turbines [17,18].In addition to its high melting point (T m = 1911 K) and thermal conductivity (>70 W K m), NiAl consists a low density (5.95 g cm 3 ) and excellent oxidation resistance.However, apart from these advantages, NiAl exhibits only low fracture toughness and ductility at room temperature as well as insufficient strength at high temperatures and poor creep resistance [19].One possible way to overcome these disadvantages is by alloying refractory metals such as Cr, Mo, Re, V and W [20].With this, a second reinforcing phase is introduced which leads to the formation of two-phase eutectic systems.During the directional solidification (DS) of these eutectic systems, different fibrous and lamellar structures of the reinforcing solid solution phase (disordered A2 ) embedded in a NiAl -matrix (ordered B2 ) evolve.At this point it should be mentioned, that the components, phases and systems are highlighted in the text by a serif free font.To distinguish between the compound NiAl and the phase NiAl , the phases are written in italic.The systems are indicated by a "-" between the components.The most promising high-temperature structural materials are the eutectic composites of the systems NiAl-Cr and NiAl-Mo.These systems exhibit also a greater improvement in the overall mechanical properties compared to the systems NiAl-W and NiAl-Re [21], with only small volume fractions of the reinforcing phase.Moreover, the here chosen systems reveal a better creep resistance compared to NiAl-V [22].Due to this, several experimental studies of the directionally solidified ternary alloys NiAl-Mo [23][24][25][26][27][28] and NiAl-Cr [28][29][30][31] have been conducted in the past.
Accordingly, the two-phase alloys of the quaternary NiAl-(Cr,Mo) system are regarded as the most promising candidates for the usage in high-temperature applications since they might possess a superior mix of properties compared to the two-phase alloys of NiAl-Cr and NiAl-Mo [31,32].Peng et al. [33] investigate the thermodynamic properties of the quaternary system and show the existence of a eutectic trough connecting the invariant reactions of the two isopleth cuts NiAl-Mo and NiAl-Cr.This trough can be seen in the three-dimensional phase-diagram of the pseudo-ternary system NiAl-Cr-Mo, see Figure 1a.The depicted phase diagram is calculated based on the CALPHAD database of Peng et al. [33] by treating Ni and Al as a single component NiAl.In Figure 1b, the temperature profile of the eutectic trough (solid line) is plotted over the Mo amount x Mo of the reinforcing phase (Cr,Mo).Similar to the invariant reactions in the isopleth cuts NiAl-Mo and NiAl-Cr, a three-phase region Liq + NiAl + (Cr,Mo) exists beneath the entire eutectic trough [8,34].The solidus line of the three-phase region is plotted in dotted lines.Figure 1c shows the eutectic trough (green line) in the liquidus projection of the pseudo-ternary system NiAl-Cr-Mo surrounded by the phase diagrams of the isopleth cuts of the systems NiAl-Cr and NiAl-Mo and of the binary system Cr-Mo.The orange and blue lines in the liquidus projection represent the equilibrium concentrations of the NiAl and the (Cr,Mo) phase calculated for concentrations along the eutectic trough and an undercooling of 2 K beneath the solidus line, respectively.
The microstructure evolution during the directional solidification reaction Liq ⇆ NiAl + (Cr,Mo) along the eutectic trough has been experimentally investigated in [21,35].Depending on the compositions of NiAl-xCr-yMo with x = 0 . . .34 at.% and y = 10.3 . . .0 at.%different arrangements of the reinforcing phase A2 (Cr,Mo) evolve embedded in the B2 NiAl -matrix.During solidification of the binary eutectic reaction of the isopleth alloy system NiAl-Cr, circular Cr-rich fibers grow in a hexagonal arrangement, embedded in a NiAl -matrix.The first transition to lamellae takes place when decreasing the Cr concentration.A lamellar structure is stable for Cr contents from 33 to 10.3 at.%.A further decrease in Cr below 10.2 at.% and a corresponding increase of Mo up to 10.3 at.% causes a morphology transition to a mixture of fibers and lamellae.A further morphology transition along the eutectic trough was found for compositions with a Cr concentration below 1.8 at.%.For these compositions Mo-rich fibers evolve embedded in the NiAl -matrix [20].Due to the highly anisotropic interfacial energies in NiAl-10Mo, the shape of the evolved Mo-rich fibers in transverse section is rather square than circular [25].Subsequent to DS the phases arrange with respect to growth direction and an increase in mechanical properties can be achieved [23,25].Further experimental studies of the quaternary system NiAl-(Cr,Mo) have been conducted for example by Whittenberger et al. [36,37] and Chen et al. [38] Additionally to these experimental investigations, numerical phase-field simulation studies for the ternary systems at their eutectic reactions NiAl-10Mo [8] and NiAl-34Cr [12,13,34,39,40] have been performed.A quantitative comparison between experimental samples and large-scale phase-field simulations of the isopleth alloy system NiAl-Cr is reported in [34] for different solidification velocities.The influence of varying solidification velocities on the adjustment process of the microstructures is investigated in [41].For the opposing isopleth alloy system NiAl-Mo the growth of Mo-rich fibers of rectangular cross section surrounded by the NiAl -matrix phase is simulated in [8].
However, phase-field simulations of the quaternary system NiAl-(Cr,Mo) are missing in the literature.Therefore, the microstructure evolution during the directional solidification process of a defined composition at the eutectic trough will be simulated in this work.As an example composition the eutectic NiAl-31Cr-3Mo system is chosen for the conducted simulation studies.The positions of this composition and of the corresponding equilibrium concentrations of the NiAl and the (Cr,Mo) phase are additionally marked in the liquidus projection of Figure 1c by corresponding colored points.At this composition the threephase region, depicted in Figure 1b has only a width of approximately 10 K, which is beneficial for the modeling of the material system for the subsequently planned phase-field simulations.Experimental micrographs of NiAl-31Cr-3Mo [36] show the adjustment of a lamellar microstructure, although based on [42] a fibrous microstructure is expected for the obtained reinforcing phase fraction of around 1 3. Also first pre-studies employing phase-field simulations indicate the growth of a merely fibrous microstructure for the composition NiAl-31Cr-3Mo.Following [43,44], the growth of a lamellar microstructure could be attributed to the anisotropic properties of the surface energies.This is another point which makes the directional solidification process of NiAl-31Cr-3Mo particularly interesting for a simulation study with the phase-field method.

Phase-Field Method
For the simulations of diffusion driven phase transitions in the quaternary system Ni-Al-Cr-Mo, a thermodynamically consistent phase-field model, based on a grand potential functional and an Allen-Cahn-type variation, is used [9,45,46].For the conditions of consideration, the system consists of three phases and four components.The N = 3 order parameters φ α, describe the local phase fractions of the three phases participating in the eutectic reaction Liq ⇆ NiAl + (Cr,Mo) in NiAl-31Cr-3Mo.To differentiate the phases α, β, . . .from their indices, the indices are marked by a ◻ symbol.The development of the N phase-fields φ α are described by the time evolution equations To couple the different timescales of the evolution equations, the relaxation parameter τ is introduced [46] and to control the inherent width of the diffuse interface of the phasefield model, the parameter ε is used.The right hand side of Equation ( 1) can be divided into three parts.In the first part rhs 1,α , the shape of the diffuse interface between the phases is modeled by two superimposing interfacial energy density contributions: the gradient energy a(φ, ∇φ) and the obstacle potential ω(φ).The formulation of the gradient energy density is expressed as where γ α β denotes the interfacial energy parameter and q α β = φ α∇φ β − φ β∇φ α defines the normal vector with q = q α β q α β as the corresponding unit vector of the interface [43].The function a c ( q) is describing the anisotropy in the interfacial energies between the phases.Following [43,44], an anisotropy formulation of the form: with qi as the Cartesian component of the unit normal vector q in the direction i, is used, to model the growth of a lamellar structure within the three-dimensional microstructures.Equation ( 3) generates a two-fold anisotropy in the solid-solid interface normals, that are perpendicular to the growth direction y [44].This leads to an anisotropy of strength δ α β in the xz plane and ± δ α β 2 in the xy and yz planes [43], respectively.Using this most simple anisotropy formulation proposed in [44], the main focus of this work remains on the investigation and validation of the generated quaternary model for NiAl-31Cr-3Mo instead of an intense investigation of the anisotropic behavior in the system.For the same reason, the influence of the elastic strain energy between two phases and the influence of the growth orientations on the phase morphology will not be considered in this work.
The formulation of the obstacle potentials ω(φ) is given in [9] and includes next to the interfacial energies γ α β the higher order term γ αβδ .This term is introduced to reflect the correct equilibrium angle conditions at the triple junctions [9,47], by suppressing the appearance of an additional phase at the phase boundaries.
The second part of the right hand side rhs 2,α from Equation ( 1) is describing the driving force for the phase transitions.This force is defined by the differences of the grand potentials ψ β, which are stored in the grand potential function ψ(φ, µ, T).By using of the projection of the thermodynamic energies instead of the energies itself, the interface energies and the interface thickness are decoupled in the grand potential approach.This is an advantage compared to free energy models.To ensure the thermodynamic consistency of the modeled system, the grand potentials are derived from the Gibbs energies of the different phases [34], which are incorporated from thermodynamic CALPHAD databases.To optimize the computational effort, the formulations of the Gibbs energies are approximated by a parabolic approach of the form: with the temperature dependent coefficients A i α(T), B j α(T) and C α(T) and the bulk concentration vector c.As shown in [8], using these approximated parabolic approaches can reduce the computational effort by more than an order of magnitude compared to the original CALPHAD formulations.
The third part of the right hand side from Equation ( 1) includes the Lagrange multiplier Λ, which is introduced to fulfill the constraint ∑ N α=1 ∂φ α ∂t = 0. To ensure mass conservation during the simulations, the time evolution of the chemical potentials is formulated and coupled to the evolution equation of the phase fields.The chemical potential vector µ consists of a parameter µ i for each component (i = Ni, Al, Cr, Mo) and its time evolution is derived from the mass balance of the K = 4 concentrations and from Fick's law as Thereby, the vectors c α(µ, T) contain the concentrations 1 c, 2 c, . . ., K c of the K components for the regarded phase α and the mobility term M(φ, µ, T) [46] incorporates the diffusion coefficient matrix D of the involved phases.The function h α is introduced in the form h α = φ 2 α (3 − 2 φ α) to interpolate between the different phases and an anti-trapping current J at [46,48,49] is applied to balance the effects of artificially enlarged interfaces from the phase-field modeling.
Due to the multiple times faster diffusion of heat compared to mass, the evolution of the temperature is described by a moving analytic approach of the form Starting from an initially imprinted temperature field, with the base temperature T 0 , the temperature T evolves, with the gradient ∇T and the velocity v ∇T in the growth direction y.In 3D simulations the growth direction is changed from y to z.

Analysis Method
The subsequently performed two-dimensional simulations of the system NiAl-31Cr-3Mo are validated by investigating the relation of the undercooling ∆T, the growth velocity at the solidification front v F and the spacing λ of the evolving lamellae pairs.This relation was derived by Jackson and Hunt for binary eutectic systems in [16] and is used in several publications [50][51][52][53][54] to validate phase-field simulations.For this reason, a short summary of the foundations of this theory should be presented.A more detailed discussion of the Jackson-Hunt-relationship and its suitability for validating phase-field simulations of rod structures is given in [55].
Based on the Gibbs-Thomson equations the undercooling ∆T at the solidification front can be given as: ∆T = ∆T con + ∆T cur + ∆T kin (7) with ∆T con describing the influence of the concentration, ∆T cur represents the influence of the solid-liquid interface curvatures and ∆T kin refers to kinetic effects.Following Equation ( 7) and the argumentation from Jackson and Hunt, the analytical relation of the undercooling ∆T, the growth front velocity v F and the lamellar spacing λ, can be formulated as: The parameters A, B and C are constant but specific for each observed alloy system.The shape of Equation ( 8) is describing the so called Jackson-Hunt curve and the spacing related to the minimum of the curve is labeled as λ ext .This extreme position refers to the minimal undercooling at the solidification front and follows the relationship Jackson and Hunt's theory, originally developed for binary systems, was extended by McCartney et al. [56] to three component, two-phase eutectics and by Choudhury et al. [57] to three component, three-phase eutectic growth in 2D.

Model
In this section, the modeling of the Gibbs energy formulations for the phases NiAl , (Cr,Mo) and Liq is described.The Gibbs energy formulations are derived based on the thermodynamic CALPHAD database of Peng et al. [33] for the full quaternary system Ni-Al-Cr-Mo and they are approximated by reproducing the workflow of the semi-automated framework from Noubary et al. [8].Therefore, first the equilibrium concentrations of the solid phases are calculated based on the thermodynamic data [33] for temperatures of T 1 = 1700 K and T 2 = 1710 K, respectively.For the liquid, the concentration of the eutectic reaction at approximately 31 at.%Cr is used.At T 2 , the CALPHAD database shows a three phase domain consisting all mentioned phases beneath the eutectic reaction (see Figure 1c).As the modeling of this three-phase region would require the implementation of additional phases, it should not be depicted in the modeled system.Thus the calculations at T 2 have been made by only regarding the two solid phases.The results of the equilibrium calculations are summarized in Table 1.
Table 1.Equilibrium concentrations of NiAl and (Cr,Mo) for T 1 = 1700 K and T 2 = 1710 K, respectively, and concentration of the eutectic reaction at the liquidus temperature T liq.= 1716.56K calculated based on [33].All entries are rounded after the second decimal place.

Temp.
Al The later derived Gibbs energy formulations are valid for concentrations in the vicinity of these pre-calculated equilibrium conditions.In a second step, the Gibbs energies for all three phases are calculated in a mesh in the vicinity of their equilibrium concentrations for the temperatures T 1 and T 2 .Based on these data, the formulations for the Gibbs energies are approximated for each temperature with a least squares method and second order polynomials, as introduced in Equation (4).Next, small adaptions of the parameters B j α(T) and C α(T) of Equation (4) are made to improve the agreement of the equilibrium concentrations.After these adjustments, the average difference between the modified formulations and the originally generated formulations is less than 0.05%.The maximum deviation with 0.15% is found for the phase (Cr,Mo) at a temperature of 1700 K.In a final step, the two corresponding functions are combined via a linear interpolation to generate a unique concentration and temperature dependent formulation for each phase.The resultant approximated formulations of the Gibbs energies are given in Equations ( 10)- (12) for the dimensionless unit J sim cells 3 .The formulations for the quaternary alloy system NiAl-31Cr-3Mo are modeled for the full quaternary system Ni-Al-Cr-Mo and depend on the Within these formulations, the comparatively high values in Equation (11) for the phase NiAl are notable, which describe the quadratic dependency on Mo c.These high values are a result of the relatively small amount of Mo in the phase, which increases the complexity of the modeling.Similar to the modeling of stoichiometric phases [14,58], a strong curvature is applied to the polynomial in the direction of the Mo concentration to prevent a wide variance of Mo c.With this, the stability of the thermodynamic system is ensured.
Due to the simplifications made in the description of the thermodynamic energies, their accuracy has to be checked.Following the validation from Noubary et al. [15], the subplots in Figure 2 show the equilibrium concentrations calculated from the CALPHAD database (solid lines) and from the Equations ( 10)-( 12) (symbols), respectively.For all four concentrations an excellent accordance can be found not only in the predefined temperature range between T 1 and T 2 (gray area) but also up to eutectic temperature T eut , marked with a red dotted line, and also down to 15 K below T 1 .It can be expected that also for smaller temperatures a good accordance can be achieved.However, as such low temperatures are not expected to appear within the simulations a further validation of the equilibrium concentrations is omitted.Within the investigated temperature range from 1685 K to 1716 K the average deviation between the equilibrium concentrations of the solids calculated from the derived Gibbs energy functions and the equilibria calculated from the CALPHAD database is ∼0.1%.The maximum deviation of ∼0.9% is found for the Mo concentration of the equilibrium concentration of NiAl at the lowest temperature of 1685 K. Despite the simplifications made in the modeling of the thermodynamic energies, this good agreement shows the suitability of the generated functions for representing the thermodynamic system and thus the applicability of the approach presented for their approximation.For the upcoming computations within the multi-physics phase-field framework PACE3D version 2.5 [59,60], the concentration of Ni is replaced by the constraint Ni c = 1 − Al c − Cr c − Mo c.With this, the Equations ( 10)-( 12) can be rearranged, which leads to a reduction of the concentrations vectors c α to K − 1 independent entries and hence to an increase of the computational efficiency of the model.All valid system and numerical parameters are summarized in Table A1 in the Appendix A in dimensionless and physical units.

Experimental Setup
The desired eutectic composition of NiAl-31Cr-3Mo (at.%) was synthesized by using the elements Ni, Al, Cr and Mo of 99.97, 99.9, 99.96 and 99% purity, respectively.The alloy was manufactured by using an arc-melting device provided by Edmund Bühler GmbH (Germany).The device chamber was evacuated multiple times and purged with Ar to ensure an oxygen-lean operation atmosphere.Furthermore, a lump of Zr metal was melted prior to the melting process, in order to remove the residual oxygen in the chamber.To achieve sufficient homogeneity, the produced buttons were flipped and re-melted at least five times.For the directional solidification the buttons were cast into a rod-shaped Cu mold.The rods were 12 mm and 170 mm in diameter and length, respectively.No specific loss of a certain element during arc-melting was noticed and an overall mass loss lower than 0.3 wt.% was measured.Therefore, the changes in alloy composition by evaporation are considered negligible.The worldwide unique crucible-free zone melting device used in this work is described elsewhere in detail [61].The experiments were performed in an Ar atmosphere with a continuous Ar flow rate of 1.5 L min.During the DS process, the sample rotates about its axis and moves upward (lift) through the induction coil that is used to heat up and locally melt the material, as shown in Figure 3b.The rotation speed R was kept constant at 10 rpm and the process velocity v textrmF varied from 42 mm h (v 42 ) to 102 mm h (v 102 ).For microstructural analysis, samples were cut by electro-discharge machining parallel and perpendicular to the growth direction (GD) and were prepared by standard metallographic procedure.Scanning electron microscopy was performed using backscattered electron imaging (SEM-BSE) on a Zeiss EVO50 microscope by Carl Zeiss AG (Germany) at 20 kV and a spot size of 500.
In order to visualize the DS process, Figure 3 shows the microstructure at room temperature before Figure 3a and subsequent to directional solidification Figure 3c.In all micrographs, dark contrast corresponds to the NiAl matrix and bright contrast to (Cr,Mo) solid solution.In the as-cast state Figure 3a, a random crystal orientation distribution of the colonies can be seen, whereas an alignment of the microstructure parallel to GD was obtained at v 42 in Figure 3c.Additionally, the heating process and the formation of the liquid zone are displayed in Figure 3b.

Simulation Setup
To investigate the microstructure evolution of NiAl-31Cr-3Mo within phase-field simulations, two-and three-dimensional setups similar to [10,34,62] are used.Starting from an initial arrangement of nuclei on one side of the simulation domain, the two solid phases (Cr,Mo) and NiAl evolve, controlled by the applied temperature gradient (Equation ( 6)) in y-direction (z-direction in 3D) and by the imposed pulling velocity v F , into a liquid.For the realization of a constant flux of the melt into the domain a Dirichlet boundary condition is applied at this liquid end of the domain.At the opposite solid end a Neumann boundary condition is used.Finally, an infinite domain, perpendicular to the solidification front, is modeled by periodic boundary conditions.
To validate the generated thermodynamic model for NiAl-31Cr-3Mo, 2D Jackson-Hunt studies are conducted with the velocities 42 mm h (v 42 ) to 102 mm h (v 102 ).For this, an initial setup with one lamellae pair of the solid phases is set at the bottom of the domain (y ≤ 20 voxel cells).While all 2D simulations have a height of 300 voxel cells (6 µm) in y-direction, the width of the domains in x-direction is varied from 28 up to 84 voxel cells, which corresponds to a physical domain width of 840 up to 2520 nm.
For the subsequently performed 3D large-scale phase-field simulations, similar to [9,10,41], an initial random Voronoi tessellation with λ ≪ λ ext is used for the initial setup of the solid phases.To ensure a microstructure growth unaffected from the boundary conditions [2], the 3D simulations are conducted in a domain with a base size perpendicular to the growth direction z of 400 × 400 voxel cells (8 µm × 8 µm) and a height of also 300 voxel cells (6 µm).Exemplary illustrations of the 2D and 3D simulation setups are given in [63] and [34], respectively.
The 2D simulations are running for 5 million time steps on 3 cores on a local computer in the institutes network and the large-scale 3D simulations are performed for 8.5 million time steps on 768 cores on the SuperMUC-NG system at LRZ Munich [64].Compared to the experiments, a three times larger thermal gradient ∇T of 60 K mm is used for all simulations.Following the theory of Jackson and Hunt [16], the evolving lamellar spacing and the undercooling at the solidification front are not affected by the thermal gradient ∇T.However, with increasing amount of the thermal gradient the required time steps to reach a convergence of the microstructure in the simulations can be reduced [41].

Results and Discussion
In the first part of this section, the results of the two-dimensional phase-field simulations are discussed.Starting from an initial setting of one lamellar pair of solid phases, the phases NiAl and (Cr,Mo) evolve depending on the investigated domain width.If the considered domain width is too large, the widths of the single lamellae start to vary periodically along the growth direction.Since no nucleation mechanism [63] is used in the simulations, the microstructure is not able to refine and hence, the oscillation can lead to an overgrowth of one of the solid phases.If the domain width is too small, a similar elimination of one lamella can occur.This leads to a stable range of the domain width from 48 to 76 cells or 960 to 1520 nm for the growth velocity v 42 and from 32 to 60 cells or 640 to 1200 nm for the velocity v 102 .In these ranges, still an oscillation of the lamellar width can occur, but with ongoing simulation time a stable growth of one lamellae pair is achieved.During these oscillations the lamellar widths adjust leading in the final microstructure to a volume fraction of the (Cr,Mo)-lamellae of 32.5 vol.% for v 42 and 32.9 vol.% for v 102 .
To compute the undercoolings of these stable simulations, first the growth heights of the solidification front is measured by detecting the solid-liquid-interface position in y-direction after 5 million time steps for each simulation.Depending on these positions y, the temperature at the solidification front is calculated by using Equation (6) and hence the undercooling at this position is determined.The results of these calculations are shown in the plots of the undercooling-spacing relationship in Figure 4 for the growth velocities v 42 and v 102 , respectively.In these studies, the spacing of the simulated lamellar structures is similar to the conducted domain widths of the simulations.Both plots show the expected curve of the undercooling with a minima λ ext at 1.09 µm for and 0.8 µm for v 102 , respectively.
For further studies and the correlation with experimental results, three-dimensional phase-field simulations are performed.The resulting microstructure after 8.5 million time steps for the growth velocity v 102 is displayed in Figure 5a.The microstructure contains blue (Cr,Mo)-lamellae within an orange NiAl -matrix.Both phases are individually illustrated in the sub-figures Figure 5b,c, respectively.The sub-figures on the right side of Figure 5 visualize the solidification front of the evolved microstructures after 2.5, 5.5 and 8.5 million time steps.Starting from a mixed structure of mainly fibers and short lamellae at 2.5 million time steps in Figure 5d, a lamellar structure with single isolated fibers has evolved after 5.5 million time steps in Figure 5e.In the last 3 million time steps, this lamellar structure is established by the consolidation of longer lamellae.The final microstructure in Figure 5f shows a well-aligned mainly lamellar structure.By continuing the simulation, it can be assumed, that the remaining elongated fibers in the upper left corner of the microstructures will also connect to a further lamella.The formation of lamellae within the simulation is mainly driven by merging events of the (Cr,Mo)-phase.Only occasionally splitting events and no overgrowth or nucleation events are detected in the microstructure.These observations are in accordance to [41], where the growth of lamellae during the deceleration of a two-phase eutectic microstructure is also driven my merging events.
The microstructure after 8.5 million time steps exhibits an average lamellar spacing of ∼1018 nm for the growth velocity v 102 .For the growth velocity v 42 an average lamellar spacing of ∼1353 nm is measured.To determine these spacings, eight cuts with a distance of fifty cells are made through the 3D simulation domains perpendicular to the aligned lamellae.Based on these cuts, the widths w of the NiAl -and (Cr,Mo)-sections are measured and their maximum, minimum and average values are stored.The sum of the widths w NiAl and w (Cr,Mo) being the lamellar spacing λ.The determined lamellar spacings of the 3D simulations are approximately 25% larger compared to the detected values λ ext at minimal undercooling conditions of the 2D simulations.These deviations are in accordance with three-dimensional phase-field simulations of the ternary system Ni-Al-Cr [34,63].
The subsequently calculated undercoolings at the solidification front for the 3D simulations are added to the plots of the undercooling-spacing relationship for the twodimensional simulations in Figure 4 with a single mark at the average lamellar spacing of the corresponding simulation.Even if the plot in Figure 4a for the growth velocity v 42 shows a visually large deviation between the single measured undercooling value of the 3D simulation and the undercooling curve of the 2D simulations, a qualitative comparison sheds light on the fact that the actual deviation is less than 0.09%.For the growth velocity v 102 the deviation is only ∼0.01%.
The solidification front microstructures after 8.5 million time steps for both velocities are presented in Figure 6 accompanied by their corresponding experimental micrographs.After this calculation time, the phase components of the microstructures have adjusted with a maximum change of less than 1% within the last one million time steps and the growth velocities as well as undercoolings have adjusted.A visual comparison of the simulation results and the experimental micrographs at the same magnification illustrates an overall agreement between the microstructures.In both structures the growth of (Cr,Mo)-lamellae embedded in a NiAl -matrix is observed.However, while the evolved lamellae in the simulated microstructures are all aligned and grow parallel to each other, the experiments perform a variety in growth and in direction of the lamellae, due to anisotropy and coupled growth of interfacial orientational variants.To further validate the agreement between the simulated and experimentally achieved microstructures, significant features like the lamellar spacing λ and the widths w of both phases are measured.For the experiments, the lamellar spacing and the respective phase widths were determined using ImageJ software.The lamellar spacing here corresponds to the distance between centers of (Cr,Mo)-lamellae.The averaged results of both measurement are summarized in Table 2.While the measured lamellar spacings for the growth velocity v 102 are in good agreement with a deviation of approximately 2%, for v 42 a deviation of approximately 9% is noticed.A reason for this deviation can be found in the limited number of the evolving lamellae in the simulation.Evaluating the average number of lamellae pairs within the simulation domain, results in 5.9 lamellae pairs for the simulation and 6.4 for a corresponding experimental domain.Hence, this lead to the conclusion that the spacings just deviate by a half lamellae pair.As discussed in [63] for the growth of two-dimensional lamellar growth, several stable settings with different numbers of lamellae pairs can occur in the same simulation domain, if the domain is larger than three times λ ext .A tendency to larger structures with fewer pairs of lamellae was additionally reported in [63], which also explains the larger spacings in 3D simulations compared to the detected values λ ext from 2D simulations.As the three-dimensional microstructures already adjust to a stable configuration, the evolution of an additional lamellae pair is not expected in the simulations until the microstructure is forced to rearrange by a disturbance, for example by a variation of the growth velocity [41].However, a comparison of the extreme points λ ext from the two-dimensional studies of the undercooling-spacing relationship exhibits that the growth of an additional pair of lamellae in the 3D simulations still leads to a stable microstructure.Due to the periodic boundary conditions in the simulation domain, only an integer number of lamellae pairs can occur.To establish a better comparison to the experimental results, an increase of the simulation domain is necessary.For the sake of integrity, the measured lamellar spacings from the experiments are also added to the plots in Figure 4.As the exact undercooling at the solidification front is not known for the experiments, their lamellar spacings are labeled with a dotted and dashed line, respectively.For these data points also a good accordance with extreme points λ ext is documented.By considering the evolved widths w of the phases, the largest deviation of ∼14% is found for the NiAl -matrix and the velocity v 102 .For the velocity v 42 , the deviations for the matrix are only 3%.The average measured values for the widths of the (Cr,Mo)-lamellae are in good agreement, especially by considering the measurement scattering in the experiments.

Phase-field simulation Experiment
Both, experiments and simulations show a slight increase in the volume fractions of the (Cr,Mo)-lamellae with increasing velocities.The volume fractions of approximately 33 vol.% for the (Cr,Mo)-lamellae in the 3D simulations are in accordance with the results of the previously discussed 2D simulations.However, the simulated volume fractions are nine to ten percent smaller compared to the experiments.This deviation in the volume fractions can be ascribed to the used thermodynamic CALPHAD database [33], since remarks in [35] state that an adjustment of the solute partitioning is needed in this database.The theoretical value of the volume fraction calculated from [33] is 34.5 vol.%.As the simulations are modeled based on the thermodynamic data, smaller volume fractions compared to the experiments are expected.The deviation of the simulated volume fractions from the thermodynamic value is between four and fife percent depending on the applied growth velocity, which is in the expected range of accuracy.

Conclusions
In this work, an extended phase-field approach for the simulation of microstructures in quaternary material systems based on a grand potential functional formulation is introduced.The Gibbs energy functions of the example system NiAl-31Cr-3Mo are modeled on the basis of thermodynamic CALPHAD databases and the equilibrium conditions of the system.Subsequently, the directional solidification process for this material system is simulated in 2D and 3D with two different velocities using the multi-physics phase-field framework PACE3D.The 2D simulations are analyzed by deriving their undercoolingspacing relationships and the 3D simulations are validated by comparison of the resulting microstructures with corresponding experimental micrographs.Based in these comparisons, the following observation are made: (i) The generated Gibbs energy functions for the modeling of the quaternary system NiAl-31Cr-3Mo show an excellent behavior in reproducing the equilibrium concentrations of the solids, given in the thermodynamic CALPHAD database.Within the defined temperature range for the modeling, the deviation between the CALPHAD data and the rebuilt equilibrium concentrations is ∼0.1%.Also for temperatures up to eutectic temperature T eut.and down to 15 K below predefined temperature range a very good accordance with a maximum deviation of ∼0.9% is detected.(ii) The undercooling-spacing relationships of the two-dimensional phase-field simulations follow the expected curve for temperature dependent solidification described in the theory of Jackson and Hunt.For both velocities a defined minimum λ ext is found.(iii) The establishment of a lamellar microstructure in 3D can be ascribed to the used twofold anisotropy formulation, as simulations without anisotropy lead to the growth of a merely fibrous structure.Although this anisotropy formulation is not motivated based on the investigated material system, the results demonstrate the capability of the approach to validate the derived model of NiAl-31Cr-3Mo to describe the microstructure evolution during directional solidification.To increase the variety within the microstructure with the formation of eutectic two-phase cell boundaries, a more general interfacial anisotropy formulation can be used including crystallographic orientations in 3D to enable the description of eutectic colony boundaries.(iv) The evolved microstructures of the 3D large-scale phase-field simulations for the quaternary system NiAl-31Cr-3Mo are in qualitative and quantitative accordance with experimental micrographs.The deviations in the lamellar spacings are ∼9% for v 42 and ∼2% for v 102 , respectively.To reduce the deviations and the impact of a single missing lamellae pair, larger simulation domains are necessary.Furthermore disturbances need to be implemented into the solidification process together with the activation of a nucleation mechanism.
Based on these findings, it can be concluded, that the presented phase-field approach is suitable to model the eutectic phase transition of quaternary systems based on thermodynamic driving forces.Further, the generated thermodynamic model for the system NiAl-31Cr-3Mo has proven to be qualified to investigate the directional solidification process of the considered quaternary system NiAl-31Cr-3Mo, based on a validation by reproducing the experimentally achieved microstructures.
Table A1.Summary of numerical and material parameters in dimensionless simulation parameters and in their physical units for the quaternary system NiAl-31Cr-3Mo.

xFigure 1 .
Figure 1.Phase diagrams of NiAl-Cr-Mo [33]: (a) Three-dimensional phase diagram of the pseudoternary system NiAl-Cr-Mo.(b) Temperature profiles of three phase region beneath eutectic trough.(c) Liquidus projection of the pseudo-ternary system NiAl-Cr-Mo (green) showing the eutectic trough and the subsequently calculated equilibrium concentrations of the NiAl (orange) and the (Cr,Mo) (blue) phase, surrounded by the phase diagrams of the systems NiAl-Cr, Cr-Mo and NiAl-Mo.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Figure 2 .
Figure 2. Comparison of equilibrium concentrations calculated from the CALPHAD database (solid lines) and from the Equations (10)-(12) (symbols) for the phases NiAl (orange) and (Cr,Mo) (blue).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Figure 3 .
Figure 3. SEM-BSE micrographs of NiAl-31Cr-3Mo in (a) as-cast and (c) as-DS condition, micrographs parallel to growth direction (GD indicated by arrow) subsequent to directional solidification at v 42 = 42 mm h.(b) Photograph of the zone melting process during operation.

vFigure 4 .Figure 5 .
Figure 4. Undercooling-spacing relationships from two-dimensional phase-field simulations of the eutectic NiAl-31Cr-3Mo system for the growth front velocities v 42 (a) and v 102 (b).The measured undercooling from the 3D simulations is added with a single rectangular mark, respectively.

Figure 6 .
Figure 6.Solidification front from 3D phase-field simulations and from corresponding experiments with the growth front velocities v 42 = 102 mm h (a,b) and v 102 = 42 mm h (c,d).

Table 2 .
Measured lamellar spacings and widths as well as phase fractions from 3D simulations and experiments.