When materials crystallize or solidify, often due to changes in temperature or concentration, they usually form complex, branchlike patterns known as dendritic structures, similar to those seen in snowflakes or frost on windows. The phase-field method excels in modeling complex interfaces in science and engineering due to its simplicity and flexibility, making it a preferred choice for modeling dendritic solidification processes. This approach's usefulness in portraying the evolution of dendritic structures has been demonstrated in seminal works by Kobayashi [14] and further explored in numerous subsequent studies in [10][11][12][13]15,16,40]. Utilizing the phase-field approach for dendritic solidification modeling entails the introduction of a scalar phase-field variable, which serves to distinguish between solid and liquid states, forming the basis for a detailed free energy model of the system. The system's total free energy is commonly postulated to comprise three fundamental components: configurational entropy, which dictates spatial
which helps identify each phase. To address the discontinuity in ÿ, it is assumed that there is a thin, smooth transition region at the interface, characterized by a width on the order of ÿ(ÿ). Moreover, we define another scalar function ÿ (ý, ý), ý ∈ Ω to denote the temperature field across the domain, which is expressed in the dimensionless form.
Hence, the slightly-refined model for flow-coupled anisotropic phase-field dendritic solidification, which adheres to the modified energy law, as developed in [4,39] (with the original formulation presented in [2,25]), is described as follows.
1ÿ 2
)
)
(2.5)
In the above context, ÿ represents the average fluid velocity field, ÿ denotes the pressure, ÿ indicates the kinematic viscosity, ÿ defines the scaled characteristic interfacial width, ÿ(ÿ) = (1ÿ 2 )(1 + ÿ), ÿ represents the diffusion rate of the temperature, ÿ and ÿ are associated with the "latent heat" (cf. [14]), the function ÿ (ÿ) = ý ′ (ÿ) with ý (ÿ) = (ÿ 2 -1) 2
, describes a double-well type potential, ÿ(∇ÿ) is the anisotropic coefficient, ÿ represents a positive mobility parameter and is set to be a constant (cf. [13,26]), ý(ÿ) = (1ÿ 2 ) 2 specifies that the transfer of latent heat is confined to the liquid-solid interface (cf. [12]), and the term ÿ(ÿ) corresponds to the variational derivative of ÿ(∇ÿ). The last term -ÿÿÿ(ÿ)ÿ, within the momentum equation (2.2) contributes additional viscous dissipative stress, acting as a dispersed momentum sink in the diffuse interface region (cf. [2,25]). Notably, this term is absent in the homogeneous solid and liquid regions, exclusively affecting the interface between these two phases.
Remark 2.1. In the original model developed by Beckermann et al. in [2], the momentum equation is given as
However, the formulation (2.6) does not satisfy the energy dissipation law, a fact that becomes evident when performing an ÿ 2 -inner product of (2.6) with the velocity field ÿ. Thus, in our previous work in [4,39], we introduced a minor yet crucial modification to the momentum equation, see (2.2). This amendment ensures compliance with energy dissipation laws and remarkably replicates the dendritic morphologies observed in Beckermann's original model, see the details in [4,39].
Generally, the definition of ÿ(∇ÿ) reads as ÿ(∇ÿ) = 1 + ÿ 4 cos(ÿΘ),
where ÿ represents the count of anisotropy folds, 0 f ÿ 4 f 1 is a parameter indicating the magnitude of anisotropy, and Θ = arctan(
). We also have ÿ(∇ÿ) g 1ÿ 4 .
When ÿ = 4 (four-fold anisotropy), the following format of ÿ(∇ÿ) is also frequently adopted:
) for 2D case;
) for 3D case.
(2.8)
For ÿ = 4, by taking the variational derivative of ÿ(∇ÿ), we derive ÿ(ÿ) as
) for 3D case.
(2.9)
Without losing generality, the Dirichlet-type boundary conditions are used for ÿ, and no-flux boundary conditions are used for ÿ and ÿ , specified as follows:
where ÿ is the normal vector pointing outward from the domain boundary. Periodic-type boundary conditions are commonly applied to all variables when simulating the model without incorporating the flow field, see [4,10,11,15,16,40]. Additionally, the initial conditions for the system (2.2)-(2.5) are specified as follows:
(2.11)
The system (2.2)-(2.5) represents a binary mixture where phase transitions occur between the liquid and solid phases. For the no-flow case (by setting ÿ = 0), the derivation of the PDE system (2.4)-(2.5) utilizes the variational formulation (VF) method, as ex- tensively described in [12]. This VF method can be viewed as a hybrid variational approach. It incorporates Allen-Cahn dynamics-an ÿ 2 -gradient flow mechanism-to formulate the ÿ-equation from the system's total free energy ý(ÿ, ÿ, ÿ ). Meanwhile, the temperature equation for ÿ employs the ÿ -1 -gradient flow, often referred to as Cahn-Hilliard dynamics, which is applied to the enthalpy in the dimensionless form, denoted as ÿ -ÿ (ÿ)
, where ÿ ′ (ÿ) = ý(ÿ). These derivations are extensively detailed in [12]. Meanwhile, the parameters used in (2.4)-(2.5) are essentially a re-scaling of those used in [12].
When incorporating the flow field ÿ, it becomes essential to restrict the fluid momentum equation to the liquid phase exclusively.
The model (2.2)-(2.5) adopts a straightforward approach, a weight factor
, which is used to enforce a limitation, ensuring that the momentum of the fluid is considered only in the liquid phase. Specifically, inside the solid phase with ÿ equal to 1, both the momentum equation (2.2) and the incompressibility condition (2.3) become nullified. Within the liquid phase with ÿ equal to -1,
3) reduce to the classical Navier-Stokes equations. This methodology was initially introduced in [2] and has since been widely adopted in phase-field modeling of complex systems, such as the binary mixture of liquid crystal and viscous fluid, see [22,36,37].
We now show that the model for dendritic growth, (2.2)-(2.5), admits the energy dissipative law, which is achieved by conducting a standard energy estimate process, described as follows. Before presenting the proof, we introduce several notations that will be consistently used in this article. The ÿ 2 -inner product of any two functions ÿ(ý) and ÿ(ý
Additionally, the ÿ 2 -norm of ÿ(ý) is given by ‖ÿ‖ = (ÿ, ÿ)
1 2 . Lemma 2.1. The system (2.2)-(2.5) satisfies the following energy law:
where
(
Proof. First, by taking the ÿ 2 -inner product of (2.2) with ÿ, and using the divergence-free condition given in (2.3) and integration by parts, we get
Second, by performing the ÿ 2 -inner product on (2.4) with ÿ ý , and applying integration by parts, we derive
Third, by computing the ÿ 2 -inner product of (2.5) with ÿ ÿÿ ÿ and using integration by parts, we derive
.16) Combining the three obtained equalities (2.14)-(2.16), and noting that the term I vanishes due to the divergence-free condition (2.3) and the boundary condition 1-ÿ 2 ÿ ⋅ ÿ| ÿΩ = 0, terms II and III cancel each other out, and the term IV also vanishes by the same reasoning as term I, we arrive at (2.12).
A few remarks are given as follows to clarify the details related to the energy law (2.12) within the model (2.2)-(2.5).
Remark 2.2. For a closed thermodynamical system, it's established that energy dissipation must occur. However, the final two terms on the right side of (2.12) may be non-negative, which casts doubt on the certainty of energy dissipation. Thus, the energy law as presented in (2.12) holds "physical" merit primarily under ideal conditions where the variable ÿ consistently stays within the [-1, 1]
Journal of Computational Physics 524 (2025) 113737 range, facilitating the diminution of total free energy. Yet, from a mathematical perspective, proving that the solution to this PDE system adheres to the maximum principle poses a significant challenge.
Similarly, designing a numerical scheme that guarantees the numerical solution for ÿ remains bounded within [-1, 1], thereby adhering to the discrete maximum principle, presents an equally daunting challenge. Therefore, rather than employing advanced techniques to strictly control the range of ÿ in the numerical solution, the "clamping" method is utilized in the design of the discrete numerical scheme to address the issue of potential non-dissipative energy. This approach involves constraining ÿ within a slightly narrower range of [-1 + ÿ, 1ÿ], with a small positive control parameter ÿ ≪ 1, as detailed in (3.35), effectively sidestepping the complexities associated with ensuring maximum bounds.
Remark 2.3. The system exhibits a notable characteristic termed the "zero-energy-contribution" property, as demonstrated in Lemma 2.1. It is found that the application of both the divergence-free condition and boundary condition results in the elimination of terms I and IV, while terms II and III are observed to counteract each other. This leads to the insight that, upon the execution of specific inner products, these terms do not affect the energy law, making it fitting to categorize them under the "zero-energycontribution" label.
The implementation of this characteristic in the design of algorithms for other coupled models is explored in the authors' previous studies, cf. [31,33,35]. However, adapting this strategy to the dendritic model discussed in this paper encounters significant obstacles, particularly due to the unique nonlinear terms and their coupling mechanisms between variables. Therefore, despite the established methodology, the distinct complexities inherent to this specific model complicate the direct application of this approach.
Now, we are ready to design an efficient fully discrete scheme using the finite element method for solving the model of flowcoupled, melt-convective dendritic solidification. We expect that the scheme can possess several "desired" attributes: unconditional energy stability, linearity, second-order temporal accuracy, and fully decoupled structure. Achieving all these favorable characteristics within a single algorithm is far from trivial, given the inherent complexities of the model. These complexities include various elements such as the nonlinear and/or coupled terms related to ÿ (ÿ), ÿ(∇ÿ), ý(ÿ), and 1-ÿ 2 , among others. Each of these factors contributes to the model's overall complexity and represents significant challenges in the development of robust and effective computational methodologies.
To develop the "desired" scheme, we start by recasting the PDE system into a version that is more amenable to algorithmic implementation, yet remains equivalent to the original formulation.
We reconfigure the PDE system (2.2)-(2.5) by integrating several nonlocal type auxiliary variables into the framework. It's imperative to emphasize that this adjusted PDE system maintains its equivalence to the original equations, serving as the cornerstone of our methodology.
To begin, we introduce a new auxiliary variable of the nonlocal type, denoted by ý (ý), aimed at converting the nonlinear portion of the energy potential into a quadratic form. The variable ý (ý) is defined in the following manner:
where ÿ 1 and ÿ 2 are two predetermined positive constants with 0 < ÿ 1 < (1ÿ 4 ) 2 . Since ý (ÿ) is a quartic polynomial with a leading term that is positive, it ensures that, upon the extraction of a quadratic component related to ÿ 2 , this potential can consistently provide a lower bound. Furthermore, considering that ÿ 2 (∇ÿ) g (1ÿ 4 ) 2 , the radicand in (3.1) is guaranteed to be bounded from below. To enhance positivity further, we incorporate the constant ý into the radicand. Utilizing the auxiliary nonlocal variable ý (ý), and transforming the nonlinear segment of the energy potential to be ý 2 , is known as the SAV method (cf. [20,24,28,38,39]).
Using the new variable ý , we reformulate (2.4) into an equivalent form:
where
It is immediately apparent that there is an equivalence between (3.2) and (2.4). It is because integrating the second equation from (3.2) and considering the initial condition of ý | (ý=0) allows the complete restoration of the formulation in (2.4). This step highlights the inherent equivalence of the two expressions.
Next, we define another nonlocal auxiliary variable, ý(ý), and construct a particular ODE for its temporal evolution, which is expressed as follows:
The distinct feature of the above ODE lies in the fact that every term on the right-hand side of (3.4) is composed of ÿ 2 -inner product terms, and their sum is exactly zero, which is confirmed by the proof in Lemma 2.1, namely, ý ý = 0. Therefore, since ý| (ý=0) = 1, (3.4) directly leads to the solution ý(ý) = 1.
Third, we introduce the third auxiliary variable, ý(ý), also of the nonlocal type, and design a specific ODE for its time evolution, that is formulated as
From the proof of Lemma 2.1, we can see that ý ý = 0. Thus, using ý| (ý=0) = 1, the solution of the ODE system (3.5) is also trivial, namely, ý(ý) = 1.
For the model described by equations (2.4)-(2.5), we utilize (3.2) to replace (2.4). Furthermore, the nonlocal variables ý and ý, along with their respective ODEs, (3.4) and (3.5), are used to alter particular terms in the system. This comprehensive combination culminates in a revised system, detailed as follows:
)
)
)
with the boundary conditions specified in (2.10) and the following initial conditions:
(3.13)
The modified system (3.6)-(3.12), incorporating three auxiliary nonlocal variables ý , ý, ý, retains its equivalence to the original model (2.2)-(2.5). The terms labeled "Q-reform" and "R-reform" essentially remain unchanged since ý = 1 and ý = 1. Consequently, these specific modifications involving ý , ý, ý to formulate an equivalent system serve as the foundational "algorithm-friendly" version, crucial for the development of the desired scheme aimed in this study. Indeed, it is also worth noting that reformulation is a strategy frequently utilized in the design of numerical algorithms, providing effective solutions to a variety of complex challenges. Examples include the modification of advective terms to enhance stability (cf. [1,8]), altering the divergence-free condition to circumvent solving the pressure Poisson equation (cf. [9,23]), and redefining energy potentials as a quadratic form in the context of IEQ/SAV methods (cf. [21,27,28]). The flexibility provided by the "reformulation" approach is critical in developing the "desired" scheme to solve the highly coupled and nonlinear dendritic phase-field model in this article.
Given that the system (3.6)-(3.12) is equivalent to the original system, it is expected to adhere to the same energy law in Lemma 2.1.
However, establishing the energy law for this "algorithm-friendly" version is primarily aimed at guiding the discrete scheme; thus, it remains essential to explicitly demonstrate how the energy law is preserved within the restructured system (3.6)-(3.12). This process is detailed as follows.
Lemma 3.1. The system, (3.6)-(3.12), upholds the following energy law:
where
(3.15)
Proof. By taking the ÿ 2 -inner product of (3.6) with ÿ, and using integration by parts and (3.7), we derive
By taking the ÿ 2 -inner product of (3.8) with ÿ ý and using integration by parts, we get
By taking the ÿ 2 -inner product of (3.9) with ÿ ÿÿ ÿ and using integration by parts, we obtain ÿ 2ÿÿ
.18) By multiplying (3.10) with 2ý , we get ý ýý |ý | 2 = ýý (ÿ, ÿ ý ) ⏟ ÿÿÿĀ ÿÿÿ⏟ VI 1 . (3.19) By multiplying (3.11) with ý, we get
By multiplying (3.12) with ý, we get
.21) By combining (3.16)-(3.21) and recognizing that terms grouped under the same Roman numeral labels, such as I 1 and I 2 , cancel each other, we derive the energy dissipation law (3.14). □ Remark 3.1. The difference between Lemma 2.1 and Lemma 3.1 in deriving the energy law emphasizes the rationale for the specific modification from the original system, (2.2)-(2.5), to its new equivalent version, (3.6)-(3.12). This modification is primarily centered around the different methods of cancelling the nonlinear terms. Using the latent heat coupling term ý(ÿ)ÿ from the original equation (2.4) as an illustrative example, it is significant to note that during the derivation of Lemma 2.1, term ýý from (2.15) is counteracted by term ýýý from (2.16
). This relationship necessitates the careful discretization of ý(ÿ)ÿ and ý(ÿ)ÿ ý to ensure that these terms effectively negate each other, often resulting in the development of coupled-type numerical algorithms.
However, shown in Lemma 3.1, the process of deriving the energy law for the modified system (3.6)-(3.12) employs a distinctly different approach to cancellation. The terms III 1 in (3.17) and V 1 in (3.18) are not required to nullify each other. Instead, term III 2 in (3.20) effectively cancels out III 1 , and similarly, term V 2 in (3.20) negates V 1 . This variation allows for different discretization methods to be applied to ýý(ÿ)ÿ and ýý(ÿ)ÿ ý when developing numerical schemes, facilitating the creation of a fully decoupled type scheme.
In this subsection, we design a fully discrete finite element scheme for solving the modified model (3.6)-(3.12), which is the equivalent system of the original model (2.2)-(2.5). The scheme is expected to be the "desired" type, that is, linear, fully decoupled, unconditionally energy stable, and second-order accurate in time. We denote ÿý > 0 as a time step size and ý ÿ = ÿÿý for 0 f ÿ f ý with ÿ = ýÿý. Let ÿ ÿ be the numerical approximation to the function ÿ(⋅, ý)| ý=ý ÿ . Some finite dimensional discrete subspaces are introduced as follows. Assuming that the polygonal/polyhedral domain Ω is discretized by a conforming and shape regular triangulation/tetrahedron mesh ℎ that is composed by open disjoint elements ÿ such that Ω = ⋃ ÿ∈ ℎ ÿ. We use ý to denote the space of polynomials of total degree at most ý and define the following finite element spaces:
where
Besides, we assume the pair of spaces (ý ℎ , ÿ ℎ ) satisfy the inf-sup condition [7]:
where the constant ÿ only depends on Ω. A well known inf-sup stable pair (ý ℎ , ÿ ℎ ) is the Taylor-Hood element [7].
The semi-discrete formulations of the system (3.6)-(3.12) in the weak form read as: find (ÿ, ÿ , ÿ, ÿ
)
)
for (Φ, Θ, ÿ, ÿ) ∈ ý ℎ × ý ℎ × ý ℎ × ÿ ℎ . We now construct the fully discrete scheme as follows. Given that the values of ÿ ÿ ℎ , ÿ ÿ ℎ , (ÿ
1 are known, we update ÿ ÿ+1 ℎ , ÿ ÿ+1 ℎ , ÿÿ+1 ℎ , (ÿ ℎ ÿ ℎ ) ÿ+1 , ÿ ÿ+1 ℎ , ý ÿ+1 , ý ÿ+1 , ý ÿ+1 , through the following three-step process:
Journal of Computational Physics 524 (2025) 113737
)
(3.36)
where
Several remarks are given to explain the scheme (3.31)-(3.38).
Remark 3.2. In the design of the above discretization approach, we skillfully combine explicit and implicit treatments to deal with nonlinear terms. For those terms involving nonlocal variables ý and ý, the implicit treatment is used to discretize ý and ý, while explicit treatment is used to discretize others (such as ý ÿ+1 ý * ÿ * ℎ in (3.31)). This hybrid discretization approach is instrumental in achieving a fully decoupled scheme.
Moreover, we employ the SAV method to convert the nonlinear portion of the free energy into a quadratic form (i.e., ý 2 ), introducing a novel nonlocal variable ý to facilitate the linearization of nonlinear terms. But our approach diverges from the SAV approach by applying a purely explicit treatment to the resultant nonlinear term ÿý (see ÿ * ý * in (3.31)), whereas the SAV method usually discretizes it as ÿ * ý ÿ+1 . To underscore this distinction, we refer to our methodology as the explicit-SAV method. Remark 3.3. We develop a projection-type approach to facilitate the decoupling of velocity field and pressure computations. The projection method has undergone extensive analysis for traditional Navier-Stokes equations that do not include a weight factor, see the detailed analysis in [3,19]. Following a comparable way, when solving the fluid momentum equation with the weight factor, the intermediate velocity ÿÿ+1 ℎ is first obtained in step 2, and then it is projected onto the divergence-free space. It's important to highlight that the velocity form resulting from the projection step consistently encompasses this weight factor due to the specific form of the divergence-free condition.
Additionally, in the context of traditional Navier-Stokes equations, it has been established that the projection-type approach can achieve second-order temporal accuracy for the velocity field while being limited to first-order accuracy for the pressure. The diminished precision in pressure is associated with the imposition of homogeneous Neumann-type boundary conditions to the pressure, thoroughly discussed in [6]. The proposed scheme numerically exhibits the same behavior, as detailed by numerical tests in Section 4.
Remark 3.4. The final velocity field with the weight factor, ÿ ÿ+1 ÿ ÿ+1 ℎ , can be demonstrated to satisfy the divergence-free condition in a discrete manner. Namely, by taking the ÿ 2 -inner product of (3.37) with ∇ÿ ℎ , ÿ ℎ ∈ ÿ ℎ , we get Remark 3.6. Due to the term ÿ * , which comprises a significant number of explicitly treated nonlinear and second-order spatial terms, the stabilization terms associated with ÿ 1 and ÿ 2 in (3.31) are crucial in ensuring the ÿ 1 -stability of ÿ. Explicit treatment of nonlinear terms is well-known for potentially causing instability. To counter this, using linear terms as stabilizers, which match the spatial order of the nonlinear terms, is a widely recognized approach, underlying the incorporation of the ÿ 1 term. Similarly, the ÿ 2 term is used to offset the explicitly handled term ÿ (ÿ) within ÿ * . For an extensive review of linear stabilization techniques in the context of gradient flow models, especially those characterized by considerable stiffness, we refer to [24].
We now proceed to demonstrate that the scheme
where
(3.44) We rewrite (3.38) as ÿ ÿ+1 ℎ ÿÿ+1 ℎÿ ÿ+1 ℎ ÿ ÿ+1 ℎ = 2ÿý 3 ∇(ÿ ÿ+1 ℎ ÿ ÿ ℎ ).
Using (3.46), we derive
where the last equality is derived using the following identity
(3.48)
We reformulate the projection step (3.38) as
(3.49)
Taking the ÿ 2 -inner product of the above equation with itself, we derive
Hence, by multiplying 2ÿý 2 ∕3 to the above equation, we derive
We rewrite (3.38) again as
(3.52)
By taking the ÿ 2 -inner product of (3.52) with 3 2 ÿ ÿ+1 ℎ ÿ ÿ+1 ℎ and using (3.41), we derive
We combine (3.44), (3.47), (3.51), and (3.53) to obtain
(3.67)
It is very easy to solve (3.66) since all terms contained in ÿ 1 and ÿ 2 are already obtained previously.
One important issue is the solvability of (3.66), which can be illustrated by showing
3 2ÿýÿ 2 ≠ 0 as follows. We set Φ ℎ = 3 2ÿý ÿ ÿ+1 ℎ2 in the second equation of (3.64) to deduce -(ÿ * ý * , 3ÿ ÿ+1 ℎ2 2ÿý ) -ÿ ÿ (ý * ÿ * ℎ , 3ÿ ÿ+1 ℎ2 2ÿý ) = 9ÿ 4ÿý 2 ‖ÿ ÿ+1 ℎ2 ‖ 2 + 3ÿ 1 2ÿý ‖∇ÿ ÿ+1 ℎ2 ‖ 2 + 3ÿ 2 ÿ 2 2ÿý ‖ÿ ÿ+1 ℎ2 ‖ 2 g 0. (3.68) By multiplying the second equation in (3.65) with 3 ÿý ý ÿ+1 2 , we get ý ÿ+1 2 (ÿ * , ÿ * ý ) = 3 ÿý |ý ÿ+1 2 | 2 g 0. (3.69) By setting Φ ℎ = ÿ ÿÿ ÿ ÿ+1 ℎ2 in the fourth equation of (3.64), we get
The combination of (3.68), (3.69) and (3.70) gives the fact of -ÿ 2 g 0. Thus (3.66) is always solvable. After ý ÿ+1 is obtained, we proceed to update ÿ ÿ+1 ℎ , ÿ ÿ+1 ℎ , ý ÿ+1 according to (3.63).
Step 2: we split ÿÿ+1 into a linear combination form, using the nonlocal variable ý ÿ+1 , which is expressed as follows:
.71) We solve ÿÿ+1 ℎ1 and ÿÿ+1 ℎ2 as follows. By substituting ÿÿ+1 ℎ with (3.71) in (3.35), and subsequently dividing the resulting equation based on the nonlocal variable ý ÿ+1 , we derive the following two sub-equations: We continue to solve the nonlocal variable ý ÿ+1 . By using the splitting form given in (3.71), we formulate (3.36) into the following form:
where
) . accuracy comparison between the developed second-order ESAV scheme and the first-order decoupled scheme developed in [39] for the anisotropic case ÿ 4 = 0.25.
The initial conditions are specified in (4.1), and other model parameters are detailed in (4.2).
size of ℎ = 2ÿ 512 , we show the ÿ 2 -numerical errors for all variables at ý = 1 both for the isotropic case (ÿ 4 = 0) and the anisotropic case (ÿ 4 = 0.25). The convergence plots reveal that our ESAV scheme achieves a second-order convergence rate for the variables ÿ, ÿ , and ÿ, and a first-order convergence rate for the pressure ÿ .
We continue to assess the convergence order regarding spatial discretizations. Similar to the tests of temporal error, to isolate the impact of spatial discretization on the error, we need to adopt a sufficiently small time step, ÿý = 1e-5, to ensure the temporal discretization error remains negligible. Fig. 4.2 depicts the computed numerical errors for different mesh sizes ℎ, specifically in the anisotropic case with ÿ 4 = 0.25. The results indicate that the velocity in the ÿ 1 norm, along with the pressure, phase-field variable, and temperature in the ÿ 2 norm, exhibit a second-order convergence rate for ℎ. The velocity in the ÿ 2 norm shows a third-order convergence rate for ℎ. These obtained convergence rates align with the anticipated theoretical accuracy for (ÿ, ÿ ) using the ÿ 2∕ÿ 1 elements, and for ÿ, ÿ using the ÿ 1 elements.
Using the same example, we illustrate the advantage of the second-order ESAV scheme developed in this article in terms of accuracy by comparing it with the first-order decoupled scheme developed in [39]. It can be observed that when the adopted time step size is relatively large, such as ÿý = 0.01, the accuracy difference between the first-order and second-order schemes is not significant, with the error generated by the second-order scheme being slightly smaller than that of the first-order scheme. As the time step size decreases, the error curve of the scheme in [39] follows a first-order convergence rate, whereas the proposed ESAV scheme achieves a second-order convergence rate. For instance, when the time step size is 1.25e-3, the error obtained using the scheme in [39] is 9.6e-5, while the error achieved by the second-order ESAV scheme is only 3.9e-6, which is approximately only 4% of the former. This demonstrates the significant accuracy advantage of the second-order scheme over the first-order scheme in [39]. In addition to accuracy, computational efficiency is a key advantage of the proposed second-order ESAV scheme. While the scheme in [39] is In the ÿ-direction, the boundary conditions for the variables ÿ, ÿ , and ÿ = (ÿ, ÿ, ý) are set as follows, ÿ ÿ ÿ| ÿ=0,ÿ = 0, ÿ ÿ ÿ | ÿ=0,ÿ = 0, ÿ| ÿ=0,ÿ = 0, ÿ| ÿ=0,ÿ = 0, ý| ÿ=0,ÿ = ý 0 .
(4.8)
We first perform simulations for the no-flow case, varying the latent heat parameters to ÿ = 1, 1.25, 1.5. Snapshots of the isosurface {ÿ = 0} before the dendritic branches touch the boundary are shown in Fig. 4.7(a), where the dendrites exhibit a maple-like shape. As ÿ increases, the branches become noticeably thinner. We then couple a flow field by setting ý 0 = -0.01 in (4.8) and re-conduct simulations. Snapshots of the isosurface {ÿ = 0} before the fastest-growing tip touches the boundary are shown in Fig. 4.7(b). The flow field significantly alters the crystal structure; similar to the 2D simulation, the branch oriented towards the flow field grows the fastest.
This article designs the first, efficient fully discrete numerical scheme for solving the hydrodynamically coupled phase-field dendritic solidification model. The scheme is anchored on the innovative explicit-SAV method, which is instrumental in the linearization of nonlinear potentials and the establishment of a decoupled structure. This method capitalizes on a distinctive property-"zero-energy contribution"-exhibited by the coupled nonlinear terms. In addition, the approach involves the definition of multiple auxiliary variables, a strategic move that aids in the reformulation of the governing system. This systematic restructuring is pivotal in achieving the desired computational efficiency and accuracy. The scheme's decoupled architecture permits the independent calculation of all variables at every time step, ensuring an operationally efficient framework. The robustness of the approach is also reinforced by its unconditional energy stability and the effectiveness has been demonstrably proven through a series of comprehensive numerical tests conducted in both 2D and 3D, underscoring its applicability and reliability in practical scenarios. To the authors' knowledge, this article's numerical algorithm represents a pioneering achievement-it is the first of the "desired" type fully discrete algorithm for solving this highly nonlinear and coupled phase-field dendritic model.
Journal of Computational Physics 524 (2025) 113737