Dendritic growth is a common phenomenon during casting of metals and has a significant impact on microstructure and mechanical property of finished products. As a widely used numerical approach, phase field model can formulate dendritic growth quantitatively and help us better understand the formation of dendritic morphology. Nowadays, it has no difficulty modeling single and polycrystalline dendritic growth of pure metals and alloys, both in 2D and 3D [[1], [2], [3], [4], [5], [6], [7]]. However, lots of previous studies only involved thermal and solute diffusion without melt flow and motion of solids. This is far from the actual situation. To analyze the influence of melt flow on growth, some simulations coupled with fluid dynamics were performed to model fixed dendrite growing in forced or natural convection [[8], [9], [10], [11], [12], [13]]. Nevertheless, these simulations only apply to dendrites attached to mould during casting but not to those moving freely in melt. Consequently, a more realistic model should not ignore movement of solids. There are three issues in the presence of motion of solid.

Firstly, it requires a suitable FSI model to solve flow field both in liquid and solid. Tong et al. [8] forced the velocity in solid to become zero by adding “distributed momentum sink” into momentum equation so that solid couldn’t move at all. Yamaguchi et al. [14,15] employed material point method to model motion, rotation, and even deformation of both single dendrite and polycrystalline structure [14,15]. However, deformation-induced liquid flow was not really solved, but given by “zero-gradient extension”. Consequently, any point in liquid would have the same velocity as the closest solid [14]. LBM (Lattice Boltzmann model) is a promising numerical approach to solve flow problems. Rojas et al. [16,17] employed it to simulate motion and growth of multiple dendrites. In their study, velocity in liquid and each dendrite was solved individually. Total force and torque acting on every dendrite were computed by regarding them as rigid bodies. An alternative is two-phase flow model. Solid is artificially given a viscosity large enough to resist deformation caused by flow [18,19].

Secondly, locations of dendrites should be updated according to computed flow field. This could be accomplished by embedding advection terms into phase field model [14,15,19] directly or solving advection equations individually [16,17].

Thirdly, an appropriate way is necessary to handle orientation evolution of dendrite and its change due to rotation. Otherwise, dendrite might grow into unphysical distorted shapes (Fig. 2(c) in Ref. [14]). Rojas et al. [17] used multiple phase field equations for multiple dendrites and solved them separately. Obviously, more dendrites mean more equations to be solved. A better strategy is to introduce orientation field into phase field model [14,15,19]. In this way, no matter how many dendrites there are, only one phase field variable is required. Treatment of orientation change due to rotation is straightforward according to the rotating rate of dendrite.

Two-phase flow model is our choice to compute flow field. In this approach, two phases can be treated as a whole rather than separately. Different from Refs. [18] and [19] which also chose this approach, we describe dendritic growth by vector-valued phase field model [[20], [21], [22]], which is computationally convenient to handle orientation periodicity (see Section 2.1.2) and collision between solids. Locations of dendrites are updated by solving advection equations individually so that phase field model still leads to symmetric and positive definite stiff matrices. After rigorous numerical tests, we simulated growth, settlement, rotation, and collision of a single dendrite and multiple dendrites in Al-Cu alloy.

The most challenging aspect in this work is the flow field computation with a steep viscosity variation across s-l interface. In this study, two-phase flow is formulated by N-S equations. Unfortunately, lots of studies were focused on huge viscosity jump for Stokes flow [[23], [24], [25], [26]] rather than N-S flow. To ensure numerical stability, an exponential expression is constructed to make the viscosity vary across s-l interface in a relatively steady way. We develop a preconditioner combined with scaling strategy, taking into account the balance between momentum and continuity equation. Two-level time stepping scheme is designed to reduce run time as far as possible (Section 2.4). Adaptive mesh refinement and parallel computation based on MPI protocol are adopted to further accelerate computation. Numerical tests show our model can withstand a tremendous viscosity ratio (10¹⁰) between solid and liquid. The numerical approaches we develop to compute flow field in this study would be of great interest in other fields like simulation of the earth’s convective mantle [27], ice sheet [28], and polymer processing techniques [23], which also involve huge viscosity variation. Our code is elaborated based on the finite element library, DEAL.II [48].

Symbols in this paper are listed in Appendix D for ease of lookup.

2.1.1. Scalar-valued phase field model

The conventional polycrystalline phase field model takes scalar-valued form and involves following variables: phase field ϕ, temperature T, solute concentration c, and orientation field α [22]. ϕ varies smoothly from -1 in the liquid to 1 in the solid through the diffuse interface. The s-l interface is defined by ϕ = 0. In this study, for pure metal, T is solved and c = 0 all the time; for alloy, T is prescribed as a constant in the whole domain and c is really computed. This is feasible as thermal diffusivity is usually several orders of magnitude larger than solute diffusivity (Section II. in Ref. [3]). α indicates the preferred growth direction of a dendrite (Fig. 1 in Ref. [22]). It ranges from -π/4 to π/4 because cubic metal with 4-fold symmetry is considered in this study. In liquid, α is meaningless but here is set as 0.

The scalar-valued phase field equation reads [15,22]

The governing equation of orientation field is [15,22]

The governing equations of temperature [1,22] and solute [3,29] are expressed as follows:

In the case of pure metal, Eqs. (1), (2), (3) are solved (D_T: thermal diffusivity); in the case of alloy, Eqs. (1), (2) and (4) are solved. Solute diffusivity ranges from D_cl in liquid to D_cs in solid. In this paper D_cs = 0 as it’s always several orders of magnitude lower than D_cl [3,29].

In Eq. (1), the diffuse interface thickness W is

where a_s (anisotropy) takes the following form in solid:

(ε₄: anisotropic strength). (x₁, x₂) is the location in local frame associated with the orientation of each dendrite (Fig. 1 in Ref. [22]). In liquid, a_s = 0. W₀ (spatial scale) is determined by the ratio W₀/d₀ (d₀: capillary length). This ratio is artificially chosen to facilitate computation in large scale. Coupling constant λ is

(a₁ is a constant, see Table 1). Capillary length is

where c_p and L are the volumetric specific heat and the latent heat. Γ (Gibbs-Thomson coefficient) and T_l-s (freezing range) are

(γ: surface tension, T_m: melting point of solvent, m: liquidus slop, c_∞: the initial concentration of the alloy, k: solute partition coefficient). When interface kinetics is negligible in cases like low-speed solidification, for pure metal (Eq. (90) in Ref. [1]) and alloy (Section IV.C.3 in Ref. [3]), the phase field relaxation time is

Time scale τ₀ is given by:

(a₂ is a constant, see Table 1). θ (dimensionless undercooling of alloy) is expressed with the help of Fig. 1:

θ is associated with Ω (supersaturation):

and (see Fig. 1)

U_T is dimensionless temperature in pure metal, and U_c is dimensionless concentration in alloy:

In Eq. (2), τ_α is the orientation relaxation time (see Table 1) and mobility function P is given by

where s, η, μ_c, β in Eqs. (1), (2) and (19) are constants given in Table 1.

Homogeneous Neumann boundary conditions are imposed on the whole boundary: N∙∇ϕ, N∙∇α, N∙∇T, N∙∇c = 0. N is the outward unit normal to the boundary. Usually, initial ϕ is given by

where d_s is the signed distance from (X₁, X₂) to s-l interface, with d_s < 0 inside and d_s > 0 outside the solid respectively. (X₁, X₂) is location in the global frame (Fig. 1 in Ref. [22]).

2.1.2. Vector-valued phase field model

When dendrites touch each other and grain boundaries begin to form, the orientation periodicity has to be taken into consideration. For example, two grains with α = ±3π/16 should have a misorientation of π/8 (±π/4 represent the same orientation). Without specific treatment, in scalar-valued model, the misorientation would be treated as 3π/8, resulting in incorrect boundary behavior (Fig. 3 in Ref. [22]).

To avoid this, it’s necessary to modify ∇α in Eq. (1) and two additional diffusion equations should be solved (Appendix B.3 in Ref. [2]). This is not easy to implement when finite element method is employed (Section 3.1 in Ref. [22]). Another option to deal with the periodicity is to transform (ϕ, α) into an ordering vector F = (F₁, F₂) [[20], [21], [22]] (this vector was denoted by “(p, q)” in previous studies [[20], [21], [22]]. In this paper, “p” is the pressure):

where α₄ = 4α. Fig. 13(c) illustrates F in an intuitional way.

From Eq. (21), it’s straightforward to acquire the evolution equation,

which will pose problems during finite element discretization (see Eqs. (A.16)-(A.20) in Ref. [22]), however. To overcome this issue, evolution equations are finally split into (see Ref. [22])

and

where A₁ and A₂ stand for the right hand side of Eqs. (1) and (2).

In vector-valued phase field model, during each time step, Eqs. (23) and (24) are solved in turn along the solid line in Fig. 2, followed by Eq. (3) for pure metal or Eq. (4) for alloy. In this study, vector-valued model is adopted to handle the contact between dendrites. Unless otherwise stated, (ϕ, α) shown below is converted from F. Fig. 13(b) gives a view of F distribution near the grain boundary and s-l interface.

In this study, two-phase flow model is used to describe motion of dendrites and flow in melt. Solid is regarded as fluid with extremely large viscosity to have enough resistance to deformation exerted by flow. Combination of two-phase flow with phase field method is straightforward [[30], [31], [32]] by associating viscosity and density with ϕ.

Momentum equation

and continuity equation

are N-S equations to describe flow field (v: velocity, p: pressure, ρf: body force, f: external force acting on per unit mass such as gravitational acceleration, ρ: density).

Solid and liquid together with their interface are treated as incompressible Newtonian fluid with stress tensor σ given by

I is the unit diagonal tensor. Operator D(∙) is defined as follows:

In previous studies on two-phase flow [[30], [31], [32]], dynamic viscosity μ is a linear function of ϕ:

f_s = (1 + ϕ)/2 and f_l = (1 - ϕ)/2 are local solid and liquid fraction. μ_s and μ_l denote viscosity in solid and liquid. A better choice is

which is equivalent to

To make solid rigid enough, we set μ_s = 10¹⁰μ_l throughout this study. The linear expression Eq. (29) is not feasible any more in this case. As displayed in Fig. 3, linear expression Eq. (29) leads to a sharp viscosity increase just after ϕ > -1, which tends to result in an ill-conditioned stiff matrix after discretization. Clearly, exponential expression Eq. (31) corresponds to a relatively steady increase. Eq. (31) combined with the scaling method in Appendix B.2 is able to handle such huge viscosity contrast.

Density of Al is related to local solid fraction:

Because we focus on solidification of Al-4 wt.%Cu in Section 4, μ_l is simply approximated by μ_l = ν_Al(l)ρ_Al(l) (ν_Al(l): kinematic viscosity of liquid Al, ρ_Al(l): density of liquid Al).

Density is associated with concentration and local solid fraction by

where

and the same applies to ρ_Cu. w_Al and w_Cu denote local mass fraction of Al and Cu.

Boundaries in this study include Γ₁ (Dirichlet boundary), Γ₂ (“do-nothing” boundary [[33], [34], [35]]). On Γ₁, both components of velocity are prescribed directly. If v = 0, it is “no-slip” condition. In Section 4, Dirichlet condition was imposed on the entire boundary. In this case, pressure should be fixed at one point to obtain a certain pressure field (Fig. 4.1 in Ref. [36]). We constrain the pressure at the left bottom corner as zero in in Section 4. On Γ₂,

Eq. (35) is meant to constrain mean pressure across this boundary as p_m (Eq. (28b) in Ref. [33]).

Plasticity can be introduced into solid without significant modifications of our current model. Ideal plasticity leads to an equivalent viscosity dependent on ∇v [37,38]:

“:” is the double dot product between two second-order tensors T₁ and T₂:

and σ_y denotes solid’s yield strength; dim is the dimension of the problem. Plasticity will be taken into account in our future work.

The advection equation describes both translational and rotational motion of dendrites with a given flow field:

φ stands for variables in the phase field model including F₁, F₂, T or c. The discretization and solver for Eq. (38) are illustrated in Appendix C. Rotation of a dendrite will change its orientation. The angular rate is

During a step dt, orientation should change ωdt, and F rotates in the following way:

Flow field computation is far more time-consuming than solving other equations due to the former’s complexity (Appendix B). In this study, we design a two-level adaptive time stepping scheme where flow field calculation has a large step size while other computations reserve a fine one. Fig. 4 displays this scheme in a flow chart.

Suppose all variables have been known at t = t_n (n: the number of last outer time step). The first task in current step is to determine dt (inner step size) and n_in (inner loops). Δt (outer step size) is n_indt. To keep the accuracy of advection equation (Eq. (38)), dt is constrained by (according to Section 6 in Ref. [39])

for each cell where ϕ_c > -1 (ϕ_c, h_c, v_c: phase field, local cell size, and velocity at cell center). Eq. (41) aims to prevent solid from moving across 0.1h_c during dt. If dt is only limited by Eq. (41), when dendrites are almost motionless, dt approaches infinity. According to our test, to guarantee stability for solving Eqs. (23) and (24), dt should also be limited by

dt should be as large as possible as long as it satisfies Eqs. (41) and (42). n_in should ensure that in all cells where ϕ_c > -1,

Since we use a variable viscosity depending on ϕ, Eq. (43) is meant to keep solid from moving across a cell during Δt. Similarly, n_in should be as large as possible but limited by n_in ≤ 50 to ensure Δt is a finite value.

Flow field is updated to t_n+1 = t_n + Δt in the outer loop and it is used to advect other variables in every inner loop. In each inner loop, for all variables to be advected, advection equation (Eq. (38)) will end up with the same stiff matrix but different right-hand sides (Eq. (C.6)). As a result, we assemble the stiff matrix of Eq. (38) and apply LU decomposition to it in the outer loop (“advection (1)” in Fig. 4). In the inner loop, we just assemble the right-hand side for different variables and utilize the same decomposition to solve Eq. (38) repeatedly (“advection (2)” in Fig. 4). Then we employ Eq. (40) to “rotate F”. Finally, we’ll solve Eqs. (23), (24) and Eq. (3) or (4) (“growth” in Fig. 4). After all inner loops are finished, we get all solutions at t_n+1.

As stated in Section 3.3 of Ref. [22], to prevent numerical oscillation when solving Eqs. (23) and (24), we maintain cell size as h_min near s-l interface (-0.99 < ϕ_m < 0.99, ϕ_m is the mean value of ϕ on all nodes in a cell). h_min corresponds to the maximum cell level L_max. L_max should be large enough so that

h₀ denotes the size of primitive coarse grid without any refinement. In regions of large velocity variation but far from interface, we limit cell size as h_v = 2^-3h₀. Figs. 12(d) and 13(a) show the mesh from different perspectives.

Polycrystalline growth without flow and motion simulated by vector-valued phase field model has been validated with care in Ref. [22]. This section is focused on testing other parts in our model. By default, in Section 3, f = 0, d₀/W₀ = 0.554 (consequently, according to Eq. (12), D = W₀²/τ₀), c_∞ = 0 (pure metal).

Our FSI model was validated by comparing the computed flow field with the calculation by commercial software Fluent [49] in this section. Forced flow in a curved channel (Fig. 5(a)) was simulated in this section.

Initially, this channel was full of stationary liquid. Liquid of the same properties continued flowing into the domain from bottom. As shown in Fig. 5(b, d), the inflow velocity at position (X₁, 0) on the inflow boundary was prescribed as (0, |sin(2πX₁/l)|) mm/s (l is the length of bottom). For Fluent, no-slip boundary condition was imposed on the concave-convex sidewalls directly. For our FSI model, the same channel was formed by attaching solid on the no-slip straight sidewalls (Fig. 5(b)). “Do-nothing” (Section 2.2) boundary condition was applied on the top to allow for outflow. Only Eqs. (25) and (26) were solved. ϕ was maintained as Fig. 5(b) all the time. In this test, W₀ = 10^-6 m.

Fig. 5 suggests our simulation agrees well with the computation by Fluent, which implies that Eq. (31) is equivalent to imposing no-slip condition between liquid and solid. In addition, the velocity distribution along the line in Fig. 5 (X₁ = 0.25 mm) is plotted in Fig. 6. Generally, these velocity profiles coincide well for the two computations. For Fluent, solid region was not included at all, so there was no velocity value inside the convex wall. For our FSI model, we give an artificially large viscosity in solid to mimic the convex wall. The velocity was almost zero, indicating that μ_s = 10¹⁰μ_l is rigid enough for solid to resist deformation caused by flow.

3.2.1. Diagonal translation of a circular grain

In this section, the methods to solve Eq. (38) (Appendix C) were examined in a manner analogous to that in Ref. [14,39] by simulating motion of a growing circular grain under prescribed flow field and constant undercooling.

Initially, the domain was filled with pure melt and a circular grain as displayed in Fig. 7(a) (α = 0). Dimensionless temperature of the melt was U_T = -0.05 all the time. Eqs. (23), (24), (38) and (40) were implemented. Thermal diffusion (Eq. (3)) and N-S equations (Eqs. (25) and (26)) were not solved. Flow field was uniformly imposed in the whole domain, driving the circular grain to translate diagonally from one corner to the other; then, flow direction was inverted and after the same period of time, the grain returned to its initial location. With such flow field, there must be substance entering domain through boundaries. To ensure it was liquid phase that entered domain, we set (g^F1, g^F2) = (0, 0) (see Appendix C) on inflow boundaries when updating locations of solid using Eq. (38).

In this test, we set ε₄ = 0 so that the grain grew as a circle without anisotropy. With such diagonal flow field, the circle’s center should have reached (0.75, 0.75) l at t = 128 τ₀ and returned to (0.25, 0.25) l at t = 256 τ₀, which agreed well with Fig. 7(a). Since uniform flow field should only result in translation without deformation (strain rate (∇v + (∇v)^T)/2 = 0), the grain was supposed to grow up as if it was motionless. To validate this, we also simulated the case without flow (Fig. 7(b)). The final s-l interfaces in two cases overlapped well (Fig. 7(c)).

3.2.2. Rotation of a growing dendrite

Dendrite’s rotation is accompanied by orientation change, which can be effectively handled by Eq. (40). Its accuracy was examined in the same manner as in Fig. 2 in Ref. [14].

Initially, a solid seed of radius r₀ = 5 W₀, U_T = 0 and α = 0 was placed at the center of a domain of l × l (l = 1000 W₀). The rest of the domain was full of undercooled pure melt (U_T = -0.65). In this test, ε₄ = 0.05. Eqs. (3), (23), (24), (38) and (40) were implemented except flow field computation (Eqs. (25) and (26)). After motionless growth for 50τ₀, a constant rotational flow field

(angular rate ω = π/2500 rad/τ₀) was imposed in the entire domain, driving the dendrite to rotate along with the temperature field (Fig. 8(c)). Again, the strategy in Appendix C (see Eq. (C.5)) guarantees melt with U_T = -0.65 entering the domain, so that temperature far from thermal boundary layer (see Fig. 8(c)) was not affected by inflow, enabling comparison with the motionless growth without any flow.

In Fig. 8(a), we mark what orientation the dendrite should have at each moment according to the prescribed angular rate, confirming the accuracy of Eq. (40) to deal with rotation of dendrite. The applied rotational flow field was supposed to drive melt to rotate with dendrite around center together as one rigid body (strain rate (∇v + (∇v)^T)/2 = 0). As a result, the discrepancy between the growth of rotational and motionless dendrite should be negligible, which is validated by the almost coincided morphology shape in Fig. 8(b). The preserved 4-fold symmetry demonstrates that the orientation change due to rotation of dendrite has been well treated. Otherwise, the dendrite will grow into a distorted shape like Fig. 2(c) in Ref. [14].

In this section, we modeled growth and settlement of a single dendrite and multiple dendrites in a domain (1 mm × 2 mm) full of Al-Cu alloy. The dendritic growth (Eqs. (4), (23) and (24)), flow field computation (Eqs. (25) and (26)), translation and rotation (Eqs. (38) and (40)) were implemented based on the time stepping scheme in Fig. 4 and parameters in Table 2. Solidification happened at constant temperature where Ω = 0.35 (see Fig. 1 and Eq. (16)). For numerical stability, gravity was imposed after t = 0.056 s when solid seed (initial radius r₀ = 5 W₀) had grown up a little bit.

Initially, the solid seed (α = 0) located at (0, 1.8) mm. Because solid has a larger density than liquid, it grew and settled down as shown in Fig. 9(a). Clearly, the downward tip grew fastest because it was on the upstream side [8]. Fig. 10(a) shows symmetric vortices on both sides of the dendrite. Similar flow pattern has been observed in simulation of a sinking cubic block [23,24]. The solute distribution is displayed in Fig. 10(b). During solidification, solute will be rejected from solid and pile up ahead of s-l interface if k < 1 [44]. The excess solute rejected from downward tip was swept away by fluid, leading to steep concentration gradient around the downward tip (Fig. 10(c)) and speeding up its growth. The growth velocity of downward tip increased gradually as a function of time, which can be attributed to continuously increasing drop velocity (Fig. 9(b)). When it approached the bottom, both growth and drop slowed down rapidly.

Multiple dendrites should be simulated for a more realistic model. In this simulation, five solid seeds were distributed in melt at t = 0 (Fig. 11(a)). As recorded in Supplementary material (growth and settlement of five dendrites. Left panel: phase field; right panel: solute concentration. See Appendix E), “soft collision” and “hard collision” took place throughout this process.

The flow pattern and solute distribution in this simulation were much more complex than those in Section 4.1. As shown in Fig. 11(b), dendrites 1 and 5 had already touched the bottom while others continued settling, surrounded by Cu-rich inter-dendrite flow. Dendrite 4 was undergoing rapid rotation because of the strong rotational flow around it. Fig. 11(c) suggests solute boundary layers around dendrites 3 and 4 began to collide (“soft collision”). Soon afterwards, “hard collision” occurred (Figs. 12 and 13). Dendrites 1 and 2 had contacted and grain boundary was forming between them. Fig. 13 shows the local details where these two dendrites touched. F points to the direction of 4α and it has a norm of (1 + ϕ)/2 (Fig. 13(b, c)). This simulation convinces that our model is capable of handling growth, motion and rotation of multiple dendrites, diffusion and advection of solute, impingement among dendrites, contact between dendrite and wall, and formation of grain boundary. The calculation was parallelized among 12 processors (Fig. 12(d)).

We measured run time of the following components (Fig. 14):

(1) Growth (assembling): assembling linear systems for Eqs. (4), (23) and (24).

(2) Flow field (assembling): assembling stiff matrix, right hand side in Eq. (B.36) and matrices to form the preconditioner.

(3) Advection (assembling): assembling stiff matrix arising from Eq. (C.6).

(4) Growth (solving): solving linear systems for Eqs. (4), (23) and (24).

(5) Flow field (solving): solving v and p.

(6) Advection (decomposition): decomposing stiff matrix for Eq. (C.6).

(7) Refinement: setting refinement indicators, refining and coarsening mesh, and transferring solutions between old and new mesh.

Fig. 14 implies that solving advection equation individually doesn’t introduce significant additional cost. Our vector-valued phase field model should have been able to simulate the coarsening of grains after complete solidification (Fig. 8 in Ref. [22]). However, even though parallelized among 12 processors (Fig. 12(d)) in a Dell Precision T7500 workstation and accelerated by two-level stepping scheme (Fig. 4), as shown in Fig. 14, flow field computation was still extremely time-consuming, especially in the presence of huge number of unknowns. There were 386,623 nodes in Fig. 12(d) (3,478,513 unknowns of v and p). As a result, we didn’t continue computation thereafter. In our future work, we can overcome this challenge by a less rigorous convergence criterion for flow field and large-scale parallelization on a cluster.

In this study, we incorporate two-phase flow model into vector-valued phase field model to simulate solidification coupled with melt flow and motion of dendrites. Solid is provided with a tremendous viscosity 10¹⁰ times as large as that in melt. We employ a series of numerical procedures to ensure numerical stability and their effectiveness was confirmed by comparing our result with the calculation by Fluent. Motion is realized by solving advection equation. Orientation change due to rotation has been taken into account. These were examined by modeling solid growing in imposed flow field. Two-level time stepping scheme, adaptive mesh refinement, and parallel computation were used to accelerate calculation. We simulated a single dendrite growing and settling in Al-Cu alloy. It was observed that downward arm grew fastest, which can be well explained in theory. We also simulated settlement of multiple dendrites. “Soft collision” of solute and “hard collision” between dendrites were observed in this process. This study is a step stone to model solidification coupled with deformation and deformation-induced inter-dendritic flow. Introducing plasticity into the two-phase flow model is straightforward without significant modification. We will perform large-scale calculation on a cluster based on the extended model in the future. Experimental validation such as in situ observation is also under consideration to make our model more convincing. Vector-valued 3D phase field model can be transformed from the well-established scalar-valued 3D model [5,45,46], which enables the extension to 3D situation.

This work was financially supported by the National Key Research and Development Program (No. 2018YFA0702900), the National Natural Science Foundation of China (Nos. 51774265 and 51701225), the National Science and Technology Major Project of China (No. 2019ZX06004010), the Key Program of the Chinese Academy of Sciences (No. ZDRW-CN-2017-1), and the Program of CAS Interdisciplinary Innovation Team, and Youth Innovation Promotion Association, CAS.

ψ^v and ψ^p denote the test functions for v and p. Multiplying Eq. (25) with ψ^v and integrating these terms by parts yield

(∙, ∙)_S and <∙, ∙>_∂S are integrations over domain S and boundary ∂S. For two scalars, vectors and second-order tensors respectively, we define

and the same applies to <∙, ∙>_∂S. Substituting Eq. (27) into Eq. (A.1) yields

In Eq. (A.5), the integration on the whole boundary could be subdivided into two parts according to the boundary conditions mentioned in Section 2.2:

N∙σ is the stress on the boundary. On Γ₁ (Dirichlet boundary), ψ^v = 0. As a result,

On Γ₂ (“do-nothing”), Eqs. (27) and (35) lead to the following:

In this study, we set p_m = 0 to allow free inflow and outflow [34]. In this case, on the whole boundary

Then, discretize Eq. (A.5) in time:

v and p are solutions at t_n+1 = t_n + Δt and v^old is the known velocity at t_n. The explicit treatment for the last two terms in Eq. (A.10) can keep the stiff matrix symmetric and avoid Newton’s iteration for the non-linear convective term (v∙(∇v)).

Multiplying Eq. (26) with -ψ^p leads to

We seek approximations for v, p

We adopt Taylor-Hood finite element Q₂² × Q₁ [26] to discretize the flow variables v and p. Q₂² implies v uses a continuous quadratic shape function ψ_j^v for two components, and Q₁ indicates p has a continuous linear shape function ψ_j^p. Z^v and Z^p are solution vectors of N_v and N_p degrees of freedom.

Substitute Eqs. (A.12) and (A.13) into Eqs. (A.10) and (A.11). Eqs. (A.10) and (A.11) are supposed to hold for ψ^v = ψ_i^v (i = 1, …, N_v) and ψ^p = ψ_i^p (i = 1, …, N_p). Then we get a linear system as

where

M^pp and R^p arise from the possible pressure constraint at one point (see Section 2.2). In fact, almost all elements in them are zero. Consequently, linear system (A.14) is very close to the well-known “saddle point structure”:

With a zero block at the bottom right, “saddle point structure” is cumbersome to solve. Solvers and preconditioners robust to positive definite stiff matrix, such as Conjugate Gradient method with Jacobi or SSOR (symmetric successive over-relaxation) preconditioner, usually don’t work in this case (see “Linear solvers and preconditioners” in Ref. [50]). To make matters worse, in this study, the huge viscosity contrast in two phases makes the situation more challenging. This section illustrates the preconditioning strategy we developed to overcome these challenges.

A popular block upper triangular preconditioner [[24], [25], [26],47] for structure like Eq. (A.18) takes the following form:

with

S is known as Schur complement:

When applied as right preconditioner, P is considered optimal (Section 2.3 in Ref. [25] and Section 3 in Ref. [47]):

GMRES-type (Generalized Minimal Residual) solver is believed to converge rapidly for Eq. (B.4) [26,47]. After we get Y, the original solution is recovered by

Eq. (B.2) is of very little practical use since S^-1 = (M^pv(M^vv)^-1M^vp)^-1 is hard to evaluate. Alternatively, we shall seek approximation to S^-1 since P only acts as a preconditioner. How to approximate S^-1 depends on M^vv, which takes different form for different momentum equation [[24], [25], [26],47]. The procedures to find appropriate approximation in this study are introduced as follows.

Eq. (4.8) in Ref. [47r gives the following approximation:

Q^v and Q^p are mass matrices, given by

H^v and H^p are defined as

Eq. (B.6) leads to

According to Section 4 in Ref. [47],

L^p is Laplace matrix for pressure

M^vv can be approximated by H^v (see “Linear solver and preconditioning issues” in Ref. [51]). Eqs. (B.11) and (B.12) give approximation to S:

Then S^-1 can be approximated by

Eqs. (B.10) and (B.15) yield

where

are weighted Laplace and mass matrices for p. Finally, approximation to P^-1 reads

Because of the tremendous viscosity variation, elements in M^vv (Eq. (A.15)) and Q_w^p (Eq. (B.18)) can vary by orders of magnitude, which may have a bad impact on effective representation of (M^vv)^-1 and (Q_w^p)^-1, and deteriorate preconditioner’s performance.

Define d^v and d^p, two scaling vectors containing the diagonal elements of M^vv and Q_w^p:

with a diagonal matrix Λ^p:

The scaled linear system equivalent to Eq. (B.4) is

Z is recovered by

and

In Eq. (B.23)

S_I^-1 is approximated by

where

The row/column scaling strategy shown above guarantees a symmetric stiff matrix.

When assembling L_w^p, to guarantee it’s positive definite, we constrain the pressure of the left bottom point in the domain as zero. Now that M_I^vv, L_w^p and (Q_w^p)_I are symmetric and positive definite, their inverse matrices can be represented by Cholesky decomposition. This pays off because the same decomposition can be used during every outmost iteration. Note that L_w^p is decomposed before scaling (see Eq. (B.33)) because it doesn’t have the magnitude issue as M_I^vv and (Q_w^p)_I (see Eqs. (A.15), (B.17) and (B.18)).

Section B.2 explains the strategy to balance elements inside M^vv and Q_w^p individually. The balance between two blocks M^vv and M^vp should not be ignored, either. Linear system Eq. (B.23) arises from Eq. (25) (momentum conservation) and Eq. (26) (incompressibility constraint) without prescribing the weight between them. If Eqs. (25) and (26) are unbalanced, M^vv and M^vp will also be affected. Since we choose iterative solver as the outmost solver, the solver may concentrate too much on the residual of one of the two equations (see “The scaling of discretized equations” in Ref. [52]).

To avoid this, we multiply Eq. (26) by a factor s_p [26]:

Accordingly, linear system Eq. (B.23) changes into the following equivalent one:

Z is recovered by

and Eqs. (B.25) and (B.26). In Eq. (B.36)

To keep balance between M_I^vv and s_pM_I^vp in Eq. (B.36), s_p is determined by

||∙||_F is the Frobenius norm of the matrix.

Eq. (B.36) is the linear system we actually solved. As the same with Ref. [26], FGMRES (flexible GMRES) serves as the outmost solver. It iterates until

||∙||₂ is the l₂ norm. The tolerance is Tol = 10^-10 throughout this study.

Appendix C. Discretization and solver of advection equation

It’s straightforward to write the weak form of Eq. (38) as

where ψ^φ denotes the test function. Using the streamline diffusion stabilized method [53r, we replace ψ^φ with ψ^φ+δv∙∇ψ^φ:

where δ is the stabilizing parameter as stated in Ref. [53].

On the inflow part ∂S- of the boundary, Eq. (38) should be augmented by boundary condition [53]:

g^φ is a function of the location on ∂S-. It prescribes the properties of the substance entering domain such as temperature (g^T) or concentration (g^c), solid fraction and orientation (g^F1, g^F2). The inflow portion ∂S- of the boundary is featured by

The outflow part must not have any boundary conditions [53]. To incorporate inflow boundary values into weak form, we rewrite Eq. (C.2) as [53]:

Time discretization yields

As stated in Section 2.3, φstands for a series of variables. According to Eq. (C.6), all variables share the same stiff matrix with different right hand sides. As velocity is updated in each outer loop, within which the stiff matrix for advection equation keeps unchanged, we just assemble one stiff matrix and apply LU decomposition to it in the outer loop. This decomposition could be utilized for all variables in each inner loop.

ϕ Phase field

α Orientation field

T Temperature

c Concentration

α₄ 4α

a_s Anisotropy (Eq. (6))

τ₀ Time scale (Eq. (12))

W₀ Spatial scale

τ_ϕ Phasefield relaxation time (Eq. (11))

τ_α Orientation relaxation time

W Diffuse interface thickness (Eq. (5))

γ Coupling constant (Eq. (7))

a₁, a₂ Constants in phase field model

d_0T, d_0c Thermal and solute capillary length (Eq. (8))

d₀ d₀ = d_0T (pure metal) or d_0c (alloy)

ε₄ Anisotropic strength

U_T Dimensionless temperature for pure metal (Eq. (17))

U_c Dimensionless concentration for alloy (Eq. (18))

θ Dimensionless undercooling for alloy (Eq. (14))

Ω Supersaturation (Eq. (16))

Y Defined by Eq. (17)

P Mobility function (Eq. (19))

k Solute partition coefficient

D_T Thermal diffusivity

D_cl, D_cs Solute diffusivity in liquid and solid

D D = D_T (pure metal) or Dc_l (alloy)

c_∞ Initial concentration of the alloy

c_l, c_s Equilibrium concentrations in liquid and solid at temperature T (Fig. 1)

T_m Melting point of solvent

m Liquidus slope

Γ Gibbs-Thomson coefficient (Eq. (9))

γ Surface tension

L Volumetric latent heat

c_p Volumetric specific heat

T_l-s Freezing range (Eq. (10), namely l₂ in Fig. 1)

s, η, μ_c, β Constants in Eqs. (1), (2) and (19)

F=(F₁, F₂) Vector-valued phase field (Eq. (21))

A₁, A₂ Right hand side of Eqs. (1) and (2)

d_s Signed distance used to apply initial phase field (Eq. (20))

v, p Velocity and pressure

f External force acting on per unit mass

σ Stress tensor

Γ₁, Γ₂ Dirichlet and “do-nothing” boundary for flow field

I Unit diagonal tensor

D(∙) An operator defined by Eq. (28)

μ Dynamic viscosity

f_s, f_l Local solid fraction: f_s = (1 +ϕ)/2, liquid fraction: f_l = (1 - ϕ)/2

μ_s, μ_l Dynamic viscosity in solid and liquid

ρ Density (Eq. (33))

ρ_Al, ρ_Cu Density of Al and Cu (Eq. (32))

w_Al, w_Cu Local mass fraction of Al and Cu

ρ_Al(l), ρ_Al(s) Density of liquid and solid Al

ρ_Cu(l),ρ_Cu(s) Density of liquid and solid Cu

N Outward unit normal to the boundary

p_m Mean pressure on “do-nothing” boundary

d^v, d^p Scaling vectors (Eqs. (B.20) and (B.21))

Λ^p Diagonal matrix constructed by dp (Eq. (B.22))

s_p The factor to balance momentum and continuity equation (see Eqs. (B.35) and (B.40))

M^vv,M^vp,M^pv,M^pp Blocks in linear system Eq. (A.14)

M_I^vv,M_I^vp,M_I^pv,M_I^pp Scaled M^vv, M^vp, M^pv, M^pp (Eqs. (27), (28), (29))

Z^v, Z^p Solution vectors of v and p (Eqs. (A.12) and (A.13))

Z_I^v, Z_I^p Scaled Z^v, Z^p (Eqs. (B.25) and (B.26))

R^v, R^p Right hand sides in linear system Eq. (A.14)

R_I^v,R_I^p Scaled R^v, R^p (Eqs. (B.30) and (B.31))

P Defined by Eq. (B.1), preconditioner in linear system Eq. (B.4)

P^-1 Inverse of P (Eq. (B.2)), approximated by Eq. (B.19)

P_I^-1 Inverse of scaled P, approximated by Eq. (B.32)

P_II^-1 Used in the final linear system Eq. (B.36), approximated by Eq. (B.38)

S Schur complement (Eq. (B.3)), approximated by Eq. (B.14)

S^-1 Approximated by Eq. (B.16)

S_I^-1 Scaled S^-1 used to form P_I^-1, approximated by Eq. (B.33)

S_II^-1 Used to form P_II^-1, approximated by Eq. (B.39)

Y Solution vector in linear system Eq. (B.4)

Y_I Solution vector in scaled linear system Eq. (B.23)

Y_II Solution vector in the final system Eq. (B.36)

Q^v, Q^p Mass matrices for v, p (Eqs. (B.7) and (B.8))

L^p Laplace matrix for p (Eq. (B.13))

L_w^p, Q_w^p Weighted Laplace and mass matrices for p (Eqs. (B.17) and (B.18))

(Q_w^p)_I Scaled Q_w^p (Eq. (B.34))

H^v, H^p Matrices used to seek approximation to S (Eqs. (B.9) and (B.10))

Tol Tolerance to judge convergence of flow field solver

φ The variable to be advected (see Eq. (38))

g^φ The boundary condition on ∂S- for variable φ

ω Angular rate (Eq. (39))

N_v, N_p The number of degrees of freedom in finite element space of v and p

Ψ^v, Ψ^p Test functions for v, p

Ψ_j^v, Ψ_j^p Shape functions for the jth degree of freedom in finite element space of v and p

(∙,∙)_S, <∙, ∙>_∂S The integration over domain S or its boundary ∂S (Eqs. (2), (3), (4))

v_c, ϕ_c Velocity, phase field at cell center

ϕ_m Mean value of ϕ on cell nodes

h_c Cell size

h₀ The size of primitive coarse grid

h_min Cell size in regions where -0.99 < ϕ_m < 0.99

h_v Cell size in regions of great velocity variation but far from s-l interface

L_max The maximum cell level (see Eq. (44))

dt, Δt, n_in Inner and outer step size, inner loops (Fig. 4)

t_n, t_n+1 The beginning and end of the (n+1)th outer time step (Fig. 4)

(X₁, X₂), (x₁, x₂) Positions in global and local frame (Fig. 1 in Ref. [22])

t Time

S Computational domain

∂S Boundary of S

∂S_- Inflow part of the boundary

Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.jmst.2020.05.005.