3.3 Numerical implementation
3.3.1 Numerical tests
It is necessary to verify the proposed MPF model as well as the numerical procedure. The simulation of dendritic growth is a
stan-dard verification tool in this field.
The evolution of a small nucleus placed in an undercooled melt is simulated. The form of the nucleus changes from round to den-dritic. This shape depends on conditions and anisotropy in surface tension and the kinetic coefficient. After some time the growing dendritic tip achieves its stationary state, which is determined by the growth velocity and the tip radius. Here we deal with anisotropy in surface tension, γ = γ0α(~n). The so called four-fold symmetry form is taken:
α(~n) = 1−34+ 44 (∇xφ)4+ (∇yφ)4
|∇φ|~ 4 , (3.22) or it is rewritten in terms of the orientation angle:
α(θ) = 1 +4cos 4θ. (3.23) Here θ is the angle between the interface and x-axis. 4 is the anisotropy parameter. The surface tension in the 3D case takes a similar form:
α(~n) = 1−34+ 44(∇xφ)4+ (∇yφ)4 + (∇zφ)4
|∇φ|~ 4 . (3.24) The corresponding expressions for capillary lengthd0 in the 2D and 3D cases are:
d0(~n) = 1−34−124(∇xφ)4+ (∇yφ)4−8(∇xφ)2(∇yφ)2
|∇φ|~ 4 ,
d0(~n) = 1−34 −124 (∇xφ)4+ (∇yφ)4+ (∇zφ)4
|∇φ|~ 4
−4(∇xφ)2(∇yφ)2+ (∇xφ)2(∇zφ)2 + (∇yφ)2(∇zφ)2
|∇φ|~ 4
!
.
Figure 3.2: The 3D dendritic crystal.
St= 0.45,δx= 0.8,∗4= 0.05,D= 4, size:
250×250×250
Figure 3.3: The 3D dendritic crystal.
St = 0.45, δx = 0.6, ∗4 = 0.0081, D = 12, size: 250×250×250
Table 3.1: Dendritic growth. δx= 0.4,τ = 1,
St λ D
√ W d0
ξ v˜ v˜- theory 0.55 0.05 4.70 2 0.277 0.0167 0.0170 0.55 0.05 7.05 3 0.185 0.0178 0.0170 0.45 0.05 9.40 4 0.139 0.00564 0.00545 0.30 0.05 23.5 10 0.055 0.000685 0.00068 0.55 0.02 4.70 2 0.277 0.00694 0.00685
λ= a2 a1
τ DT
ξ⊥2L, Ca= a2ξ⊥2γξ2
τ DT . (3.25)
An example of a 3D dendritic crystal is shown in Fig. 3.2. The color corresponds to the local temperature. The cold regions on the tips indicate the influence of the Gibbs-Thomson effect. Only one eighth part of the whole crystal was simulated.
The simulation results for 2D cases are collected in Table.3.1.
The results show good agreement with the theoretical predictions [26].
0 0.01 0.02 0.03 0.04 0.05
10 100 1000 10000
40 32 24 16 8
ε
Nit
αdt
Figure 3.4: Value of the diabatic pa-rameter and time step vs. time. St = 0.30see Table 3.1. Dashed lines corre-spond to the maximum and minimum allowed values of thediabatic parame-ter . The ratio of the new to the old time step is fixed: αdt= 1.05.
Figure 3.5: Time step and diabatic parameter for the 3D dendritic crys-tal. The dashed line corresponds to the maximal possible value of the time step for the central-forward finite difference method. St = 0.45, δx = 0.8, ∗4 = 0.05, D= 2, size:250×250×250
The evolution of the time step and corresponding value of the diabatic parameter are shown in Fig. 3.4. It is seen that in the initial stage the time step grows, but at some point it decreases and then stabilizes. During the period of decrease, the crystal shape changes from round to dendritic. The birth of dendrites is accompanied by an increase in the temperature gradient near the dendritic tips, which causes an increase in the local diabatic parameter.
The corresponding plot for the 3D case is shown in Fig. 3.5.
The dashed line indicates the maximal possible value of the time step for the central-forward finite difference method. The two times larger time step may not seem like a very significant achievement, especially if we remember that it requires a six times larger mem-ory and ∼ 3 time more floating point operations. However, the computation of the heat transport takes only ≈ 15% of the total
computation time, making an ≈80% increase in performance.
The dendritic crystal shown in Fig. 3.3 has a very small anisotropy. The exhibited morphology corresponds to the dublon type, according to the morphology diagram given by E. Brener et al. [9]. This case was computed with and without symmetry as-sumptions. The simulation results are equal.
3.3.2 Equilibrium shape
The next example is the calculation of the equilibrium shape of a crystal. Knowledge of the character of anisotropy in surface tension comes from analysis of the form of extremely slow growing crystals.
A similar procedure is reproduced in the phase-field simulation.
Solidification with small undercooling and slow growth kinetics is considered. Latent heat is not taken into account.
The anisotropy of the surface energy is taken again in the form of the four-fold symmetry, Eq.(3.22). The corresponding equilib-rium shape can be described by the expression [54]:
x = α(θ) cos(θ)−α0(θ) sin(θ),
y = α(θ) sin(θ) +α0(θ) cos(θ). (3.26) Simulations start from the round nucleus,r0 = 70 l.u. and they are made with the following physical parameters:
Γ = 0.69 l.u., µ= 0.090 l.u.. (3.27) The equilibrium form is achieved after some time, see Fig. 3.6. It is seen that the initial nucleus melts in some regions and grows in others.
A comparison of the phase-field simulation and the analytical solution is shown in Fig. 3.7. This result allows us to conclude that
50 150 250 350
50 150 250 350
Figure 3.6: Evolution of a crystal with the anisotropy = 0.2. Isolines corre-spond to the different time: t= 0, 2000, ..., 12000l.u.
the proposed modification of the phase-field equation reproduces the anisotropy effects with high accuracy, even in the case of strong anisotropy.
−1 1
−1 0 1
0
−1 0 1
−1 0 1
−1 0
1
−1 0 1
−1 0 1
−1 0 1
Figure 3.7: Equilibrium forms of crystals for different values of anisotropy:
= 0.05, 0.1, 0.2, 0.5. The analytical form is marked by a line and the phase-field simulation is plotted with a filled area.
Chapter 4
Phase-field models for binary alloys
The previous chapter demonstrates application of the phase-field model for simulation of the dendritic growth in a pure, undercooled melt. It has been extended to the model for solidification in binary alloys by Wheeler, Boettinger and McFadden [69], the WBM model.
The WBM model that has been used widely [69, 12, 67, 34] is de-rived in a thermodynamically-consistent manner. In this model, any point within the interfacial region is assumed to be a mix-ture of a solid and a liquid, both with the same composition. The phase field parameters can be determined not only through the sharp interface limit condition, but also through the finite inter-face thickness condition [29]. It has been shown that the model can correctly reproduce the solute trapping phenomena at a high interface velocity [12, 29].
A careful study of the WBM model shows that there is limita-tion on the width of the transilimita-tion region [29, 30]. The existence of this limitation pushed Kim, Kim and Suzuki to develop the so
48
called KKS model for binary alloys [30]. This model is derived in a thermodynamically-consistent manner, as well. The difference be-tween the WBM and the KKS model is in the definition of the free energy density for the interfacial region.
Differential equations of the KKS model were solved by the finite difference method. The first attempts to simulate GeSi crystals demonstrated that stability restrictions of the FD method on the time step make simulation of the system impossible. Application of the MRLK scheme to the concentration equation allows one to increase the time step by 10-100 times, in comparison with the FD scheme.
This chapter has the following structure. In the first two sec-tions, a short overview on the theory of phase transition in binary alloys is given. The two following sections contain overviews of the WBM and the KKS models and the numerical procedure for the solving of the phase-field equations is further considered. In the last section, we present an application of the developed technique to the NiCu alloy.