Dual-scale phase-field simulation of Mg-Al alloy solidification

Phase-field simulations of the nucleation and growth of primary α-Mg phase as well as secondary, β-phase of a Mg-Al alloy are presented. The nucleation model for α- and β-Mg phases is based on the “free growth model” by Greer et al.. After the α-Mg phase solidification we study a divorced eutectic growth of α- and β-Mg phases in a zoomed in melt channel between α-phase dendrites. The simulated cooling curves and final microstructures of α-grains are compared with experiments. In order to further enhance the resolution of the interdendritic region a high-performance computing approach has been used allowing significant simulation speed gain when using supercomputing facilities.


Introduction
Magnesium alloys are among the most promising structural materials for lightweight applications due to their low density compared to aluminum and steel. In Mg technology the most applied alloying element for technical alloys like AZ91, AM50 (cast) or AZ31 (wrought) is aluminum. The mechanical properties of such Mg-Al alloys are primarily determined by their micro-structure which typically consists of two main phases: the Mg-rich hexagonal close packed (HCP) α-phase and a near stoichiometric Mg 17 Al 12 β-phase, and further minor secondary phases due to third alloying elements. From a scientific point of view it remains a challenge to simulate the dendritic solidification structure of the primary HCP-α-phase. In particular the anisotropy function of the interface energy between liquid and the α-phase which determines the growth structure of the dendrite, remains under debate [1,2]. A second question is the morphology of interdendritic eutectic, which is decisive for the corrosion properties of the alloy [3]. In the present work we apply the phase-field method [4,5,6,7] to simulate the microstructure evolution of Mg-Al alloys during solidification. We apply the anisotropy function as given in [8] which produces a dendrite trunk normal to the basal plane of the HCP crystal. While several authors reported 2-dimensional and 3-dimensional phase-field simulations eg. [8,9,2] of Mg-alloy solidification, a new challenge is to reveal the eutectic structure, which is assumed to happen in a divorced mode, but is at least one order of magnitude finer than the primary dendritic structure.

Experimental details
The binary Mg-Al test alloys were prepared by permanent mold indirect chill casting as described in [10]. The alloys were melted in an electric resistance furnace under a protective atmosphere

Phase-field model
In this study we apply the multi-phase-field method on a meso-scale with locally linearized phase diagram information [8]. In this method each grain of a distinct crystallographic phase is attributed by its own phase-field variable φ p ( x, t) ∈ [0, 1]. Three distinct crystallographic phases are considered: liquid melt, Mg-rich α-phase and a near stoichiometric M g 17 Al 12 βphase, where each solid grain is attributed by an identifier of the crystallographic phase and its orientation in space (for details see [5,8]). For the interface between melt and the α-phase an anisotropic interface energy σ pq = σ pq ( n) as a function of the interface normal n with respect to the crystallographic orientation of each grain is used [8]. The interface energy between the melt and β-phase as well as between the solid α-and β-phases assumed isotropic. In order to consider the effect of the external cooling rate on nucleation and growth of the solid phases we solve the heat balance in average over the whole domain with n xyz discrete elements: where T is the temperature of the system,Q ext is the heat extraction rate, ρ and c p are the mass density and the specific heat which are taken constant throughout the simulation.ψ pq are the local transformation rates between pairs of phases which sum up to the total transformation rate of individual phasesφ p = N q =pψ pq in the interface or a junction point with N phases (for details see [11]). L pq is the latent heat associated with the individual transformations between phases p and q. In the present study we treat all possible transformations between liquid ↔ α, liquid ↔ β and α ↔ β phases.  This particular version alters the point of maximum solubility in the α from 11.6% to 7.75% to match the area of the original liquid ↔ α two-phase area in the phase diagram. Thus the solidus slope (m α↔liq. ) was adjusted from the original -18.43 K/% to -27.58 K/%. Also C α↔β and T α↔β were adjusted to retain m α↔β . Consequently the original α ↔ β intersection-point (1229.84K ; 34.83%Al) was adjusted to (1301.84K ; 34.20%Al). All of the used values for the characterization of the phase diagram are summarized in table 1.  Table 2. Nucleation model parameters particle density for α-phase nucleation K α 0.2 · 10 14 m −3 particle density for β-phase nucleation K β 1.5 · 10 16 m −3 a of the normal size-distribution (mean) −1.5 · 10 −7 m b of the normal size-distribution (standard deviation) 1.0 · 10 −7 m

Nucleation model
In order to capture the dependence of the nucleation density on the cooling rate we employ the so-called "free growth model" by Greer and coworkers [13]. This model is based on a given distribution of seed particles of different sizes. Assuming perfect wetting conditions between the primary α-phase and the seed particles, the hemispheric cap of the α-phase on the seed particles defines the critical undercooling for free growth of the nucleating phase which is inversely proportional to the size of the seed particles [12]. In our simulations we start from a random distribution of the particles with a given density and size distribution over the entire simulation domain initially filled with liquid melt. The utilized normal distribution is given by Each particle is assigned an effective size 'd' which determines the critical undercooling for its activation. The resulting number of particles is typically higher than the actual number of the activated particles (active nucleation sites). Activation happens if a grid cell, which contains a particle, reaches it's critical undercooling, i.e. if the difference between actual temperature and melting temperature dependent on the local alloy concentration exceeds the local critical undercooling. A similar model has also been applied by Eiken et al. [8]. The parameters of the utilized nucleation model are summarized in table 2. Note, that only the tail of a normal distribution is used, thus a is negative, but only seeds with a positive size are created. We also assume a significantly higher nucleation density of the β-phase in the melt channels between the primary α-phase dendrites.

Simulation results
The simulations are performed on two different scales. Although, the primary solidification in three-dimensional domains can be investigated in order to reveal the effect of the cooling rate on the nucleation density of α-phase dendrites (see Fig. 3), it is computationally demanding to study the entire solidification process in 3D over the temperature range of more than 200K even with high cooling rates. In contrast, in the 2D simulations it is possible to investigate the entire solidification process starting from the primary α-phase dendritic solidification and continuing with the eutectic solidification of the secondary β-phase upon cooling below the eutectic point. It shall be noted, that on a scale of primary α-phase dendrites the eutectic structure in the interdendritic region cannot be resolved and as the result the secondary β-phase builds a bulk structure between the primary dendrites as can be seen in Fig. 4. In order to investigate eutectic growth in detail we zoom into the channel between two primary dendrites upon reaching the eutectic temperature. Here we select the scale of the simulations by a factor of 5 finer than the scale of primary α-phase dendritic solidification simulations. The heat extraction rate is assumed constant, which relates to an applied cooling rate (neglecting release of latent heat) CR =Q ρcp . We note that this cooling rate is different from the cooling rate evaluated from the experiments, due to the difference in definition. We will compare both entities by comparing  Figure 3. Results of a 3D simulation of primary α-phase nucleation and growth.
the solidification time (see below). All simulations of primary solidification are performed in a domain with periodic boundary conditions initially filled with homogeneous melt. The simulations of divorced eutectic growth start from a concentration profile taken around the melt channel of a 2D simulation of primary α-phase growth. Then a second particle distribution for β-phase nucleation is generated using the parameters given in Table 2.   Fig. 1 Simulations of consecutive α-and β-phase nucleation and their divorced eutectic growth were carried out in 2D. The applied cooling rates were varied from 5K/s to 25K/s. The resulting micro-structures are shown in Fig. 4. The solidified microstructure consists of a primary α-phase dendrites (blue) surrounded by the thin channels of β-phase (orange). The volume fraction of the β-phase is about 5% in all simulations. A clear trend of increasing α-phase grain density with increasing cooling rate can be observed which is in agreement with experimental observations (see Fig. 5). This behavior is achieved naturally by the implementation of a free growth model proposed by Greer et al. [13] which allows to reproduce the recoalescence effect responsible for the increased nucleation density at higher cooling rates. The corresponding simulated cooling curves are shown in Fig. 6.  Figure 6. Simulated cooling curves.
The solidification experiment with water cooling takes approximately 20 seconds from fully liquid to fully solid which matches our simulations with an applied cooling rate of 25 K/s. Thus the simulated microstructure shown in Fig. 4(e) corresponds to the experimentally obtained microstructure shown in Fig. 5(c). Comparison between simulation and experiment shows good agreement in grain density and connectivity. The characteristic grain-size in the simulation is about 150 µm which is in correspondence with the experiment for the same estimated cooling rate. Due to the prohibitive computational expense of slower solidification conditions, the lowest simulated cooling rate was 5 K/s and takes about 110 seconds of simulated time to fully solidify. The experimental results for air or furnace controlled cooling were much slower taking 600 and 1700 seconds respectively. However a clear trend of coarser grains with decreasing cooling rate can be observed.
Another important factor is the closed shell β-phase formation which is responsible for increasing the corrosion-resistance of Mg-Al-alloys. Thus good correspondence of our simulation results with experimental observations allows for computationally assisted design of corrosion resistant Mg-alloys. Of course there is a trade off between fine grain structure (better mechanical properties) and low connectivity of primary dendrites (better corrosion properties). The balance between mechanical and corrosion properties can be further tuned and investigated using the proposed method.
In order to investigate in detail the β-phase nucleation and growth kinetics and the formation of a divorced eutectic microstructure with the tertiary α-phase nucleation in the residual melt these simulations were performed in 3D. The results of a divorced eutectic solidification simulation are shown in Fig. 7.
The simulations starts with a residual melt channel with the eutectic composition of 31% and nearly eutectic α-phase at the bottom and the top of the simulation volume. This system is cooled down starting at 712.0 K, slightly above the eutectic temperature. The α-phase is still slowly growing on this scale. As soon as the undercooling is sufficient, nucleation automatically takes place and generates a large number of β-grains in the residual melt. The high nucleation density of 3 · 10 13 particles per m 3 was chosen because it is expected for the particles to be aggregated in the residual melt. With an abundant supply of nucleation sites and our chosen cooling rate of 5 K/s we get around 50 nuclei within the channel. A strong tendency of β-phase growing towards and at the α-liquid-interface can be observed. This is because Al is segregated strongly by the growing α-phase into the melt close to the interface, since the solubility for Al in α is only 11.6% at the eutectic temperature of 709.4 K. This leads to a higher concentration of Al in the melt close to the slowly growing α-phase. The β-phase, consisting of about 40% Al grows preferably in the area of the increased concentration. Due to this circumstances a rapid, full enclosure of the α-phase by the β-phase can be observed. Thus, further growth of the α-phase is prohibited. Continuing β-phase growth depletes the concentration in the residual melt. It is impossible for the melt at that state to fully solidify as β-phase due to the lack of Al. This leaves the system in a state where a tertiary α-phase nucleation is necessary to fully solidify. Investigations of this behavior will be presented elsewhere.

Outlook
We have presented 2D simulations of consequent α-and β-phase evolution, shown in Fig. 4. 3D simulations of the same accuracy will require about two orders of magnitude more memory and computational power, which is currently not feasible on a single compute node.  Fig. 7. The domain is subdivided into hundreds of blocks shown to the right with the rank-dependent color-coding.
Therefore, distributed-memory parallelism based on block-structured domain decomposition has been implemented in OpenPhase, which will be used for further 3D simulations on large grids. Fig. 8 shows a benchmark-simulation of eutectic solidification in α-channels. The computational domain of 300 × 300 × 300 grid points is subdivided into hundreds of blocks that are assigned to 128 processors. As the phase-field values can only change in the interface regions, OpenPhase only calculates the updates of the phase-fields in these regions. While this approach drastically  improves the performance, it leads to a very inhomogeneous distribution of the computational load and requires sophisticated load balancing techniques. An over-decomposition of the computational domain allows a graph-partitioner to assign blocks to processors, so that in the worst case the load imbalance is bounded by the computation time on the most expensive block. Therefore, expensive blocks are divided into multiple parts, when the implemented work load model predicts a benefit for the wall clock time for the next time steps. The system is rebalanced at runtime, whenever it is needed. Figure 9 shows performance results on the supercomputer JuRoPA at Forschungszentrum Jülich with up to 1024 processors. Excellent weak scaling is observed: it can be seen that the computation time does not change significantly when the total work load per process is kept constant. This allows the realization of large scale simulations utilizing a parallel cluster without increasing the wall time. For primary solidification perfect weak scaling results were achieved with up to 2048 processes on grids with up to 800×800×1600 grid points. With this new possibilities of large-scale simulations, we hope to set up a full 3D solidification simulation analogous to the 2D calculations shown in Fig. 4. Also, due to an increased resolution due to higher count of grid points, a simulation can be carried out resolving both: the primary as well as the eutectic solidification on their different length-scales. For a realistic setup a grid of about 600 × 600 × 600 grid points is needed. Thus the time needed for one timestep is estimated to be about 2 seconds with 4096 processes. A calculation of 25 seconds real time would take approximately 11.5 days. [14] 8. Summary In this paper we presented an approach to simulate the solidification of Mg-Al alloys. A phase diagram based on Thermocalc TCS 4.0 database TCBIN was linearized and adjusted for consistency. The "free-growth" model was implemented into OpenPhase for realistic nucleation behavior of both α-and β-phases. 3D simulations of the first stage of primary solidification were carried out. A full solidification cycle was simulated in 2D including a primary α-as well as secondary β-phase nucleation and growth in one simulation, with close to 200 K of cooling in between the two events. Resulting cooling curves as well as simulated microstructures were compared to experiments. Remarkable agreement has been obtained for characteristic grain size and β-phase closed shell formation. Due to the insuficient resolution a dual-scale approach was employed to model an interdendritic region of eutectic melt prior to final solidification. Divorced eutectic behavior was reproduced showing a need for tertiary nucleation of α-phase due to full shielding of primary dendrites after a short period of time. The possibility for future highly parallelized computations of large scale 3D domains with high resolution at different scales was