Full-field simulation of solidification and forming of polycrystals

The phase-field method has emerged as the method of choice for simulation of microstructure evolution and phase-transformations in material science. It has wide applications in solidification and solid state transformations in general. Recently, the method has been generalized to treat large deformation and damage in solids. A through process full-field simulation will be presented starting from solidification and ending with the evolution of damage during large deformation. Aspects of numerical discretization, efficient numerical integration and massive parallelization will be discussed.


Introduction
Due to their low density, magnesium alloys are of high importance for modern lightweight structures and are widely used for various applications in automotive and aerospace industries as well as in the consumer electronics.The most applied alloying element in Mg alloys is aluminum which is used for AZ91, AM50 (cast) or AZ31 (wrought) technical alloys.The mechanical properties of Mg-Al alloys are primarily determined by their microstructure which 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.
An exemplary SEM-micrograph of the later simulated Mg-5%Al is reproduced in fig. 1.The addition of Al in small concentrations has a positive effect on the corrosion properties of the Mg-Al alloys, mainly due to the formation of a percolating β-phase network [3] which is less electrochemically active compared to the α-phase.Experiments, as in [4] revealed that the formation of a closed shell of β-phase, preventing α-phase dendrites from building networks, is a good strategy to significantly increase the corrosion resistance of cast Mg-Al alloys.Furthermore, it severely hinders the galvanic corrosion if Mg-alloys are brought into contact with more noble metals.Thus, tailoring the Mg-Al alloy microstructure is of great importance for controlling its corrosion resistance [5].In this work, we apply the phasefield method [6,7,8] to simulate the solidification microstructure of Mg-Al alloys.We outline the strategy for the full solidification cycle simulation using the phase-field method followed by a damage simulation.
The information about the as cast morphology and the distribution of the secondary β-phase in particular shall point towards improved materials and process design, dark: α-phase, bright: β-phase also in terms of corrosion.The proposed approach enables metallurgical process control by correlation of metallurgical parameters to metallic phase fractions and mechanical properties.This will significantly contribute to the development of Integrated-Computational-Materials-Engineering (ICME) by adding service-life aspects to computational materials design.

Phase-Field Model
In this study, we apply the multi-phase-field method on a meso-scale with locally linearized phase diagram information [9].In this method each grain of a distinct crystallographic phase is attributed by its own phase-field variable ߶ (, ‫)ݐ‬ ∈ [0,1].
Three distinct crystallographic phases are considered: liquid melt, Mg-rich α-phase and a near stoichiometric Mg 17 Al 12 β-phase, where each solid grain is attributed by an identifier of the crystallographic phase and its orientation in space.
The governing equations for the evolution of phase fields and the concentration are given below (for details see [7,9]). with Here ߶ is the phase field parameter; k, l and m are the phase field indexes; μ kl is the interface mobility between phases k and l; N is the number of non-vanishing phase fields/grains in a given point; I is the generalized curvature term; η is the interface thickness; h is the shape function related to the phase field contour; Δ݃ is the thermodynamic driving force; ‫ܦ‬ is the diffusion coefficient of Al in phase k; ‫ܒ‬ is the antitrapping current.
The domain size of our simulations is set to 300µm in every direction with a grid spacing of 1µm resulting in 300³ cells.The interface thickness η is set to the length of 5 grid-cells.
In order to consider the effect of the external cooling rate on nucleation and growth of solid phases we solve the heat balance in average over the whole domain with n xyz discrete elements: Here, ܶ is the temperature of the system, ܳ ̇௫௧ is the heat extraction rate, ߩ and ܿ are the specific heat and the mass density which are taken as constant throughout the simulation.
߰ ̇ are the local transformation rates between pairs of phases which sum up to the total transformation rate of individual phases ߶ ̇ = ∑ ߰ ̇ ே ஷ in the interface or a junction point with N phases (for details see [11]).‫ܮ‬ is the latent heat associated with the individual transformations between phases k and l.In the present study, we treat all possible transformations between liquid and α, liquid and β and α and β phases.

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 [12].This model is based on a given distribution of seed particles of different size.Assuming perfect wetting condition 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 [13].
In our simulations, we start from a random distribution of 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 its 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 in [9,10].For further details and used parameters see [14].

Mechanical framework
Using standard notation in continuum mechanics, the deformation mapping maps a material point ξ (position vector) belonging to the reference configuration to its position in the actual deformed configuration x(ξ).Locally, this mapping can be described by the deformation gradient: with u being the displacement vector.The deformation gradient has the polar decomposition with the rotation tensor R and the right stretch tensor U.
For small deformation, where డ‫ܠ‬ డ ≪ can be assumed, and for defining a predictor strain increment in the iterative solution procedure as outlined below, we neglect the nonlinear terms, approximating డ డ ~డ డ‫ܠ‬ .Hence, this assumption results in the definition of the engineering strain Δε and the rotation increments ΔR: Although the given description of deformation does not rely on tracing back the current configuration to the original configuration, it is helpful to integrate the strain (either total, elastic or plastic) by a proper strain measure.
For the sake of additivity, we use the logarithmic strain: The basic idea is to formulate an iterative scheme that, first, by iterative solution of the small strain mechanical problem finds the mechanical equilibrium : And, second, at the same time applies the geometry change in an Eulerian setting.The geometry update is given by the previously defined rotation increment (eq.4) and by the displacement u itself that can be applied to any conserved scalar variable ρ by the compressible transport equation: Rotation increments Δ‫܀‬ are explicitly applied to any tensorial variables τ In detail, the scheme is structured as follows.Assuming that at time t the system starts in mechanical equilibrium, the reference state is identified with the current state: The system will be driven out of mechanical equilibrium by: 1) changing the boundary conditions, e.g., by applying external load, 2) changing the material's constants, e.g., by diffusion of species or temperature and by the impact of damage, 3) plastic relaxation ઽ or 4) phase transformation or grain boundary motion associated with volume change and the eigenstrain ઽ * .In the following, we will consider in detail the effect of damage on the elasticity tensor C, plastic relaxation, and phase transformation summing up to stress-free strain increment.The present mechanical state is no longer divergence free: Here the superscript () denotes the zeroth iteration.An elastic corrector strain ઽ ଵ for the first iteration is calculated defining a corrector displacement ‫ܝ‬ ଵ : In the next stage, the intermediate displacement ‫ܝ‬ ଵ is applied by a geometry update according to eq. ( 12) and the intermediate rotation Δ‫܀‬ ଵ is applied to update the Cauchy stress with the intermediate stress increment With the new stress increment, the stress-free strain Δઽ has to be updated according to the plastic relaxation and the stiffness tensor ۱ according to the damage evolution as described in the damage section.Eq. ( 12) describes an update of the actual stress state of a material's point that is accumulated from the incremental stresses and rotations based on a so-called co-rotational objective time derivative in line with [17].With the new stress predictor Δો ଵ at hand, the iteration procedure is repeated.Convergence is reached after n iterations when the corrector strain vanishes ઽ → 0 as well as the rotation increment Δ‫܀‬ → .Then the stress tensor ો becomes divergence-free.
The presented approach can be coupled with any plasticity framework in a straightforward manner.Here, we apply a simple damage independent viscoplastic model in order to demonstrate the effect of damage on elastic softening during plastic yielding.The influence of damage instead will be considered by differentiation of nominal and effective stress space (see [18]) where the effective stress is the driving force for creep.With a strain rate constant A and a stress exponent k we define the plastic strain rate as Thus, accumulation of large plastic strain is straight forward, simplifying its coupling with plastic strain dependent models, like the damage model.
In order to model failure of the modeled microstructure, a damage model published in [18] is used.It assumes the stiffness tensors being dependent on the equivalent plastic strain: Where the equivalent plastic strain ‫‬ is calculated via integration of the equivalent plastic strain rate:

Result
The governing equations (eq.1-3) are solved with an explicit Euler integration scheme.Our simulations start with a cubical fully liquid domain with an edge length of 300µm size at 900K (above the melting temperature).
From this stage, a fixed heat extraction is applied at each time step.After a few simulated time steps, the system is undercooled and nucleation occurs.The nucleated grains of primary α-phase grow in a dendritic fashion and fill up most of the volume.The nucleation density is highly dependent on the heat extraction rate applied (see [15]).
After the growth slows down, the heat extraction rate becomes again increasingly more dominant than the latent heat production and the system cools down further.When the domain gets cooled below the eutectic temperature, nucleation and growth of the secondary βphase occurs.The resulting microstructure is shown in Fig. 2. The eutectic solidification happens on a different scale and is usually treated in a multi-scale approach (see [14]).In this work we can assume the rest of the melt at this point in time to decompose into secondary β-phase and tertiary α-phase and form a fine-grained eutectic structure.Because of the scale difference of our primary grain structure and the interdendritic eutectic, we treat the eutectic as one phase with continuous properties rather than a mix of two phases.
Having the solidified microstructure, a tensile test like calculation is performed.In general, all through the mechanical simulation, the phase-field evolution can be considered to ensure the proper interface profile.It has however reduced dynamics due to the low temperatures of mechanical testing.The system is strained along the yaxis until failure is observed.By failure we mean the formation of a large region of maximal damage ‫ܦ‬ = 0.95 and a stress drop as seen in fig. 3.In the current simulation, this happens after the applied strain of approximately 10%.
The stiffness coefficients tensor used here is taken from [16] and can be seen in tab.1.The plasticity model is parametrized by k = 5 and A = 1e-13 Pa -5 .The damage parameters are set as ‫‬ = 0.08 and ‫‬ ௫ = 0.10.The 10% deformed system is shown in fig. 4. Each of the grains is modelled to be individually oriented.Due to the stiffness anisotropy, the different grain orientations result in a stress field that starts plastic strain evolution and finally creates damage and softens the material.Since there is only a slight variation in the stiffness coefficients, the resulting damage does not localise in the grain boundaries.Instead, after there is a maximally damaged material, the damaged region developes in an almost spherical manner.However, the stress-strain curve presented in fig. 3 reproduces qualitatively realistic behaviour (compare Song et al. [19]).

Parallelism
Large scale phase-field simulations, such as the simulation of Mg-Al alloy solidification mentioned above, require significant amounts of memory and computation time, thus exceeding the capabilities of single node computers.Distributed-memory parallelism using MPI allows us to overcome these limitation and makes simulations of large computational domains feasible.A simple domain decomposition that divides the domain into equally sized sub-domains can yield an insufficient performance due to load-imbalances caused by an inhomogeneous distribution of the computational load in the domain.
The cause for these inhomogeneous load distributions is the active parameter tracking as proposed by Vedantam et al. [1].
Here, only the interactions between locally existing phase fields in each grid point are calculated, which significantly reduces the total computational complexity when dealing with a large number of different phase fields.This practice, however, can lead to inhomogeneous load distributions as the computational complexity in a grid point is dependent on the number of locally active phase fields in this grid point.
In order to fully utilize the computational resources of modern parallel computer architectures, the loadbalancing techniques presented by Tegeler et al. [2] are applied.This involves an over-decomposition of the computational domain, i.e. more sub-domains than MPIprocesses as seen in fig.6, so that a greedy graphpartitioning algorithm can be used to assign sub-domains to processes with the goal of minimizing the loadimbalance and communication cost.
Even if the computational load is balanced between processes, waiting times can occur at synchronization points, i.e., halo exchange between the sub-domains if the computational load is not balanced between all synchronization points.Therefore, multiple synchronization points within one time step are avoided by increasing the size of the halo to allow multiple stencil operations per time step.
In order to improve the worst case load-imbalance guaranteed by the graph-partitioning algorithm, we use bisection of computationally costly sub-domains into two smaller sub-domains of reduced computational load along a cutting plane, which splits the computational load into two approximately equal parts.
As the bisection of a sub-domain introduces additional computations on the halo as well as communication, an optimal sub-domain load is determined at runtime taking into account the hardware properties and the specific load distribution and only subdomains of higher load are split into smaller domains.
For the simulation of Mg-Al alloy solidification as seen in fig. 2 a total of 6 compute nodes with 16 cores and 128GB of RAM each were used.
In order to benchmark the scaling properties of the parallelization, a set performance tests using the simulation of Mg-Al alloy solidification with the multiscale approach mentioned above have been conducted.In a channel between two primary α-phase dendrites a secondary β-phase nucleates and grows much faster than the primary α-phase.This is depicted in fig.6 together with the domain decomposition generated by the presented load balancing techniques.
It shows that the sub-domains near the border of the channel have been split into smaller sub-domains as the computational load is high in this region, due to interfaces between αand β-phase and the liquid.
The results of the scaling benchmark are shown in fig.7 and it can be seen that if the computational domain per core is kept constant some fluctuation in the weak scaling performance are observable, but no significant decline in performance as the largest time per step on the largest domain is even slightly smaller compared to smallest domain in the set of benchmark runs.

Conclusion
The presented work incorporates modelling of two physical processes occurring in a magnesium system: the solidification process that results in a typical Mg microstructure and a subsequent mechanical testing process quantifying the quality of the resulting alloy microstructure.The used approach follows the spirit of computational execution of experimental work and defines a possible workflow corresponding to modelling of a material processing chain.Both parts are executed within the phase-field framework profiting on its morphology describing potential accompanied by mathematical simplicity.Due to specification of realistic model parameters the observed results are comparable to experimental observations.With this prove of concept for the used work procedure at hands, the authors plan to use precise mechanical data for a specified system and a crystal structure informed plasticity to improve the mechanical model gaining a one to one relation between the experimental findings and the simulation predictions.

Figure 3 .
Figure 3. Stress-strain curve obtained during the simulation of damage in the Mg-5wt.%Al system.

Figure 5 .
Figure 5. Damage evolution (black surface) during the last calculation steps in zoomed view.

Figure 6 .
Figure 6.Interdendritic Mg-Al alloy solidification with the domain decomposition shown on the right-hand side.Different colours indicate different process ranks.

Figure 7 .
Figure 7. Combined weak (dashed lines) and strong (solid lines) scaling performance for the interdendritic Mg-Al alloy solidification as seen in fig.6.