Simulating microgalvanic corrosion in alloys using the PRISMS phase-field framework

In this prospective paper, we first review the existing simulation tools to simulate microgalvanic corrosion during free immersion. Then, we describe a recently developed application that employs PRISMS-PF, an open-source, high-performance phase-field modeling framework. The model employed in the application accounts for the electrochemical reaction at the metal/electrolyte interface and ionic migration in the electrolyte to determine the evolution of the corrosion front. We present the implementation details for the application and discuss its features such as super-linear parallel scaling performance for a sufficiently large system. Finally, we demonstrate the capability of the application by simulating corrosion of the matrix phase of an alloy near a secondary phase particle in two and three dimensions.


Introduction
Corrosion of metals is a long-standing economic and engineering problem. The total annual direct cost of corrosion was estimated to be over $2.5 trillion in 2013. [1] It also limits rapid development of new metallic materials due to a long time scale over which the degradation occurs. For example, Mg alloys are lightweight and have excellent strength-toweight ratio and exhibit versatile processability via rolling, extrusion, and casting when alloyed with Al. [2] However, such alloys suffer from localized corrosion triggered by the microgalvanic coupling between the Mg 17 Al 12 intermetallic particles and the Mg matrix. [3] While different forms of corrosion may take place on Mg alloys, [4] galvanic corrosion, either through contact with a dissimilar metal or by alloying, is a common occurrence because of the electronegativity of Mg. [2] Secondphase precipitates also lead to microgalvanic or local coupling effects. [5][6][7] In addition, localized corrosion is often observed. [8] The conditions associated with each form of corrosion have not been fully determined, and how different alloy chemistries and precipitate/grain microstructures lead to such a range of corrosion behavior is yet to be quantitatively understood. Thus, it is clear that we must further develop the science behind their corrosion behavior before Mg alloys can be widely used as a structural material in a corrosive environment. [2] Only when we have a better mechanistic understanding underlying the corrosion process, we will be able to conduct science-driven alloy design and optimization of Mg alloys for a range of applications. However, due to the long lifetime of structural materials and slow progression of corrosion under non-accelerated conditions, experiments can take many months or sometimes years to produce data and resulting insights. Therefore, it is critical that modeling and simulation of corrosion are integrated into the material design cycle for Mg alloy development.
Recent advances in the modeling of corrosion employing the phase-field, [9][10][11][12][13][14] level set, [15] and sharp interface [16,17] methods allow us to predict how corrosion proceeds in different materials and conditions, taking into account the different electrochemical properties of phases and grain boundaries present in the system. However, in order to bring such capability to the broader community and to advance the science of corrosion, it is important to establish an open-source software that is flexible, easy to use, and highly efficient. Existing computational tools that are publicly available freely or for a fee to simulate the evolution of the corrosion front include the COMSOL Corrosion module [18] and the HOGNOSE [19] application based on the MOOSE [20] framework. In COMSOL Multiphysics®, which is a commercial software, the moving corrosion interface is tracked by using either the Arbitrary Lagrangian-Eulerian method or the level set method. The evolution of this interface is governed by electrochemical reaction kinetics and ionic transport within the electrolyte, which includes diffusion and migration. However, because it is a commercial software with a high license fee, COMSOL Multiphysics® is not available to all researchers and engineers in the material science and engineering community. Furthermore, for the same reason, it is frequently unavailable on many high-performance-computing clusters, which limits its application to large 3D problems. Additionally, meshing of complex microstructures within COMSOL Multiphysics® requires the users to have advanced knowledge of meshing packages. [21,22] In contrast to COMSOL Multiphysics®, the HOGNOSE application is developed based on the open-source MOOSE framework, a finite element parallel computational framework targeted at the solution of coupled, nonlinear partial differential equations. The MOOSE framework is highly modular and features a plug-in infrastructure that simplifies definitions of physics, material properties, multiphysics coupling, and postprocessing. The HOGNOSE application simulates early-stage irradiation-enhanced corrosion of Zircaloy-4 in pressurized water reactors using the phase-field method. The implemented model describes the evolution of the metal-oxide interface in fuel cladding, where corrosion is enhanced by the release of iron from secondary phase particles. However, galvanic or microgalvanic effects are not considered to be contributing factors to corrosion in this scenario. Therefore, to our knowledge, there is not yet a publicly available application or software tool to simulate galvanic or microgalvanic corrosion with MOOSE or any other open-source frameworks.
To enable the wider materials science community to easily and efficiently simulate galvanic or microgalvanic corrosion and to integrate the Materials Genome Initiative [23] (MGI) approach into their research, we offer a new corrosion application within PRISMS-PF titled corrosion_microgalvanic. [24] Below, we discuss the model equations implemented in the application, the corresponding implementation details in PRISMS-PF, and the application's features such as parallel performance and 2D and 3D capabilities. PRISMS-PF is an open-source high-performance framework for phase-field simulation of microstructure evolution. [25,26] This framework has been developed within the Predictive Integrated Structural Materials Science (PRISMS) Center. [26] The code is built on the deal.II finite element library [27] and exhibits excellent computational performance, which is enabled by the use of a matrix-free variant of the finite element method (FEM). This method allows for an efficient time integration, in addition to mesh refinement capabilities and multi-level parallelism (distributed memory parallelization, task-based threading, and vectorization). The framework is designed with a modular, application-centric code structure for ease of use and contains 28 built-in applications. Performance for a solidification benchmark problem [28] has been shown to meet or exceed that of other phase-field frameworks (including MOOSE). [25] Some functionalities of the framework include built-in methods for nucleation, order-parameter remapping for grain growth simulations, and an iterative solver for nonlinear partial differential equations. Recent applications of the PRISMS-PF framework include simulations of precipitate morphology and composition, [29] pitting corrosion of microstructures from additive manufacturing, [30] and dendritic solidification in Mg alloys. [31]

Model equations
The model utilized in this application considers migration of ions in the electrolyte, electrochemical reaction at the metal/ electrolyte interface, and the motion of the interface due to the reaction. Below, we discuss the equations used to describe these phenomena, which are based on the corrosion model described by Chadwick et al. [9] with an addition of an order parameter to represent the cathodic phase. Additional details of the model can be found in Ref. 9. All the variables and symbols used in the model equations below are listed in Table 3.

Phase-field equations
In this application, we employ the phase-field method as a tool to track the moving metal-electrolyte interface as well as the interfaces between various phases of the metal. Field variables known as order parameters are used to describe the geometry of all phases. Each order parameter has a value of 1 within the corresponding phase and 0 outside that phase. At the interface between the phase and another phase, the order parameter smoothly transitions between 0 and 1. In this formulation, we use ψ and ϕ j to represent the order parameter for the electrolyte phase and the jth metal phase, respectively. Furthermore, the free energy functional, F , of the system can be written in terms of the order parameters and their gradients as [32] where ǫ is the gradient energy coefficient and represents the system volume. Note that the quantities F and ǫ 2 are normalized by the height of the double well, W, which has units of J/ m 3 . The bulk free energy f 0 is given by [32] We note that we only consider two phases in the metal for this application: ϕ 1 indicates the anodic phase and ϕ 2 indicates the cathodic phase. Nonetheless, Eqs. 1 and 2 are kept in their general form to allow the consideration of more phases. The evolution of the anodic metal/electrolyte interface is governed by Cahn-Hilliard equations with a source term, [33,34] where M ϕ 1 and M ψ are the mobilities for phases ϕ 1 and ψ , respectively, and v is the normal component of the velocity of the anodic metal/electrolyte interface, which is related to the anodic reaction current density (RCD), i rxn,1 , and an interpolation function for the anodic phase, ξ 1 (that is equal to 0 in the cathodic phase and 1 in the anodic phase) via (3) where V m is the molar volume of the corroding metal, z m is the charge number of the corresponding metallic cation in the electrolyte, and F is Faraday's constant. We note that Eq. 5 is a slightly modified version of the expressions provided by Refs. 9, 15, and 35 for the interface velocity. The equations used to calculate i rxn,1 and ξ 1 are described in the next subsection. It is important to note that the terms proportional to v on the right-hand side of Eqs. 3 and 4 are not due to the bulk motion. Rather, these terms are added as source terms to model the motion of the electrolyte/anodic phase interface due to the corrosion reaction at the interface. Note that multiplication by the magnitude of the gradient of ψ confines the region where these terms are nonzero. In Eq. 3, the Cahn-Hilliard dynamics, rather than the Allen-Cahn dynamics, was chosen to describe the evolution of ϕ 1 to ensure that the change in mass of anodic phase is due only to the corrosion reaction. Following Ref. 9, we use the same approach for the electrolyte phase (Eq. 4) but with the opposite sign of the source term to ensure the direction of the interfacial displacement is correct. We do not consider any deposition on the cathodic phase, and thus the evolution of ϕ 2 is given by the Cahn-Hilliard equation, [33] where M ϕ 2 is the mobility for the phase ϕ 2 . While corrosion does not affect the cathodic phase, it is important to evolve this equation so that the order parameters adjust accordingly to the free energy functional in Eq. 1 as other order parameters evolve; without doing so can lead to undesired numerical artifacts and instability. In this work, we do not consider the deposition of the corrosion product, although it is a potential valuable extension, as discussed later. Finally, all the mobilities are set to be equal and are dependent on i rxn,1 and ξ 1 . These mobilities are given by an expression similar to that provided by Chadwick et al. [9] : where 2δ = 2 √ 2 ǫ 2 is the equilibrium interfacial thickness. The mobility coefficients in Eq. 7 are defined in such way that the dimensionless Peclet number, 2δv/M , is unity to ensure that the moving interface can maintain the profile required by the free energy functional. [9] Electrochemical equations The ionic transport in the electrolyte is caused by migration, diffusion, and convection. We assume that there is continuous replenishment of the ions in the electrolyte through sufficiently large fluid flow, i.e., the well-stirred limit. Under these assumptions, no ionic-concentration gradients exist in the electrolyte, due to which the diffusion of ions in the electrolyte does not need to be considered. Additionally, fluid flow does not need to be explicitly considered for in this limit. Such a setup is consistent with other studies reported in the literature [13,36,37] and is valid for systems in contact with flowing electrolyte like structural components operating in rivers and seas. Consequently, the current density in the electrolyte, i e , is given by [38] where κ e is the ionic conductivity of the electrolyte and e is the electrostatic potential in the electrolyte. Current continuity requires that in the electrolyte. Furthermore, at the metal/electrolyte interface, the normal component of i e must be equal to the reaction current density, i rxn : where n m/e is the unit normal vector of the metal/electrolyte interface pointing out of the metal. Natural (zero-normal-derivative) boundary conditions are imposed at all the other boundaries in the system. Equations 8-10 are combined using the Smoothed Boundary Method (SBM) [9,[39][40][41] to obtain an expression for e within the electrolyte. The SBM applies mathematical identities in calculus to transform a partial differential equation to an alternative form that implicitly includes the boundary conditions. Here, since these equations are only valid within the electrolyte or its boundary, the SBM is applied with ψ as the domain parameter. The resulting equation for e is given by Note that even though i rxn is calculated throughout the system, its value is only meaningful in the interfacial regions and affects the electrolyte potential only where |∇ψ| is nonzero, i.e., at the metal/electrolyte interface.
To calculate i rxn in Eqs. 10 and 11, we must first obtain the cathodic and anodic RCDs independently. In general, RCD is given by the product of the corrosion current and an exponential function containing the overpotential (which includes the corrosion potential). The anodic RCD is obtained by a modified Tafel relation as [38] where i corr,1 is the corrosion current density for the anodic phase, η 1 is the anodic overpotential given by m − e − E corr,1 ; m and E corr,1 represent the applied electrostatic potential of the metal and the corrosion potential of the anodic phase, respectively; and A 1 is the corresponding Tafel slope. Figure 1 shows a typical Tafel plot and how the corrosion parameters ( i corr , E corr , and A ) are obtained. In Eq. 12, the current density smoothly approaches i max as i 0,1 e η 1 /A 1 becomes large, and thus i max sets the maximum current density. This constraint is added to account for any kinetic limitations such as limiting ionic transport in the electrolyte or charge transfer current density that may exist in the system due to surface passivation or the deposition of the corrosion product. [42] We assume no such limitation on the cathode, and thus, the cathodic RCD is simply given by the standard Tafel relation [38] where η 2 = � m − � e − E corr,2 , and the cathodic corrosion current density, the corrosion potential, and the Tafel slope for the cathodic phase are denoted by i corr,2 , E corr,2 , and A 2 , respectively. Note that we assume that the metallic conductivity is large enough to ensure that m is nearly constant throughout the metal. Additionally, we assume corrosion to advance under an unbiased condition, i.e., free immersion, and thus we set m = 0 . Furthermore, the use of natural (zero-normalderivative) boundary conditions for e on all the borders of the system in contact with the electrolyte ensures that the net current along the metal/electrolyte interfaces is zero. This can be shown by applying the divergence theorem on the volume (area in 2D) occupied by the electrolyte along with Eq. 8. Finally, the reaction current density in Eqs. 10 and 11 is then calculated from the cathodic and anodic RCD as a linear combination of the two: where ξ j is an interpolation function that varies from 0 to 1 , which indicates the region of the corresponding phase. In this work, we defined the interpolation function as where ζ is a small number added to the denominator to avoid division by zero. Extensive testing indicates that ζ must be small enough but not too small, ≈ 10 −3 , in order to ensure smooth behavior near the interfacial region (see Fig. 7 in the appendix).

PRISMS-PF implementation
The governing equations of the model are implemented into the PRISMS-PF framework, [24] which employs a matrix-free variant of FEM. We provide a list of variables/parameters names used in the model equations and the corresponding names used in the code in Table I. An explicit forward Euler scheme with constant step size is used for time integration. Natural (zero-normal-derivative) boundary conditions are applied for all fields at all the boundaries of the computational domain, but these conditions can be modified to periodic or Dirichlet Figure 1. A schematic of a typical Tafel plot. The magenta dashed lines represent the linear approximations to the anodic and cathodic polarization curves. E corr and i corr are obtained from the intersection of these lines and the Tafel slopes are obtained as the slopes of these lines. [43] Table I. A list of variables and parameters names used in the model equations and the corresponding names used in the PRISMS-PF Microgalvanic Corrosion code (corrosion_microgalvanic). [24] Variable/parameter name in the model equations Variable/parameter name in the code EcorrCathodic ζ lthresh boundary conditions. The equations are discretized in space with a rectangular adaptive mesh with first-order elements in two or three dimensions. The minimum linear element size is chosen such that the interfaces are well resolved, i.e., an interface contains at least four nodes across the region. The highest level of mesh refinement is selected dynamically for all regions where either of ϕ 1 , ϕ 2 , or ψ lie between the values 10 −4 and 1 − 10 −4 . The smallest and largest element size are set to ~ 3 × 10 −8 m and ~ 1.23 × 10 −7 m, respectively. The time step is set as 0.01 s, and the mesh adaption is carried out every 2000 time steps (i.e., 20 s). The order and domain parameters are bounded to the range [0,1] at the start of every time step to avoid negative values of the order/domain parameters near the interfaces. Moreover, all the interfaces are equilibrated for 100 s without corrosion to ensure that ϕ 1 , ϕ 2 , and ψ are consistent with the governing equation in the static condition based on the initial condition. The equilibrium electrostatic potential field, e , is computed at every time step using the integrated nonlinear solver functionality. These default settings can be modified by appropriately changing the following files in a manner described in Ref. 25 and the PRISMS-PF user manual [44] : parameters.prm: This is a parsed text file that contains information about the material parameters such as i corr,j , E corr,j , A j , κ e , V m , and z m ; geometry-related parameters like dimension and size; mesh-related parameters such as element order, adaptivity criteria and frequency, and level; solver parameters like tolerance criterion and value, time step, number of solver iterations, initial guess for e ; and boundary conditions. All these parameters can be modified by the user according to their requirement. We recommend a combination of literature survey and experimental measurements, which include polarization measurements, to determine the material parameters for different systems. The rescaled gradient energy coefficient, ǫ 2 , which controls the interfacial thickness, is also set in this file. The value of this parameter should be chosen so that the interfacial thickness is much smaller (~ 10%) than the length scale of the smallest feature of interest in the problem. Accordingly, the size of the mesh at the highest refinement level should be set to ensure that there are 3 ∼ 5 elements across each of the interfaces and the time step should be small enough to ensure the stability of the time stepping algorithm. We also note that the solver parameters and miscellaneous parameters such as ζ may need to be tuned if model parameters are modified. As with all PRISMS-PF applications, the microgalvanic corrosion application does not need to be re-compiled if only parameters.prm is modified.
ICs_and_BCs.cc: This is a code file that has the information about the initial state of the system. Users can modify the initial geometry of the system by appropriately setting the fields ϕ 1 , ϕ 2 , and ψ . We note that the PRISMS-PF also supports a file input for setting the initial conditions of the system.
Other files: The governing equations of the model are set in equations.cc, and the main class of the application and its member functions and variables (including model constants) are declared in customPDE.h.
The application outputs data in the Visualization Toolkit Unstructured Points Data (VTU) format. We recommend using a tool such as VisIt or Paraview to visualize and analyze the output files. Both these software packages support Python scripting, which enables automation of visualizations and postsimulation analysis.

Parallel performance
The parallel performance of the application was tested for a 3D system with ~ 4.8 × 10 7 degrees of freedom, and for 500 time steps. The system consisted of a cylindrical cathodic phase ( β ) (~ 21% volume fraction in the metal) surrounded by the anodic phase ( α ), and both phases are in contact with an electrolyte, as shown in Fig. 2(a). The symmetry of the setup allows one to perform the simulations in only one quarter of the system. The result from the strong scaling test, performed on the "ICX-normal" queue (Intel Xeon Platinum 8380, 2.3 GHz, 80 cores, 256 GB DDR4 RAM) of the TACC Stampede2 cluster, is shown in Fig. 2(b). As can be seen, the application exhibits super-linear parallel performance even up to 3200 cores, which is enabled by an increase in the available cache as the number of cores is increased. Thus, we conclude that the application is well suited for investigating scientifically and technologically relevant large systems with an excellent parallel efficiency. In the next section, we demonstrate the utility of the application in investigating microgalvanic corrosion in Mg alloys.

Demonstration examples
To demonstrate the utility of the application, we simulated the microgalvanic corrosion effect in an Mg alloy. We employed the properties of AM30 alloy (Mg-3wt%Al-0.4wt%Mn) for the α phase and of Mg 17 Al 12 for the β phase, as in Refs. 36  We first present results for a 2D system, where the anodic ( α ) and cathodic ( β ) phases are present in alternating lamellae. By employing natural boundary conditions, we can represent this configuration with a simulation domain that includes only half of the α phase and half of the β phase of one full period, as shown in Fig. 3. We demonstrate that the application captures the evolution of the corroding interface and the electrolyte potential and use it to illustrate the effect of the area fraction of β phase on the corrosion rate. Finally, we demonstrate the ability of the application to simulate the evolution of the corroding interface in 3D by using the 3D system shown in Fig. 2(a). Figure 3 shows the evolution of the anodic metal/electrolyte interface along with the evolution of the electrolyte potential, e for a system with 20% area fraction of β (defined as a fraction of the entire domain width). All other model parameters are provided in Table II. For the corrosion parameters, we adopt the values determined by Beura et al., [37] which were based on the polarization curve measured by Deshpande. [36] As can be seen, the anodic region corrodes at a faster rate near the anodic metal/cathodic metal interface than away from the cathodic region. The nonuniformity of the corrosion rate occurs due to the spatial variation  --of the reaction current density (RCD), which is higher near the anodic metal/cathodic metal interface throughout the duration of the corrosion simulation, as shown in Fig. 4. (Note that the shape of the α-phase interface near the junction with the β-phase interface is slightly rounded due to the diffuse interface effect. This rounding can be reduced by employing numerical parameters that provide a thinner interface.) There appear to be two distinct regions, one near the anodic metal/cathodic metal interface, which is dominated by the accelerated corrosion, and the other characterized by almost constant corrosion rate, leading to nearly uniform corrosion. This indicates that there are at least two length scales associated with the corrosion process even when only migration is considered: a short length scale over which the corrosion current decays quickly, and another length scale (which is greater than the simulation domain size for the simulation shown in Fig. 3) over which the corrosion current persists.

2D simulations
Next, we investigated the effect of the initial area fraction of β on the corrosion rate by selecting 5 values of the β fraction: 5%, 10%, 15%, 20%, and 25%. Figure 5(a) shows the comparison of the evolution of the corrosion RCD integrated along the anodic metal/electrolyte interface (hereafter, integrated current density, or ICD) for all the cases, and the location of the corroding interface is compared at t = 2480 s in Fig. 5(b). Two trends can be observed in these results. First, the ICD increases with the initial β fraction for a fixed time, including at t = 0 . The fact that increasing the β region exposed to the electrolyte increases the ICD indicates that the corrosion process here is limited by the cathodic reaction. This was confirmed by conducting additional simulations and analysis covering a larger range of volume fractions to identify regimes that are limited by the cathodic reaction and the anodic reaction, as well as the crossover regime. Second, the ICD increases with time for a given initial β fraction. This occurs because the corrosion of the α phase exposes more β phase to the electrolyte, leading to a greater β/electrolyte interfacial length (equivalent to area in 3D). This increase in β/electrolyte interfacial length also reduces the cathodic-reaction limitation and thus accelerates corrosion. This result illustrates the importance of microstructure; specifically, when the cathodic phase is elongated (e.g., in eutectic microstructure), the increase in the exposed surface of the cathodic phase will lead to an increase in corrosion rate; thus, any corrosion will accelerate the subsequent corrosion in this limit.

3D simulation
The 3D capability of the application is demonstrated in Fig. 6, where the evolution of the metal/electrolyte interface and the ICD are shown for the 3D system from Fig. 2. The morphological evolution of the corroding surface is similar to that of the 2D simulations, i.e., a fastest corrosion is observed close to the α/β interface, leading to a circular trench near the β (cathodic) phase. However more uniform corrosion proceeds  away from the interface, again suggesting an existence of two lengths scales discussed above. Furthermore, the overall corrosion rate increases with time (as shown in Fig. 6(b)), again suggesting that growth in the exposed surface area of β , as well as that of α due to its nonplanar evolution, leads to further acceleration of corrosion. Since such evolution is strongly dependent on the microstructure, the design of alloys for improved corrosion resistance should consider not only the alloy composition but also the processing and the resulting microstructure. We note that small fluctuations in Fig. 6(b) result from the limited mesh resolution used for the simulation to reduce the computational cost. They can be eliminated by increasing the mesh resolution and reducing the time step size; however, based on the 2D studies, the overall evolution is not affected by the fluctuations.

Conclusions
One of the critical foundations required in realizing the visions outlined in the Material Genome Initiative [23] is the software infrastructure. There is a growing number of open-source software, as well as commercial software, but for any given application, the choices are still limited. In this prospective article, we focused on simulation of corrosion. We first provided an overview of the existing tools for simulating corrosion. We then summarized the features and capabilities of the open-source corrosion_microgalvanic [24] application developed within the PRISMS-PF framework. The application is capable of simulating both 2D and 3D systems and exhibits super-linear parallel performance for sufficiently large simulations. To show the utility of this application, we employed it to investigate the effect of the β (cathodic) phase fraction on the corrosion rate, where the microstructure was assumed to be of lamellae (in 2D) or a cylinder embedded in the α matrix (in 3D). It was found that the corrosion rate in the systems examined was limited by the cathodic reaction, and thus increasing the β/electrolyte interfacial length in 2D or area in 3D led to an increase in the corrosion rate. Similarly, an increase in the corrosion rate with time was observed for a given simulation as the α phase corrodes and exposes more of the β phase to the electrolyte. The distribution of the corrosion current, as well as the morphology of the corroded surface, suggests the existence of two length scales, one that is small and is associated with a larger corrosion current near the α/β interface, and the other on a larger scale, greater than the simulation domain size we examined. Furthermore, it was found that both the 2D and 3D systems exhibit similar corrosion behavior. Detailed analyses of these results, as well as an expanded parametric study to develop a quantitative understanding of the effect of microstructure on microgalvanic corrosion in Mg alloys, are beyond the scope of this paper and will be left for a future publication. While corrosion_microgalvanic application currently does not account for diffusion of ions in the electrolyte phase, it is trivial to add it in the same manner as the original corrosion application. Moreover, new features can be added to enhance the applicability of the corrosion-related problems. For example, deposition of corrosion products from cathodic reactions, which has an important role in self-limitation of corrosion, could be implemented by the same scheme used to evolve an interface during corrosion. These additional developments are also left for future work.

Appendix
Effect of ζ Figure 7 shows the effect of ζ on ξ 1 for the 2D simulation with 20% β phase fraction. As can be seen, reducing ζ value causes irregularities in ξ 1 , which affects i rxn . This effect necessitates a careful selection of the ζ value; it must be small enough but not too small. Our extensive testing showed that ζ = 10 −3 results in smooth evolution of ξ 1 .
See Table III.

Code availability
The application files can be found in the PRISMS-PF GitHub repository (https:// github. com/ prisms-center/ phase Field).

Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.

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.