The second category of the simulation is devoted to macroscopic solidification problems.
The first example is the heat flow and the interface shape in a horizontal Bridgman crystal growth, which is a popular process for growing compound semiconductor crystals. A typ-ical configuration is illustrated in Fig. 11a; the crucible is neglected. The material is kept in a molten state on the left (with the temperature TH higher than the melting point Tm), while on the right, the heat is removed for solidification (with a cooler temperature at TC).
Due to extremely low growth rate, the interface movement is neglected. For a stationary state with constant TH(= 1) and TC(= 0) the position and the morphology of the interface
FIG. 10. (a) Final adaptive mesh (27,472 cells with 7 levels of grids;δ = 5 × 10−4) for 18 spheres in a heated square; (b) the calculated velocity and thermal fields. Spacing of isotherm contours is 0.1.
are coupled with the heat transfer and melt flow. This problem has been solved extensively as a benchmark problem for various Rayleigh numbers (Pr= 0.71) by Newton’s method [7] and multigrid schemes [34]. Therefore, it is a good candidate for comparison. The first refined mesh (1922 cells with four levels andδ = 2 × 10−3) based on the interface refinement is shown on the left-hand side of Fig. 11b, while the body-fitted coordinate is shown on the right (a coarse grid); Ra= 2 × 104. The structured grid for comparison is finer being 320× 160 (in the fifth level). The comparison of the calculated isotherms is shown on the left figure of Fig. 11c, where the refinement is based on both the interface and velocity errors; the flow pattern is illustrated by the dashed lines. As shown, the agreement is excellent for both isotherms and the interface shape. Furthermore, the calculated Nusselt number (Nu= 1.980) is in good agreement with the one by the structured grid (1.977);
the error distribution is similar to Fig. 6b. Again, the error contours for velocity are shown
FIG. 11. (a) Schematic of the horizontal Bridgman configuration; (b) initial adaptive mesh (left) having 1922 cells with 4 levels andδ = 2 × 10−3and body-fitted mesh (right); (c) comparison of isotherms (flow is indicated by the dash-dotted lines) and the final mesh and errors (right); spacing for isotherms= 0.05.
on the right figure by the dashed lines; all the errors in each stages of refinement are added. Again, at the beginning, the largest errors (−1.56 and 0.63) are distributed near the boundary. After refinements, the maximum error in the order of 0.5 is distributed uniformly slightly away from the initial refined zone; one can find the small clusters of error contour at|e| = 0.1 there. The final mesh (4010 cells) on the right of Fig. 11c indicates that the refinement indeed follows the error estimators. The total number of cells on the final mesh (4010) is one order smaller than that by the structured mesh (320× 160). As mentioned previously at the beginning of the calculation, a larger interface thickness (6δ) needs to be used. During mesh refinement, one can further reduce this value to reach the sharp interface limit. However, in some cases, ifδ is reduced too fast, the new interface may be out of the
FIG. 12. (a) Schematic of one-dimensional directional solidification; (b) adaptive moving mesh (only half is down) for calculation; the infinite domain is simulated by a finite domain (L= 4); δ = 1 × 10−2.
refined region. Then, the calculation can be greatly retarded and sometimes diverged. On the other hand, without mesh refinement, gradually reducing the interface thickness during iterations is useless.
In the previous case, there is no heat of fusion released due to the steady-state assumption.
To further examine the effect of the heat of fusion, we also perform calculation for a one-dimensional solidification problem as shown in Fig. 12a, where its analytical solution on the interface position d(t), the solid (Ts), and liquid (Tl) temperatures can be derived [38, 39]
(TC= 0, TH= 1, and Tm= 0.5) as d(t) = 2γ√
t; Ts = Tm
erf(γ )erf(x/√
4t); (26)
Tl = 1 +Tm− 1
erf(γ )erf(x/√ 4t), whereγ is a constant and can be calculated implicitly by
Tm= erf (γ )
1+
√π
St γ exp(γ2)erfc (γ )
. (27)
A sample moving mesh is shown in Fig. 12b (1826 cells). In this case, the solid phase does not exist at t= 0−, but this does not cause any trouble in calculation. In order to take a bigger time step for integration, the refined zone cannot be too small. Otherwise, the interface can easily move out of the refined zone leading to the failure of refinement. If the interface
FIG. 13. (a) Comparison of calculated isotherms (dashed lines) at different times with the analytical results (solid lines); St= 1 and δ = 1 × 10−2; (b) comparison of calculated interface position (symbols) as a function of square root of time at different Stefan numbers with the analytical results (solid lines);δ = 1 × 10−2.
speed is low, a much thinner refined zone can be used. As will be illustrated shortly, this is the case for the dendritic growth using the phase-field model. The calculated isotherms at different times are shown in Fig. 13a. As shown, the agreement with the analytical solution is excellent. The smaller thermal gradient in the melt side of the interface at longer times is due to the release of the heat of fusion. The positions of the interface at different Stefan numbers (St) are also plotted against with the square root of time, and they happen to be linear (as shown by the first formula of Eq. (26)). Again, the agreement is also quite satisfactory.
Extension to the solidification in a crucible (Pr= 0.71 and Ra = 105), as shown in Fig. 14a, is also straightforward. Initially, the melt is kept at TH= 1. Due to zero thermal gradients, the melt has no convection. When t= 0+, we set the walls to be cold at TC= 0, except the upper one, which is still assumed to be adiabatic. In addition to the solidification at the walls (with a new phase formation), melt convection also starts. The calculated result showing isotherms and flow pattern (dashed-line) at t= 0.05 is shown in Fig. 14b (on the left), and the mesh is shown on the right (about 8000 cells). The larger isotherm spacing in the melt is also due to the heat of fusion as well as the melt convection (St= 1 and δ = 5 × 10−3). Again, in this case, we have purposely using a bigger refined zone in order to use a larger time step for integration. The refinement on the velocity error is also performed, as illustrated in Fig. 14c. Again, more cells are allocated near the interface, where the momentum boundary develops (about 500 cells are added). We have put the isotherms with and without this refinement on the left figure of Fig. 14c for comparison, but the difference is not significant. Because the upper boundary is adiabatic, the isotherms at the top are perpendicular to the top wall. At the end of solidification, the melt convection also becomes weaker due to the smaller space for convection. The thermal gradients become much smaller as well. The calculation can continue until the temperature of the whole domain approaches TC. The disappearance of the melt and the interface do not cause any problems during calculation.