Numerical Simulation: Revealing the Secrets of "Temperature and Flow" in Crystal Growth

Numerical Simulation: 

Revealing the Secrets of "Temperature and Flow" in Crystal Growth​

The growth of large-sized crystals is a highly systematic engineering endeavor. Different growth methods—such as the Czochralski method, flux method, vertical gradient freeze method—all require precise temperature control and a stable melt environment over time scales ranging from several hours to hundreds of hours to allow crystals to grow slowly along the intended direction.

Although the growth process may appear "slow and controllable," furnaces operate at high temperatures of 1000–2000 K, and their enclosed interiors make direct measurement challenging. As a result, key information such as melt flow, heat distribution, radiative losses, and interface morphology often cannot be observed directly. Numerical simulation thus becomes an essential tool for understanding these "invisible processes," allowing us to uncover the secrets of "temperature and flow" in crystal growth.

I. Macroscopic Physical Mechanisms in Crystal Growth[1,2]

Inside a crystal growth furnace, multiple material regions coexist, including the crystal, melt, crucible, insulation materials, heaters, and gas passages. Heat transfer varies across these regions, yet they interact with each other, forming a complex multiphysics system.

1. Heat Conduction in Solids (Crystal, Crucible, Insulation, etc.)[3]

Heat in solids is primarily transferred through conduction, governed by the basic equation:

ρ C_p (∂T/∂t ) = ∇·(k ∇T) + Q

T (Temperature): Fundamental quantity describing the thermal state within the material.

ρ (Density): Mass per unit volume of the material.

C_p (Specific Heat Capacity): Indicates how quickly a material heats up when absorbing heat.

k (Thermal Conductivity): Indicates how well a material conducts heat.

• High k (graphite, metals)→ Rapid heat diffusion → Uniform temperature.
• Low k (ceramics, fiber wool)→ Good thermal insulation.

Q (Heat Source): Energy generated by resistive heating, induction heating, etc.

The thermal conductivity of solids influences the shape of the crystal’s solid–liquid interface (concave/convex), which in turn affects growth stability.

2. Heat Conduction and Convection in the Melt​

The melt not only conducts heat but also flows, so heat can be transferred both by conduction and by being carried along with the flow:

ρ C_p (∂T/∂t + u·∇T) = ∇·(k ∇T) + Q

u (Velocity Field): Describes the direction and speed of melt flow.

u·∇T: Represents the transport of the temperature field due to convective effects.

∇T (Temperature Gradient): The direction and magnitude of temperature change.

While melt conducts heat slowly, convection is strong, so temperature changes rely more on fluid flow. Melt flow originates from:

Buoyancy-Driven Convection: Temperature difference → Density difference    ρ(T) = ρ[1 − β(T − T₀)]

• β (Thermal Expansion Coefficient): Proportional change in volume upon heating.

Marangoni Effect: Surface tension varies with temperature γ(T).

• Higher temperature → Lower surface tension.
• Directional flow is induced in the wetting layer.

Crystal/Crucible Rotation: Forced convection.

This means that melt behavior is almost entirely determined by the temperature field, while the melt behavior, in turn, affects the temperature field—a typical strongly coupled system. Critical issues related to crystal quality, such as compositional segregation and interface irregularities, are closely linked to the flow inside the melt.

3. Flow of Melt and Gas (Navier–Stokes Equations)[4]

The flow inside the melt follows the classic fluid dynamics equations:

ρ(∂u/∂t + u·∇u) = −∇p + ∇·μ(∇u + (∇u)^T) +ρ₀ g [1-β(T-T₀)] + F

Where:

−∇p (Pressure Term): Force driving fluid flow.

∇·μ(∇u + (∇u)ᵀ) (Viscous Term): Describes "internal friction" within the fluid—how "viscous" it is. High viscosity hinders flow.

ρ₀g[1-β(T-T₀)] (Buoyancy Term): Buoyancy under the Boussinesq approximation, related to gravity.

F (External Force Term): Represents forces such as Marangoni force, Lorentz force, and shear forces from crystal and crucible rotation.

In crystal growth, the most important driving force typically comes from: temperature difference → density difference → buoyancy-driven convection.

4. Radiative Heat Transfer, Non-Negligible at High Temperatures[5]

When temperatures exceed 1000 K, radiation becomes the dominant mode of heat transfer. Surface-to-environment radiation is often simplified as:

q_rad = ε σ (T⁴ − T_env⁴)

ε (Emissivity): The surface's ability to emit infrared radiation.

σ: Stefan–Boltzmann constant.

T_env: Ambient temperature.

This implies: the higher the temperature, the stronger the radiation—and the relationship follows T⁴, meaning a slight temperature increase can lead to a significant rise in radiative heat loss.

This equation applies only to surface-to-environment radiation. For mutual radiation among multiple surfaces inside a crystal growth furnace, a surface-to-surface radiation model is required, involving view factors or equivalent radiation equations. If the medium participates in internal radiation (e.g., semi-transparent crystals like LBO, BBO, YAG, YVO₄), more complex radiation transfer equations (RTE) are needed, though engineering simulations often employ approximate models like P1 or DOM.

5. Electromagnetic Field Equations for Induction Heating[6]

Many high-temperature crystal growth furnaces use induction heating, where the crucible or metal components are heated by eddy currents (Joule heating) generated in an alternating magnetic field. The electromagnetic field is calculated using Maxwell’s equations:

​∇ × ( (1/μ) ∇ × A ) + σ ∂A/∂t = J_coil

A (Magnetic Vector Potential): Its curl represents the magnetic field B=∇×A.

μ (Magnetic Permeability): Determines magnetic field penetration.

σ (Electrical Conductivity): Reflects the material’s ability to generate eddy currents.

J_coil :Coil Current Density.

Eddy current dissipation (Joule heating) → heats the crucible and melt. The eddy current heat source can be written as: Q_eddy = σ |∂A/∂t|². This heat source is directly incorporated into the heat transfer equation, coupling the thermal and electromagnetic fields. Moreover, eddy current distribution strongly depends on frequency—higher frequencies lead to stronger skin effects, concentrating heat near the material surface.

6. Boundary Conditions and Initial Values​

The model must specify boundary conditions for each physical domain:

Temperature boundaries  account for heat input (boundary heat sources) and heat loss (convection, radiation, or fixed heat flux).

Fluid boundaries​ at solid walls adopt a "no-slip" condition; free surfaces must consider surface tension effects and convective heat transfer.

Electromagnetic boundaries​ assign excitations at coil current or magnetic potential inputs, with outer boundaries typically treated as far-field or magnetically insulated.

Initial values: Quasi-steady-state calculations can use empirical values or results from a previous steady state as initial guesses; strict transient simulations should start from a reasonable initial temperature field and a stationary/existing flow field to ensure physical reliability in time evolution.

II. How Are These Equations Solved?​

These equations are nonlinear, strongly coupled, and involve multiple physical domains, making direct analytical solutions impossible. Numerical simulation software (e.g., COMSOL, Fluent) or dedicated crystal growth simulation tools (e.g., FEMAG, CGSim) discretize the equations by dividing the geometry into small grids. Common numerical methods include Finite Element Method (FEM), Finite Volume Method (FVM), and Finite Difference Method (FDM). The software iteratively solves for temperature, velocity, and electromagnetic fields at each grid point.

The solving process can be simplified into three steps:

1. The geometry is divided into many small grids (physical quantities are assumed approximately uniform within each grid), as shown in Figure 1.
2. Each partial differential equation is discretized into algebraic equations on the grid.
3. The solver iterates continuously until the temperature field, velocity field, and electromagnetic field gradually converge to a stable solution.

Figure 1. Finite Element Mesh Generation for a 2D Axisymmetric Model​

III. What Key Results Can Simulations Reveal?​

By solving the equations described above, we can obtain information about the furnace interior that is entirely unobservable through direct measurement. Key results are illustrated below.

Figure 2. Temperature Distribution

 in the Crystal-Melt Region and Interface Morphology (FEMAG Results)​

Figure 2 shows the temperature distribution within the crystal-melt region and the crucible, along with the morphology of the solid-liquid interface. The isotherms reveal the temperature difference between the crystal and melt sides. The black curve between the crystal and melt represents the solid-liquid interface. Where the isotherms are closely spaced, it indicates a steep thermal gradient, which can induce thermal stress. The convex or concave shape of the interface directly indicates local heat accumulation or dissipation.

Figure 3. Flow Streamlines and Vortex Structures Inside the Melt (FEMAG Results)​

Figure 3 illustrates the flow streamlines and vortex structures within the melt. The color indicates the local velocity magnitude (redder color indicates higher velocity), and the arrows show the flow direction. In a low-viscosity melt, a symmetric primary vortex is observed, featuring a descending flow along the centerline, followed by ascending flow along the sides, and a surface return flow towards the center. This recirculation is driven by buoyancy due to temperature differences and the Marangoni effect (surface tension gradients). If the surface return flow is significant, it can cause local depression at the crystal interface.

Figure 4. Electromagnetic Field Distribution Generated by the Induction Coil [8]

Figure 4 shows the magnetic flux distribution generated by the induction coil during operation. The redder the color, the higher the magnetic flux density. The higher magnetic flux density is concentrated between the coil turns and near the crucible. The alternating magnetic field induces eddy currents within the conductive crucible. These eddy currents generate heat via the Joule effect, with their distribution showing a characteristic "skin effect": higher frequencies confine the current and heating to a shallower layer near the surface. In engineering design, the magnetic field distribution directly guides the selection of coil turns, spacing, and operating frequency to ensure uniform heating and avoid localized hot spots or undesirable flow disturbances in the melt.

Figure 5. Gas Flow Circulation Within the Furnace​

Figure 5 illustrates the distribution of gas flow streamlines in a local region of the furnace. The density of the streamlines indicates the magnitude of the local gas velocity—higher density corresponds to greater flow speed—while the arrows denote the flow direction. Multiple convective cells are observed in the gas flow region, predominantly located above the heater and near the insulation walls, suggesting significant flow velocity and enhanced cooling effects in these areas. The overall and local gas flow patterns within the crystal growth furnace determine the convective heat transfer distribution throughout the chamber, thereby influencing the melt surface temperature and the local environment at the solid-liquid interface.

IV. Summary​

Crystal growth is a complex system driven by the interplay of heat conduction, convection, radiation, and electromagnetic heating, with different physical laws governing each region. Numerical simulation provides a window into understanding these physical processes. It allows us not only to predict crystal quality but also to optimize furnace design, improve process stability, reduce trial-and-error costs, and enhance yield. As computational power advances, these "secrets" of crystal growth are becoming increasingly clear, paving the way for more controllable and designable crystal growth processes.

Reference:

1. Min, Nai-ben. Fundamentals of the Physics of Crystal Growth. Shanghai Scientific and Technical Publishers, 1980.
2. Chernov, A. A.; Wu, Zi-qin; Hong, Yong-yan; Gao, Chen. Modern Crystallography, Vol. 3: Crystal Growth.
3. Dupret, F.; Ryckmans, Y.; Wouters, P.; Crochet, M. J. Numerical Calculation of the Global Heat Transfer in a Czochralski Furnace. Journal of Crystal Growth 1986, 79 (1), 84–91. https://doi.org/10.1016/0022-0248(86)90419-7.
4. Yeckel, A.; Derby, J. J. Computer Modelling of Bulk Crystal Growth. In Bulk Crystal Growth of Electronic, Optical & Optoelectronic Materials; Capper, P., Ed.; John Wiley & Sons, Ltd: Chichester, UK, 2010; pp 73–119. https://doi.org/10.1002/9780470012086.ch3.
5. Howell, J. R.; Mengüç, M. P.; Daun, K. J.; Siegel, R. Thermal Radiation Heat Transfer, Seventh edition.; CRC Press / Taylor & Francis Group: Boca Raton, 2021.
6. [6]Richard E. Haimbaugh. Theory of Heating by Induction, Practical Induction Heat Treating, 2nd ed., ASM International, 2015, p 9-23, https://doi.org/10.31399/asm.tb.piht2.t55050009.
7. Enders-Seidlitz, A.; Pal, J.; Dadzis, K. Development and Validation of a Thermal Simulation for the Czochralski Crystal Growth Process Using Model Experiments. Journal of Crystal Growth 2022, 593, 126750. https://doi.org/10.1016/j.jcrysgro.2022.126750.

[8]Tavakoli, M. H.; Wilke, H. Numerical Study of Heat Transport and Fluid Flow of Melt and Gas during the Seeding Process of Sapphire Czochralski Crystal Growth. Crystal Growth & Design 2007, 7 (4), 644–651. https://doi.org/10.1021/cg060383y.