Influence of static magnetic fields and solutal buoyancy on silicon dissolution into germanium melt

(1)

Influence of static magnetic fields and solutal

buoyancy on silicon dissolution into

germanium melt

by

Anton Kidess

BSc., Friedrich-Alexander Universität Erlangen-Nürnberg, 2006

A Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of

MASTER OF APPLIED SCIENCES in the Department of Mechanical Engineering

©Anton Kidess, 2009.

University of Victoria, Department of Mechanical Engineering

(2)

Supervisory committee

Influence of static magnetic fields and solutal buoyancy on silicon dissolution into germanium melt

by

Anton Kidess

BSc., Friedrich-Alexander Universität Erlangen-Nürnberg, 2006

Dr. Sadik Dost, Supervisor

(Department of Mechanical Engineering) Dr. Henning Struchtrup, Member

(Department of Mechanical Engineering) Dr. Dante Canil, Member

(3)

Abstract

Supervisory Comittee

Dr. Sadik Dost, Supervisor

(Department of Mechanical Engineering) Dr. Henning Struchtrup, Member

(Department of Mechanical Engineering) Dr. Dante Canil, Member

(Department of Earth and Ocean Sciences)

Abstract

Elemental semiconductors like silicon and germanium have been used since the beginning of the electronics industry. Silicon has dominated research and production and thus silicon based devices can be produced at the lowest cost using the most mature technology. While dopants can be used to tailor the electric properties of the semiconductor within certain lim-its, more flexibility is gained using compound semiconductors such as silicon-germanium. The electric properties of a compound semiconductor are highly dependant on the com-position, which in turn is influenced by the dissolution reaction and flow characteristics during the growth process.

Liquid phase diffusion (LPD) is a solution growth technique that has been proposed to grow silicon-germanium seed crystals for other growth techniques. The dissolution of silicon is a limiting factor for the growth rate in LPD and also Bridgman growth tech-niques. Investigation of the dissolution process is aimed at increasing the growth rate while still maintaining maximum uniformity of the crystal composition. To accomplish this, a static magnetic field was utilized in experiments done by Armour [AD07]. The experimental results showed that a top seeded configuration without magnetic fields leads to a diffusion driven process and homogeneous dissolution, while the addition of a strong

(4)

0.8 Tesla magnetic field resulted in non-uniform and slightly increased dissolution. This work is complementary to the experimental investigation and aims to help understand the influence of magnetic fields on silicon dissolution.

For this work, an OpenFOAM magnetohydrodynamics application including heat and species transport and three different magnetic force models has been developed and vali-dated. The simulations done show that an isothermal state is reached within 90 seconds if no temperature gradient is imposed. Additional simulations with a temperature gradient helped to rule out a possible thermal leak in the experimental system, confirming that it must have been close to isothermal. Since the solutal expansion coefficient of SixGe1−x

has not been measured properly to the Author’s knowledge, two possible values for the expansion coefficient have been considered. It has been found that the exact value of the solutal expansion coefficient does not have a great influence on the results of this work.

(5)

List of Tables

1.1 Macro-scale models . . . 4

1.2 Comparison of CFD codes . . . 6

2.1 Coefficients for the generic transport equation . . . 9

2.2 Thermophysical properties of the Si-Ge system . . . 17

(9)

List of Figures

1.1 Dissolution experiment results after 20 minutes . . . 5

2.1 Dissolution experiment crucible (all units in mm) . . . 8

2.2 Simulation domain . . . 8

3.1 3D control volume . . . 22

3.2 2D view of a control volume . . . 23

3.3 Blocking strategies . . . 29

3.4 Skewed grid . . . 29

3.5 Final grid shape . . . 30

3.6 High resolution grid used for OpenFOAM simulations . . . 31

4.1 Solutal buoyancy benchmark results at non-dimensional times 0.75, 1.35 and 3.15 (concentration contour left, temperature contour + streamlines right) . 34 4.2 Horizontal Bridgman benchmark with no additional equation for the mag-netic force, OpenFOAM . . . 35

4.3 Horizontal Bridgman benchmark with no additional equation for the mag-netic force, FASTEST . . . 36

4.4 Horizontal Bridgman benchmark with a variable magnetic field, OpenFOAM 36 4.5 Horizontal Bridgman benchmark with induced electric field, OpenFOAM . 37 4.6 Cylinder benchmark without additional equations for the magnetic force, OpenFOAM . . . 38

4.7 Cylinder benchmark with a variable magnetic field, OpenFOAM . . . 39

4.8 Cylinder benchmark with induced electric field, OpenFOAM . . . 39

4.9 Cylinder benchmark with induced electric field, FASTEST-3D results for two different time marching schemes . . . 40

(10)

LIST OF FIGURES

5.2 Concentration isolines for an isothermal simulation without magnetic field; obtained with the Crank-Nicholson time marching scheme (SixGe1−x) . . . 43

5.3 Concentration isolines after 20 minutes of dissolution, reduced expansion coefficient (SixGe1−x) . . . 44

5.4 Velocity vectors in a plane cut normal to x, reduced expansion coefficient . 45 5.5 Velocity in a plane 1mm from the dissolution interface, B=0.8T . . . 45 5.6 Silicon concentration at the crucible centerline and 1cm offset with and

without magnetic field, reduced expansion coefficient . . . 46 5.7 Concentration isolines after 20 minutes of dissolution, high expansion

coef-ficient (SixGe1−x) . . . 47

5.8 Velocity vectors in a plane cut normal to x . . . 48 5.9 Velocity in a plane 1mm from the dissolution interface, B=0.8T . . . 48 5.10 Relative difference of concentration between reduced and high solutal

ex-pansion coefficient results . . . 49 5.11 Relative difference of concentration between reduced and high solutal

ex-pansion coefficient results with 0.8 Tesla magnetic field . . . 49 5.12 Temperature profile . . . 51 5.13 Concentration isolines for reduced solutal expansion coefficient with a

tem-perature gradient (SixGe1−x) . . . 51

5.14 Velocity vectors for reduced solutal expansion coefficient with a temperature gradient . . . 52 5.15 Concentration at centerline and offset locations with βC = 0.17 and a

tem-perature gradient . . . 53 5.16 Concentration isolines for high expansion coefficient with temperature

gra-dient (SixGe1−x) . . . 54

5.17 Velocity vectors for reduced expansion coefficient with temperature gradient 54 5.18 Concentration at centerline and offset locations with βC = 0.53 and a

(11)

Nomenclature

B Magnetic field

βC Solutal expansion coefficient

βT Thermal expansion coefficient

C Concentration of silicon cp Heat capacity D Diffusion coefficient J Electric current E Electric field η Magnetic diffusivity g Gravity

λ Thermal conduction coefficient FL Lorentz force

µ Dynamic viscosity µm Magnetic permeability

p Pressure

ΦE Electric scalar potential

ρ Density

(12)

Nomenclature σ Electric conductivity

T Temperature

(13)

Acknowledgements

I would like to thank Dr. Sadik Dost for the opportunity and financial support to pursue this research. I am grateful to Professor Franz Durst for supporting my wish to study abroad. Many thanks also for the continued support of Neil Armour and Alexander Raufeisen, both always finding time to answer my questions. Finally, I would like to thank Nicholas Audet and Blagovest Levitcharsky from 5NPlus Inc for the opportunity and support during a three month graduate internship in Montreal.

(14)

1 Introduction

1.1 Silicon-germanium

Germanium was the material used to produce the first transistor in 1947. Just a few years earlier, germanium diodes started replacing vacuum tubes in electronic devices. Ger-manium continued to be the leading material in the emerging electronics industries until the 1960’s, when high purity silicon rapidly became the dominant semiconductor material [Sei97]. With decreasing commercial interest in germanium, much of the research and devel-opment was directed towards silicon. Fuelled by high frequency applications that emerged in the 1990’s, especially personal communication services such as GSM and CDMA, the need for “faster” semiconductors arose. A promising candidate which has received lots of attention is silicon-germanium (Si-Ge), since it combines high mobility germanium to-gether with mature silicon manufacturing technology (Si-Ge is compatible with existing technology). Si-Ge can also be used as a substrate for Gallium-Arsenide, which can be used to make even faster chips. Other applications include photo-detection, photovoltaics and thermal imaging.

Crystal growth techniques from the liquid phase can be categorized in two groups: melt growth and solution growth. With melt growth, the semiconductor material is first heated past its melting point, then cooled to form a crystal. Often when trying to grow single crystals a seed is required. The crystal seed serves as a nucleation site to aid solidification and is a small single crystal usually grown by other means. Melt growth techniques are used in most commercial applications since they produce high quality material in large quantities. From the melt growth techniques, the Czochralski crystal pulling method is by far the most popular for commercial applications.

In most of the applications mentioned above, Si-Ge single crystals with a specific com-position are needed. Growing Si-Ge with melt growth techniques is problematic due to the large miscibility gap and the very different physical properties of silicon and germanium. The high temperature gradient applied is another issue, since it leads to large stress in the growing crystal and to varying radial compositions since the composition is highly

(15)

temper-CHAPTER 1. INTRODUCTION

ature dependent. Silicon is highly reactive at elevated temperatures, which promotes the inclusion of impurities into the crystal.

Solution growth techniques work at a temperature that is slightly higher than the lowest melting temperature of the constituents. The component with the lowest melting tem-perature is used as the solvent in which the source material slowly dissolves before the mixture re-solidifies on a substrate. Since the temperature is significantly lower than in melt growths techniques, the growing crystal will be subjected to less thermal stress. How-ever, since the growth rate is limited by the dissolution speed, solution growth is usually much slower than melt growth.

A solution growth process is liquid phase diffusion (LPD), which has been proposed as a method to grow seeds for growth techniques like the travelling heater method (THM) and liquid phase electro-epitaxy (LPEE) . To grow silicon-germanium crystals using the LPD technique, a solid silicon source is placed on top of a germanium melt and a solid germanium substrate. A temperature gradient is applied such that the silicon source starts dissolving into the germanium melt and resolidifies as silicon-germanium on the lattice matched germanium substrate.

In many practical situations the buoyancy induced motion in the melt is not desired and efforts have been made to control the flow. One possibility to actively control growth parameters is to apply rotating or static magnetic fields. Static fields have been used since the 1960’s to reduce turbulence in Czochralski melts. Rotating fields can impose steady flows that dominate over the natural convective flow, while static fields are used to dampen flows. Weak magnetic fields are used to compensate fluctuations in the flow, and strong magnetic fields are used to suppress convection entirely. More information on the use of magnetic fields in crystal growth can be found in the review paper by Series and Hurle [SH91].

1.2 Modelling of crystal growth

Heat and mass transfer in the growth system have a great effect on the quality of the grown crystal. Composition uniformity, impurity concentration and crystal structure are all affected by the flow of the solution during the growth process. In order to grow high quality crystals a good understanding of the flow in the melt is highly beneficiary. Accompanying experiments with computer simulations is necessary not only because of the high cost associated with experiments, but also because the experimental study of transient flows

(16)

CHAPTER 1. INTRODUCTION

in the melt is difficult due to very high temperatures, closed systems and opaque liquids. Ideally non-contact measurement techniques are used, which include visualization using X-rays and ultrasound. Both require tracer particles that have to match the density of the solution very well to avoid sedimentation or floating. Additionally, experimental side effects such as quenching are not well resolved.

Crystal growth phenomena occur on various scales and the numerical model has to be selected to appropriately resolve the scale. Atomistic scale phenomena include nucleation, early stages of crystal growth and growth of nano crystals. Simulations on this scale are usually conducted using Monte Carlo or molecular dynamics type of simulations. On the macroscopic scale phenomena like convection, turbulence and heat transfer to participating media are investigated using 2D or 3D time dependent models. The models commonly used are finite difference, volume and elements, all of which are Eulerian continuum models. Another category are Lagrangian techniques such as smoothed particle hydrodynamics (SPH). An overview of key features is given in table 1.1.

SPH is a novel technique developed for astrological simulations but also well suited for crystal growth applications. Instead of using a mesh, SPH introduces a number of particles into the domain and determines the influence of a particle on every other particle within a “smoothing radius”. Since particle motions are tracked directly, interface tracking, e.g. in phase change applications, is simplified. Also, adaptive resolution reduces to increasing or decreasing the number of particles in a certain region. Typically the number of particles to achieve a certain accuracy is lower than the number of grid points needed with a Eulerian technique. However, computational neighbour search is a complex task and needs highly optimized algorithms to be feasible. Currently, the finite volume method is still more efficient for most applications [RD07].

1.3 Scope of this work

The interactions of the flow and the magnetic field are complicated, especially when con-sidering secondary effects such as induced electric fields or bending of isomagnetic lines. Numerical calculations are used to overcome the problems encountered in understanding experimental flows and optimize the process to enhance yield and quality of the grown crystals.

Experiments investigating the effect of a vertical static magnetic field on the silicon dissolution have been done by Neil Armour [AD07]. Somewhat counter-intuitive Armour

(17)

CHAPTER 1. INTRODUCTION Table 1.1: Macro-scale models

Model Characteristics

Finite difference easy to implement, but not very flexible Finite volume very flexible, conservative, but sensitive to grid

distortions

Finite element very flexible, but not conservative

Smoothed particle hydrodynamics grid-less, simplifies interface tracking e. g. for solidification simulations, adaptive resolution simple, performance has to be improved

finds a static magnetic field enhances the silicon dissolution rate close to the crucible walls. Figure 1.1 shows the obtained crystals after 20 minutes of top-seeded dissolution with and without a magnetic field, where the different interface shapes are clearly visible. The simulations carried out in this work are complementary to the experiments done by Armour and aim to solve the question about the exact influence of the magnetic field.

Previous investigations of the flow with the commercial finite volume based CFD package Ansys CFX have only been successful for the thermal field without magnetic forces [AD07]. While CFX is well suited for general purpose fluid flow problems, problems arise with strong magnetic fields and the numerical challenges involved solving Magnetohydrodynamic flows. Even though velocity magnitudes become very small when strong fields are applied, they are not easy to investigate numerically due to very strong body forces weakening diagonal dominance of the system and thin magnetic boundary layers that are difficult to resolve. The CFX simulations could only be made convergent with a scaled gravity and magnetic field reducing the body forces. However, the results obtained might be unphysical. The problems encountered with CFX will likely occur in all other commercial general purpose CFD packages (since the source is not available the numerics can’t be optimized for the problem at hand), so in this work simulations will be done with open source codes optimized for use in crystal growth with magnetic fields.

(18)

CHAPTER 1. INTRODUCTION

Figure 1.1: Dissolution experiment results after 20 minutes

(a) No magnetic field (b) Magnetic field of 0.8 Tesla

1.4 Multiphysics software

Numerous CFD codes are available under an open source license1 _{and vary widely in their}

capabilities, documentation and community involvement. Development of these codes is mostly driven by the academic institutions where they originated, but in some cases have gained enough momentum to sustain a code contributing community.

For this work, a set of essential requirements and desirable features was created and used to evaluate available CFD codes. Essential requirements are working solvers for the momentum, concentration and temperature fields as well as the ability to solve conjugate heat transfer problems. Unlike smaller additions like solutal buoyancy, adding those fea-tures to a code is considered far beyond the scope of this work. Desirable feafea-tures are compatibility with a decent pre- and post-processing tool, ease-of-use, active development, decent documentation and an active community.

A comparison of three codes that satisfy the essential requirements and at least a subset of the desirable features is presented in table 1.2. OpenFOAM (FOAM = field operation and manipulation) and FASTEST 3D (Flow Analysis by Solving Transport Equations Simulating Turbulence 3-Dimensional) were utilized for this work because the author has more experience with the finite volume method and ELMER is based on the finite elements. FASTEST has been used extensively for simulations of the Czochralski process at the

(19)

CHAPTER 1. INTRODUCTION Table 1.2: Comparison of CFD codes

Feature OpenFOAM ELMER FASTEST-3D

License GPL GPL individually licensed to

academic institutions

NS, temp, concentration equations + + +

Used for crystal growth problems - + ++

Discretization FV FE FV

Documentation ++ +

-Community ++ +

-Extensibility, flexibility ++ +

-fluid mechanics institute of the university of Erlangen-Nuremberg (Lehrstuhl fuer Stroemu-ngsmechanik, Universitaet Erlangen-Nuernberg). It is a finite volume CFD code working with multiblock structured meshes and is available under an open source academic license. OpenFOAM is a general framework to solve partial differential equations developed at the Imperial College in London, England. Although most of the standard solvers are targeted towards fluid dynamics problems, it also includes some financial and structural solvers as well. OpenFOAM is written in C++ and makes great use of object-oriented design principles, which enables the development of elegant “applications” that directly resemble the underlying partial differential equations. It is very flexible thanks to a wide selection of linear solvers and interpolation techniques, and boasts a large user community and available commercial consulting.

In the early stages of this work, FASTEST was solely used to simulate the dissolution process due to its maturity with respect to crystal growth problems. However, throughout this work problems with the magnetic field implementation became apparent (see section 4) that could not be solved. Instead, an OpenFOAM solver was developed and utilized to simulate the flow in the melt with boundary conditions based on the non-magnetic results obtained with FASTEST and CFX.

(20)

2 Problem Formulation

2.1 Geometry

The main focus of this work was put on the top seeded configuration, where the silicon is placed on top of the germanium melt, because the bottom seeded configuration was shown to be a lot less stable in the experiments done by Armour [ADL07]. A drawing of the dissolution experiment setup is shown in figure 2.1. To simplify the simulation, only the source, solution and the left, right and bottom quartz walls are included in the simulation domain. The solution domain is shown in figure 2.2. The solid domains were not included in the simulations done with OpenFOAM. Instead, only the melt was considered with a high resolution, and the boundary conditions were adapted accordingly.

2.2 Governing equations

The mass, momentum, energy and species transport equations for a viscous, incompressible liquid are displayed in equations 2.1 to 2.4. In three dimensions there are six equations for the six variables u1, u2, u3, p, T and C, which are the velocity, pressure, temperature

and concentration, respectively. Depending on the model used for the Lorentz force FL,

additional variables are introduced. Se and R are energy sources / sinks and reaction

terms, respectively. Since no reactions are modelled within the domain and no heat is generated or absorbed they will not be considered from now on.

∂ui ∂xi = 0 (2.1) ∂ρui ∂t + ∂ρuiuj ∂xj = ∂ ∂xj µ∂ui ∂xj ! − ∂p ∂xi + ρg i+ FL,i (2.2) ∂ρcpT ∂t + ∂ρcpujT ∂xj = ∂ ∂xj λ∂T ∂xj ! + Se (2.3)

(21)

CHAPTER 2. PROBLEM FORMULATION

Figure 2.1: Dissolution experiment crucible (all units in mm)

Figure 2.2: Simulation domain

Quartz, 2mm thick 25 mm diameter

30 mm 10 mm Silicon source (solid)

Germanium melt (liquid) g

B

y

x z

(22)

CHAPTER 2. PROBLEM FORMULATION Table 2.1: Coefficients for the generic transport equation

Quantity Φ ΓΦ QΦ Mass 1 0 0 Momentum Uj µ −∂P ∂xj + ρgj+ FL,i Thermal energy cpT λ Se Species concentration C ρD R ∂ρC ∂t + ∂ρujC ∂xj = ∂ ∂xj ρD∂C ∂xj ! + R (2.4)

The six transport equations can be written as a generalized transport equation (equation 2.5). An overview of the coefficients can be found in table 2.1.

∂ρΦ ∂t + ∂ρujΦ ∂xj = ∂ ∂xj ΓΦ ∂Φ ∂xj ! + QΦ (2.5)

2.3 Boussinesq approximation

An assumption that significantly simplifies solving the momentum balance equation 2.2 is that the mixture density is constant in the solution domain. If buoyancy is a driving force this assumption will lead to wrong results because density gradients due to varying thermal and solutal conditions are not considered. A possibility to resolve this problem without introducing a variable density is to consider a non-constant density only in the body force term of equation 2.2. This is known as the Boussinesq approximation and yields good results if the change in density is small and does not affect the mixture variables significantly.

Expression 2.6 shows a Taylor expansion of the density as a function of temperature and concentration where higher order terms were ignored. By introducing the thermal and solutal expansion coefficients, we can rewrite the body force in the momentum equations as shown in equation 2.7. For constant gravitational acceleration, the term ρ0g can be

(23)

CHAPTER 2. PROBLEM FORMULATION ρ(T, C) = ρ0+ ∂ρ ∂T _T 0,C0 (T − T0) + ∂ρ ∂C _T 0,C0 (C − C0) (2.6) ρg= ρ0g+ g [βT(T − T0) + βC(C − C0)] (2.7)

Tabulated values of the solutal expansion coefficient have to be used with care. Unlike the thermal buoyancy coefficient, which is commonly given in 1/K but can be used for values in °C, the solutal expansion coefficient differs greatly for molar concentrations and concentrations in weight%. FASTEST expresses concentrations in wt%, while the devel-oped OpenFOAM application works with molar%. If the solutal expansion coefficient is not given with respect to the proper unit, a conversion will have to be done.

Another possible pitfall is that the solutal expansion coefficient is usually given in terms of weight or molar percentages, while computer codes typically use values between zero and one (this is certainly the case for FASTEST and OpenFOAM). The resulting body force differs by a factor of 100 if the conversion is not carried out correctly. While the thermal expansion coefficient is measured reasonably well at least for elemental silicon and germanium, no solutal expansion measurements have been done for liquid silicon-germanium solutions to the Author’s knowledge.

FASTEST already implements the Boussinesq approximation for temperature variations. However, the effect of concentration variations on the mixture density was not taken into account, so the program code had to be extended to include solutal buoyancy. This modi-fication was validated using a square cavity test case based on work done by Bergman and Hyun [BH96]. The same benchmark is used to validate the OpenFOAM implementation. The results are presented in chapter 4.

2.4 Magnetic field source term

2.4.1 MHD approximation

A set of assumptions which will greatly simplify the Maxwell equations is known as the MHD approximation (see Chandrasekhar [Cha61] and Hughes [HY66], pp 142-148).

1. |V |2 _c2_{. Clearly in this work all velocities will be well below the speed of light and}

(24)

2. The induced magnetic field is assumed to be much smaller than the applied external field. With this assumption the magnetic field will be the same in all reference frames. 3. B = µmH for conductors and moderate electric fields.

4. For conductors, the displacement current is much smaller than the conduction current J and hence can be neglected.

5. The electric field is small compared to the magnetic field.

6. In liquid conductors, the space charge and its time derivative can be neglected. Also the conductivity is assumed constant.

2.4.2 Equations

Using the MHD approximation described in the previous section the Maxwell equations describing the influence of fluid motions on the electromagnetic field can be expressed as shown in 2.8-2.11. Here J is the electric current, B is the magnetic field, E is the electric field and µm is the magnetic permeability.

∇ · J = 0 (2.8)

∇ · B = 0 (2.9)

∇ × E = −∂B

∂t (2.10)

∇ × B = µmJ (2.11)

The inverse effect is described by the Lorentz force given by equation 2.12, with the simplifications mentioned in section 2.4.1.

FL= J × B (2.12)

The simplest way to determine the magnetic field force is to use Ohm’s law (2.13)

J = σ (E + u × B) (2.13)

to replace the current density in the Lorentz equation (2.12) and neglect the induced electric field, leaving equation 2.14.

(25)

Writing the Lorentz force in this form introduces no additional variables and therefore no additional boundary conditions have to be specified. It is the magnetic field implementation that requires the least amount of computing power.

Two more sophisticated approaches will be described in the following two sections. All three models (no additional equation, variable magnetic field and induced electric field) are tested with two benchmark cases by Hadid and Henry [HH96, HH97].

Variable magnetic field

Chandrasekhar and Hughes use equation 2.11 to eliminate the current density and express the Lorentz force solely in terms of the magnetic flux (equation 2.15).

FL= 1 µm(∇ × B) × B = − 1 µm B ×(∇ × B) = 1 µm(B · ∇)B − 1 2µm ∇B2 (2.15) Using this form of the Lorentz force, an additional equation is introduced and solved for B. The third Maxwell equation describes the change of the magnetic field with time and is used as a starting point to derive an additional equation. Replacing the induced electric field in equation 2.10 with Ohms law 2.13 leads to equation 2.16.

∇ ×(J

σ − u × B) = − ∂B

∂t (2.16)

The fourth Maxwell equation 2.11 is used to replace the current density J, which gives the desired equation for the time varying magnetic field (equation 2.17).

∂B

∂t = ∇ × (u × B) − ∇ × ( ∇ × B

µmσ ) (2.17)

If the magnetic permeability is assumed to be constant, we can use the identity from equation 2.18 to simplify the ordinary differential equation and obtain equation 2.19.

∇ ×(∇ × B) = ∇(∇ · B) − ∇2B= −∇2B (2.18) ∂B ∂t = 1 µmσ ∇2B+ ∇ × (u × B) (2.19)

(26)

CHAPTER 2. PROBLEM FORMULATION form of the magnetic field equation 2.21:

∇ ×(u × B) = (B · ∇)U − B(∇ · u) + u(∇ · B) − (u · ∇)B = −(B · ∇)u − B(u · ∇) (2.20)

∂B ∂t = η∇

2_{B −}_{(B · ∇)u − B(u · ∇)} _(2.21)

The resulting equation is very similar to the momentum equations and can be solved by similar means. This similarity is most obvious in the magnetic diffusivity η = 1

µσ,

which has the same units as the diffusivity of momentum (kinematic viscosity). To solve this equation, and also ensure the magnetic field stays divergence free (2.9, analog to the continuity equation for incompressible fluids), an artificial magnetic field pressure is introduced. The pressure term will tend to zero for a converged solution and is only needed to solve the equation by the PISO algorithm common for the momentum equation. Additional boundary conditions for the magnetic field and the magnetic field pressure will have to be specified using this approach.

The solution approach outlined in this section is common in MHD flows of plasmas, but requires knowledge of the magnetic permeability material property. A reasonable assumption for conducting liquids is to use the permeability of vacuum ([HY66], 141f), however it can be avoided using a different approach to express the Lorentz force. Also, since four additional equations will be solved, this model will require the most computing power.

Induced electric field

Using Ohm’s law (equation 2.13) to replace the current density in the Lorentz equation 2.12 leads to expression 2.22.

FL= σ(E + u × B) × B (2.22)

The induced electric field presents an additional variable and thus an additional equation is needed to fully specify the problem. This equation, given in terms of the electric scalar potential φE, is obtained by taking the divergence of Ampère’s law (equation 2.11). The

(27)

−∂ρE

∂t = µm∇ · J = 0 (2.23)

Following the MHD approximation, the overall charge is constant and the left hand side of the equation is zero. Using the charge conservation equation in Ohm’s law (2.13) leads to a Laplace type equation 2.24.

∇2_φ

E = ∆φE = ∇ · (u × B) (2.24)

The relationship between the induced electric field and the electric scalar potential is described by equation 2.25.

E = −∇φE (2.25)

With this model, an additional boundary condition for the scalar potential is needed. Since only one additional equation is required, the additional computation cost is moderate.

2.5 Boundary and initial conditions

Boundary conditions for the velocity, temperature and concentration variables have to be set to fully specify and solve the transport equations. Since concentration and velocity are not computed outside the fluid region, conditions only have to be set on the boundaries of the solution zone.

2.5.1 Boundary conditions for the momentum equation

For the velocity, non-slip conditions are used at the crucible walls and the source-solution interface bounding the fluid region (equation 2.26).

u _wall= 0 (2.26) ∂P ∂n _wall= 0 (2.27)

A problem with the pressure solution is that in the closed crucible all domain bound-aries will be assigned a Neumann type, zero-gradient boundary condition (equation 2.27).

(28)

Fortunately, in such systems the absolute value of the pressure is not of interest, and a solution can be obtained by setting a fixed pressure in an arbitrary point.

2.5.2 Boundary conditions for the concentration equation

At the crucible walls a no mass-flux condition is applied, which is equivalent to applying a zero concentration gradient at those boundaries (equation 2.28). At the dissolution interface the saturation condition 2.29 is applied. The condition is derived from the silicon-germanium phase diagram [LRS01], with C denoting the concentration of silicon in molar% and T being the temperature in Kelvin.

∂C ∂n _wall= 0 (2.28) C _{interf ace}= 1 − q 1412 + 1600 395 − T − 40 √ 395 √ 395 (2.29)

2.5.3 Boundary conditions for the temperature equation

It is assumed that the heat flux and temperature are continuous across the inner boundaries of the solution domain (equations 2.30 and 2.31), and Dirichlet type boundary conditions are used at the outer walls of the quartz crucible (equation 2.32).

Tliquid= Tsolid (2.30) λ∂T ∂n _liquid= λ ∂T ∂n _solid (2.31) T _{crucible wall}= 1373 K (2.32) ∂T ∂n|top−wall= −5 000 K_/_m _{←→ A · λ}∂T ∂n|top−wall = −420 W (2.33)

The initial velocities are set to zero, as well as the initial silicon concentration in the melt. The domain is initialized to 1073 K. For FASTEST, the boundary condition at the crucible walls are set to 1373 K. For the simulations with OpenFOAM the temperature at the lower and side domain boundaries was linearly raised to 1373 K within 90 seconds.

(29)

For the top domain boundary a gradient condition as in equation 2.33 was used, with heat transfer rates of zero and 420 Watt.

2.5.4 Boundary conditions for the induced electric field equation

The quartz ampule at the border of the domain is an insulator, so the current density normal to the quartz wall will be equal zero (equation 2.34).

J · n= 0 (2.34)

Applying this condition to Ohm’s law (equation 2.13) leads to a boundary condition for the induced electric field. Knowing that the velocity will be zero at a no-slip, non-permeable boundary, the boundary condition reduces to a zero normal gradient of the electric scalar potential (equation 2.35).

∂φE

∂n = (u × B) · n = 0 (2.35)

2.5.5 Boundary conditions for the variable magnetic field equation

On surfaces normal to the magnetic field flux vector, a Dirichlet condition is set (equation 2.36). On all other surfaces, a von Neumann condition is chosen (equation 2.37).

B|axial = B0 (2.36)

∂B

∂n|sidewall= 0 (2.37)

As for the physical pressure, von Neumann conditions for the artificial magnetic pressure are set on all walls. To solve the variable magnetic field equations with these boundary conditions, the magnetic field pressure is set to a constant value at an arbitrary location in the domain.

2.6 Physical parameters

The physical properties of the silicon source, the germanium solution and the quartz cru-cible are listed in table 2.2. The solution is treated as a binary mixture with a small amount of solute, so the mixture properties will be close to those of liquid germanium. Much of

(30)

the available material property data has been collected in the 1960’s and early 70’s, before interest in germanium faded and silicon became the dominant material used in semicon-ductor technology. Not all of the available data is reliable [LRS01], since measurements are problematic due to difficulties associated with the high melting temperatures of silicon and germanium. In particular, the solutal expansion and the binary diffusion coefficients are not well measured. If temperature dependent data is available, values at a temperature of 1373 K are used; otherwise data collected at the melting point is used.

Table 2.2: Thermophysical properties of the Si-Ge system

Property Material Symbol Value Unit Reference

Density Solution ρ 5490 kg_/_m3 [NH92]

Source 2301

Quartz 2200

Dynamic viscosity Solution µ 7.4 × 10−4 kg_/_m·s [NH92]

Thermal conductivity Solution λ 43.0 W_/_m·K [YSFN02]

Source 23.7

Quartz 2.0

Specific heat Solution cp 387 J_/_kg·K [NH92]

Source 976

Quartz 1200

Thermal expansion coefficient Solution βT 0.94 × 10−4 K−1 [NH92]

Solutal expansion coefficient Solution βC 0.012 1_/_wt% [RCF85]

0.0053 1_/_mol% [Yil06]

0.0017 from Pb-Sn ratio

Electric conductivity σ 1.405 × 106 Ω−1m−1 [SVZ96]

Diffusion coefficient Solution D 2.5 × 10−8 m2

(31)

3 Numerical Solution

To investigate the flow and concentration in the melt the governing equations and boundary conditions are discretized and solved numerically. The pre-processors used to mesh the solution domain are Ansys ICEM and the OpenFOAM utility blockMesh. The ICEM grid was used together with FASTEST 3D and the blockMesh grid together with OpenFOAM to solve the equations introduced in the previous chapter.

3.1 Solution properties

3.1.1 Boundedness

An important criterion for the usability of a numerical solver is its boundedness. Certain bounds are implied for the solution variables, e.g. a density greater than zero or a con-centration between zero and one, and must be maintained by the solver. Unfortunately, boundedness can only be guaranteed by first order schemes. With higher order schemes unbounded solutions can occur, but this problem can often be resolved by using finer grids.

3.1.2 Consistency

A consistent numerical scheme guarantees that the difference between the exact derivative and the discretized differentiation operator becomes zero in the limit of vanishing grid spacing ∆x or time step ∆t. This difference is called truncation error, because it arises through truncating the Taylor series used to deduce the differentiation operator. It does not provide information about the global quality of the solution, but how large the local error is made at every time step. For its evaluation we assume we have an exact approximation at the time step ti−1/ the node xi−1and then evaluate the error made to obtain the function

value at the time step t / node x.

The performance of a scheme can be classified by the dependence of the truncation error on the grid spacing; an n-th order method will have a truncation error that decreases proportional to (∆x)n_{. If the truncation error is proportional to a ratio of ∆x and ∆t}

(32)

CHAPTER 3. NUMERICAL SOLUTION

special care has to be taken to insure consistency of the solution. An example for an additional condition to ensure consistency is the Courant-Friedrichs-Lewy condition (CFL, [CFL28]; see also 3.2.1).

3.1.3 Stability

Another important quality criterion is the stability of a solution method. Stable methods ensure that approximation or round-off errors (which are like noise around the real solution) are not amplified but damped and thus stability is an important condition for convergence.

3.1.4 Convergence

If the difference between the discretized solution and the exact solution vanishes with a grid spacing that goes to zero, the solution method is convergent. This is similar to the definition of consistency, however now the global error is evaluated instead of the local truncation error. According to Lax’s equivalence theorem, convergence follows stability and consistency for a properly posed initial value problem.

3.2 The finite volume method

FASTEST and OpenFOAM use the finite volume method to discretize the governing partial differential equation. The numerical scheme can be formulated for the generalized transport equation 3.1 in its integral form:

ˆ V ∂ρΦ ∂t dV + ˆ V ∂ρuiΦ ∂xi dV = ˆ V ∂ ∂xi ΓΦ ∂Φ ∂xi dV + ˆ V QΦdV (3.1)

Using Gauss Theorem, the volume integrals can be converted to surface integrals: ˆ V ∂ρΦ ∂t dV | {z } I + ˆ S ρuiΦdS | {z } II = ˆ S ΓΦ ∂Φ ∂xi dS | {z } III + ˆ V QΦdV | {z } IV (3.2) I: Unsteady term II: Convection term III: Diffusion term IV: Source term

(33)

As can be seen in equation 3.2, three different approximations have to be made to solve the equation: The time derivative, the volume and the spacial derivative. In the following sections 3.2.1 - 3.2.4 approximations to solve the four terms will be introduced.

3.2.1 Unsteady term

Computation of the unsteady term involves computation of the CV volume and the time derivative (equation 3.3). Computation of the volume is trivial for cuboid elements.

ˆ V ∂ρΦ ∂t dV ≈ _∂ρΦ ∂t ∆V P (3.3) The time derivative can be approximated using various schemes. Implicit schemes are inherently stable and therefore can be used with arbitrary time steps. Like explicit methods larger time steps will lead to less accurate solutions, but always follow the accurate solution while explicit methods might start oscillating. Explicit methods have to obey the Courant-Friedrichs-Lewy condition. The CFL condition depends on the form of the PDE examined. For the hyperbolic equation 3.4 the CFL condition is shown in equation 3.5 as it was originally derived [CFL28]. ∂2u ∂t2 = ∂2u ∂x2 (3.4) ∆t ∆x < C (3.5)

The CFL condition for solving the transport equations is shown in equation 3.6. The value of the expression is often referred to as the Courant number Co. It assures that no fluid particle can travel further than a mesh step ∆x in a time period ∆t. The derivation of this condition can be found in the referenced literature [FP02].

Co= u · ∆t

∆x <1 (3.6)

Generally it takes more programming effort to implement implicit methods, but this is not a concern when readily available software packages are used. Even when using implicit methods sometimes the CFL criterion has to be monitored, e.g. when using a segregated PISO algorithm to couple pressure and velocity in CFD codes (see e.g. [Iss86]). Automatically adjusted time steps based on a Courant number limit help to achieve the best possible performance.

(34)

A first order scheme (backward Euler / implicit Euler method) is shown in equation 3.7, where n, n+1 denote the time step and the index P denotes that the derivative is evaluated at the cell centre of a control volume.

_∂ρΦ ∂t n+1 P ≈ (ρΦ) n+1 P −(ρΦ)nP ∆t (3.7)

Equation 3.8 shows a second order implicit scheme.

_∂ρΦ ∂t n+1 P ≈ 3(ρΦ) n+1 P −4(ρΦ)nP + (ρΦ)n−1P 2∆t (3.8)

A hybrid second-order scheme which mixes implicit and explicit time marching is known as the Crank-Nicholson scheme. Like other implicit schemes it is unconditionally stable but not bounded.

3.2.2 Convection term

A surface integral has to be calculated to compute the convection and diffusion terms of equation 3.2. One of the advantages of the finite volume method is that this equation is valid for the solution domain as well as each control volume, so the integral can be written as the sum of fluxes through the control volume faces (equation 3.9).

ˆ s ρuiΦdSi = X f ˆ sf ρuiΦdSi,f (3.9)

A simple, second order accurate approximation of the integral over a face is the midpoint rule (3.10), where it is assumed that the integrand is constant over the cell face surface.

ˆ

sf

ρuiΦdSi,f ≈(ρui)fΦfSf (3.10)

The value of Φ is defined at the centre of the control volumes and therefore has to be interpolated between two cell centres to obtain the value at the face centre needed in equation 3.10. It is advisable to use a second order approximation, since that would keep the overall error second order. Using a linear interpolation between two neighbouring nodes is simple and retains the second order error. Equation 3.11 exemplarily shows the interpolation for the eastbound face of a control volume (see also figures 3.1 and 3.2). This is analog to the central difference scheme in finite differences and is thus called CDS.

(35)

CHAPTER 3. NUMERICAL SOLUTION Figure 3.1: 3D control volume

Φe= (1 − λe) ΦP + λeΦE, λe=

xe− xP

xE− xP (3.11)

Unfortunately, the CDS approximation produces oscillatory results in the presence of strong convection, i. e. for Péclet numbers larger than two. Refining the grid and thus lowering the cell Péclet number is a remedy, however this will lead to quickly growing equation systems in 3D simulations. Another approach is to use a first order method, since they can guarantee boundedness unconditionally (see also section 3.1.1). FASTEST implements the upwind discretization scheme (UDS) according to equation 3.12.

Φe= max(ux,0)ΦP + max(−ux,0)ΦE (3.12)

A disadvantage of UDS besides being only first order accurate is that it is numerically diffusive, which will lead to smoothing of rapid changes of a variable and smearing out of numerical errors. A straightforward way to reduce the negative effects of UDS while still inhibiting oscillations is mixing both methods controlled by a parameter β. This is known as “flux blending” and is shown in equation 3.13.

Φ = βΦCDS+ (1 − β)ΦU DS (3.13) 3.2.3 Diffusion term

As shown in section 3.2.2 for the convective term, the integral can be expressed as the sum of fluxes through the control volume faces (equation 3.14). Again the midpoint rule is used

(36)

CHAPTER 3. NUMERICAL SOLUTION Figure 3.2: 2D view of a control volume

to approximate the integrand over the cell face, which is demonstrated for the east facing surface in equation 3.15. ˆ S ΓΦ ∂Φ ∂xi dS=X f ˆ Sf ΓΦ ∂Φ ∂xi dSi,f (3.14) ˆ Sf ΓΦ ∂Φ ∂xi dSi,f ≈ ΓΦ ∂Φ ∂x e∆S x,e (3.15)

The diffusion term does not cause convergence problems, so the CDS approximation (equation 3.16) can be used without any restrictions. It is noteworthy that the CDS ap-proximation will work best for uniform grids. For non-uniform grids, some of the accuracy is lost, but it still outperforms UDS.

_∂Φ ∂x1 e ≈ ΦE−ΦP xE− xP (3.16) 3.2.4 Source term

A simple second-order approximation for the source term is replacing the integral by a product of the cell volume (which is easy to compute especially for cuboid elements) and the value of the source term at the cell centre, as shown in equation 3.17. The pressure term in the momentum source requires special care (see section 3.3). Similarly, for source terms arising from an applied magnetic field a special treatment may be needed and is described in section 2.4.

ˆ

V

(37)

3.2.5 Algebraic equation system

Assembling all four discretized terms and expressing the equation system in matrix form will lead to a sparse coefficient matrix A, the vector of unknowns Φ and the right hand side of the equation Q (equation 3.20). Discretizing a three dimensional domain leads to a matrix A with seven non-zero diagonals: The main diagonal with coefficients for ΦP, and

the sub-diagonals with coefficients for the surrounding points ΦE, ΦW, ΦN, ΦS, ΦT and

ΦB. The dimension of the matrix is equal to the number of control volumes in the domain.

The assembly will exemplarily be demonstrated for terms II and III of the generalized transport equation 3.2. IIf − III f = ˙mcfΦU DSF + γ ˙mcf(ΦCDSF −ΦU DSF ) | {z } QDF f + ˙md f(ΦP −ΦF) = ˙mc fmax( ˙mcf,0)ΦP + ˙mcfmax( ˙−m c f,0)ΦF + ˙mdfΦP − ˙mdfΦF + QDFf = Af_PΦP + AFΦF + QDFf (3.18)

Summing equation 3.18 over all faces of a control volume leads to the computational molecule shown in equation 3.19. Doing so for all control volumes leads to the equation system 3.20, which is preferably solved with multigrid or conjugate gradient methods due to their superior performance.

X

f

IIf − IIIf = (AeP + AwP + AnP + AsP + AtP + AbP)ΦP + AEΦE

+AWΦW + ANΦN + ASΦS+ ATΦT + ABΦB+ QP (3.19)

A ·Φ = Q (3.20)

3.2.6 Treatment of boundary conditions

Two types of physical boundary conditions and a virtual boundary condition are needed to simulate silicon dissolution: Dirichlet boundary conditions at walls for velocity, tem-perature and species, a von-Neumann condition for temtem-perature and connective conditions between virtual block boundaries in multiblock meshed domains. The latter is merely

(38)

an instruction to copy node values from one block to another and will therefore not be mentioned further on.

Dirichlet boundary conditions are equivalent to setting a constant value for the value of Φ at the boundary location. Therefore the corresponding term can be moved to the right hand side of the equation (3.21).

APΦP + AEΦE + AWΦW + ANΦN + ASΦS+ ATΦT = QP − ABΦB (3.21)

Von Neumann boundary conditions (3.22) are used when setting a constant flux at a boundary. It is equivalent to copying close boundary values to the boundary nodes.

∂Φ ∂xi = 0 → ΦP −ΦN xi,P − xi,N = 0 → ΦP = ΦN (3.22) (AP + AB)ΦP + AEΦE + AWΦW + ANΦN + ASΦS+ ATΦT = QP (3.23)

3.3 Pressure coupling

For incompressible flows, solving the momentum equations is not straight forward. The three momentum equations have to be used to solve for pressure and velocity, which add up to four variables. This problem does not arise for compressible flows because the con-tinuity equation can be solved for the additional density variable and the pressure can be connected to the density through an equation of state. For incompressible flows the conti-nuity equation reduces to ∇ · u = 0, which contains no information about the pressure and only poses a scalar constraint on the velocity vector. It is possible to solve the momentum equation using the constraint by computing the pressure and velocity simultaneously, how-ever the computational cost of this fully coupled approach will be multiple times higher than that of the segregated approach described in the following.

Although derivation of a pressure equation would be possible with the differential form of the parent equations, it is crucial that all terms in the derived pressure equations are dis-cretized in the same way as in the parent equations and therefore the disdis-cretized equations will be used.

ρ δt(u

∗

i − uni) = H(u∗i) − ∆ipn+ Si (3.24)

Equation 3.24 shows the base form of the discretized momentum equation, where un i

(39)

CHAPTER 3. NUMERICAL SOLUTION and pn_{are the “guessed” velocity and pressure, ∆}

i is the discretized gradient operator, Si

contains the source terms independent of u and H contains all other terms. In practise the guessed velocity and momentum fields will either be a pre-set initial condition or the values from the previous time-step.

To increase stability, it is helpful to separate the diagonal elements of H and treat them implicitly (equation 3.25). ρ δt(u ∗∗ i − uni) = H 0_(u∗ i) + A0u∗∗i −∆ipn+ Si (3.25)

Rearranging gives an equation which can be solved for u∗

i in the first predictor step

(equation 3.26). (ρ δt− A0)u ∗ i = H 0_(u∗ i) − ∆pni + Si+ ρ δtu n i (3.26)

The new velocity field that is calculated in the predictor step will not satisfy the con-tinuity equation, so it will have to be corrected with a pressure field that corresponds to the continuity constraint. To derive the pressure equation, the predictor equation 3.26 is subtracted from the velocity update equation 3.27, leaving equation 3.28.

(ρ δt − A0)u ∗∗ i − ρ δtu n i = H 0_(u∗ i) − ∆p ∗ i + Si (3.27) (ρ δt − A0)(u ∗∗ i − u ∗ i) = −∆i(p∗− pn) (3.28)

Applying the divergence operation to the velocity increment equation 3.28 leads to a Poisson type equation (3.29) from which the new pressure field can be determined. Since the new velocity u∗∗

i is supposed to fulfil the continuity equation ∆iui= 0, the corresponding

term can be dropped. The updated pressure field can now be used in the velocity increment equation 3.28 to obtain the corrected velocity.

∆i[(

ρ δt − A0)

−1_∆

i](p∗− pn) = ∆iu∗i −∆iu∗∗i (3.29)

Unfortunately, the convergence of this algorithm is poor. In fact, it will not converge at all unless the computed pressure and velocities are heavily underrelaxed, which is done when using the SIMPLE (semi-implicit method for pressure-linked equations) algorithm. Guidance values are αU = 1 − αP, e.g. αP = 0.2 and αU = 0.8 [FP02].

(40)

underrelaxation by applying a second and possibly third pressure corrector step analogous to the first corrector. However, this is based on the assumption that the momentum discretization can safely be frozen through the series of correctors, which is only valid for small time steps. Although this poses a limitation similar to explicit methods (Courant number limit), the PISO method is very useful for low velocity applications (e.g. flows with magnetic dampening). A detailed comparison of SIMPLE and PISO type schemes can be found in [Bar98].

FASTEST uses the SIMPLE algorithm, while the OpenFOAM solver is based on a PISO loop.

3.4 Grids and mesh quality

3.4.1 Grid types

The discretized equations will be evaluated at points defined by a numerical mesh rep-resenting the geometric domain. The mesh can either be structured or unstructured. A structured grid is equivalent to a Cartesian grid - it consists of points which can be uniquely identified by a set of indices i, j, k and grid lines that can be numbered consecutively. Block-structured grids are a subtype of Block-structured grids which divide the solution domain into two levels: A very coarse division which contains large segments of the domain (blocks) and a fine structured mesh within these coarse blocks. The same solvers that are used for simple structured grids can be applied to block-structured grids.

Unstructured grids are very flexible in placing their grid points. Commonly tetraeders or triangles are used, but the control volumes may have any shape and the number of neighbouring nodes is not restricted. Unstructured grids have advantages in geometric flexibility, local grid refinement and dynamic moving of nodes, but are more demanding compared to block-structured grids regarding memory consumption, speed and effort of parallel processing. The resulting equation system no longer has a regular, diagonal struc-ture and many efficient solvers are not applicable to unstrucstruc-tured grids. Unstrucstruc-tured grid generation usually requires the use of complex algorithms.

Grid generation is not trivial and requires some experience. The quality of a structured grid, which is measured in quantities like control volume angles and control volume edge ratio, has a large impact on the convergence of a simulation because the finite volume method is sensitive to grid distortions. Grid distortions appear when mapping the hexag-onal elements to round shapes, e.g. when meshing a pipe. An overview of properties for

(41)

CHAPTER 3. NUMERICAL SOLUTION Table 3.1: Grid types

Structured grid Unstructured grid

Memory matrix has a regular structure and can be stored very compressed

irregular, non-diagonal structure

CPU simple addressing, efficient solvers complex addressing, matrix reordering to reduce bandwidth,

less efficient solvers Geometric flexibility hexahedral elements are

problematic with round shapes

can adapt to any given geometry

Adaptive grids single nodes cannot be added or removed

nodes can be easily added and removed

Moving grids grid movement can lead to unfavourable elements

element shape is more flexible

Parallel processing parallel processing is simple grid decomposition and boundary information exchange is more

involved

both grid types is given in table 3.1.

3.4.2 Grid generation

FASTEST 3D works with block-structured (multiblock) meshes that can be generated by Ansys ICEM, a commercial grid generator. Blocking is an important step during grid generation that ultimately determines the maximum grid quality (see regular vs. O-Grid) and the maximum parallelization efficiency, since parallelization in FASTEST is done by distributing blocks to the available CPUs.

OpenFOAM includes mesh conversion utilities for meshes in various formats. Flu-ent meshes seem to work best and can be exported from ICEM. However the included blockMesh utility works well for primitive shapes and was found more convenient to use. blockMesh is a command-line tool that uses an input file defining key geometric locations, connections between those points and blocking information to generate a mesh.

A straightforward approach to block the cylindrical domain of the dissolution model is to create a bounding box around the domain and map it to fit the round shape, as depicted in figure 3.3a. This approach works reasonably well for cubic domains, but will produce highly skewed elements in the corners of the bounding box when mapped to a round geometry (figure 3.4).

(42)

(a) Bounding box (b) O-grid

Figure 3.3: Blocking strategies

(43)

CHAPTER 3. NUMERICAL SOLUTION Figure 3.5: Final grid shape

A better way to create the blocking is to employ an O-Grid using five blocks: a cuboid block in the centre of the round domain and four trapezoidal blocks around the centre cuboid (figure 3.3b). The trapezoidal blocks produce much favourable elements when mapped to the round geometry. The size of the cuboid can be varied to obtain the best possible mesh quality.

In the cylindrical magnetic force benchmark (see 4.2.2), OpenFOAM turned out to be less sensitive to grid distortions than FASTEST. Grid optimization is a bit easier using ICEM due to the immediate visual feedback it provides. However OpenFOAM is more tolerant towards grid quality so less optimization was needed and using blockMesh was sufficient. FASTEST and OpenFOAM provide two different ways to work with non-optimal grids: FASTEST applies underrelaxation to the equations if convergence is problematic, while OpenFOAM can apply non-orthogonal correctors and automatic time step adjustion. It was found that underrelaxation comes with a greater performance penalty than the non-orthogonal correctors.

The grid used with FASTEST including the source and quartz domains is shown in figure 3.5.

For the simulations done with OpenFOAM the quartz and source domains were not considered and the boundary conditions have to be adapted accordingly. Since multiple

(44)

(a) side view (b) top view

Figure 3.6: High resolution grid used for OpenFOAM simulations

species are only present within the melt, dropping the solid domains only affects the solution of the temperature and induced electric field. The new temperature boundary conditions are based on calculations done with FASTEST and the CFX results by Armour [AD07]. Since the particles are stationary in the source material, no electric field is induced by the magnetic field and the boundary conditions used are still the von Neumann condition given in 2.5.4.

The grid utilized for the OpenFOAM simulations is displayed in figure 3.6. With 220 000 control volumes it is finer than the FASTEST grid. To increase the resolution of the Hartmann layer, the grid spacing was refined towards towards the top and bottom surface.

3.5 Stability considerations with strong magnetic forces

Considering the Lorentz force in its simplest form 3.30, the identity 3.31 can be used to obtain the vector 3.32. The Ai terms in the expanded vector denote additional terms that

(45)

come in when considering the induced electric field, but are not important for the stability analysis. FL= σE((u × B) × B) (3.30) (u × B) × B = (B · u)B − B2_u _(3.31) Fx Fy Fz = σE     Ax+ u3· Bx· Bz−u1· Bz2 − u1· By2+ u2· Bx· By Ay+ u1· By· Bx−u2· B_x2− u2· B_z2+ u3· By· Bz Az+ u2· Bz· By−u3 · B_y2− u3· B_x2+ u1· Bz· Bx     (3.32)

The Lorentz force is a body force that acts on the liquid, and hence is included into the momentum balance equations (2.2) on the right hand side. The quadratic coefficients linked to the solution variable will always be negative and cause convergence problems if kept in the source term. To improve convergence, Kumar et al. [KDD07] suggested an implicit treatment by moving the quadratic expressions to the left hand side of the discretized momentum equations, which increases diagonal dominance of the coefficient matrix. In OpenFOAM, this can be done by using the fvm::Sp() operator as shown in listing 3.1.

Listing 3.1: OpenFOAM implicit source term

- 1/ rho * s i g m a * (( fvc :: g r a d ( p h i E )) ^ B ) + 1/ rho * s i g m a * (( B & U )* B )

(46)

4 Validation

4.1 Solutal buoyancy implementation

The source code modification to support solutal buoyancy (see section 2.3) was validated using a test case presented by Bergman and Hyun [BH96]. Bergman and Hyun conducted two dimensional simulations of the flow, temperature and concentration fields of liquid Pb-Sn exposed to horizontal gradients of temperature and species. From the range of parameters investigated by them, the results for a simulation on a grid of 65 x 65 nodes, a Rayleigh number of Ra = 100 and a buoyancy ratio N = −10 were selected for comparison. The results obtained both with FASTEST and OpenFOAM are in good comparison with the values found in the reference article. Concentration contours as well as streamlines at dimensional times 0.75, 1.35 and 3.15 are shown figure 4.1. Information on the non-dimensionalization can be found in the reference article.

4.2 Magnetic field implementation

Three different implementations for the magnetic field described in detail in section 2.4 have been validated using two studies by Hadid and Henry as reference. In the follow-ing, the setup found in reference [HH96] will be referred to as the cavity benchmark and reference [HH97] as the cylinder benchmark. All three implementations were tested with OpenFOAM. With FASTEST only the simple form without an additional equation and the induced electric field were investigated.

The cavity benchmark was evaluated at Ha = 100 and Gr = 1 · 104_{. The cylinder}

benchmark parameters were Ha = 150 and Gr = 5·104_{. The configuration of the magnetic}

field aligned with gravity turned out to be the most problematic to reproduce accurately, so the results presented in the following sections will be based on that configuration.

(47)

CHAPTER 4. VALIDATION

Figure 4.1: Solutal buoyancy benchmark results at non-dimensional times 0.75, 1.35 and 3.15 (concentration contour left, temperature contour + streamlines right)

(48)

Figure 4.2: Horizontal Bridgman benchmark with no additional equation for the magnetic force, OpenFOAM

4.2.1 Cavity benchmark

Compared to the reference paper, a slightly higher grid resolution of 35x19x19 was used.

No additional equation

Without any additional equation, transversal flow is inhibited everywhere except at the very ends of the pipe. The longitudinal flow is constant over the width of the channel and varies linearly over its height. This only very roughly matches the benchmark results. A plot of the velocity vectors is shown in figure 4.2. The FASTEST results in figure 4.3 are similar.

Variable magnetic field

Benchmark results using this approach were better than the previous approach, but also did not show flow patterns that satisfy the reference paper results. The global flow pattern is similar, but the curvature of the velocity profile (x and z components) is far less pronounced. A plot is shown in figure 4.4.

(49)

Figure 4.3: Horizontal Bridgman benchmark with no additional equation for the magnetic force, FASTEST

Influence of static magnetic fields and solutal buoyancy on silicon dissolution into germanium melt