Effect of Mass-balance Error on Evolution of Solid-liquid Interface during Solidification and Melting in Cellular Automaton Simulations

Abstract

To quantitatively predict the microstructure using a cellular automaton (CA) model, problems such as the calculation accuracy of the interface curvature and occurrence of mass-balance errors must be overcome. In this study, we focus on mass-balance errors, develop a unique CA model that considers the solid–liquid interface movement associated with solidification and melting, and investigate the effects of mass-balance errors on microstructural evolution. Mass-balance errors occur in CA models based on a two-domain approach that separately defines the solute concentration distribution in solid and liquid phases. We identified the causes of mass-balance errors and quantitatively estimated them, and then conducted a simulation of the microstructural evolution of Al-4.5 mass%Cu alloy during isothermal undercooling and continuous cooling with and without mass-balance error corrections. During isothermal undercooling, if the mass-balance error is not corrected, the solute concentration increases during the coarsening of the dendrite arms, and the increase is greater with increasing undercooling, which causes the solid–liquid interface to melt despite the isothermal process. During continuous cooling, if the mass-balance error is not corrected, the solute concentration decreases during arm coarsening in the middle of solidification. Under the simulation conditions, the absolute percentage error of the solute concentration relative to the initial composition exhibited maximum values of 1.31% and 14.9% during isothermal undercooling and continuous cooling, respectively. Under both conditions, the values are almost zero if the mass-balance error is corrected. Therefore, mass-balance errors can be appropriately corrected in the proposed CA model based on a two-domain approach.

1. Introduction

Microstructure formation during alloy solidification, which affects the formation of casting defects, such as segregation, porosity, and cavities, is an important process that must be elucidated.¹⁾ The primary phases of the microstructure formed during alloy solidification are dendrites.²⁾ Therefore, numerical simulation models based on phase-field (PF)^3,4,5,6,7,8) and cellular automaton (CA)^{9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30)} methods have been investigated to predict the morphology of dendritic growth. Currently, the PF method is widely used as a highly accurate simulation method for examining the evolution of microstructures during solidification.^31,32,33) It is particularly attractive for simulating complex shapes such as dendrites without using advanced numerical methods. This is because the solid–liquid interface is expressed based on a diffuse interface. However, this diffuse interface must be discretized with fine computational grids, which increases the computational cost. Moreover, the PF method requires asymptotic analysis in the derivation of parameters to satisfy the sharp interface concept in the diffuse interface, and it is not easy to understand the derivation.^5,6,7) The CA method incurs lower computational cost and is easier to understand than the PF method. However, the solid–liquid interface is expressed based on a sharp interface, which renders the calculation of the interface curvature challenging.

Microstructural evolution during solidification is a free-boundary problem. Therefore, the solid-liquid interface must satisfy the Gibbs–Thomson relationship. In the PF model, the interface parameters unique to the PF method, known as PF parameters, are derived via asymptotic analysis to accurately satisfy the Gibbs–Thomson relationship in the solid–liquid interface region as a diffuse interface. Thus, a quantitative PF simulation can be performed using the appropriate parameters.^5,6,7,8) In the CA model, the Gibbs–Thomson relationship is directly applied to a solid–liquid interface based on a sharp interface. Therefore, the interface curvature must be calculated, as it directly affects the accuracy of CA simulations. In CA simulations, the computational domain is generally partitioned into cubes in a three-dimensional (3D) system (squares in a two-dimensional (2D) system), which are known as “cells.” In single-phase solidification, the solid fraction, f_S, is defined as a parameter that provides the state of the cells, where f_S = 1 denotes a solid cell, f_S = 0 is a liquid cell, and 0 < f_S < 1 is an interface cell. One or two interface cells are located between a solid cell and liquid cell, and the interface curvature is calculated based on the solid fraction of the interface cell, surrounding solid cells, and liquid cell. However, numerical errors occur in the calculation of the interface curvature because only a few interface cells exist between solid and liquid cells.^9,10,11) To reduce the effects of mesh anisotropy and achieve quantitative simulations, various solid–liquid interface evolution algorithms^12,13,16) and numerical methods for curvature calculation^14,15) have been proposed. In recent years, the height function method, which is used in computational fluid dynamics, has been reported to demonstrate a high accuracy in calculating the interface curvature in 2D CA simulations.^14,15) Wei et al. proposed a quantitative CA model using the height function method and confirmed that the error due to mesh anisotropy is as small as that of the PF model based on a simulation of the equilibrium shape. Furthermore, 2D simulations of the growth of dendrites, lamellar eutectics, and irregular eutectics were performed using a quantitative model.¹⁴⁾

As mentioned above, an important issue in the CA model is the calculation of interface curvature with high accuracy. Additionally, in the CA model based on a sharp interface, the discontinuous distribution of solute concentration at the solid–liquid interface caused by solute partitioning must be addressed. This discontinuity was treated as the solute concentration in the interface cell. Depending on the manner in which it is addressed, mass-balance errors are expected to occur, implying that solute conservation in the system will not be satisfied. To the best of our knowledge, the treatment of the solute concentration in the interface cell can be categorized into two approaches: the two-domain approach^{16,19,20,21,22,23,24,25,30)} and the one-domain approach.^26,27,28,29) The former defines the solute concentrations in the solid and liquid phases separately based on the solid–liquid interface, which satisfies the local equilibrium assumption, whereas the latter connects the solute concentrations of the solid and liquid phases continuously and smoothly by defining the average solute concentrations of the solid and liquid phases based on the solid fraction of each cell.

Based on the sharp interface, the position of the solid–liquid interface was determined using the solid fraction of the interface cell. With the solid–liquid interface as the boundary, the solute concentration of the interface cell is differentiated between the solid- and liquid-phase sides. However, the interface position within an interface cell cannot be easily and precisely determined using numerical methods. Therefore, in the two-domain approach, the interface cell is regarded as a liquid phase,^{16,19,20,21,22,23,24,30)} and the solute concentration of the interface cell is set to the equilibrium concentration of the liquid phase based on the local equilibrium assumption. This approach can facilitate the management of the solute concentration in the interface cell and avoid significant changes in the solute concentration caused by solute diffusion between the interface and liquid cells, which is effective when the solid fraction of the interface cell is approximately unity. However, this approach calculates the solute diffusion without considering the solid fraction in the interfacial cells, which implies that the solute concentration is not conserved.^9,25,30) Meanwhile, in the one-domain approach, the solute concentration in the solid phase is divided by the equilibrium partition coefficient in the interface cell, thereby reducing the difference in concentration between the solid and liquid phases and resulting in a continuous solute concentration distribution between them.^27,28,29) Therefore, the solute diffusion equation does not include an explicit solute partitioning term and the solute concentrations are conserved. However, in this approach, the equilibrium partition coefficient must be constant; for example, in the case of binary alloys, the liquidus and solidus lines must be straight. In alloy systems with curved liquidus or solidus lines, microsegregation cannot be predicted quantitatively. Therefore, from the viewpoint of industrial applications of CA simulations, such as quantitative prediction of microsegregation, a CA model based on the two-domain approach is necessary, and the mass-balance error generated by the two-domain approach must be elucidated. However, the mass-balance error generated by this approach has not been sufficiently examined.

Moreover, most CA simulations focus on microstructure formation during the solidification process and tend to disregard the movement of the solid–liquid interface associated with melting. This implies that the Gibbs–Thomson relationship at negative curvature is disregarded. To quantitatively predict the dendritic morphology, the CA model must be able to manage the movement of the solid–liquid interface that accompanies melting and solidification. In the CA model based on the one-domain approach, several simulations that consider both solidification and melting have been reported.^27,28,29) This is likely because the movement of the interface is easier to manage because discontinuities in the solute concentration distribution are avoided. However, in CA models based on the two-domain approach, few simulations that consider melting have been reported. Therefore, there has been little discussion regarding the mass-balance error caused by the movement of the solid-liquid interface due to melting.

In this study, we developed a CA model that considers the movement of the solid–liquid interface during solidification and melting using a two-domain approach. Based on the developed model, we clarified the cause of the mass-balance error generated by the movement of the solid–liquid interface during solidification and melting and performed 3D simulations of dendritic growth during isothermal undercooling and continuous cooling under calculation conditions with and without corrections to the mass-balance error. Based on the obtained simulation results, we investigated the effect of the mass-balance error on microstructural evolution.

2. CA Model for Dendritic Growth

In this study, a CA model based on a two-domain approach that considers the movement of the solid–liquid interface during solidification and melting was developed. In the CA simulations, the computational domain was partitioned into cubic cells in three dimensions and three cell states were defined: liquid, solid, and interface cells. The temperature, solute concentration, solid fraction, and interfacial curvature of each cell are determined. These values are typically obtained by numerically solving the governing equations presented in this section using the finite difference method.

Temperature T was calculated using the energy equation when the heat flow to the mold must be considered for a larger system. In this study, to use a relatively small system, the temperature was specified unidirectionally along the z-axis under a constrained temperature gradient G and constant cooling rate R_c, as expressed in Eq. (1).

  

$$T(z)=T_0+Gz-R_{c}t,$$(1)

where T₀ is the initial temperature at z = 0, and t is the time.

The solute concentration, C, was calculated using Eq. (2) for the liquid and solid cells, and Eq. (3) for the interfacial cell. The second term on the right side of Eq. (3) represents the solute partition. The solute concentrations in the solid and liquid cells are given by C only; no variables for the solute concentrations that differ for each phase are used.

  

$$\frac{\partial C}{\partial t}=\nabla \cdot (D_{\text{i}}\nabla C)$$(2)

  

$$\frac{\partial C}{\partial t}=\nabla \cdot (D_{\text{L}}\nabla C)+C_L^{*}(1-k)\frac{\partial f_{\text{S}}}{\partial t},$$(3)

where i = L and S represent the liquid and solid phases, respectively; D_L and D_S are the solute diffusion coefficients in the liquid and solid phases, respectively; and k is the equilibrium partition coefficient.

In the two-domain approach, the interface cell is regarded as liquid, and the solute concentration is given by the equilibrium solute concentration of the liquid phase, $C_L^{*}$, expressed in Eq. (4). Therefore, the solute diffusion in the liquid cell adjacent to the interface cell was calculated using Eq. (2). In the solid cell adjacent to the interface cell, the equilibrium solute concentration of the solid phase, $C_{\mathrm{S}}^{*}$ = $k{C}_L^{*}$, was calculated based on the local equilibrium assumption. Assuming that the solute concentration of the solid phase in the interface cell is $C_{\mathrm{S}}^{*}$, solute diffusion can be calculated using Eq. (2).

  

$$C_{\text{L}}^{*}=\frac{T-T_{\text{M}}}{m_{\text{L}}}+\frac{\mathrm{\Gamma }\cdot {\kappa }_{\mathrm{w}\mathrm{m}\mathrm{c}}}{m_{\text{L}}},$$(4)

where T_M is the melting point of the solvent element, m_L is the liquidus slope, Γ is the Gibbs–Thomson coefficient, and κ_wmc is the weighted mean curvature,^19,22) which includes the interfacial energy anisotropy for cubic alloys, and is expressed as shown in Eq. (5).

  

$$\begin{array}{l}{\kappa }_{wmc}=(3{\epsilon }_4-1)\cdot ({\partial }_x{n}_x+{\partial }_y{n}_y+{\partial }_z{n}_z)\\ \hspace{1em}\hspace{1em}-48{\epsilon }_4\cdot (n_x^2\cdot {\partial }_x{n}_x+n_y^2\cdot {\partial }_y{n}_y+n_z^2\cdot {\partial }_z{n}_z)\\ \hspace{1em}\hspace{1em}+12{\epsilon }_4\cdot Q\cdot ({\partial }_x{n}_x+{\partial }_y{n}_y+{\partial }_z{n}_z)\\ \hspace{1em}\hspace{1em}+12{\epsilon }_4\cdot (n_x\cdot {\partial }_{x}Q+n_y\cdot {\partial }_{y}Q+n_z\cdot {\partial }_{z}Q)\end{array}$$(5)

  

$$n_x=\frac{{\partial }_x{f}_{\text{S}}}{\left|\nabla f_S\right|},\mathrm{\hspace{0.17em} } n_y=\frac{{\partial }_y{f}_{\text{S}}}{\left|\nabla f_{\text{S}}\right|},\mathrm{\hspace{0.17em} } n_z=\frac{{\partial }_z{f}_{\text{S}}}{\left|\nabla f_{\text{S}}\right|},\mathrm{\hspace{0.17em} } Q=n_x^4+n_y^4+n_z^4$$

Here, ε₄ is a parameter that provides the anisotropic strength; and n_x, n_y, and n_z are the x-, y-, and z-components of the normal vector of the interface, respectively.

The change in the solid fraction of the interface cell at the n-th time step Δf_S⁽ⁿ⁾ is expressed by Eq. (6).

  

$$\mathrm{\Delta }{f}_{\text{S}}^{(n)}=\frac{C_{\text{L}}^{*(n)}-C_{\text{L}}^{\text{diff}}}{C_{\text{L}}^{*(n)}(1-k^{(n)})}$$(6)

where $C_{\mathrm{L}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ is the solute concentration in the interface cell obtained by calculating solute diffusion. The solid fraction of the interface cell was calculated by integrating Δf_S at each time step. When the solid fraction of the interface cell becomes 1 at the n^th time step, the interface cell is converted to a solid cell, the liquid cells in the Moore neighborhood are converted to interface cells, and the solid–liquid interface shifts toward the liquid phase (solidification). However, when the solid fraction of the interface cell becomes zero at the n^th time step, the interface cell is converted to a liquid cell, the solid cells in the Moore neighborhood are converted to interface cells, and the solid–liquid interface shifts toward the solid phase (melting). To calculate the solid fraction, the limiting condition in Eq. (7) is imposed to prevent the solid fraction from exceeding 1 owing to solidification, whereas the limiting condition in Eq. (8) is imposed to prevent the solid fraction from becoming smaller than zero owing to melting.

  

$$\mathrm{\Delta }{f}_{\text{S}}^{(N)}=min\left(\mathrm{\Delta }{f}_{\text{S}}^{(n)},\mathrm{\hspace{0.17em} } 1-f_{\text{S}}^{(N-1)}\right)$$(7)

  

$$\mathrm{\Delta }{f}_{\text{S}}^{(N)}=-min\left(\left|\mathrm{\Delta }{f}_{\text{S}}^{(n)}\right|,\mathrm{\hspace{0.17em} } f_{\text{S}}^{(N-1)}\right),$$(8)

where N denotes the time step at which the limit of the solid fraction during solidification and melting was exceeded.

3. Mass Balance Error Caused by Model Based on Two-domain Approach

3.1. Estimation of Mass Balance Error

In the CA model based on the two-domain approach, mass-balance errors arise first in the solute diffusion between the interface cell and adjacent cells (liquid and solid cells), second in the solute partition in the interface cell, and third in the calculation of the correction for the change in solid fraction (Eqs. (7) and (8)).

First, we explain the mass-balance error that occurs when solute diffusion is calculated using the finite difference method. When the solute concentration in the interface cell is higher than that in the adjacent cell, the solute diffuses from the interface cell to the adjacent cell. In the opposite case, the solute diffuses from the cell adjacent to the interface cell. Here, we assume that the liquid and solid cells are adjacent to both sides of the interface cell in one dimension (1D). Thus, the relationship between the solute concentrations of the interface cell and adjacent cell can be classified into the following four cases, with the solute concentrations of the solid, liquid, and interface cells denoted as C(S), C(L), and C(I), respectively:

Case-1: C(I) > C(L) and C(I) > C(S)

Case-2: C(I) > C(L) and C(I) < C(S)

Case-3: C(I) < C(L) and C(I) > C(S)

Case-4: C(I) < C(L) and C(I) > C(S)

In Case-1, the solute diffuses from the interface cell to the solid and liquid cells; in Case-4, it diffuses from the solid and liquid cells to the interface cell; and in Case-2 (Case-3), it diffuses from the interface cell to the liquid cell (solid cell), and from the solid cell (liquid cell) to the interface cell. Figure 1 shows a schematic illustration of the solute diffusion around the interface cell in Case-1 and Case-4. For simplicity, we assume a 1D solute concentration distribution in which the interface cell is partitioned into liquid and solid regions based on the solid fraction, and the liquid and solid cells are adjacent to each other on the liquid and solid sides of the interface cell, respectively. When the solute concentration in the interface cell at the n−1th time step, that is, C⁽ⁿ⁻¹⁾ = $C_L^{*(n-1)}$ (solid black line), diffuses into the adjacent liquid cell, the solute concentration becomes $C_{\mathrm{L}\mathrm{S}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ (solid blue line) and then diffuses into the solid cell, and the solute concentration in the interface cell at the n^th time step becomes $C_{\mathrm{L}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ (solid red line). Here, $C_L^{*(n-1)}$ is the equilibrium solute concentration in the liquid phase at the n−1^th time step. $C_{\mathrm{L}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ and $C_{\mathrm{L}\mathrm{S}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ are the transient solute concentrations in the interface cell after solute diffusion. In Case-1, the solute corresponding to the orange area in the liquid and solid cells shows a mass-balance error that increases the average solute concentration within the system, whereas in Case-4, the solute corresponding to the green area shows a mass-balance error that decreases the average solute concentration within the system. In Case-2, solute diffusion into the liquid cell increases the amount of solute, whereas diffusion out of the solid cell decreases the amount of solute. In Case-3, solute diffusion into the solid cell increases the amount of solute, whereas diffusion out of the liquid cell decreases the amount of solute. Therefore, solute diffusion in the interface cell increases and decreases the solute equivalent to the sum of the orange and green areas, respectively. The total mass-balance error during the solute diffusion, $\mathrm{\Delta }{C}_{\mathrm{e}\mathrm{r}\mathrm{r}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$, is expressed by Eq. (9). Notably, $\mathrm{\Delta }{C}_{\mathrm{e}\mathrm{r}\mathrm{r}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ is the same for all four cases.

  

$$\begin{array}{l}\mathrm{\Delta }{C}_{\text{err}}^{\text{diff}}=(C_{\text{L}}^{*(n-1)}-C_{\text{LS}}^{\text{diff}})\cdot f_{\text{S}}^{(n-1)}+(C_{\text{LS}}^{\text{diff}}-C_L^{\text{diff}})\cdot f_{\text{S}}^{(n-1)}\\ \hspace{1em}\hspace{1em}\hspace{1em}=(C_{\text{L}}^{*(n-1)}-C_L^{\text{diff}})\cdot f_{\text{S}}^{(n-1)}\end{array}$$(9)

Fig. 1. Schematic illustration of solute diffusion around the interface cell. (a) Case-1 and (b) Case-4. (Online version in color.)

Next, we explain the mass-balance error that occurs when calculating solute partitioning in the interface cell. Figure 2 shows a schematic illustration for explaining the mass-balance error that occurs when calculating the solute partition, which is performed after the solute-diffusion calculation. In the CA model based on the two-domain approach, solute partitioning is considered in the second term on the right side of Eq. (3), which is equivalent to the change in solid fraction (Eq. (6)) in the simulations. Based on Fig. 2 and Eq. (6), the relationship between the change in the solid fraction and the difference in the solute concentration corresponds to the area of the blue frame being equal to that of the red frame before solute partition. However, during solidification (Δf_S > 0), the solute in the green region shown in Fig. 2(a) decreases, whereas that in the orange region, which corresponds to the difference in the equilibrium solute concentration in the solid phase that occurs during one time-step update, increases. During melting (Δf_S < 0), the amount of solute in the orange region shown in Fig. 2(b) increases. Therefore, the total mass-balance error in solute partition $\mathrm{\Delta }{C}_{\mathrm{e}\mathrm{r}\mathrm{r}}^{\mathrm{p}\mathrm{a}\mathrm{r}}$ is expressed as shown in Eq. (10), for both solidification and melting.

  

$$\mathrm{\Delta }{C}_{\text{err}}^{\text{par}}=(C_L^{\text{diff}}-C_{\text{L}}^{*(n)})\cdot f_{\text{S}}^{(n-1)}+(C_{\text{S}}^{*(n)}-C_{\text{S}}^{*(n-1)})\cdot f_{\text{S}}^{(n-1)}$$(10)

Fig. 2. Schematic illustration showing the mass-balance error that occurs when calculating solute partition, which is performed after the solute diffusion calculation. (a) Solidification, (b) melting. (Online version in color.)

Finally, we discuss the mass-balance errors that occur when converting the cells from interface cells to solid or liquid cells. This error occurs when the value of the solid-fraction change is corrected using the limiting conditions in Eqs. (7) and (8) for solidification and melting, respectively. Figure 3 shows a schematic illustration for explaining the mass-balance error that occurs when the value of the solid-fraction change is corrected using the limiting conditions in Eqs. (7) and (8). The limiting condition for solidification (Eq. (7)) is that when the solid-fraction change calculated using Eq. (6) exceeds 1 − f_S⁽ⁿ⁻¹⁾, the solid-fraction change is corrected to a smaller value than that calculated using Eq. (6). As shown in Fig. 3(a), this corresponds to the solid-fraction change calculated using a solute concentration higher than the solute concentration after calculating solute diffusion. This difference in solute concentration is defined as $\mathrm{\Delta }{C}_{\mathrm{L}}^{\mathrm{s}\mathrm{o}\mathrm{l}}$ (> 0). The corrected solid-fraction change Δf_S^sol is expressed as shown in Eq. (11).

  

$$\mathrm{\Delta }{f}_{\text{S}}^{\text{sol}}=\frac{C_{\text{L}}^{*(n)}-(C_L^{\text{diff}}+\mathrm{\Delta }{C}_{\text{L}}^{\text{sol}})}{C_{\text{L}}^{*(n)}(1-k^{(n)})}$$(11)

Fig. 3. Schematic illustration of the mass-balance error that occurs when the interface cell is converted into solid or liquid cells. (a) Solidification, and (b) melting. (Online version in color.)

Substituting Eq. (6) into Eq. (11) and rearranging, $\mathrm{\Delta }{C}_{\mathrm{L}}^{\mathrm{s}\mathrm{o}\mathrm{l}}$ is expressed as shown in Eq. (12).

  

$$\mathrm{\Delta }{C}_{\text{L}}^{\text{sol}}=C_{\text{L}}^{\text{*(}n)}(1-k^{(n)})(\mathrm{\Delta }{f}_{\text{S}}^{(n)}+f_{\text{S}}^{(n-1)}-1)$$(12)

The limiting condition for melting (Eq. (8)) is that when the solid-fraction change calculated using Eq. (6) becomes smaller than − f_S⁽ⁿ⁻¹⁾, the solid-fraction change is corrected to a larger value than that calculated using Eq. (6). As shown in Fig. 3(b), this corresponds to the solid-fraction change calculated using a solute concentration lower than the solute concentration after calculating solute diffusion. The difference in the solute concentration was defined as $\mathrm{\Delta }{C}_{\mathrm{L}}^{\mathrm{m}\mathrm{e}\mathrm{l}}$ (< 0). The corrected solid-fraction change Δf_S^mel is expressed as shown in Eq. (13).

  

$$\mathrm{\Delta }{f}_{\text{S}}^{\text{mel}}=\frac{C_{\text{L}}^{*(n)}-(C_{\text{L}}^{\text{diff}}+\mathrm{\Delta }{C}_{\text{L}}^{\text{mel}})}{C_{\text{L}}^{*(n)}(1-k^{(n)})}$$(13)

Substituting Eq. (6) into Eq. (13) and rearranging, $\mathrm{\Delta }{C}_{\mathrm{L}}^{\mathrm{m}\mathrm{e}\mathrm{l}}$ is expressed as shown in Eq. (14).

  

$$\mathrm{\Delta }{C}_{\text{L}}^{\text{mel}}=C_{\text{L}}^{*\text{(}n)}(1-k^{(n)})(\mathrm{\Delta }{f}_{\text{S}}^{(n)}+f_{\text{S}}^{(n-1)})$$(14)

Therefore, the mass-balance error that occurs when converting a cell from an interface cell to a solid or liquid cell, $\mathrm{\Delta }{C}_{\mathrm{e}\mathrm{r}\mathrm{r}}^{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}}$, is expressed as shown in Eq. (15) and corresponds to the sky blue area in Fig. 3.

  

$$\mathrm{\Delta }{C}_{\text{err}}^{\text{conv}}=\mathrm{\Delta }{C}_L^{corr}\cdot f_{\text{S}}^{(n-1)}$$(15)

where $\mathrm{\Delta }{C}_{\mathrm{L}}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}$ is the solute-concentration change corrected when converting the interface cell into a solid or liquid cell, that is, $\mathrm{\Delta }{C}_{\mathrm{L}}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}$ = $\mathrm{\Delta }{C}_{\mathrm{L}}^{\mathrm{s}\mathrm{o}\mathrm{l}}$ when Δf_S⁽ⁿ⁾ is corrected to Δf_S^sol, $\mathrm{\Delta }{C}_{\mathrm{L}}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}$ = $\mathrm{\Delta }{C}_{\mathrm{L}}^{\mathrm{m}\mathrm{e}\mathrm{l}}$ when Δf_S⁽ⁿ⁾ is corrected to Δf_S^mel, and $\mathrm{\Delta }{C}_{\mathrm{L}}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}$ = 0 when Δf_S⁽ⁿ⁾ is not corrected.

Consequently, the total mass-balance error ΔC_err in the CA model based on the two-domain approach is expressed as the sum of mass-balance errors due to solute diffusion and partitioning in the interface cell and the cell-state conversion, which is shown in Eq. (16).

  

$$\begin{array}{l}\mathrm{\Delta }{C}_{\text{err}}=\mathrm{\Delta }{C}_{\text{err}}^{diff}+\mathrm{\Delta }{C}_{\text{err}}^{par}+\mathrm{\Delta }{C}_{\text{err}}^{conv}\\ \hspace{1em}\hspace{1em}\hspace{0.5em}=\left\{C_L^{*(n-1)}-C_{\text{S}}^{*(n-1)}-(C_{\text{L}}^{*(n)}-C_{\text{S}}^{*(n)})+\mathrm{\Delta }{C}_{\mathrm{L}}^{corr}\right\}\cdot f_{\text{S}}^{(n-1)}\\ \hspace{1em}\hspace{1em}\hspace{0.5em}=\left\{C_L^{*(n-1)}(1-k^{(n-1)})-C_{\text{L}}^{*(n)}(1-k^{(n)})+\mathrm{\Delta }{C}_{\mathrm{L}}^{corr}\right\}\cdot f_{\text{S}}^{(n-1)}\end{array}$$(16)

3.2. Correction of Mass-balance Error

As mentioned in the previous section, the CA model based on the two-domain approach generates a mass-balance error calculated using Eq. (16), at each time step. Mass balance is satisfied by appropriately correcting this error at each time step. Notably, the mass-balance error estimated using Eq. (16) is described only by the equilibrium concentration and the solid fraction of the interface cell. In other words, because the mass-balance error can be calculated using only the values at the interface cell, it can be easily calculated even when extended to 2D or 3D systems. In this study, the ΔC_err calculated for each interface cell is corrected by equally distributing it to the liquid cells inside a circle with a radius of 5Δx from the cell center. If no liquid cell exists within a circle, ΔC_err is corrected by adding it to the solute concentration of the liquid phase in the interface cell. This correction procedure is similar to simulating solute-concentration fluctuations using artificial random noise, which is typically used for dendritic growth.^31,32,33)

4. Validation of Effect of Mass Balance Error

In this study, we numerically simulated the single dendritic growth of an Al-4.5 mass%Cu alloy during isothermal undercooling and continuous cooling with and without mass-balance error correction and investigated the effect of the mass-balance error on the microstructural evolution. Owing to mesh anisotropy, it is difficult for the CA model to accurately simulate dendrite morphologies with different tip growth directions. This is because of the accuracy of the numerical calculation of interface curvature. However, it should be noted that mesh anisotropy does not affect the occurrence of mass-balance errors and their corrections. Therefore, in this study, we discuss the correction of mass-balance errors in simulations where the dendrite tip growth direction, which corresponds to the (100) direction, coincides with the coordinate axis.

4.1. Calculation Conditions

The calculation conditions and physical properties of the Al–Cu alloy³⁴⁾ used in the calculations are listed in Table 1. The simulation was performed using a 3D system, and the calculation domain was partitioned into n_x × n_y × n_z cells with a cell size of Δx = 1.0 μm, where n_x, n_y, and n_z represent the numbers of cells in the x-, y-, and z-directions, respectively. Zero Neumann boundary conditions were applied in all the directions. Dendritic growth began in the initial solid phase, specified at the center of the xy-plane at z = 0. Figure 4 shows the liquidus and solidus lines for the range of temperatures and Cu compositions used in this study. The equilibrium phase diagram of the Al–Cu alloy was obtained using the thermodynamic calculation software Pandat.³⁵⁾ The equilibrium partition coefficient k and equilibrium solute concentration in the liquid phase $C_{\mathrm{L}}^{\mathrm{e}\mathrm{q}}$ were obtained as a function of temperature T by fitting the curved liquidus and solidus lines via polynomial approximations.

  

$$\begin{array}{l}k=-4.55422\times {10}^{-9}{T}^3+1.53427\times {10}^{-5}{T}^2\\ \hspace{1em}-1.67869\times {10}^{-3}T+6.13362\end{array}$$(17)

  

$$\begin{array}{l}{C}_{\mathrm{L}}^{\mathrm{e}\mathrm{q}}=-2.14691\times {10}^{-9}{T}^3-1.01624\times {10}^{-3}{T}^2\\ \hspace{1em}\hspace{1em}+1.49142T+504.993,\end{array}$$(18)

where temperature is expressed in units of Kelvin. The equilibrium solute concentration in the liquid phase of the interface cell at T, $C_L^{*}$(T), was obtained by substituting the temperature T_GT while considering the Gibbs–Thomson relationship into $C_{\mathrm{L}}^{\mathrm{e}\mathrm{q}}$(T), that is, $C_L^{*}$(T) = $C_{\mathrm{L}}^{\mathrm{e}\mathrm{q}}$(T_GT), where T_GT = T + Γκ_wmc. Because the purpose of this study is to discuss the effect of mass-balance errors on microstructural evolution, we do not discuss the curvature-calculation error, which is another important issue in the CA model. However, the interface curvature based on the weighted mean curvature was calculated using a finite-difference scheme based on bilinear interpolation, which has been reported to reduce the curvature-calculation error.¹³⁾ Furthermore, the correction of the mass-balance error is not affected by the curvature calculation.

Table 1. Calculation parameters and material properties of Al–Cu alloy used in the simulation.

	Symbol	Value	Units
Diffusion coefficient in liquid phase	DL	3.0 × 10−9	m2 s−1
Diffusion coefficient in solid phase	DS	3.0 × 10−13	m2 s−1
Gibbs–Thomson coefficient	Γ	1.6 × 10−7	K m
Anisotropic strength	ε4	0.02	–
CA cell number in x-direction	nx	128	–
CA cell number in y-direction	ny	256	–
CA cell number in z-direction	nz	256	–
CA cell size	Δx	1.0	μm

Fig. 4. Liquidus and solidus curves of primary phase of Al–Cu alloy calculated using thermodynamic calculation software. (Online version in color.)

In the simulation of isothermal undercooling, the effect of mass-balance error on microstructural evolution was investigated under constant values of k and $C_{\mathrm{L}}^{\mathrm{e}\mathrm{q}}$. The computational domain was set to four constant undercooling degrees, that is, ΔT = 1.0, 1.5, 2.0, and 2.5 K, and the calculations were performed up to 2000000, 1500000, 1000000, and 500000 timesteps, respectively.

In the simulation of continuous cooling, the effect of the mass-balance error on the microstructural evolution was investigated under k and $C_{\mathrm{L}}^{\mathrm{e}\mathrm{q}}$ varying with temperature. Simulations were performed under four calculation conditions: two cooling rates R_c = 0.5 Ks⁻¹ and 1.0 Ks⁻¹ for G = 5.0 Kmm⁻¹, and two cooling rates R_c = 1.0 Ks⁻¹ and 2.0 Ks⁻¹ for G = 10.0 Kmm⁻¹. T₀ was set to 922.1 K, and the calculation was performed until the temperature of the CA cell at the top of the calculation domain (n_z = 256) reached the eutectic temperature of the Al–Cu alloy (821.4 K).

5. Results and Discussion

5.1. Isothermal Solidification under Constant Undercooling

Figure 5 shows snapshots of the dendritic growth morphology obtained with and without correction for each undercooling. At ΔT = 1.0 and 1.5 K, only coarse dendrite primary arms were observed with almost no development of secondary arms, whereas at ΔT = 2.0 and 2.5 K, a dendritic morphology was observed with the development of secondary arms. Although no clear difference in morphology was observed between the simulation results with and without the correction of the mass-balance error, the morphology of the secondary arms at ΔT = 2.0 and 2.5 K differed slightly. Because the necessity of correcting mass-balance errors is ambiguous based on merely comparing the dendrite morphology, we calculated the percentage error in the solute concentration, PE(%), using Eq. (19) and verified the effect of the solute concentration on the growth process.

  

$$\text{PE}(\%)=\frac{C_{ave}-C_0}{C_0}\times 100$$(19)

where C₀ is the initial solute concentration (4.5 mass%Cu) and C_ave is the average solute concentration over the entire computational domain.

Fig. 5. The dendritic growth morphologies obtained with and without correction were simulated under isothermal undercooling conditions. (a) ΔT = 1.0 K, (b) ΔT = 1.5 K, (c) ΔT = 2.0 K, and (d) ΔT = 2.5 K. (Online version in color.)

Figure 6 shows the relationship between PE(%) and elapsed time with and without correction for the mass-balance error. In the simulation without correction, the PE(%) increased with time, and the value of PE(%) was larger in the simulation with higher undercooling. The increase in PE(%) can be explained using Eq. (16). During isothermal undercooling, the equilibrium solute concentration in the liquid phase and the equilibrium partition coefficient did not change over time. Here, we assumed that the interface curvature changed minimally during the time step Δt. Therefore, the mass-balance error for each time step is $\mathrm{\Delta }{C}_{\mathrm{L}}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}$·f_S⁽ⁿ⁻¹⁾, and because $C_{\mathrm{L}}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}$ = $C_{\mathrm{L}}^{\mathrm{s}\mathrm{o}\mathrm{l}}$ during solidification, the mass-balance error increases. The arms may melt during isothermal undercooling because of coarsening and coalescence. However, the solid-liquid interface moves in a direction that reduces the interfacial area. Therefore, the interface cells are converted into solid cells, $C_{\mathrm{L}}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}$ = $C_{\mathrm{L}}^{\mathrm{s}\mathrm{o}\mathrm{l}}$, and the mass-balance error increases during isothermal undercooling. Specifically, the PE(%) values for 500000 time steps at ΔT = 1.0, 1.5, 2.0, and 2.5 K were 0.0519%, 0.142%, 0.166%, and 1.31%, respectively. For ΔT = 1.0, 1.5, and 2.0 K, the mass-balance error continued to increase even after 500000 time steps, and PE(%) was 0.127% at 2000000 time steps for ΔT = 1.0 K, 0.327% at 1500000 time steps for ΔT = 1.5 K, and 0.297% at 1000000 time steps for ΔT = 2.0 K. Compared with these results, the mass-balance error at ΔT = 2.5 K is extremely large. Because PE(%) is considered to be related to the dendritic morphology, we investigated the interfacial area density, S_v, for each condition. S_v is defined as the total interfacial area divided by the total volume of the system,^36,37,38) where a larger S_v indicates a more complex solid–liquid interface morphology. The maximum values of S_v at ΔT = 1.0, 1.5, 2.0, and 2.5 K were 0.010, 0.015, 0.025, and 0.047 μm⁻¹, respectively, and the S_v at ΔT = 2.5 K was significantly higher than those under other conditions. This corresponded to the highest number of interface cells in the simulation at ΔT = 2.5 K. Therefore, the mass-balance error that occurs when an interface cell is converted into a solid cell, $\mathrm{\Delta }{C}_{\mathrm{L}}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}$, is the main factor that increases PE(%). Additionally, the mass-balance error was the largest in the simulation with ΔT = 2.5 K, which indicated the highest S_v value. In contrast, in the simulation with correction, PE(%) did not increase; moreover, it became almost constant at 100000–200000 time steps. Thus, one can confirm that the PE(%) value for all conditions remained extremely low and that almost no mass-balance error occurred. This proves that the mass-balance error in the conversion of the interface cell was appropriately corrected. Although the value is slightly negative, it remains within the range of the numerical calculation error.

Fig. 6. Relationship between percentage error in solute concentration, PE (%), and elapsed time for mass-balance error in simulations of isothermal undercooling. (a) Without and (b) with correction. (Online version in color.)

Figure 7 shows the relationship between the volume fraction of solids, g_S, and the elapsed time under each undercooling. In all simulations, the volume fraction of solids increased rapidly in the first 100000–200000 time steps. Subsequently, in the simulations with correction, the volume fraction of solids gradually increased until it converged to a certain value. The convergence value is the solid fraction, f_S^eq = ($C_{\mathrm{L}}^{\mathrm{e}\mathrm{q}}$−C₀)/$C_{\mathrm{L}}^{\mathrm{e}\mathrm{q}}$ (1−k), calculated based on the lever rule for each undercooling, that is, f_S^eq = 0.091, 0.131, 0.168, and 0.203 at ΔT = 1.0, 1.5, 2.0, and 2.5 K, respectively. f_S^eq was not reached within the total time-step range set in this study. This is because of the lower growth rate at lower undercooling. We believe that the solid fraction will reach f_S^eq if we continue the calculations. However, in the simulation without correction, the solid volume fraction increased rapidly and then decreased, which was more pronounced under higher undercooling conditions. The decrease in the solid volume fraction was due to an increase in the solute concentration caused by the mass-balance error, which induced melting. As shown in Fig. 6, in the simulation without correction, the mass-balance error continued to increase with time. Correspondingly, the solid fraction decreased. In particular, in the simulation at ΔT = 2.5 K, both the decrease in the volume fraction of solids and the mass-balance error were the largest, and one can conclude that the decrease in the volume fraction of solids depends on the amount of mass-balance error. During isothermal undercooling, the solid fraction did not decrease thermodynamically. However, because melting of the solid–liquid interface was considered in this model, we were able to confirm the inappropriate phenomena caused by the mass-balance error. Therefore, the mass-balance error must be corrected to appropriately simulate microstructural evolution under isothermal undercooling.

Fig. 7. Relationship between volume fraction of solid and elapsed time in simulations of isothermal undercooling. (a) ΔT = 1.0 K, (b) ΔT = 1.5 K, (c) ΔT = 2.0 K, and (d) ΔT = 2.5 K. (Online version in color.)

5.2. Solidification under Continuous Cooling

Figure 8 shows snapshots of the dendritic growth morphology obtained with and without correction for each temperature gradient and cooling rate. Under the same temperature gradient, the higher the cooling rate, the higher is the drawing rate, V = R_c/G. Meanwhile, under the same cooling rate, the smaller the temperature gradient, the higher is the drawing rate. Therefore, the higher the drawing rate, the higher was the growth rate, thus resulting in a complex dendrite morphology with developed side–arms. During continuous cooling, similar to during isothermal undercooling, we could not readily clarify the effect of mass-balance errors on the dendrite morphology. Therefore, we performed a quantitative evaluation using PE(%) and S_v.

Fig. 8. The dendritic growth morphologies obtained with and without correction simulated under continuous cooling conditions. (a) G = 5.0 Kmm⁻¹, R_c = 0.5 Ks⁻¹, (b) G = 5.0 Kmm⁻¹, R_c = 1.0 Ks⁻¹ (c) G = 10.0 Kmm⁻¹, R_c = 1.0 Ks⁻¹, and (d) G = 10.0 Kmm⁻¹, R_c = 2.0 Ks⁻¹. (Online version in color.)

Figure 9 shows the relationship between PE(%) and the elapsed time obtained in simulations with and without correction under each condition. The horizontal axis represents the dimensionless time, t/t_end, where t_end is the time at which the eutectic temperature is reached. In the simulation without correction, the higher the V and R_c, the larger was the negative value of the mass-balance error. The increase in PE(%) with the negative values can be explained using Eq. (16), which is similar to the simulation results for isothermal undercooling. In the continuous cooling of an alloy with k < 1, the equilibrium solute concentration in the liquid phase increases as the temperature decreases. As the equilibrium partition coefficient does not change rapidly in one time step, the second term in the parentheses on the right-hand side of Eq. (16) is larger than that of the first term, and the mass-balance error value becomes negative. The effect of the third term is less significant than those of the first and second terms in the continuous-cooling process because this mass-balance error occurs only when f_S > 1 or f_S < 0 in the interface cell. Specifically, the mass-balance error at G = 10.0 Kmm⁻¹, R_c = 2.0 Ks⁻¹ was the largest among the four conditions, with a PE(%) of −14.9% at t/t_end = 1, whereas the mass-balance error at G = 5.0 Kmm⁻¹, R_c = 0.5 Ks⁻¹ was the smallest, with a PE(%) of −3.01% at t/t_end = 1. The absolute values of these PE(%) were significantly higher than those obtained during isothermal undercooling. This is because the main cause of the mass-balance error during isothermal undercooling is the change in the cell state (the third term of Eq. (16)), whereas it is the solute diffusion and distribution during continuous cooling (the first and second terms of Eq. (16)). Calculations of solute diffusion and partitioning for the interface cell were performed at each time step. Therefore, during continuous cooling, a large mass-balance error occurred unless it was corrected. However, in the simulation with correction, the absolute value of PE(%) was less than 0.05%, which is extremely low compared with the simulation result without correction. Additionally, the behavior of PE(%) did not show a clear dependence on the solidification conditions, such as the temperature gradient and cooling rate.

Fig. 9. Relationship between percentage error in solute concentration, PE (%), and dimensionless time, t/t_end, for mass-balance error in simulations of the continuous cooling process. (a) Without and (b) with correction. (Online version in color.)

Figure 10 shows the relationship between the interfacial area density, S_v, and volume fraction of solids, g_S, obtained from the simulation under each condition. Regardless of whether correction was applied, under all conditions, S_v increased in the early stages owing to the development of dendrite arms and then decreased owing to the coarsening of the arms. S_v reached its maximum when g_S was in the range of 0.5 to 0.6, and in all conditions, the maximum value of S_v in the simulation without correction was slightly larger. However, this difference in S_v was marginal and could not be clearly confirmed from the dendritic morphology shown in Fig. 8.

Fig. 10. Relationship between interfacial area density, S_v, and volume fraction of solids, g_S, based on simulations of continuous cooling. (a) G = 5.0 Kmm⁻¹ and (b) G = 10.0 Kmm⁻¹. (Online version in color.)

Figure 11 shows the relationship between PE(%) and g_S for the simulations without correction under each condition. In the simulation without correction, the mass-balance error increased significantly to a negative value as g_S increased, beginning when g_S reached 0.3–0.4. This indicates that the mass-balance error occurred more significantly during the coarsening of the dendrite arms than during their development. During continuous cooling, a dendrite skeleton develops in the early stages of solidification, and then solidification gradually progresses because of the coarsening of the arms. In the simulation, when t/t_end was approximately 0.07, g_S was approximately 0.4, indicating that most of the time steps were in the coarsening process. Therefore, mass-balance errors occurred primarily during the calculation of solute diffusion and partitioning when the growth rate of the solid–liquid interface was low. This consistently explains the absence of significant difference in the dendrite morphology with and without correction. Thus, when predicting only the dendrite morphology in the early stages of solidification, the correction for the mass-balance error may not necessarily affect the results significantly. However, when predicting the microsegregation that occurs until the end of solidification, correction is necessary.

Fig. 11. Relationship between percent error in solute concentration, PE (%), and volume fraction of solids, g_S, based on continuous cooling simulations. (Online version in color.)

Figure 12 shows the solute concentration distributions in the liquid phase on the yz-cross section at g_S = 0.5 for the simulations with and without correction under G = 10.0 Kmm⁻¹ and R_c = 2.0 Ks⁻¹. Because of unidirectional solidification under a temperature gradient, the solute concentration in the liquid phase decreased toward the +z direction in both simulations. The color distribution showed that the simulation with correction had a higher solute concentration in the liquid phase than that without correction. The PE(%) at g_S = 0.5 in the simulation without correction was −1.75%, and the average solute concentration in the system decreased. The maximum solute concentrations in the liquid phase at g_S = 0.5 in the simulations with and without correction were 8.820 and 8.674 mass%Cu, respectively, and the average solute concentrations in the liquid phase were 8.158 and 8.007 mass%Cu, respectively. In contrast, the maximum solute concentrations in the solid phase at g_S = 0.5 in the simulations with and without correction were 1.151 and 1.128 mass%Cu, respectively, and the average solute concentrations in the solid phase were 0.830 and 0.824 mass%Cu, respectively. The maximum and average values of the solute concentration in the solid phase did not differ significantly with or without correction. In this model, the solute concentration in the solid phase did not artificially increase or decrease, because the mass-balance error was corrected by the solute concentration in the liquid phase. Therefore, the mass-balance error due to solute diffusion occurred mainly between the interface and liquid cells, possibly because the diffusion coefficient of the solid phase was several orders of magnitude smaller than that of the liquid phase.

Fig. 12. Solute concentration distribution in the liquid phase at yz-cross section of g_S = 0.5 during continuous cooling under G = 10.0 Kmm⁻¹, R_c = 2.0 Ks⁻¹. (a) Without and (b) with correction. (Online version in color.)

Finally, we investigated the microsegregation behavior in simulations with and without correction for the mass-balance error, that is, the change in the average solute concentration in the liquid phase with the solidification process. Figure 13 shows the relationship between the solid volume fraction and the average solute concentration in the liquid phase, C_L/C₀, obtained by simulations with and without correction under the conditions of G = 10.0 Kmm⁻¹ and R_c = 2.0 Ks⁻¹. For comparison, the relationship between the solid fraction and solute concentration in the liquid phase based on equilibrium (lever rule) and non-equilibrium solidification (Scheil) is shown. In the Al–Cu alloy, the relationship between the solid fraction and the average solute concentration in the liquid phase, with or without correction, is closer to the Scheil relationship than to the lever rule because the diffusion coefficient of solute Cu in the solid phase is small. Comparing the simulations with and without correction, the solute concentration in the simulation with correction was higher and closer to the Scheil relationship. This is because the mass-balance error was negative without correction, as shown in Fig. 11. Therefore, in the quantitative CA simulation using the two-domain approach, the microsegregation behavior is significantly affected if the mass-balance error is not corrected. In particular, in multicomponent alloys, the solute concentration in the liquid phase is directly related to the solidification path; therefore, it is expected to have a significant effect on the quantitative predictions of the composition and size of inclusions and precipitates. Consequently, we confirmed that the CA model developed in this study can satisfy solute conservation and quantitatively predict microsegregation behavior under various solidification conditions and alloy systems.

Fig. 13. Relationship between the solid volume fraction and average liquid solute concentration of simulations with and without correction (G = 10.0 Kmm⁻¹and R_c = 2.0 Ks⁻¹). (Online version in color.)

In addition, to confirm the effect of cell size on the correction of mass-balance error, we performed simulations using different cell sizes (1.5 and 2.0 μm) and obtained similar conclusions with the simulations using a cell size of 1.0 μm. In the simulations with correction, the PE(%) was almost zero regardless of the cell size, and the smaller the cell size, the smaller the error value. The absolute values of PE(%) under isothermal undercooling and continuous cooling were <0.005% and 0.2%, respectively. In the simulations without correction, PE(%) increased to the positive side under isothermal undercooling conditions and decreased to the negative side under continuous cooling conditions; the trend did not change with cell size.

6. Conclusions

In this study, we developed a CA model based on a two-domain approach that can reduce mass-balance errors. In this model, the solid–liquid interface movement associated with solidification and melting was considered to accurately simulate the solid–liquid interface movement that satisfies the Gibbs–Thomson relationship. Using this model, we performed simulations of dendrite growth during isothermal undercooling and continuous cooling, and investigated the effect of mass-balance error on microstructural evolution.

During isothermal undercooling, no correction to the mass-balance error resulted in an increase in the solute concentration when the dendrite arms coarsened, with a greater increase observed under more significant undercooling. Consequently, the solid–liquid interface melted and the solid volume fraction decreased. However, during continuous cooling, no correction to the mass-balance error resulted in a decrease in the solute concentration when the dendrite arms coarsened. The mass-balance error from continuous cooling solidification was larger than that from isothermal undercooling solidification. This is because, during isothermal undercooling, the mass-balance error that occurs when calculating the conversion of an interface cell into a solid-phase cell is dominant, whereas during continuous cooling, the mass-balance error that occurs when calculating solute diffusion and partitioning becomes dominant.

The CA model developed in this study reduced the mass-balance error in the simulations under various solidification conditions to almost zero. In particular, during continuous cooling, the mass-balance error was reduced to almost zero even when a curved liquidus corresponding to the actual equilibrium phase diagram and equilibrium partition coefficient as a function of temperature were used. Under all calculation conditions, the effect of the mass-balance error on the dendrite morphology in the early stage of solidification was insignificant. However, this mass-balance error significantly affected the quantitative prediction of microsegregation that occurred in the final solidification stage. Therefore, the mass-balance error must be corrected for a CA model based on the two-domain approach, and we believe that if a highly accurate method for calculating the interface curvature is implemented in the proposed CA model, it will be widely used as a practical quantitative model in the future.

Statement for Conflict of Interest

The authors declare no conflict of interest regarding this manuscript.

Acknowledgments

This study was supported by JSPS KAKENHI (grant numbers JP 22H01839 and 23K23107) and the 33rd ISIJ Research Promotion Grant.

References

• 1)   J. A.  Dantzig and  M.  Rappaz: Solidification, 2nd ed., EPFL Press, (2016), 539. ISBN978-2-940222-97-1
• 2)   R.  Trivedi and  W.  Kurz: Int. Mater. Rev., 39 (1994), 49. https://doi.org/10.1179/imr.1994.39.2.49
• 3)   R.  Kobayashi: Physica D, 63 (1993), 410. https://doi.org/10.1016/0167-2789(93)90120-P
• 4)   R.  Kobayashi: Exp. Math., 3 (1994), 59. https://doi.org/10.1080/10586458.1994.10504577
• 5)   A.  Karma: Phys. Rev. Let., 87 (2001), 115701. https://doi.org/10.1103/PhysRevLett.87.115701
• 6)   B.  Echebarria,  R.  Folch,  A.  Karma and  M.  Plapp: Phys. Rev. E, 70 (2004), 061604. https://doi.org/10.1103/PhysRevE.70.061604
• 7)   M.  Ohno and  K.  Matsuura: Phys. Rev. E, 79 (2009), 031603. https://doi.org/10.1103/PhysRevE.79.031603
• 8)   M.  Ohno: Phys. Rev. E, 86 (2012), 051603. https://doi.org/10.1103/PhysRevE.86.051603
• 9)   K.  Reuther and  M.  Rettenmayr: Comput. Mater. Sci., 95 (2014), 213. https://doi.org/10.1016/j.commatsci.2014.07.037
• 10)   L.  Beltran-Sanchez and  D. M.  Stefanescu: Metall. Mater. Trans. A, 34 (2003), 367. https://doi.org/10.1007/s11661-003-0338-z
• 11)   L.  Beltran-Sanchez and  D. M.  Stefanescu: Metall. Mater. Trans. A, 35 (2004), 2471. https://doi.org/10.1007/s11661-006-0227-3
• 12)   W.  Wang,  P. D.  Lee and  M.  McLean: Acta Mater., 51 (2003), 2971. https://doi.org/10.1016/S1359-6454(03)00110-1
• 13)   L.  Wei,  X.  Lin,  M.  Wang and  W.  Huang: Mater. Disco., 3 (2016), 17. https://doi.org/10.1016/j.md.2016.06.001
• 14)   L.  Wei,  Y. Q.  Cao,  X.  Lin,  M.  Wang and  W.  Huang: Comput. Mater. Sci., 156 (2019), 157. https://doi.org/10.1016/j.commatsci.2018.09.019
• 15)   K.  Reuther and  M.  Rettenmayr: Comput. Mater. Sci., 166 (2019), 143. https://doi.org/10.1016/j.commatsci.2019.04.050
• 16)   M. F.  Zhu and  D. M.  Stefanescu: Acta Mater., 55 (2007), 1741. https://doi.org/10.1016/j.actamat.2006.10.037
• 17)   M. F.  Zhu and  C. P.  Hong: ISIJ Int., 41 (2001), 436. https://doi.org/10.2355/isijinternational.41.436
• 18)   M. F.  Zhu and  C. P.  Hong: ISIJ Int., 42 (2002), 520. https://doi.org/10.2355/isijinternational.42.520
• 19)   S.  Pan and  M.  Zhu: Acta Mater., 58 (2010), 340. https://doi.org/10.1016/j.actamat.2009.09.012
• 20)   Y.  Natsume and  K.  Ohsasa: Mater. Sci. Eng., 84 (2015), 012050. https://doi.org/10.1088/1757-899X/84/1/012050
• 21)   Y.  Natsume: J JFS, 87 (2015), 636 (in Japanese). https://doi.org/10.11279/jfes.87.636
• 22)   Y.  Natsume: ISIJ Int., 56 (2016), 1378. https://doi.org/10.2355/isijinternational.ISIJINT-2016-058
• 23)   J.  Ogawa and  Y.  Natsume: Mater. Sci. Eng., 861 (2020), 012059. https://doi.org/10.1088/1757-899X/861/1/012059
• 24)   J.  Ogawa and  Y.  Natsume: J JFS, 94 (2022), 187 (in Japanese). https://doi.org/10.11279/jfes.94.187
• 25)   S. C.  Michelic,  J. M.  Thuswaldner and  C.  Bernhard: Acta Mater., 58 (2010), 2738. https://doi.org/10.1016/j.actamat.2010.01.010
• 26)   A.  Zyska: Materials, 14 (2021), 3393. https://doi.org/10.3390/ma14123393
• 27)   Q.  Zhang,  H.  Fang,  H.  Xue,  S.  Pan,  M.  Rettenmayr and  M.  Zhu: Scientific Reports, 7 (2017), 17809. https://doi.org/10.1038/s41598-017-17857-2
• 28)   Q.  Zhang,  H.  Xue,  Q.  Tang,  S.  Pan,  M.  Rettenmayr and  M.  Zhu: Comput. Mater. Sci., 146 (2018), 204. https://doi.org/10.1016/j.commatsci.2018.01.032
• 29)   Q.  Zhang,  X.  Shen,  Z.  Wang,  S.  Zhang and  M.  Zhu: J. Mater. Sci., 56 (2021), 10393. https://doi.org/10.1007/s10853-021-05958-3
• 30)   X.  Liang,  C.  Bos,  M.  Hermans and  I.  Richardson: Metall. Mater. Trans. B, 54 (2023), 1088. https://doi.org/10.1007/s11663-023-02742-3
• 31)   T.  Takaki: ISIJ Int., 54 (2014), 437. https://doi.org/10.2355/isijinternational.54.437
• 32)   M.  Ohno: ISIJ Int., 60 (2020), 2745. https://doi.org/10.2355/isijinternational.ISIJINT-2020-174
• 33)   D.  Tourret,  H.  Liu and  J. L Lorca: Prog. Mater. Sci., 123 (2022), 100810. https://doi.org/10.1016/j.pmatsci.2021.100810
• 34)   W.  Kurz and  D. J.  Fisher: Fundamentals of Solidification, 4th ed., Trans Tech Publications Ltd, Switzerland, (1998), 294. ISBN0-87849-804-4
• 35)   S.-L.  Chen,  S.  Daniel,  F.  Zhang,  Y. A.  Chang,  X.-Y.  Yan,  F.-Y.  Xie,  R.  Schmid-Fetzer and  W. A.  Oates: Calphad, 26 (2002), 175. https://doi.org/10.1016/S0364-5916(02)00034-2
• 36)   Q.  Du,  D. G.  Eskin,  A.  Jacot and  L.  Katgerman: Acta Mater., 55 (2007), 1523. https://doi.org/10.1016/j.actamat.2006.10.035
• 37)   R.  Mendoza,  J.  Alkemper and  P. W.  Voorhees: Metall. Trans. A, 34 (2003), 481. https://doi.org/10.1007/s11661-003-0084-2
• 38)   H.  Neumann-Heyme,  K.  Eckert and  C.  Beckermann: Acta Mater., 140 (2017), 87. https://doi.org/10.1016/j.actamat.2017.08.021