Dendrite operating state in directional solidification of AlCu binary system: numerical benchmark test with the OpenPhase software

A scientific benchmark test is carried out for a multi-phase-field model with double-obstacle potential by performing three-dimensional simulations of dendritic growth under directional solidification. The effects of key numerical parameters of the multi-phase-field model such as numerical resolution and interface width on the dendrite tip operating state are studied, optimal parameter values are set, and where the operating state becomes independent of varying these parameters is elaborated. Some uncertainties in the proper choice of effective interface mobility in the thin-interface limit are discussed and a pragmatic solution is adopted. The binary alloy Al-Cu with 4 at.% Cu is chosen as the material system because it has been used in many previous experimental and numerical studies to investigate dendrite morphology under directional solidification. The recently developed sharp phase-field model by Finel and colleagues is adapted to the double-obstacle potential function and included in the benchmark test. It is shown how the sharp phase-field model helps in achieving agreeable convergence with larger discretization, thereby reducing the computational cost significantly. The benchmarks are performed using the OpenPhase software.


Introduction
Usually, metallurgical processes begin with solidification, the rate of which is determined by the alloy system and the applied cooling conditions, with the final microstructure evolving in a self-organizing manner.Nowadays, the availability of advanced computational facilities means that microstructure modeling and simulation are attracting tremendous interest from both academia and industry, and there are currently numerous models for microstructure modeling across length scales [1][2][3][4].Thanks to these models and modern computer architecture, our understanding of solidification under various conditions has evolved significantly over the past few decades.However, there remains a fundamental and long-standing problem that is the length scale at which solidification takes place [5].
Free dendritic growth in melts with lower or higher undercooling has been the subject of interest for decades.In a first attempt to model dendritic solidification, Ivantsov [6] attempted to solve the diffusion problem by describing both the thermal and solutal diffusion fields in the bulk liquid ahead of the moving interface with a paraboloid of revolution assuming an isothermal and iso-solutal interface.By dropping the Gibbs-Thompson effect (i.e., that curved interfaces must be at a lower temperature than the actual melting temperature), it was possible to find the paraboloid of revolution moving with constant velocity as a solution of the solidification problem.Unfortunately, this leads to loss of length-scale information, with only the product of the dendrite tip's curvature radius R and solidification speed V given as a function of the interface undercooling.
Later, Langer et al proposed replacing the conventional maximum-velocity condition with a marginal stability criterion involving VR 2 equaling a constant [7].Based on Ivantsov's solution combined with this marginal stability condition, a series of models of free dendritic growth-such as the Lipton-Glicksman-Kurz (LGK) model-was proposed to obtain individual relations for the dendrite operating state [8][9][10].Trivedi et al extended the marginal-stability model to rapid solidification [11], but although the models of free dendritic multi-component alloy thermodynamics, transformation strains, and fluid flow), then one must be aware of the correct bounds for selecting numerical parameters.We offer this benchmark to the community to check their own implementations.

Review of theories of free dendritic growth in binary alloys
There are multiple reviews of models of the evolution of basic free dendritic growth [3,4,32], and all those reviewed models are aimed at an accurate definition of the dendrite operating state, i.e., the relationship among far-field undercooling, tip velocity, and tip radius.In doing so, they predict the tip velocity, tip radius, and tip undercooling of an unbranched single needle-shaped dendrite growing into a melt with bulk composition c ∞ and total undercooling ΔT defined as where ΔT t , ΔT c , and ΔT r are the thermal, solutal, and capillary undercooling contributions, respectively.
As mentioned in section 1, Iantsov found a solution to the diffusion problem by modeling the moving boundary layer as a paraboloid of revolution assuming an isothermal and iso-solutal interface.Ivantsov's solution is given by where Ω is the dimensionless supersaturation, E 1 (Pe) is the exponential integral function, and Iv(Pe) is the Ivantsov function [6].
Later, Lipton et al [7] proposed the LGK model by combining Ivantsov's solution with marginal stability theory.According to this model, the total undercooling ahead of the dendrite tip is given by where L is the latent heat of fusion, C p is the heat capacity, m l is the liquidus slope, C 0 is the nominal composition, k is the equilibrium partition coefficient, τ is the Gibbs-Thomson coefficient, R is the radius of the spherical interface, and P c and P t are the solutal and thermal Peclet numbers, respectively.The second equation needed to identify the dendrite operating state uniquely is formulated from the condition of marginal stability as given by Dd VR

2
, 5 0 2 where σ * is the newly introduced interface stability parameter in this model, D is the coefficient of diffusion in the liquid, and d 0 is the capillary length.Furthermore, σ * = 0.025 in this model is identified as the stability constant and is independent of the material system.The dendrite operating state is obtained by solving equations 1-5 with an iterative scheme such as the Newton-Raphson one.The LGK model is suitable for solidification under local equilibrium conditions with low undercooling.
Although marginal stability theory predicts the correct scaling 5, it misses an important aspect, i.e., interface energy anisotropy.The weak direction or the interface stiffness σ * at the dendrite tip leads to a deviation from the parabolic shape assumed by Ivantsov and thus a small protrusion at the dendrite tip [33], and this triggers the growth of the dendrite in the preferred direction.Theoretically, this problem is solved by microscopic solvability theory, which considers interface energy anisotropy, and the scaling behavior is given by an expansion for low Peclet numbers [34,35], i.e., where σ * is again the stability parameter introduced in marginal-stability-based models, but in microscopic solvability theory it is treated as a function of crystal anisotropy.For lower Peclet numbers and crystal anisotropy, the relationship between σ * and the interface energy anisotropy ò 4 is given by , where σ 0 is a constant pre-factor.A full solution is only possible by numerical simulation, and to date the phase-field approach is the only numerical technique able to cope with this problem.

Thin-interface phase-field model with double-obstacle potential
In OpenPhase, the multi-phase-field model is implemented with a double-obstacle potential function.In the multi-phase-field formulation, each grain α is distinct from the others with respect to its orientation or phase (or both) and is represented uniquely by the phase-field parameter f α , although we restrict ourselves to the twophase system of solid and liquid.See elsewhere for the complete details of the model formulation [18].Here, we provide the key equations governing the evolutions of the phase fields f and composition c, as well as those for the driving force, interface stiffness, interface anisotropy, and interface mobility taken from references [18] and [33].The interface stiffness γ * is defined by equation 10 with the average interface energy γ 0 and the interface stiffness anisotropy ò, which in the case of fourfold symmetry is related to the interface energy anisotropy as ò = 15ò 4 .Δg is the driving force for transformation defined by equation 9: In equations 7-11, f is the phase-field parameter of phase α and is defined as 0 f α 1, M eff is the effective interface mobility calculated from equation 11, η is the interface width, and n x , n y , and n z are the components of the normal vector calculated from the phase-field variable with Furthermore, c is the mixture composition with c l and c s being the phase compositions in the liquid and solid phases, respectively, D l is coefficient of diffusion in the liquid, Δg is the driving force, ΔS is the entropy difference, m l is the liquidus slope, μ is the physical interface mobility, and T m and T are the melting and interface temperatures, respectively.
A numerical technique that we used to easily obtain suitable initial transients was to apply a soft driving force limited according to The allowed driving force is given by the highest interface stabilizing force g max ps h D = .After the initial transient, the driving force is no longer limited.
In addition to the model equations 8 and 9, one must specify how the solute partitioning is modeled, i.e., calculating the phase concentrations c l and c s from a given physical composition c = fc s + (1 − f)c l .This is done as an approximation by demanding that at every point in space, the partitioning follows the equilibrium partitioning, i.e., where C s eq and C l eq are the equilibrium phase concentrations of the solid and liquid phases, and m s and m l are equilibrium slopes of the solidus and liquidus lines, respectively.Finally, the temperature evolution is modeled via the frozen temperature profile with a constant cooling rate and an imposed thermal gradient.

Sharp phase-field model for double-obstacle potential
Usually, phase-field models are derived by defining a free energy functional on the space of all possible phasefield functions {f} and then applying a variational ansatz to derive the equation of motion for f [18].Depending on the exact analytical form of the potential function f, the phase-field equation is solved by a traveling-wave solution.If the double-well potential f is used, then the traveling-wave solution will be of a tangent hyperbolic form.However, if the double-obstacle potential is used, then the traveling-wave function takes the form The simple form given above is valid if neither a driving force nor a diffusing entity such as heat or solute concentration is included in the model.If a driving force is included, then the offset parameter x 0 can become time dependent.
The SPF method by Finel [23] is based on two ideas: (i) start the model derivation with a discretized description of the free energy functional and the Laplacian term; (ii) reverse the two derivation steps-first, the desired traveling wave-solution is defined, then the potential function f is determined, which generates the desired traveling-wave function as a solution in the phase-field equation, even though a discretized version of the Laplacian is used.See elsewhere for a detailed explanation [23].Here, this procedure is carried out for the cosinusoidal traveling-wave solution in equation 15, which leads to a modified version of the double-obstacle potential from equation 14.
Let the discrete version of the Laplacian operator be where the indices p, q, and r denote steps along the main directions of the numerical grid toward neighboring nodes, Δx is the width between adjacent nodes, and the coefficients S pqr are determining the a Laplacian stencil.
With this, the discrete version of the free energy functional is written as For a simple cubic grid, the shift vector  r pqr can be written as  r x p q r , , with p, q, r ä {0, − 1, 1}.The discretized version of the variational derivative gives The potential function f is now determined by demanding that the terms in parenthesis in equation 19 vanish, if f is a plane traveling wave according to equation 15 in a designated direction  n .The direction vector  n is normalized.This means that the three-dimensional (3D) phase-field function is assumed to behave as • 〉 is the standard scalar product.This gives ( ) ( ) ( ) of a plane traveling wave with direction  n at the neighboring grid points, where at the central point the traveling wave has value equal to  r ( ) Note that the threshold values u pqr f and l pqr f are particularly important for recognizing when the adjacent nodes are part of the bulk region and not the diffuse interface.Because the cosinusoidal traveling wave is non-analytic at this transition and has constant values throughout the bulk phases, this differentiation is theoretically important and has been shown in additional simulations (not reported here) to have a significant impact on the ability of f to cure numerical interface locking.

Numerical implementation
In the OpenPhase library, the phase-field model equations 7 and 8 are solved via an explicit first-order time integration scheme (forward Euler scheme).The spatial discretization of the Laplacian operator is achieved via the isotropic stencil [37] with stencil components defined by S p q r p q r p q r p q r 64 15 if , , 0, 0, 0 , where p, q, r = − 1, 0, 1.
The key advantage of the phase-field model implementation in OpenPhase is its rigorous consideration of the interface properties.Furthermore, as reported by Tegeler et al [31], maintaining the steady-state phase-field contour across the interface at all times during the simulation ensures constant width of the interface and results in accurate consideration of the interface properties that are key to this benchmark study.

Simulation setup
The numerical parameters and physical constants used in this study are given in table 1.To reduce the computational time and considering the symmetry of the dendrite, only a quarter of it simulated in the 3D domain, with zero-flux boundary conditions on all walls except for the top wall, where one-dimensional (1D) diffusion extension is introduced.The extension length is four times that of the 3D simulation domain.Initially, the simulation domain is filled with liquid of the model binary alloy Al-Cu with an initial composition c 0 of 4 at.% Cu.A small solid seed is placed at the lower corner of the simulation domain with coordinates (0, 0, 0), and the crystal growth axes are aligned with the coordinate axes.Because of the imposed interface energy anisotropy and thermal gradient, the dendrite grows in three dimensions with its primary arm growing along the z axis and its secondary and tertiary arms growing in the xy plane.The solidification conditions relate to those in a previous study of one of the present authors on primary spacing selection in directional solidification [33].
The applied temperature profile is chosen such that at the beginning of the simulation, the solid seed is at the initial temperature T 0 from table 1. Together with the chosen temperature gradient ∂T/∂z and the globally applied cooling rate  T , which are also listed in table 1, the steady-state growth velocity is already determined via  v T T z ( ) = ¶ ¶ .A necessary condition for the steady state is therefore that the velocity obtained from the front tracking should match the velocity imposed via ∂T/∂z and  T .This was checked and fulfilled for every simulation shown later.The unknown quantities of the steady state obtained from the simulation results are therefore the radius of curvature R and undercooling ΔT at the dendrite tip, which are evaluated for every time step once the dendrite morphology is stabilized.
A moving-frame technique is applied such that when the solidification front reaches one third of the simulation domain length, every value (e.g., phase field, composition, and temperature) is shifted down along the growth direction (the z axis).With this technique, the growing dendrite tip is always kept around the same position within the computational domain, which allows for longer computational time in a smaller domain to achieve the steady state.
The radius of curvature at the dendrite tip is evaluated by first extracting the f = 0.5 contour and then fitting a parabola along one of the principal directions, the x axis.Let the fitted parabola be of the form f (x) = ax 2 + bx + c, then according to differential geometry, the radius of curvature is given by R = 1/(2a) [38].

Results and discussion
4.1.Dendrite operating state obtained from marginal stability and microscopic solvability theories As explained in section 2, within a given material system, possessing knowledge of any one of the three components of the dendrite operating state namely V, R, and ΔT allows us to predict the remaining two components using free dendrite growth models.We have implemented both the LGK model and the microscopic solvability model (with the stability parameter obtained from table VIII in reference [26]) in the Python programming language to predict the dendrite operating state.The plots in figure 1 show the predicted tip radius and solute undercooling for given solidification velocity by the LGK and microscopic solvability models.For the chosen material system and a steady-state velocity of 40 μ/s, the LGK model, which is based on the marginal stability theorem, predicts a tip radius of 4.9 μm and a solutal undercooling of 3 K.On the other hand, the microscopic solvability theory simplified for three dimensions and lower Peclet numbers predicts a tip radius of 3.96 μm and an undercooling of 2.5 K.

Verification of symmetry boundary conditions
To determine whether imposing symmetry boundary conditions would affect the dendrite tip operating state, we compared the dendrite operating state by performing simulations of a full, half, and quarter dendrite.The numerical parameters used for simulating a quarter dendrite were those given as the reference simulation setup in table 1.To simulate a half dendrite, we doubled the dimensions of the simulation box in the x direction and shifted the position of the nucleation seed to (70, 0, 0).For the full dendrite, we doubled the system dimensions in both the x and y directions and shifted the position of the nucleation seed to (70, 70, 0).
The dendrite morphology obtained from each simulation is shown in figure 2. In figure 3, we superimpose the dendrite contours from each simulation to check for deviation due to the imposed symmetry condition; in this figure, the red, green, and black curves represent the dendrite contours of the full, half, and quarter dendrites, respectively.The dendrite contours are superimposed by keeping the contour of the half dendrite in place and sectioning the full dendrite by the x axis while mirroring the quarter dendrite with respect to the z axis.Here, the dendrite contour obtained from the quarter dendrite is shifted away from the center by half a grid point because of the imposed numerical conditions.As can be seen, the contours at the tip location for all three dendrites match very well, and the results of this test lead to the conclusion that the symmetry condition does not significantly influence the dendrite tip operating state.Therefore, we continued the benchmark simulations using a quarter dendrite with the nucleation seed placed at (0, 0, 0).

3D simulation results of dendritic growth under directional solidification
Before discussing the benchmark test results, we present simulation results obtained for the evolution of the dendrite in three dimensions for the reference setup given in table 1.The growth of the dendrite from the initial transient to the final steady state is shown in figure 4. The results obtained for the evolution of the phase fields show that during the initial transient state, the dendrite develops secondary dendritic arms.However, once the steady state is reached, the growth is concentrated along the primary growth axis.The secondary arms are not seen at the later time steps because of the activation of the moving frame along the growth axis.
Figure 5 shows the distribution of solute concentration in the whole system under steady-state dendritic growth.Figure 5  far-field diffusion, a 1D diffusion extension is implemented, as explained in appendix A. The variation of solute concentration in the solid phase is not significant because diffusion in the solid phase is neglected in the simulations.
Based on the steady-state simulation results, we calculated the dendrite tip velocity, radius, and undercooling, as described in appendix B. The dendrite tip velocity and contour in the steady state are shown in figures 6 and 7, respectively.As shown in figure 6, the dendrite tip grows initially with a solidification velocity that is much lower than the pre-set value, but it increases gradually until it reaches the steady-state velocity.Depending on the numerical parameters and initial undercooling, the tip velocity can also exceed the steadystate value, but it eventually returns to the steady-state velocity.For the reference setup, the dendrite tip achieves Figure 3. Superposition of dendrite contours obtained from full, half, and quarter dendrite.Red: dendrite contour of full dendrite sectioned parallel to x axis.Green: dendrite contour of half dendrite.Black: dendrite contour of quarter dendrite mirrored with respect to z axis to construct a half dendrite.To match the dendrite contour obtained from quarter dendrite simulation with dendrite contours obtained from half and full dendrite simulation, it's dendrite contour must be shifted away from the center by half a grid point due to the imposed numerical conditions.the preset solidification velocity of 40 μm/s after a simulation time of 50 s, after which it grows at a constant velocity, and this state is defined as the steady state.
Figures 7(a) and (b) show the dendrite contour in the steady state constructed by tracking the f = 0.5 interface contour over time.The green curve is the obtained dendrite trunk contour, whereas the blue curve is a parabola fitted to the tip portion of the obtained dendrite contour.The details of the curve fitting are given in appendix B.
Figure 8 shows the calculated solute undercooling and tip radius of the fitted dendrite contour.As shown here, the tip radius had a minimum before it attained the steady-state value of 4.8 μm.Similarly, the phase-field simulations with chosen grid resolution and interface width predict the solute undercooling to be 3.60 K.

Selection of interface mobility
Before starting the benchmark test, we must ensure that the simulation setup satisfies the diffusion-controlled limit of solidification.We use the interface limit of Karma's model [25] adapted to the double-obstacle potential in [18] and with the anti-trapping current (the last term in equation 8) and the effective interface mobility  [equation 11].For the chosen Al-Cu system with 4 at.%Cu and an interface width of η = 2.5 μm, we calculate an effective interface mobility of 4 × 10 −10 m 4 /(J•s).However, our studies show that even though the dendrite tip attains the steady state with this mobility, the calculated tip radius is higher than expected [12,30,33].To understand the influence of interface mobility on the tip radius, we performed multiple simulations with the SPF model adapted to the double-obstacle potential and with interface mobilities of 4 × 10 As depicted in Figure 9, when the effective interface mobility is set to its theoretical limit, the simulation reaches a steady state, yielding a dendrite tip radius of 6.5 μma value significantly greater than the anticipated tip radius.Nonetheless, with the increase of effective interface mobility, the dendrite tip sharpens, resulting in a reduced tip radius.Notably, our simulations indicate that the influence of interface mobility on dendrite tip radius variation diminishes once the ratio of M eff /M th is set to 3. A similar convergence tendency is observed for the standard PF simulations with interface widths of 3.5-5 μm.Note that the scheme applied here does not  concentrate the driving force for the growth of the solid phase to the center of the interface (f ∼ 0.5), as in Karma's original work.Also, an averaging procedure for the driving force in the normal direction through the interface is applied [39].All of this might affect the fully quantitative applicability of the scheme, and so further work is needed to clarify the observed discrepancy, e.g., a local mobility correction as proposed in reference [40].
For the present work, the results of varying the effective mobility lead to the conclusion that good convergence is achieved with a mobility that is three times that predicted from equation 11, and this was used throughout the benchmark test.

Benchmark test
To investigate how key numerical parameters of the phase-field model affect the dendrite operating state, we conducted several simulations in which we varied the grid spacing and interface width.In each simulation, the  simulation space was sampled at the sampling points described in figures 11 and 12 and table 2. The sampling points were chosen so that the interface width was always within the numerical limits suggested by previous studies [30,33].As described in section 3.1, each simulation began with a solid seed placed at (0, 0, 0), and a steady-state velocity of 40 μm/s imposed through a constant thermal gradient and cooling rate.The dimensions of the simulation domain were adjusted by increasing or decreasing the number of grid cells by multiples of Δx while maintaining an overall domain size of 70 μm ×70 μm ×300 μm in each simulation.We obtained the steady-state velocity and parabolic tip radius for each simulation as described previously, then we plotted the calculated dendrite tip radius against the varying numerical parameters.
As shown in figures 11 and 12 and table 2, benchmark simulations were performed for grid resolutions of Δx = [0.5, 0.75, 1] μm.For each grid resolution, the number of grid cells in the interface was varied in the range of 2.5 to 6, which resulted in an absolute interface width of 2-5 μm.Furthermore, for each configuration of Δx and η, phase-field simulations were carried out with the standard PF model and the SPF model adapted to the double-obstacle potential.Because of the huge computational cost associated with a grid resolution of Δx = 0.5 μm, only one simulation was performed with η = 3 μm.

Benchmark test results
The benchmark test results are summarized in table 3 and figures 13-15.As highlighted in table 3, the dendrite tip radius, velocity, and temperature results enable the complete characterization of the dendrite operating state,  facilitating the calculation of remaining key parameters like tip undercooling, solutal Peclet number, and stability constant.For the chosen set of data points, the observed variation in the calculated tip radius remains For all the simulations the variation among the obtained tip radii is around 10% (excluding the case when the interface width is already 5 μm because it is not possible to resolve a tip radius of below 5 μm with such high interface width).Similarly, the overall variation in the dendrite tip temperatures and resulting tip undercooling are also within a close range of less than 10%.In theory, enhancing the convergence of the obtained tip radii is achievable, albeit demanding a finer grid resolution, such as Δx = 0.5 μm.Unfortunately, an examination of the computational effort needed for both Δx = 1 μm and Δx = 0.5 μm reveals that achieving such a benchmark with a finer grid resolution is practically unfeasible.
The values of the solutal Peclet number, calculated using P c = VR/(2D) where V and R are steady-state tip velocity and tip radius from the phase-field simulations, fall within the range of 0.025 to 0.033, closely aligning with those reported in previous studies [26,30].Additionally, examining table 3 reveals that the values of the stability constant, determined through σ * = 2Dd 0 /(R 2 V ), surpass the stability constants applied in both marginal-stability and microscopic-solvability investigations.As referenced in [26], an anisotropy exceeding 3% leads to notable deviations between the phase-field outcomes and predictions of free dendritic growth models.
The computational effort details provided in the final column of table 3 illustrate that when utilizing the standard thin interface phase-field model, a 1-second solidification process simulation with a coarser grid resolution (e.g., Δx = 1 μm) demands approximately 5 hours.Remarkably, by halving the grid resolution, the simulation duration dramatically extends to 209 hours on the same machine.Further analysis of the computational details leads us to the conclusion that opting for a lower grid resolution substantially increases computational effort for the same interface width.However, the computational load reduces significantly through the adoption of the SPF model, enabling the use of narrower interface widths even with coarse grid resolutions.For instance, when considering a simulation with η = 3, the conventional PF model necessitates Δx = 0.75 μm, whereas SPF permits the use of Δx = 1 μm, effectively reducing the computational burden by over 50 times.
Furthermore, figures 13-15 illustrate how the obtained tip radius varies concerning the essential numerical parameters η, n, and Δx.In each figure, the solid lines represent results from the thin-interface PF model with double-obstacle potential, while the dotted lines depict those from the SPF model.Even though figures 13 and 14 contain the same information, it is worthwhile arranging the results with respect to different parameters because it is observed that even for the same interface width, the number of grid cells in the interface influences the obtained tip radius.From figure 13(a), we can observe that for grid resolutions of Δx = 1 μm and Δx = 0.75 μm, the variation in the tip radii increases as the number of grid cells in the interface increased from 2.5 to 6.To avoid interface pinning, a minimum of three and half cells was maintained in the standard PF simulations, whereas with SPF it was reduced to two and half grid cells.However, as seen in figure 13(a), the tip radius obtained from SPF for n = 2.5 shows significant difference between the chosen grid resolutions of Δx = 1 μm and 0.75 μm.Further, as n increased beyond 4, for both PF and SPF models, there was significant difference between the tip radius obtained from grid resolutions of Δx = 1 μm and Δx = 0.75 μm.Similarly, figure 13(b) shows that for a given interface width and n value, the tip radii obtained from both phase-field models agree well.For Δx = 1 μm, the number of interface grid cells were restricted to 5 as beyond that it would result in a interface width of higher than 5 μm which cannot resolve a tip radius of less than 5 μm.
Furthermore, figure 13 shows that the tip radius has local maxima at integer values of n (n = 4, 5, and 6) and fluctuates when n is changed from integer to non-integer values.The observed deviation in the tip radius between integer and non-integer numbers of grid cells in the interface is due to the applied driving force being averaged in the interface.The driving-force averaging is always performed over an integer number of grid cells, which affects the driving-force values for non-integer numbers of grid cells in the interface, and consequently a small deviation in the tip radius is observed.

Influence of interface width (η) on dendrite tip radius
The variation of tip radius with interface width η is plotted in figure 14: figure 14(a) shows how the tip radius varies with η at a fixed grid spacing, while figure 14(b) shows the same but at a fixed n value.From figure 14(a), it is clear that the variation of the tip radius with interface width is more scattered than that with n.Further, figure 14(a) shows a significant divergence between the tip radii obtained from the standard PF model with Δx = 0.75 μm and 1 μm at a given interface width.This can be attributed to the integer versus non-integer numbers of grid cells in the interface as mentioned above.A similar trend is observed with the SPF simulation results.The oscillations in the tip radius for Δx = 1 μm are higher than those for Δx = 0.75 μm.However, with the SPF model, convergence can be achieved even with interface widths below 3 μm.Furthermore, the computational cost required to achieve the steady state with the SPF model for Δx = 1 μm and η = 2.5 μm is almost a third of that required to do so with the standard PF model in the same configuration.
Furthermore, figure 14 leads to the conclusion that the calculated R tip is always higher than or equal to the interface width.Once again, this reiterates the fact that when using the phase-field method, one cannot resolve length scales smaller than the chosen numerical interface width.

Influence of grid resolution on dendrite operating state
Figure 15 shows how the tip radius varies with the chosen grid resolution.Clearly, even with the same grid resolution, the dendrite operating state can vary significantly depending on the number of interface grid cells and resulting interface width.Furthermore, figure 15 shows that as the grid becomes coarser, the deviation in the obtained tip radius increases for the same number of grid cells in the interface.This is because increased Δx results in increased interface width, which reduces the ability to resolve finer length scales.If a sufficient number of grid cells is maintained in the interface and its absolute width is below the critical limit, then the grid resolution does not significantly influence the dendrite operating state.However, as Δx is increased from 0.5 μm to 0.75 μm and then to 1 μm, the computational time decreases significantly, and the same trend is observed for the results obtained from the SPF model.

Summary and conclusions
A scientific benchmark test was carried out for a multi-phase-field model with double-obstacle potential (implemented in the OpenPhase software library) simulating 3D dendritic growth under directional solidification.As a model alloy, we used Al-Cu containing 4 at.%Cu.Considering the dendrite's symmetry, a quarter dendrite was simulated in 3D space with given tip velocity, and the steady-state dendrite tip radius was calculated by fitting a parabola to the dendrite contour.To evaluate the benchmark, multiple simulations were carried out by varying key numerical parameters such as grid resolution (Δx), number of grid cells in the interface (n), and interface width (η).From the simulation results, we draw the following conclusions.
• The results obtained for the dendrite tip radius are within a tolerance of 1 μm, which suggests the robustness and reliability of the PF implementation in the OpenPhase software library.
• The interface width determines the accuracy of the obtained dendrite tip radius, while the number of grid cells in the interface determines the computational cost required to achieve that accuracy.
• For an interface width below the smallest length scale in the simulation, the spatial discretization Δx does not affect the simulation results significantly.
• The driving-force averaging implemented in the OpenPhase library results in a fluctuating tip radius when the number of grid cells in the interface changes from integer to non-integer values.
• The effective interface mobility obtained from Karma's thin-interface asymptotics and applied to the current phase-field model with double-obstacle potential is of the correct order of magnitude but is too low for pure diffusion control.Varying the interface mobility shows that convergence is achieved when the effective interface mobility is taken as three times the theoretical value.This also highlights the need for efficient local mobility correction schemes.
• Finally, the new SPF model [23] was adapted to the double-obstacle potential function and included in the benchmark test.By using the new SPF model, the simulation grid resolution can be halved, resulting in the number of grid points being reduced by a factor of 2 3 in three dimensions.Together with the larger time step made possible by a coarser spatial resolution, an overall theoretical improvement in the calculation time by a factor of more than 10 can be expected.
Finally, note again that for the Finel SPF model, we did not re-evaluate Karma's thin-interface limit.However, it is safe to assume that the phase-field part remains unchanged because Finel's approach only remedies the numerical discretization errors.Knowing this, one should also consider those discretization errors for the solute diffusion profile.A future task will be to consider this in a similar manner as for the phase-field profile, which becomes particularly important if one goes to the low discretization width of two and a half grid cells within the interface.

Figure 1 .
Figure 1.Dendrite operating state (tip radius and tip undercooling for given steady-state velocity) calculated from models of free dendritic growth based on marginal stability theory (LGK model) and microscopic solvability theory.
(a)  shows the distribution of solute composition in the simulation domain, figure 5(b) shows the solute redistribution at the solid-liquid interface, and figures 5(c) and (d) show the distribution of solute concentration along the primary growth axis as line scans.To facilitate free dendritic growth and simulate

Figure 2 .
Figure 2. Convergence test for symmetry boundary conditions: dendrite morphologies obtained for quarter, half, and full dendrite.

Figure 4 .
Figure 4. Evolution of quarter dendrite in 3D system at different time steps: (a) start of solidification, with solid seed visible in lower left corner of box; (b) dendrite contour at 5 s; (c) dendrite morphology at 10 s, with growth along the z axis dominating that along the x and y axes; (d) dendrite morphology at 50 s, with steady state established and dendrite morphology stable thereafter.

Figure 5 .
Figure 5. Evolution of solute concentration in solid and liquid phases: (a) concentration of Cu in whole system (for the nominal composition of 4 at.%Cu, the maximum composition of 4.69 at.%Cu is observed in the liquid and the minimum composition of 0.64 at.%Cu is observed in the solid); (b) maximized view of Cu concentration distribution in dendritic growth region; (c) line scan of Cu composition along z axis at x = 0 and y = 0; (d) line scan of Cu composition along z axis at solid-liquid interface.

Figure 6 .
Figure 6.Evolution of dendrite tip velocity over time.Pink: calculated tip velocity.Cyan: preset solidification velocity, also identified as steady-state velocity.

Figure 7 .
Figure 7. (a) Dendrite contour in steady state.(b) Magnified view of dendrite tip contour highlighting parabolic fitting.

Figure 8 .
Figure 8. Dendrite tip radius and tip undercooling calculated from steady-state dendrite contour obtained from phase-field simulations.

Figure 9 .
Figure 9. Results of convergence tests for interface mobility: obtained dendrite tip radius vs. time for selecting different interface mobilities.

Figure 10 .
Figure 10.More results of convergence tests for interface mobility: (a) dendrite contours obtained for different interface mobilities; (b) magnified view of dendrite tip contour highlighting parabolic fitting.

Figure 11 .
Figure 11.Sampling points displayed as interface grid cells vs. grid spacing.

Figure 12 .
Figure 12.Sampling points displayed interface grid cells vs. interface width.

4. 5 . 1 .
Influence of number of interface grid cells (n) on dendrite tip radius Figure13shows how the tip radius varies with the number of grid cells in the interface: figure13(a)shows how the tip radius varies with n at a fixed grid spacing, while figure 13(b) shows the same but at a fixed interface width.

Figure 13 .
Figure 13.Variation of tip radius with the number of grid cells in the interface (n): (a) Variation of tip radius with n at a fixed grid spacing (b) Variation of tip radius with n at a fixed interface width.Solid lines: tip radius obtained from standard thin-interface PF model.Dotted lines: tip radius obtained from SPF model.

Figure 14 .
Figure 14.Variation of tip radius with interface width (η): (a) Variation of tip radius with η at a fixed grid spacing (b) Variation of tip radius with η at a fixed number of grid cells in the interface.Solid lines: tip radius obtained from standard thin-interface PF model.Dotted lines: tip radius obtained from SPF model.

Figure 15 .
Figure 15.Variation of tip radius with spatial resolution Δx.For Δx = 0.5, only one simulation is performed.For Δx = 0.75 and for Δx = 1.0, multiple simulations are performed by varying η and n.

Table 1 .
Numerical and physical parameters of benchmark test for directional dendritic solidification.

Table 2 .
Sampling points selected for numerical benchmark.

Table 3 .
Benchmark test results: dendrite tip radius, temperature, and undercooling, solutal Peclet number, and stability constant obtained for chosen sample points.: simulation duration in hours to run the solidification process of 1 second.All the simulations are carried out on a single node, equipped with 12 cores, and each core has 4GB of memory. *