We have presented some examples for solidification using the enthalpy model, where the phase-field variable is estimated from the temperature as well as the hyperbolic tan-gent function (Eq. (5)). This approach is suitable to macroscopic solidification, where the capillary contribution due to the curvature of the interface can be ignored. However, for microscopic solidification, such as the dendritic growth, interface curvature and thus the interfacial energy play an important role. In such a situation, the local interface tempera-ture is no longer the equilibrium melting point from the phase diagram. One has to incor-porate the Gibbs–Thompson equation [15, 40] into account. Therefore, other approaches, such as the phase-field simulation or the level set method can also be used to obtain the phase-field variable. To this kind of simulation, the present numerical method can be adopted easily as well. In fact, we have found that the solution becomes much easier due to the more diffusive interface obtained by the phase-field equation. Therefore, the final examples are devoted to phase-field simulation of dendritic growth. The model developed by Karma and Rappel [40] is also adopted, which is thermodynamically consistent. For compari-son purposes, we have purposely ignored the melt convection first. The effects of forced convection on the growth will be illustrated shortly. The governing equations for dimen-sionless temperature (θ ≡ Cs(T − T∞)/ H) and the phase-field variable [40] are written
FIG. 14. (a) Schematic of solidification in a square crucible; (b) calculated isotherms and flow patterns (left) at t= 0.05 and the mesh (right); (c) comparison of the isotherms (left) after further mesh refinement (right) based on velocity errors; the isotherms based on the velocity refinement are represented by the dashed lines; St= 1 and δ = 5 × 10−3.
as, respectively,
∂y)/(∂ϕ/∂x)], and τ0is a time constant;ε is a parameter for the interface energy anisotropy andβ0the seed orientation. Since the four-fold symmetry is considered, ifβ0= 0◦, one can take advantage of symmetry, and only a quarter of the domain is needed for diffusive growth.
However, when seed orientation is arbitrary, a full domain is required. To examine effect grid anisotropy, we have used a full domain for computation. In addition, for comparison purposes, the phase-field variableφ now ranges from 1 for solid and −1 for liquid, and φ = 0 is at the interface. Due to the slow development of the dendrite, typically the time step is usually small (∼0.01 τ0). Therefore, the refined zone can be very thin, which saves a tremendous amount of cells, especially for the simulation of low supercooling, where an extremely large domain is necessary. A sample of calculation is shown in Fig. 15a for the mesh (the first quarter), the isotherms (the second quarter), and the phase-field variable (the third quarter) for a fully developed dendrite (the fourth quarter) at the dimensionless supercooling ( ) being 0.55; ≡ Cl(Tm− T∞)/ H, where T∞is the far-field temperature.
We have purposely magnified the mesh at the interface locally to show the mesh refinement.
In this case, eight levels of grids are used and the cell number is 129036 for a full domain in Fig. 15a (only the mesh of a quarter domain is shown). In other words, the ratio of the largest to the smallest cell sizes is 28. Therefore, once the refinement can be carried out effectively, the mesh can cover different length scales for simulation. If an explicit integration scheme is used, the largest time step for integration is limited by the smallest cell size (the Courant–Friedrichs–Lewy (CFL) condition). Although a fully implicit scheme is used for the integration, due to the thin refined zone, to avoid numerical instability, the advancement of the new interface needs to be inside the refined zone. This is particularly true when side branching becomes important [41].
The comparison of our calculated dendrite tip speed for the supercooling of 0.55 with the solvability theory [40] is shown in Fig. 15b. Again, as shown the solvability limit can be reached as well and the agreement is very good. To examine the effect of growth orientation, we have used a full domain in this case. This also allows us to further examine the effect of grid anisotropy. Different seed orientations are considered. Forβ0= 0◦, the tip selection is in the x- or y-direction, and the agreement with the theory is very good. In this case, the calculation of the tip speed can be more accurate due to the search of zero contour of the phase-field variable at each time step. Although the calculated dendrite shapes with different orientations agree quite well with that ofβ0= 0◦, the calculation of the tip speed is more tedious and less accurate. Nevertheless, the error of the tip speed is still within 3%; the dendrite tip grows slightly faster forβ0= 30◦and 45◦. Indeed, the grid anisotropy slightly affects the tip speed and its effect can be further reduced if the grid is further refined. Furthermore, the calculation is tested for several supercoolings, and the agreement with reported values is very good. Furthermore, the performance is also comparable to that
FIG. 15. (a) Calculated dendrite shape (the fourth quarter), mesh (the first quarter), isotherms (the second quarter), and the phase-fields (the third quarter); a local lookup of the mesh is further illustrated by arrows indicating the multilevel and multiscale nature of the refinement; the spacing= 0.055 for isotherms and 0.2 for phase-field variable; (b) evolution of tip speeds (for three different growth orientations); the solvability theory is shown by the dashed horizontal line. The dendrite shapes of different growth orientations are overlapped for comparison.
FIG. 16. (a) Calculated dendrite growth for phase-field variable (upper left), mesh (upper right), isotherms (lower lest), and flow patterns (lower right) at different times; the contour spacing= 0.055 for isotherms;
(b) comparison of tip speed evolution at three locations with the results by Beckermann et al. [15].
described by Provatas et al. [25] using an adaptive finite element method that the CPU time is proportional to the domain size. However, as compared with the their refinement scheme using triangular cells, the present approach seems to be much easier. Detailed comparison with the previous study and theory can be found elsewhere [41]. Furthermore, the present
approximation is fully conserved due to the nature of the FVM approximation, so that the extension to the solutal calculation, which requires a highly conserved scheme, may be easier. Furthermore, the converge speed of the calculation at each time step is also faster than that by enthalpy model with a very smallδ. The interface thickness can be judged from the contours (10 contours) of the phase-field variable in Fig. 15a. In addition, according to Karma and Rappel [40], if the parameters are carefully chosen, such a thickness approaches the sharp-interface limit.
The final example is carried out for the dendritic growth under a forced convection for = 0.55. D = 4, and ε = 0.05. This example is the same as the one used by Beckermann et al. [15]; they used a uniform mesh for calculation, and they showed by simulation for the first time the effect of flow on the dendrite tip speed. The calculated dendrite shapes (phase fields) and meshes at different times are shown in the top of Fig. 16a; the corresponding isotherms and flow fields are shown in the lower figures. As shown, the tip at the upstream grows faster due to the thinner thermal boundary layer. The evolution of the tip velocities is further illustrated in Fig. 16b; the results of Beckermann et al. [15] are also put together for comparison. As shown, they are in good agreement. However, the cell numbers used here (2300, 5678, 7430 at t= 15, 70, 100τ0, respectively) are one order of magnitude smaller than that used by Beckermann et al. (512× 256), while our smallest cell size is also smaller than theirs. Due to the fully implicit scheme used here for integration, a much bigger time step being about 0.2τ0 can be used (0.01τ0 for an explicit scheme on the phase-field equation), and this reduces overall CPU time dramatically. To save space, further discussion on the dendritic growth, especially at low supercooling ( = 0.25 and 0.1) under a convective environment will be discussed elsewhere [41]. Furthermore, the grid at the root level has some effects on the result. Too coarse a grid will affect the overall growth environment.
5. CONCLUSIONS
An adaptive finite volume method using dynamic data structures, implemented in FOR-TRAN90, is developed for solidification problems. The method uses a fixed-grid and two-phase equations to treat free or moving interfaces as well as the associated heat transfer and fluid flow through a phase-field variable, which is defined by temperature fields and inter-face thickness. In addition to the refinement on the interinter-faces or boundaries, error estimators are also used for the refinement and the coarsening.
For the simulation of two-phase problems having a sharp interface, the accuracy of the fixed-grid approach using a structure-grid is usually limited by the cell size. However, through the adaptive mesh, the problem could be resolved while using much less cells for calculation. The strategy proposed here starts from a coarse grid having a thick interface.
As the problem converges at coarse levels, mesh refinement or coarsening is carried out according to geometric constraints or error estimators, followed by the interface thinning.
Through the extensive benchmarkings with the front tracking approach, we have found that the shape-interface problem can be well modeled. In addition, this approach is particularly useful for a domain having complex immersed interfaces (or with interface merging or splitting) that can be difficult to tackle by front tracking.
The present approach can be further extended to phase-field simulation. We have illus-trated successfully the simulation of dendritic growth using the phase-field equation. For the
case without convection, the calculated dendrite tip speeds at various undercoolings agree very well with the solvability theory. With a forced convection, the calculated result also agrees well with previous calculation, but uses much less cells. Therefore, we believe that further extensions to other multiphase systems, such as bubbling and sintering, may be fea-sible. One may also use level sets or other definitions of the phase-field variable to describe the moving interfaces. Due to the simplicity of the refinement scheme and the related data structures, the extension to three-dimensional problems is straightforward. However, one needs to keep this in mind that the computational cost for 2D problems is proportional to the arclength of the interface (for a thin refined zone). However, for 3D problems, the cost is proportional to the overall interfacial area. Therefore, a dramatic increase of CPU time is expected. Furthermore, although the approach is also suitable for a highly diffusive inter-face, as the interface thickness increases the computational cost increases as well. A careful decision needs to be made for the refinement inside the interface region for a cost-effective simulation.
Finally, due to the nature of the multilevel in the implementation, multigrid methods [34, 36] fit into the present scheme perfectly. As one may anticipate, the solution phase can be significantly enhanced. Although it is not crucial in this study due to the much smaller problem size, it may be highly desirable for 3D simulation.
ACKNOWLEDGMENTS
The authors are grateful for the support from the National Science Council and the National Center for High Performance Computing of the Republic of China under Grant NSC89-2214-E002-040. Fruitful discussion with Prof. G. Amberg in the phase-field simulation of dendritic growth is also appreciated. Valuable comments from the reviewer are acknowledged.
REFERENCES
1. C. J. Chang and R. A. Brown, Natural convection in steady solidification: Finite element analysis of a two-phase Rayleigh–Benard problem, J. Comput. Phys. 53, 1 (1984).
2. C. J. Kim and M. Kaviany, Numerical method for phase change problems with convection and diffusion, Int. J. Heat Mass Trans. 35, 457 (1992).
3. M. Lacroix, Predictions of natural-convection-dominated phase-change problems by the vorticity–velocity formulation of the Navier–Stokes equations, Numer. Heat Trans. 22, 79 (1992).
4. R. Viswanath and Y. Jaluria, Computation of different solution methodologies for melting and solidification problems in enclosures, Numer. Heat Trans. B 24, 77 (1993).
5. K. Tsiveriotis and R. A. Brown, Boundary-conforming mapping applied to computations of highly deformed solidification interfaces, Int. J. Numer. Meth. Fluids 14, 981 (1992).
6. C. W. Lan, Newton’s method for solving heat transfer, fluid flow and interface shapes in a floating molten zone, Int. J. Numer. Meth. Fluids 19, 41 (1994).
7. M. C. Liang and C. W. Lan, A finite-volume/Newton method for a two-phase heat flow problem using primitive variables and collocated grids, J. Comput. Phys. 127, 330 (1996).
8. V. R. Voller and C. Prakash, Fixed grid numerical modeling methodology for convection–diffusion mushy region phase-change problems, Int. J. Heat Mass Trans. 30, 1709 (1987).
9. W. D. Bennon and F. P. Incropera, Continuum model for momentum, heat and species transport in binary solid–liquid phase change systems, Int. J. Heat Mass Trans. 10, 2161 (1987).
10. M. Lacroix and V. R. Voller, Finite difference solutions of solidification phase change problems—Transformed versus fixed grids, Numer. Heat Trans. B 17, 25 (1990).
11. V. R. Voller and C. R. Swaminathan, General source-based method for solidification phase change, Numer.
Heat Trans. B 19, 175 (1991).
12. M. Yao and Henry de Groh III, Three-dimensional finite element method simulation of Bridgman crystal growth and comparison with experiments, Numer. Heat Trans. A 24, 393 (1993).
13. W. Shyy and M. Chen, Steady-state natural convection with phase change, Int. J. Heat Mass Trans. 33, 11, 2545 (1990).
14. M. Fabbri and V. R. Voller, Numerical solution of plane-front solidification with kinetic undercooling, Numer.
Heat Trans. B 27, 467 (1995).
15. C. Beckermann, H.-J. Diepers, I. Steinbach, A. Karma, and X. Tong, Modeling melt convection in phase-field simulations of solidification, J. Comput. Phys. 154, 468 (1999).
16. J. A. Sethian and J. Strain, Crystal growth and dendritic solidification, J. Comput. Phys. 98, 231 (1992).
17. S. Chen, B. Merriman, S. Osher, and P. Smereka, A simple level set method for solving Stefan problems, J. Comput. Phys. 135, 8 (1997).
18. M. Sussman, A. S. Almgren, J. B. Bell, P. Colella, L. H. Howell, and M. L. Welcome, An adaptive level set approach for incompressible two-phase flows, J. Comput. Phys. 148, 81 (1999).
19. D. A. Drew, Mathematical modeling of two-phase flow, Ann. Rev. Fluid Mech. 15, 261 (1983).
20. H. S. Udaykumar, H. C. Kan, W. Shyy, and R. T.-S. Tay, Multiphase dynamics in arbitrary geometries on fixed Cartesian grids, J. Comput. Phys. 137, 366 (1997).
21. H. S. Udaykumar, R. Mittal, and W. Shyy, Computation of solid–liquid phase fronts in the sharp interface limit on fixed grids, J. Comput. Phys. 153, 535 (1999).
22. T. Ye, R. Mittal, H. S. Udaykumar, and W. Shyy, An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries, J. Comptut. Phys. 156, 209 (1999).
23. C. S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25, 220 (1977).
24. D. Juric and G. Tryggvason, A front-tracking method for dendritic solidification, J. Comput. Phys. 123, 127 (1996).
25. N. Provatas, N. Goldenfeld, and J. Dantzig, Adaptive mesh refinement computation of solidification microstructures using dynamic data structures, J. Comput. Phys. 148, 265 (1999).
26. R. Tonhardt and G. Amberg, Dendritic growth of randomly oriented nuclei in a shear flow, J. Crystal Growth 213, 161 (2000).
27. I. Demirdzic and S. Muzaferija, Numerical method for coupled fluid flow, heat transfer and stress analysis using unstructured moving meshes with cells of arbitrary topology. Comput. Meth. Appl. Mech. Eng. 125, 235 (1995).
28. S. Muzaferija and D. Gosman, Finite-volume CFD procedure and adaptive error control strategy for grids of arbitrary topology, J. Comput. Phys. 138, 766 (1997).
29. S. R. Mathur and J. Y. Murthy, A pressure-based method for unstructured meshes, Numer. Heat Trans. B 31, 195 (1996).
30. L. H. Howell and J. B. Bell, An adaptive mesh projection method for viscous incompressible flow, SIAM J. Sci. Comput. 18, 996 (1997).
31. M. J. Berger and J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, J. Comput.
Phys. 53, 484 (1984).
32. C. M. Rhie and W. L. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA J. 21, 1527 (1983).
33. S. V. Patankar, Numerical Heat Transfer and Fluid Flow (Hemisphere, Washington, DC, 1980).
34. C. W. Lan and M. C. Liang, Multigrid method for incompressible heat flow problems with an unknown interface, J. Comput. Phys. 152, 55 (1999).
35. H. L. Stone, Iterative solution of implicit approximations of multidimensional partial differential equations, SIAM J. Numer. Anal. 5, 530 (1968).
36. Z. J. Wang, A fast nested multi-grid viscous flow solver for adaptive Cartesian/Quad grids, AIAA-96-2091, June 1996.
37. I. Demirdzic, Z. Lilek, and M. Peric, Fluid flow and heat transfer test problems for non-orthogonal grids:
benchmark solutions, Int. J. Numer. Meth. Fluids 15, 329 (1992).
38. V. R. Voller, Similarity solution for the solidification of a multicomponent alloy, Int. J. Heat Mass Trans. 40, 2869 (1997).
39. N. Palle and J. A. Dantzig, Adaptive mesh refinement scheme for solidification problems, Metall. Mater.
Trans. A 27A, 707 (1996).
40. A. Karma and W. J. Rappel, Quantitative phase-field modeling of dendritic growth in two and three dimensions, Phys. Rev. E57, 4323 (1998).
41. C. W. Lan, C. M. Hsu, and C. C. Liu, Efficient adaptive phase filed simulation of dentritic growth in a forced flow at low supercooling, J. Crystal Growth, in press.