Evolution of morphology in solidifying aluminium alloys

Evolution of morphology in

solidifying aluminium alloys

Evolution of morphology in

solidifying aluminium alloys

Proefschrift

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus Prof. dr. ir. J.T. Fokkema, voorzitter van het College voor Promoties,

in het openbaar te verdedigen op dinsdag 13 maart 2007 om 15.00 uur door

Winrik Onno DIJKSTRA

Dit proefschrift is goedgekeurd door de promotor:

Prof. ir. L. Katgerman

Samenstelling promotiecommissie:

Rector Magnificus voorzitter

Prof. ir. L. Katgerman Technische Universiteit Delft, promotor Dr. ir. C. Vuik Technische Universiteit Delft,

toegevoegd promotor Prof. dr. M. Cross University of Wales Swansea Prof. dr. R.M.M. Mattheij Technische Universiteit Eindhoven Dr. J. Bruining Technische Universiteit Delft Prof. dr. B.J. Thijsse Technische Universiteit Delft Prof. dr. J.Th.M. de Hosson Rijksuniversiteit Groningen

This thesis is part of the research programs of the Foundation Netherlands Institute for Metals Research (NIMR) and the Foundation for Fundamental Research on Matter (FOM).

Copyright c 2006 by W.O. Dijkstra

All rights reserved. No part of the material protected by this copyright notice may be reproduced or utilized in any form or by any means, electronic or mechanical, including photocopying, recording or by any information storage and retrieval system, without the prior permission of the author.

Contents

1 Introduction

1.1 Outline of thesis . . . 2

2 Scale invariant description of solidification

2.2 Constitutional equations . . . 8

2.3 Scaling of the constitutional equations . . . 10

2.3.1 Dimensionless numbers . . . 12

2.4 Landau transformation . . . 13

2.5 Observations . . . 15

3 Numerical approximation of solidification

3.2 Time–integration . . . 20

3.3 Boundary conditions . . . 21

3.4 Finite difference discretization . . . 22

3.5 Interface conditions at the liquid–solid interface . . . 26

3.6 How not to include the interface conditions . . . 30

3.7 The incomplete limit method . . . 33

4 Stability and accuracy of the solidification model

ii CONTENTS

4.2 An a posteriory estimate of conditions . . . 38

4.3 Influence of the boundary condition location . . . 39

4.4 Criteria for the solute concentration . . . 40

4.4.1 Low BCG residual domain . . . 40

4.4.2 Restrictions for long time durations . . . 49

4.4.3 Estimates of the interface solute concentration . . . 52

4.4.4 Influence of the interface condition . . . 61

4.5 Criteria for the temperature . . . 66

5 Stability related to the interface conditions

5.2 Stability of the Incomplete Limit Method . . . 77

5.3 The overshooting problem . . . 81

6 Discretization with the collocation method

6.2 Model description . . . 89

6.3 Residual method . . . 90

6.4 Using an overdetermined set of system equations . . . 91

7 A numerical 1–D simulation of solidification

7.2 Simulation of solidification . . . 95

7.3 Stability of the results . . . 96

7.4 Influence of solidification shrinkage . . . 100

8 Network model of fluid flow in semi–solid aluminum alloys

8.2 Model description . . . 107

8.2.1 Basic concept . . . 108

0.0 CONTENTS iii

8.2.3 Solidification and melting . . . 111

8.2.4 Mixing . . . 113

8.2.5 Conduction and Diffusion . . . 113

8.3 Observations and numerical results . . . 114

9 The solidification rules of the network model

9.2 Parameter analysis of the solidification rules . . . 121

10 How heat is spread in the network model

10.2 Description of the model geometry . . . 132

10.3 Heat transport in a heterogeneous mush . . . 133

10.3.1 Model parameters . . . 134

10.3.2 Physical parameters . . . 137

11 Discussion and conclusions

11.2 Numerical methods of the 1–D solidfication models . . . 146

11.3 Single time step of the 1–D solidification model . . . 146

11.4 Numerical experiments with the 1–D solidification model . . . 147

11.5 Description of solidification with solidification rules . . . 149

11.6 The network model . . . 150

11.7 Characteristics of the network model . . . 150

Bibliography

Symbols

Summary

iv CONTENTS

Acknowledgment

Chapter 1

Introduction

During the solidification of molten alloys often dendritic shaped grains will be formed. The grains grow in size while more and more liquid becomes solid. If enough grains with a large enough extend have been formed, the grains start to entangle each other. Eventually, they form a relative rigid solid skeleton. In the skeleton there exist usually many liquid pockets and channels.

The skeleton together with the liquid forms a mushy zone. While solidification proceeds the mush changes its morphology considerably. This development of the morphology is quit important for the as cast materials properties. Today, the morphological development of the mush is too complex to be simulated in all de-tail. In this thesis a new approach to the numerical simulation of the solidification and the fluid flow within the mush is investigated.

For the simulation of the fluid flow within the mush a network model is presented. The study of the overall fluid flow requires that a simulation should be capable to simulate the development of the morphology within a whole region of the mush. Such a region has usually a considerable amount of grains included. It is therefore inevitable to restrict the simulation of the local morphology around a single grain to the most essential aspects.

The main idea behind the network model is that the overall morphology of the liquid pockets, pores, and channels within the mush can be represented by a network of connected and geometrically simplified channels. The transition from morphological complex grain structures to straight channels induces the required restriction.

Network models have since a fairly long time been successfully used to simulate the fluid flow through porous media. Most networks have in common that the pores filled with the liquid do not change at all. In that sense they are static. Contrary to the static models, the network model used in the simulations of the

2 Introduction

mush is dynamic. In the model the channels will change their size to account for the ongoing solidification.

One important part of the network model is therefore the solidification within individual channels. The solidification rules actually used within the network simulation are quite simple. The reason for this is that the numerous channels within the network prohibits more complicated rules which require more compu-tational time to determine the amount of growth.

The complex processes of solidification that occur in reality within the channels need themselves ambitious numerical simulations. One can therefore be satisfied when a solidification model describes some aspects of solidification reasonable well. However, to get a first impression of the solidification as could occur within a channel an accurate one dimensional solidification model is numerically investi-gated.

As will be shown with this solidification model, it is already a challenge to solve the model equations numerically. It is therefore sensible to treat the solidification model independently from the network simulation. Furthermore, the usable set of model equations of the fairly general approach will be restricted mainly by the numerical solvability.

1.1

Outline of thesis

The chapters at the beginning of this thesis introduce the one dimensional solid-ification model. The topic of the following chapters is the network model. In the final chapter results are discussed and conclusions are given.

In Chapter 2 the constitutional equations of the one dimensional solidification model are given. For the numerical solution it is advantageous to represent these equations in a dimensionless and Landau transformed form.

The following Chapter 3 describes some Finite Difference discretisations of the model equations. Various possibilities exist to implement the interface condi-tions. Depending on the chosen implementation the numerical solution might be successful or not.

The stability and solvability of the model is numerical studied in Chapter 4. Apart from the well known numerical instabilities new ones due to the inclusion of the interface are discovered.

The effects of the different discretizations, respectively of the possible variants to include the liquid–solid interface, are investigated in Chapter 5. Results of the Incomplete Limit Method are presented here. At the end of this chapter the problem of overshooting is addressed.

1.1 Outline of thesis 3

Difference Method. The Collocation Method shows similar or even worse insta-bilities.

An example of 1–D solidification with a Finite Difference Method is given in Chapter 7. It shows that with quite a lot of care the instabilities can be avoided. The needed criteria are discussed.

Chapter 8 introduces the network model. Some typical network simulations are described. Finally the main results of fluid flow within the mush predicted by the model are enclosed within this chapter.

In Chapter 9 solidification is discussed again. This time the simplified rules used within the network model are investigated. It turns out that the rules behave well, i.e. they include all essential effects which one might expect to find within the network model.

The heat distribution within the network model is investigated in Chapter 10. Apart from the heat distribution the influence of a number of relevant parameters is studied.

Chapter 2

Scale invariant description

of solidification

2.1

Introduction

The physical basis of every model describing solidification is the assumption that across the liquid–solid interface mass, heat, solute and momentum is conserved. Expressed in mathematical form the conservation of these properties leads to four conservation equations. In the text below they are also called interface equations. Together with the fundamental equations, i.e. the assumption of incompressible media, Fick’s law, the Fourier equation, and the fundamental equation of dynam-ics, these equations describe in principle solidification completely.

Though this basis is very general, it is about the only common basis among all numerical solidification models of two phases. The equations can be brought in various forms. This reformulation allows to include additional assumptions, to reduce the mathematical complexity, and to adopt a more suitable form which might be used for specific numerical methods. In literature various derivations of such equations exist. Here a derivation similar to the one as can be found in the book of Rappaz et al. [46] is used.

As the derivation is common practice in physics there is no need to repeat the whole derivation here again. There are only one or two points that ought to be clarified in order to understand deviations from other models found in literature. The general form of mass balance at the liquid–solid interface is

vs 

_x=Iρs− vl 

_x=Iρl= vi(ρs− ρl). (2.1) It is assumed that the liquid–solid interface in the 1–D model is at position I. The liquid of density ρl is placed on the right side while the solid phase with density

6 Scale invariant description of solidification

ρs is on the left side of the interface. The interface moves with velocity vi and the liquid and solid velocities are vl, respectively vs. Since the network models showed a tendency to depend on solidification shrinkage, it was chosen explicitly to keep the individual densities of both phases. Many models found in literature already deviate here by assuming equal densities.

Another deviation is encountered when introducing the almost obvious assump-tion of a non moving solid phase (vs = 0). Equation (2.1) can then be written as vl  _x=I = vi  1 − ρs ρl  . (2.2)

In 1–D only the normal component of the general form of conservation of momen-tum is used ρs (vs  _x=I)2− vi2
− ρl (vl  _x=I)2− vi2
= σs− σl. (2.3) It will change accordingly to

(ρl+ ρs)vi2− ρl(vl 

_x=I)2= σs− σl. (2.4)

The normal stresses at the interface σl and σs are usually written in terms of the kinetic mobility and surface stresses. These terms are commonly used to introduce the Gibbs–Thompson relation. In this case the terms depend on the actual interface curvature. Equations (2.2) and (2.4) form two linear equations with the two unknowns vi and vl

_x=I. Eliminating vl  

_x=I one can readily derive a relation for the interface velocity vi. Most models in literature use this kind of procedure.

In the model investigated here and for the purpose of the network model the Gibbs–Thompson relation is neglected. In 1–D no sensible choice of interface curvature exists. Furthermore, if solidification shrinkage is explicitly included in the model, it is quite natural to keep the option of an additional inclusion of thermal shrinkage. To allow for this the solid movement vsthough small cannot be set equal to zero. It depends on the fundamental equations of dynamics used to determine stress and strain relations. One might therefore ask whether there is no other means to determine the interface velocity vi.

In statistical physics some media are known for which the Gibbs–Thompson re-lation is of minor importance for the micro-structure formation of solids. The main reason is that on small scales inter-atomic binding forces might dominate completely. In metal alloys the Gibbs–Thompson relation is certainly not to be neglected, but for the growth rate of a flat interface an (semi) empirical equation exist. As stated in the book of Saito [50], for crystal growth from the melt the growth is described by Wilson–Frenkel formula

2.1 Introduction 7

The growth rate vi depends on the kinetic coefficient K, the interface tempera-ture T and the difference of the free energy δµ between both phases. kB is the Boltzmann constant. For a small undercooling δT at the interface, i.e. a small deviation from the melting temperature and not too low temperatures, the growth rate can be approximated by

vi= βδT. (2.6)

This relation is only an ideal one, but as explained in Galenko and Zhuravlev [23] the kinetic coefficient β can be adapted to experimental observation. The ex-pression for the actual interface undercooling can then also be used to include the Gibbs–Thompson relation adequately within the model. In the model inves-tigated here, only an undercooling due to deviations of the interface temperature and concentration from the liquidus is considered.

As this is a description of a first approximation, thermal shrinkage is neglected. For obvious reasons Gibbs–Thompson relation will neither be considered. On the one hand, it is therefore assumed that the solid will not move (vs = 0), that the conservation of momentum (2.4) is fulfilled without notice and that mass conservation (2.2) can be used to determine the liquid velocity at the interface. On the other hand, Wilson–Frenkel formula (2.6) is added to the model equations in order to get the right atomistic interface behavior.

There is just one additional note to be made on the use of the Wilson–Frenkel formula. The formula depends heavily on the atomistic roughness of the interface and fails in case of a smooth interface. For metals one has usually a rough interface1_{. Because of the required roughness solidification, as investigated here,} can therefore not be diffusion driven. Even in case of melting the model remains different from models describing pure dissolving. Or in other words, diffusion is in the model of crucial importance but if diffusion is too slow the boundary layer in front of the interface will pile up rather than decelerate the growth rate. Though there are many numerical approximations and simulations of solidifica-tion described in literature only a few similar models were found. Nearly the same model equations were used in Galenko and Zhuravlev [23] in some varia-tions. The essential difference to the work presented here is the kind of numerical approximation used. The approximation of Galenko and Zhuravlev was based on (phase) averaging and the application of an appropriate finite volume method. The averaging is in the present model not considered because it smooths out the boundary layer in front of the interface. For this reason it is difficult to combine the averaging with the Wilson–Frenkel formula in a sensible way.

The pile up of solute in front of the interface itself is considered in various other works. Notable is the work by Crusius et al. [13] which is now part of the foun-dation of the software package DICTRA [14]. Nowhere clearly mentioned, this is a diffusion driven model par excellence. It is therefore quite a different model as mentioned above. The ingenious method to solve the model equations does really

8 Scale invariant description of solidification

not alter anything about this. In contrary, it is not clear how the method could be used in a non diffusion driven model.

Another model which resembles the one proposed here is described in Zhao et al. [62]. It differs as explained above in the sense that it uses a generalized Gibbs– Thompson relation. The interface velocity however is determined indirectly using the conservation of heat at the liquid–solid interface. Apart from this a standard weighted residual method is applied which will not be used here.

So far it is explained why some other models differ from the present one. In addition the relations to similar models have been mentioned. In the following sections the model equations for solidification of a binary alloy will be stated. Next chapter will be dedicated to the numerical methods used to solve the equation system. It might then become even more clear, that most models differ not only in the model equations but also in the method used to find a numerical solution.

2.2

Constitutional equations

The domain of the 1–D model is [0,X]. Initially, the left part [0,I] of the domain is solid while other part is liquid. Upon solidification the solid part will grow and the solid–liquid interface which has distance I from the left side of the domain will move to the right. The most natural representation of this situation might be a long solid slab of small diameter and length I. The slab is connected to a liquid cavity of maximal extension X.

As already stated, in the liquid and solid a number of fundamental equations describe the distribution of the heat and the solute concentration. These equations are the Fourier laws

∂Tα

∂t − aα∆Tα= 0, (2.7)

and Fick’s laws

∂Cα

∂t − Dα∆Cα= 0. (2.8)

Here the temperature Tαis either the temperature in the liquid, i.e. for α=l or the temperature in the solid (for α=s). The same holds for the solute concentration Cα. Only one single solute, i.e. copper, is incorporated for the solidification of the bi–alloy. Dαis the diffusion constant of the liquid or solid and aαis the thermal diffusivity in the liquid or solid respectively. Convective heat and mass transport have been neglected here. Convection in a freezing slab is small due to the small flow velocities involved.

2.2 Constitutional equations 9

At the solid–liquid interface I, mass, solute and heat are conserved as well. Using mass conservation (2.2) two additional interface equations can be derived out of the fundamental equations of solute and heat:

vi  k − ρs ρl  Cl  _x=I = Dl ∂Cl ∂x   _x=I+− Ds ∂Cs ∂x   _x=I−, (2.9) and vi(ρl− ρs)cplTl  _x=I− viL = κl ∂Tl ∂x    x=I+− κs ∂Ts ∂x    x=I−. (2.10)

In these last equations κs and κl are the conductivity in the solid and liquid which are not to be mistaken for the partition coefficient k (see Equation (2.13)). Furthermore, cpl is the specific heat within the liquid while, L is the latent heat

of fusion. It is indicated whether the derivative of the solute or liquid profile has to be taken.

The velocity viof the interface will be determined by the Wilson–Frenkel formula, which for a small undercooling δT at the interface can be approximated by

vi= β  Tf− Tl  _x=I + mCl  _x=I  . (2.11)

This equation couples conservation of solute and heat at the interface with the help of the kinetic coefficient β. The melting temperature of the pure solvent is Tf and the liquidus of the binary phase diagram has slope m.

Finally at the interface the temperature on both sides of the interface remains the same while the solute concentration has a discontinuity, i.e. at the solid liquid interface Ts  _x=I = Tl  _x=I, (2.12) and Cs  _x=I = kCl  _x=I. (2.13)

Last Equation (2.13) is not due to conservation of solute but is a consequence following directly out of the binary phase diagram.

These model equations imply two additional restrictions. Firstly, rapid solidifica-tion processes have to be excluded. In rapid solidificasolidifica-tion Equasolidifica-tion (2.13) is not necessarily fulfilled at the interface. This restriction will pose usually no problem. Only very extreme boundary conditions (BCs) have to be avoided.

10 Scale invariant description of solidification

This second restriction is far more difficult to meet. A sensible choice of initial solute concentration and temperature profiles avoids to encounter an impossible situation just at the very beginning. Even when the initial conditions are selected with care, the solute in front of the interface might pile up. The concentration at the interface might in consequence rise above any reasonable value. This effect, termed overshooting, will be discussed later on.

Not a real restriction but rather implicit assumption is the fact that the model is supposed to describe solidification of an adequate eutectic alloy like Al–Cu. Only solidification with solute concentrations lower than the eutectic concentration is admissible. Eutectic growth itself is not investigated here.

2.3

Scaling of the constitutional equations

With the model equations a number of variables have been introduced, i.e. the temperature profiles Tα, the solute concentration profiles Cα, and the velocity of the interface vi. Most other quantities are constants and material properties. Later on, two additional constants, i.e. the mesh size h and the time interval δt are introduced by the discretization. In order to make the solution more or less independent of the actual material and discretization, the variables will be scaled in such a way that the quantities within the equations become non dimensional. Using Xσ as the specific length and tσas the specific time duration, the following scaling will produce equations with in total 12 different non dimensional num-bers (i.e. Fohl, Fohs, FoDhl, FoDhs, Ex, He, Tr, Sp, Kih, k, Xχ, and hχ). The numbers will be discussed in next section. The scaling of the temperature, solute concentration and interface velocity are

Tα= 1 Tf Tα, (2.14) Cα= 1 Ceut Cα, (2.15)

with the eutectic concentration Ceut,

ν = tσ Xσ

vi. (2.16)

The direction within the slab (or domain) x and time t will be scaled as well:

2.3 Scaling of the constitutional equations 11

The Fourier number for heat conduction Fohα, the Fourier number for diffusion FoDhα and the P´eclet number Pehαare dimensionless numbers common for heat transfer equations: Fohα= καtσ ραcpαXσ 2, (2.19) FoDhα= Dαtσ Xσ2 , (2.20) and PeXσα= Xσvi Dα . (2.21)

Though the P´eclet number PeXσα is dimensionless it is not a constant number

which usually occurs in interface equations. It is not constant due to the fact that the interface velocity will change during solidification. However, the use of the P´eclet number can be avoided by applying the relationship

PeXσαFoDXσα= ν. (2.22)

With these dimensionless quantities above equations can be rewritten. Equations (2.7) and (2.8) will become

∂Tα

∂τ − FoXσα∆χTα= 0, (2.23)

∂Cα

∂τ − FoDXσα∆χCα= 0. (2.24)

The equations (2.9) and (2.10) can be written as  k −ρs ρlC l  _χ=I = FoDXσl ν ∂C_l ∂χ   _χ=I+  −FoDXσs ν ∂C_s ∂χ   _χ=I−  , (2.25) and  1 − ρs ρlT l  _χ=I− L ρlcplTf = FoXσl ν ∂T_l ∂χ   _χ=I+  −FoXσs ν ρs ρl cps cpl ∂T_s ∂χ   _χ=I−  . (2.26) Here I is the position of the interface in the scaled coordinate system. The position of the interface depends on time τ , i.e. I = I(τ ). In the scaled coordinate system, the total length of the domain X will become X. The Wilson–Frenkel formula reads in dimensionless form

1 = ν Xσ tσβTf +Tl  _χ=I − mCeut Tf Cl  _χ=I. (2.27)

12 Scale invariant description of solidification

2.3.1

Dimensionless numbers

Apart from the Fourier numbers FoXσα and FoDXσαless common dimensionless

numbers for the constants occurring in above equations can be defined. First, there is an expansion number Ex

Ex = ρs ρl

, (2.30)

which defines the inverse amount of solidification shrinkage. The ratio of the maximal possible heat within a given liquid volume and the latent heat of fusion is defined by the heat number He

He = L

ρlcplTf

. (2.31)

This number is related to the Stefan number2St = cpl_LδT. The relative difference between the eutectic temperature and the melting temperature of the pure solvent is given by the temperature range number Tr

Tr = −mCeut Tf

. (2.32)

The ratio of the specific heat in the liquid and solid forms another number Sp

Sp = cps

cpl

. (2.33)

This number is sometimes referred to as the non–dimensional specific heat. Finally there is a dimensionless number whose meaning is not obvious and which is here called kinetic number KiXσ

KiXσ =

Xσ tσβTf

. (2.34)

Together with partition coefficient k, which is dimensionless as well, and the dimensionless length of the domainX, the number of relevant parameters within the simulation is reduced to 12.

The specific length and time can in principle be chosen arbitrarily. As will be shown in the next chapter, the choice of the mesh size h as the specific length Xσ and the time step duration δt as the specific time tσ, allows to formulate simple criteria for stability. In this work this choice, leading to the dimensionless numbers Fohl, Fohs, FoDhl, FoDhl, and Kih, will be used if not stated otherwise. Other interesting choices of specific lengths are possible. For example, below the total length X of the domain will be used as the specific length. When changing

2_{Note that the Stefan number is based on the undercooling at the interface δT . In the}

2.4 Landau transformation 13

the specific coefficients the form of the model equations will not change but the dimensionless numbers will change accordingly. If the total length X is used, the dimensionless numbers Fohl, Fohs, FoDhl, FoDhl, and Kih will become FoXl, FoXs, FoDXl, FoDXl, and KiX. The number Xχ will become 1 and the new dimensionless parameter will become hχ =_X1h.

In next section a domain transformation is introduced, i.e. the Landau transfor-mation. This transformation is sometimes quite useful especially from a numerical point of view. As will be seen, the Landau transformation does not interfere with the scaling applied to get a dimensionless representation.

2.4

Landau transformation

In principle, the model equations can be solved directly for a suitable situation using the indicated domain. The movement of the interface from the left to the right side of the domain involves some additional numerical considerations. Firstly, the domain of the liquid and solid part will change and the meshing or the choice of the collocation points (see below) has to be adapted during the simulation. Secondly, as indicated in [13] the movement of the interface might lead in some methods to a breach in mass conservation at the interface. Especially when the boundary layer at the liquid–solid interface becomes pronounced chances seem to be high that this will happen.

In the geometrical simple situation of a 1–D domain the interface movement can be incorporated by a Landau transformation. This coordinate transformation will continuously transform the liquid and solid domain independently into the domain [0,1]. In this sense it brings, together with the introduction of dimensionless numbers, the whole model to a more standard situation, which allows to draw fairly general conclusions.

Before formulating the Landau transformation in mathematical terms, it must be realized that the Landau transformation has its own practical limitations. Firstly, in more complex domains than considered in this description, application of this transformation can become elaborate. Secondly, the transformation involves the inclusion of some additional terms within the constitutional equations. These terms might in turn have there own numerical issues.

As already mentioned, the Landau transformation is a coordinate transformation which will transform each point within liquid and solid to the domain [0,1]:

x ∈ [0, I(τ )] → ξ ∈ [0, 1], (2.35)

x ∈ [I(τ ), Xχ] → ζ ∈ [0, 1], (2.36)

with

ξ(χ, τ ) = χ

14 Scale invariant description of solidification

and

ζ(χ, τ ) = χ − I(τ ) Xχ− I(τ )

. (2.38)

The constitutional equations can be rewritten in terms of the new coordinates. It should be noted that the partial derivatives ofTs are

∂Ts ∂χ = 1 I ∂Ts ∂ξ , (2.39) and ∂2Ts ∂χ2 = 1 I2 ∂2Ts ∂ξ2 . (2.40)

The partial derivative of Tswith respect to time becomes ∂T_s ∂τ  χ =∂Ts ∂τ  ξ −νξ I ∂Ts ∂ξ , (2.41)

on the understanding that last two terms of the equation are taken at a constant ξ, and constant time τ respectively. For derivatives of the liquid part, i.e. Tl a similar transformation results. Of course, the same equations hold also for the appropriate solute concentrations profiles.

In terms of the new coordinates the constitutional equations (2.23) and (2.24) become in 1–D ∂Ts ∂τ − νξ I ∂Ts ∂ξ − Fohs I2 ∂2Ts ∂ξ2 = 0, (2.42) ∂Tl ∂τ − ν(1 − ζ) Xχ− I ∂Tl ∂ζ − Fohl (Xχ− I)2 ∂2_T l ∂ζ2 = 0, (2.43) ∂Cs ∂τ − νξ I ∂Cs ∂ξ − FoDhs I2 ∂2_C s ∂ξ2 = 0, (2.44) ∂Cl ∂τ − ν(1 − ζ) Xχ− I ∂Cl ∂ζ − FoDhl (Xχ− I)2 ∂2Cl ∂ζ2 = 0. (2.45)

The conditions (2.25) and (2.26) at the solid–liquid interface will change as well,

(k − Ex)Cl  _ζ=0= FoDhl ν(Xχ− I) ∂C_l ∂ζ   _ζ=0+  −FoDhs νI ∂C_s ∂ξ   _ξ=1−  , (2.46) (1 − Ex)Tl  _ζ=0− ·He = Fohl ν(Xχ− I) ∂T_l ∂ζ   _ζ=0+  −Fohs νI Ex · Sp ∂T_s ∂ξ   _ξ=1−  . (2.47)

2.5 Observations 15 and Cs  _ξ=1= kCl  _ζ=0. (2.50)

Note, the transformed equations remain in dimensionless form. The use of the dimensionless number Xχ is by now also clear. The position of the interface I is not a constant but will change during the simulation.

In case that the total domain length is taken as the specific length of the model, the equations will have following form:

∂Ts ∂τ − νξ fs ∂Ts ∂ξ − FoXs (fs)2 ∂2Ts ∂ξ2 = 0, (2.51) ∂Tl ∂τ − ν(1 − ζ) 1 − fs ∂Tl ∂ζ − FoXl (1 − fs)2 ∂2Tl ∂ζ2 = 0, (2.52) ∂Cs ∂τ − νξ fs ∂Cs ∂ξ − FoDXs (fs)2 ∂2Cs ∂ξ2 = 0, (2.53) ∂Cl ∂τ − ν(1 − ζ) 1 − fs ∂Cl ∂ζ − FoDXl (1 − fs)2 ∂2Cl ∂ζ2 = 0, (2.54) (k − Ex)Cl  _ζ=0= FoDXl ν(1 − fs) ∂C_l ∂ζ   _ζ=0+  −FoDXs νfs ∂C_s ∂ξ   _ξ=1−  , (2.55) (1 − Ex)Tl  _ζ=0− He = FoXl ν(1 − fs) ∂T_l ∂ζ   _ζ=0+  −FoXs νfs Ex · Sp∂Ts ∂ξ   _ξ=1−  , (2.56) 1 = νKiX+Tl  _ζ=0+ Tr ·Cl  _ζ=0, (2.57) Ts  _ξ=1=Tl  _ζ=0, (2.58) and Cs  _ξ=1= kCl  _ζ=0. (2.59)

In above equations the fraction of solid fs is used.

2.5

Observations

16 Scale invariant description of solidification 2 µm 150 µm Fohl 1831.97 0.33 Fohs 3416.21 0.61 FoDhl 0.25 4 · 10−5 FoDhs 3 · 10−5 4 · 10−9 Ex 1.12 1.12 He 1.6 · 10−4 1.6 · 10−4 Tr 0.14 0.14 Sp 0.91 0.91 Kih 5 · 10−5 4 · 10−3 k 0.17 0.17 Xχ 75.00 1.00 hχ 1.00 0.01

Table 2.1: Typical values for the 12 dimensionless numbers for 2 length scales. The constants of Al–Cu are used for the determination of the values and are taken from [32, 55]. The mesh size h is assumed to be 2 µm, the duration of the time interval 2 · 10−4 s, the length X of the whole solidification domain 150 µm, and the kinetic coefficient β is 0.2 m/sK.

solidification domain X was assumed to be about half the size of a typical grain, namely 150 µm.

Within Table 2.1 the Fourier number of the diffusion in the solid FoDhs and the kinetic number Kih are small. The small Fourier number shows that at least for the solidification of Al–Cu alloys the diffusion within the solid is of minor importance. In other words back–diffusion is comparatively small and might be completely neglected.

2.5 Observations 17

Chapter 3

Numerical approximation of

solidification

3.1

Introduction

In the last chapter a number of constitutional equations describing the 1–D so-lidification model have been given. The equations are distinct from those used by most other solidification models. The application of the Wilson–Frenkel re-lation is at least for metal alloys nearly unique. As indicated in the previous chapter, most other models include the Gibbs-Thompson relation. Other models may sometimes deviate as well from the described one in that sense that they leave out parts of minor importance or that equations are reduced in complexity. All these differences have certainly their impact on the results of the simulation. But even though the results might vary quite substantially, the constitutional equations remain similar. For this reason one would expect that the numerical treatment, if it is not the same is very similar for nearly all the models, too. Actually, a rather broad range of methods is used in solidification simulations. However, about three main routes exist which often recur.

The three main routes are averaging, weak formulations, and phase–field like approaches. Averaging is often used in simulations describing solidification in materials of somewhat bigger extend. Usually even whole billets can be treated nowadays. The idea to use averaging on much smaller scale, i.e. within microme-ters, is quite old. The main reason for application of this kind of averaging is its relative simplicity and the possibility to use readily existing finite volume codes. As an alternative for averaging one can rewrite the constitutional equations using a weak formulation. This means that a candidate test function is made to fulfill

20 Numerical approximation of solidification

the constitutional equations when integrated over the appropriate sub–domain. The weak formulation leads in quite a natural way to the application of the Finite Element Method. The conservative approach of taking the sub–domains to be either the liquid or solid parts establishes a front tracking method. Front tracking methods can in complex situations and especially when two interfaces join, become rather cumbersome.

To overcome problems with front tracking, phase–field models have been devel-oped. This is just another kind of weak formulation in which the equations are reformulated in terms of a phase–field function. This function indicates within the whole domain and in a smooth way the transition from the liquid into the solid part. Quite similar ideas where used for the development of other models like the Level Set Method.

In the numerical methods presented here to solve the 1–D solidification model equations none of the above routes was used. Instead a direct method, which couples the solid and liquid profiles directly at the interface has been investigated. The direct method belongs to the front tracking methods and it can in principle solve the equations by applying Finite Difference (FD) discretization schemes. For its simplicity the idea of a direct method might have been investigated a long time ago. It is not quite clear whether this is the case or not. A reasonable search in literature delivered no references to earlier works. One reason might be that, though the principles are clear, there exist a few drawbacks which can lead the method to fail completely. In this section possible deficiencies of the method and how to avoid them will be discussed.

Apart from its simplicity the direct method is preferred to the methods of the main route because of a few other reasons. In averaging it is difficult to know when the interface concentration and temperature are resolved accurate enough to be substituted within the Wilson–Frenkel relation.

The weak formulation requires much more consideration. Especially extensions are far from simple to include. They usually require that the constitutional equa-tions are brought into a weak form anew. Especially the inclusion of thermal shrinkage is difficult within the weak formulation. For the 1–D model described here originally an inclusion of thermal shrinkage was intended. This in turn was supposed to enhance the the 2–D network model.

3.2

Time–integration

3.3 Boundary conditions 21

form

∂Ts

∂τ = OTs. (3.1)

HereTsis used but depending on the equation to be looked at it might be replaced byTl,Cs, orCl. The exact form of the linear operator O depends on the equation, too. At the moment it is enough to know that it is a linear operator which can be discretized later on.

For the time integration of these equations the relative stable backward Euler method is applied. In the backward Euler method above equation (3.1) is ap-proximated by

T0

s−Ts= δτ · OTs0, (3.2)

where δτ is a sufficient small time step andT0

sis the new temperature at time τ + δτ . The last equation can be rewritten in the form

T0 s= A

−1_T

s, (3.3)

where A−1 is the inverse of another linear operator, namely the transfer operator A, which is defined as

A = 1 − δτ · O. (3.4)

The inverse of operator A can be determined by first discretizing the operator and then to invert the resulting matrix A. In this section the discretization will be done with the FD method as explained below.

3.3

Boundary conditions

The transfer operator A does not yet contain any boundary conditions at all1. This means that the interface conditions (2.25), (2.26), (2.29), and (2.28), which are boundary conditions of the moving liquid–solid interface, still have to be included within the system of linear equations.

Furthermore, at the outer boundaries of the liquid and solid additional boundary conditions (BC) exist. For the boundary conditions at the liquid side (at χ = Xχ) of the temperature profile a Neumann condition and a Dirichlet condition at the solid side (at χ = 0) are chosen. For the profile of the solute concentration the conditions are interchanged such that we have a Neumann condition at the solid side and a Dirichlet condition at the liquid side:

Tl  _χ=0=TBC(τ ), (3.5) Cl  _χ=X χ =CBC(τ ), (3.6)

1_{Note that the operator A}−1_{has only been formally introduced. Without inclusion of the}

22 Numerical approximation of solidification ∂Ts ∂χ   _χ=X χ = 0, (3.7) and ∂Cs ∂χ   _χ=0= 0. (3.8)

The idea behind this choice of BCs is that a natural situation resembling the one within the network models should be considered. The BCs state that the liquid temperature TBC(τ ) is conserved at the liquid side. To get solidification started energy is taken out of the system at the solid side. This can be induced by lowering the temperature at the solid boundary. Distribution of the solute is more likely to proceed by convection than by back diffusion. Therefore, the solute concentrationCBC(τ ) is specified at the boundary of the liquid. At the solid side the solute mass should remain conserved.

There are various ways to incorporate the boundary conditions within the sim-ulation. Depending on the method used to discretize the transfer operator, the discrete form A and Ts (or one of the other profiles respectively) can be changed to represent the boundary conditions directly. This is common practice when using the Finite Difference Method (FDM).

For the interface conditions (2.25) and (2.26) at the liquid–solid interface it is convenient to apply an explicit kind of discretization. The interface conditions on themselves are ordinary equations which can be discretized in the same way as the operator A. Ideally, leaving out one mesh point (or later on to retain one collocation point) and including the discrete IC instead within the matrix A will lead to a solution which respects the BC.

Here, the ICs and BCs are directly incorporated into the linear transfer operator. As will be seen later on, the exact way how the interface conditions are handled within a method turns out to be of eminent importance for the accuracy and success of the method.

3.4

Finite difference discretization

Since Fourier’s equation (2.7) or Fick’s equation (2.8) together with Neumann and Dirichlet BCs form a classical system to demonstrate the FDM, not every detail of the basic discretization of the transport operator is discussed here. The first and second derivative is represented by the schemes

3.4 Finite difference discretization 23

The distance between two points of the regular grid is hχandT represents here any of the fieldsTs,Tl,Cl, orCs. Successive values at points on the grid are indicated by Ti−1, Ti, and Ti+1. As stated previously, the discretizations adjoining the boundaries need to be changed in order to include the BCs. All first derivatives disappear at the boundaries even when a Landau transformation is involved. The second derivative with a Dirichlet BC at the solid (left) side becomes

∂2_T ∂χ2   _χ=χ i ' Ti+1− 2Ti hχ2 . (3.11)

The missing term Ti−1is the known boundary value TBC and it becomes now a source term of the form

s|χ=χi = −

TBC hχ2

. (3.12)

If a Neumann condition is needed the second derivative will be only slightly dif-ferent, ∂2_T ∂χ2   _χ=χ i ' Ti+1−Ti hχ2 . (3.13)

With homogeneous Neumann BCs this second derivative is zero and no additional source term is required. Though application of these Neumann BCs is simple it bears a small ambiguity. If the place of the boundary is assumed to be χi, the last point of the solid at which the second derivative is determined might either be taken to be the boundary itself, i.e. χi, or the point nearest to the boundary, i.e. χi+1. It is merely a question from which side, i.e. from the outside or inside of the domain, the derivative ∂T

∂χ   _χ=X

χ

is approximated. In the limit of an infinitesimal small hχboth solutions will be identical. As will be seen later on in our situation of a finite grid spacing, the two possibilities lead to a lower, respectively upper bound of the exact solution. Note, that the use of the approximation scheme (3.9) centered around the boundary corresponds to taking the second derivative at the boundary.

The reason to state this small ambiguity which is usually neglected or at best handled as an approximation error is that the introduction of an additional inter-face with its own interinter-face conditions exaggerates the encountered approximation error. Another reason is that at the solid–liquid interface itself such small changes in the approximation might lead to stunning effects. As will become clear later on it really depends on the way we discretize and how the limit to infinite small grid spacings is taken. This itself is of course no novelty. For example the CFL condition (see for example [46]) just states the same.

The above approximations of the first and second derivative permits the dis-cretization of the transfer operator A. The operator will consist of a matrix A and a source vector s,

24 Numerical approximation of solidification

The vectorT0, which later on becomes the solution vector, contains the values of the field at the individual grid points. If there are in each phase N grid points, the 2N × 2N matrix A belonging to the concentration field Cl and Cs has the form A =         1+FoDhs −FoDhs −FoDhs 1+2FoDhs −FoDhs −FoDhs 1+2FoDhs −FoDhs · · · · −Fo_Dhl 1+2FoDhl −FoDhl −Fo_Dhl 1+2FoDhl         (3.15)

and the source vector s becomes

s =         0 0 · · · · · · · · · FoDhl·CBC         (3.16)

The distance between the grid points h and the time step δt do not occur ex-plicitly within the A and s anymore. These constants of the approximation are already contained within the dimensionless Fourier-Numbers FoDhl = Dlδt/h2 and FoDhs= Dsδt/h2. Within the representation of matrix A the omitted entries are all zero. For the determination of the temperature distribution within the domain, a second similar discretization matrix is considered. Both matrices and the two source vectors only differ in the used dimensionless Fourier-Numbers and in the interchanged BCs.

The generation of the numerical solution starts with either an initial profile or a solution of an earlier simulation. This profile description which will be calledC in the case of a concentration profile is nothing else than the concatenated vectors of the solid and liquid profiles, i.e.

C =         Cs1 Cs2 · · · · · · Cl2N −1 Cl2N         . (3.17)

A vectorC0 of the same structure will be generated as the solution according to the equation

3.5 Finite difference discretization 25

Up to this point the derivation of above equations is similar to the well known FDM used to solve many diffusion and heat transport problems. The last Equa-tion (3.18) suggests that once the inverse of matrix A is determined one can propagate through time by successive matrix multiplications and vector subtrac-tions. This is however not the case in a moving interface problem. As will be shown in next section, when the interface conditions are included matrix A will dependent on the current interface velocity.

This dependence of matrix A on the interface velocity induces the need to solve a rather huge system of linear equations every time step2. This considerably slows down the generation of time advanced solutions. Usually the reason to prefer the FDM or Finite Element Method (FEM) to the collocation theory, as described later on, is that Collocation Theory generates a full transport matrix A while those generated by FDM are sparse. As linear equations with sparse matrices can be solved more efficiently, for example by an iterative method like the Conjugate Gradient Method, and because of the fact that this has to be done for every time step anew, one might be tempted to completely discard the collocation theory3_. The main reason why the FDM does not excel as one might expect at first is that to describe the rather steep pile–up in front of the solid–liquid interface usually requires a fine mesh grid for the FDM. This immediately leads to a huge matrix A while for the description with collocation a matrix of modest size is already suf-ficient. Generally it is difficult to predict when the one or the other method will excel. The reason for this will become clear in next section. The experience that was acquired with a simple implementation, i.e. without adaptive meshing and ’special tricks’ like the Fast Multipole Method, suggests that collocation in form of a Spectral Element Method with a low order polynomial basis and the FDM are about equal demanding in terms of computational efforts.

In Figure 3.1 the computational time required to solve the equation systems for typical determinations of the solute concentration profiles is shown. A Bi– Conjugate Gradient method was used without any preconditioning. This is not the fastest method that might be used with the FDM but the algorithm preferred for analyzing stability4_{. Clearly visible is a near linear dependence of the required} time from the size of transfer matrix A. Interesting is that this behavior seems not to be related to the stability as discussed in next Chapter 4.

2_{The reason is that usually the known inverse of A for a certain interface velocity cannot be}

reused for an efficient determination of A−1belonging to another interface velocity. The same holds as well for generated products like those of the LU–decomposition.

3_{Note that in special cases even more efficient methods than the Conjugate Gradient Method}

exists. For example, for a 1–D domain direct matrix–solvers can be superior (see e.g. [9]).

4_{This thesis does not investigate the fastest possible method to solve the solidification model.}

26 Numerical approximation of solidification 0 500 1000 1500 0 0.5 1 1.5 N ∆ t [s]

Figure 3.1: Computational time required to solve the linear equation system. N is the number of mesh elements per phase. The size of the transfer matrix is 2N × 2N . Shown is the time required for the BCG method without preconditioning and with a duration of the integration step of 2.7 ms.

3.5

Interface conditions at the liquid–solid

inter-face

In the previous section, which describs the discretization with the FDM, the interface conditions at the solid–liquid interface are not mentioned. The omission allowed a simplified representation of the general picture without inclusion of all details. In this section the simplified representation is amended in order to get the interesting moving boundary problem of solidification. This is basically done by specifying specific dots within the transfer matrix A.

Before introducing additional interface conditions one needs to know where the solid–liquid interface is to be found. The current interface position should be known by the current solution or the initial solution and the new interface location can be determined using the Wilson–Frenkel relation (2.11). More interesting is the question where, i.e. on which row, the interface condition should be applied within the matrix A.

3.5 Interface conditions at the liquid–solid interface 27

some additional terms within the scaled Equations (2.42) until (2.45)) which con-tain a first derivative. These terms which express the temporal stretching of the domain of one phase, has to be discretized as well.

If no Landau–transformation is used the first M grid points may fall within the solid part. The rest of the 2N grid points will then belong to the liquid. The interface conditions have to be imposed on the matrix rows around row number M . During the propagation over time the interface will move and is likely that at some point a grid point moves from one side of the interface to the other side. This in turn changes the number M and the matrix A has to be adapted accordingly. A problem arises when a grid point switches from one side of the interface to the other side. Its old profile value then corresponds still to the old phase. The generation of successive solutions is only possible if at all times all grid points have reasonable profile values assigned to them. In order to guarantee this assignment a linear interpolation (or extrapolation) is preformed to determine a new profile value whenever a grid point crosses the interface.

More complicate interpolation methods are possible. Within the discretization based on the Collocation Method a higher order polynomial interpolation was used. As far as can be seen a more accurate interpolation of the values for the crossing grid points does not really improve the solution. A reason for this is that near the solid–liquid interface a lot of grid points must exist for the description of the solute pileup. For the concentration profile the grid points are thus so near to each other that the linear interpolation is accurate enough. The temperature profile seems to be much less sensitive to inaccurate value approximations. The interface conditions itself can in principle be discretized using the same schemes as already described. The terms of this discretization can then simply be added to the right matrix entries. Actually, there exist three possible variations. Though each of them describes the same model, they lead to somewhat differ-ent results. The three variants are possible with or without application of the Landau–transformation and they are possible for the temperature and concentra-tion field. To explain the differences only the determinaconcentra-tion of the concentraconcentra-tion without Landau–transformation is demonstrated in what follows.

In the first variant the interface condition (2.25) is discretized and the condition (2.29) is directly applied to this discretization which yields

ν(k − Ex)CM = FoDhlCM +2− CM

k 

− FoDhs(CM −CM −1), (3.19) with the specific length chosen equal to the distance between the grid points. The condition (2.29) in this case is

CM +1=C M

k . (3.20)

28 Numerical approximation of solidification

interface and the first point M +1 within the liquid or on the solid–liquid interface are both supposed to lie so near the liquid–solid interface that they can be related to each other by (2.29). A similar equation for the first point M + 1 in he liquid can be written using (2.29) in its original form

CM = kCM +1, (3.21)

which is

ν(k − Ex)CM +1= FoDhl(CM +2−CM +1) − FoDhs(kCM +1−CM −1). (3.22)

With these schemes the interface conditions can directly be included within the transfer matrix A =         ··· ··· ··· −Fo_Dhs 1+2FoDhs −FoDhs −Fo_Dhs co1 0 −FoDhl −FoDhs 0 co2 −FoDhl −FoDhl 1+2FoDhl −FoDhl ··· ··· ···         , co1= 1 + FoDhs+ 1 k FoDhl+ (k − Ex) ν
, co2= 1 + FoDhl+ (k − Ex)ν + kFoDhs. (3.23)

Only the middle part of the whole matrix A, starting with row number M − 1, is explicitly shown. The whole matrix now has dimension 2N × 2N if there exist 2N grid points in total.

3.5 Interface conditions at the liquid–solid interface 29

The resulting transfer matrices of the second and third variant become

A =         ··· ··· ··· −FoDhs 1+2FoDhs −FoDhs −FoDhs co1 −kFoDhs −FoDhs 0 co2 −FoDhl −FoDhl 1+2FoDhl −FoDhl ··· ··· ···         , co1= 1 + 2FoDhs, co2= 1 + FoDhl+ (k − Ex)ν + kFoDhs, (3.24) and A =         ··· ··· ··· −FoDhs 1+2FoDhs −FoDhs −FoDhs co1 0 −FoDhl −1 kFoDhl co2 −FoDhl −Fo_Dhl 1+2FoDhl −FoDhl ··· ··· ···         , co1= 1 + FoDhs+ 1 k FoDhl+ (k − Ex) ν
, co2= 1 + 2FoDhl. (3.25)

Note that the constant co1 is on row M − 1 for the second variant while on row M for the third variant. Within the second and third variant one grid point lying either in the solid or liquid, respectively, will not be part of the solution. The value at this grid point can be determined out of the new solution by applying Relation (3.21). For the determination of the concentration profiles the source vector s of Equation (3.16) remains the same as stated for the first variant. The corresponding discretizations for the determination of the temperature dis-tribution can be retrieved from above matrices used for the solute concentration determination. Thereby following substitutions are needed:

k → 1,

FoDhl → Fohl,

FoDhs → Ex · Sp · Fohs (for ICs).

(3.26)

30 Numerical approximation of solidification

For the source vector s additional terms have to be introduced which lead for the first variant to s =               Fohs·TBC 0 · · · 0 ν · He ν · He 0 · · · 0               . (3.27)

The terms with ν are on row number M and M + 1. For the second and third variant the term at row M and at row M + 1, respectively, have to set to zero.

3.6

How not to include the interface conditions

There exist a few other methods to include the interface conditions within the transfer matrix. Some have been implemented without success. Here only one will be mentioned. It looks very attractive to directly use the position of the liquid–solid interface within the element of grid points M and M + 1. This would in a certain sense allow to eliminate the three variants discussed earlier. The Fourier law can also be applied up to the very last grid point next to the interface. This can be done from both sides and without loosing the validity of the interface condition (2.29).

Figure 3.2: Schematic drawing of the solid and liquid temperature profiles. Between the last grid point within the solid M and the first point within the liquid M + 1, an additional point I is drawn. It corresponds to the current position of the solid–liquid interface.

3.6 How not to include the interface conditions 31

solid–liquid interface is at point I which is between the last point within the solid M and the first point within liquid M +1. The solid temperature profile is linearly extended till grid point M + 1, while the liquid profile starts at point M . In this way both profiles cross each other at the solid–liquid interface. The interface is located a fraction f of the distance between two grid points (hχ) to the right of point M .

The interface condition (2.26) can be rewritten in terms of the values at the grid points M and M + 1 plus the interface value, i.e. at point I,

ν(1 − Ex)TI− νHe = Fohl(TM +1−TI)(1 − f ) − Ex Sp Fohs(TI−TM)f. (3.28) In the solution vector the interface value TI does not occur. This term must therefore be eliminated. The equation can be rewritten as an expression of TI:

TI =

νHe + Fohl(1 − f )TM +1+ Ex Sp FohsfTM ν(1 − Ex) + Fohl(1 − f ) + Ex Sp Fohsf

. (3.29)

When writing the second derivative of the Fourier equation in discrete form, equation (3.28) can be used to get a first derivative in the direction of the interface,

∂2_T ∂χ2   _χ=χ M = FohsTM −1− (1 + f )FohsTM+ f FohsTI, (3.30) and ∂2T ∂χ2   _χ=χ M +1 = (1 − f )FohlTI− (2 − f )FohlTM +1+ FohlTM +2. (3.31) With substitution ofTI by Expression (3.29) the second derivatives become

∂2T ∂χ2   _χ=χ M = ψνf He · Fohs +FohsTM −1 − (1 + f )Fohs− ψEx Sp(f Fohs)2
TM +ψf (1 − f )Fohs· FohlTM +1, (3.32) and ∂2T ∂χ2   _χ=χ M +1 = ψν(1 − f )He · Fohl ψf (1 − f )Ex Sp Fohs· FohlTM − (2 − f )Fohl− ψ(1 − f )2Fohl2
TM +1 +FohlTM +2, (3.33) where ψ = ν(1 − Ex) + (1 − f )Fohl+ f Ex Sp Fohs
−1 . (3.34)

32 Numerical approximation of solidification

velocities ν its advisable, especially for the estimation of the concentration pro-files, to multiply the equations by ψ−1_{. The transfer matrix for the temperature} determination becomes A =         ··· ··· ··· −Fo_hs 1+2Fohs −Fohs −ψ−1_Fo hs co1 co3 co3 co2 −ψ−1Fohl −Fo_hl 1+2Fohl −Fohl ··· ··· ···         , co1= 1 + ψ−1(1 + f ) − f2Ex Sp Fohs
Fohs, co2= 1 + ψ−1(2 − f ) − (1 − f )2Fohl
Fohl, co3= −f (1 − f )Ex Sp FohsFohl. (3.35) The elements co1 and co2 are on row number M and M + 1, respectively, due to the BCs of the temperature profile. As already mentioned the produced latent heat of fusion leads to corresponding terms within the source vector s

s =               Fohs·TBC 0 · · · 0 νf He · Fohs ν(1 − f )He · Fohl 0 · · · 0               . (3.36)

This method looks quite reasonable as it is, up to different material constants in the liquid and solid, build symmetrically around the solid–liquid interface. This symmetry is even retained within the transfer matrix. The transfer matrix for the determination of the concentration profiles poses no difficulty anymore. As the principles remain the same the resulting matrix and source vector are not shown here. Note that in addition relation (2.29) is fulfilled.

Though good looking, the method appeared to have a bad numerical behavior. The reason for this did not became quite apparent. In an ordinary simulation where the results are generated for a certain duration in time and where the interface is moved according to the newly determined interface velocities, the method is prone to form wiggles around the interface. This wiggles often tend to increase during the simulation. Furthermore, they are found not to depend on the Courant–Friedrichs–Lewy (CFL) computational instability.

3.7 The incomplete limit method 33

on the fraction f . This corresponds well with the appearance of wiggles because during an ordinary simulation a near random variation of the fraction f occurs.

3.7

The incomplete limit method

The three variants to include the interface conditions within the transfer matrix A demonstrate that the FDM together with the interface conditions (2.26) and (2.25) leads to quite a bit of ambiguity. The natural way to proceed is to select an excellent example situation, validate all three variants with the example and to select the best variant. This kind of selection, however, does not lead to a method which is guaranteed to describe all physical effects.

In the concrete situation, it is rather tempting to select one of the variants for which the interface concentration remains within bounds. Whether this is the natural method which behaves as expected or whether the method only suppresses numerically the physical effect of a solute pile–up in front of the solid–liquid interface is difficult to decide.

In this section another approach which might lead to some additional insight in the kind of solute pile–up in front of the interface is introduced. Contrary to selecting variants a more physical reasoning can be given. The variants already mentioned are not erroneous. They derive logically from the constitutional equations and the corresponding interface conditions. In order to derive a new method either another numerical solution technique or another form of the fundamental equations of the model has to be applied.

The later application is used here. One of the most general equations for the conservation across the solid–liquid interface serves as a starting point

d dt Z Ω (ρϕ)αdV + Z ∂Ω (ρϕ)α(vα− vi) · ndS + Z ∂Ω jα· ndS = Z Ω ˙ QαdV. (3.37)

Figure 3.3 shows the small flat volume element Ω located at the solid–liquid interface I. This element of thickness  will be moved together with the interface with velocity vi in the direction of the interface normal n. As it is a very general formula, the solid and liquid next to the interface might each have their own velocity vα. The value ϕ can correspond to the mass, enthalpy, unit volume of the solute or even momentum. The fluxes jαat each side of the interface and the source Q within the element are used to guarantee conservation. In this form this equation can be found in the book of Rappaz et al. [46].

34 Numerical approximation of solidification

Figure 3.3: Small volume element of width  located at the solid–liquid interface I. The solid is indicated in a gray shade.

same limit the whole equation (3.37) can be written as

ρlϕl(vl− vi) + jl
n − ρsϕs(vs− vi) + js
n = 0. (3.38) As soon as ϕα and jα are replaced with the right physical quantities, this last equation will together with the assumption vs= 0 lead to the well known interface conditions of the concentration profile (2.25), of the temperature profile (2.26) and of mass conservation (2.1).

The limit of an infinitesimal small volume element, though valid in a physical sense, is when used in a numerical method not always making sense. The reason is that two limits namely the one of an infinitesimal small grid spacing h and one of the element thickness  have to be considered. Both limits however depend somehow on each other. A thickness  much smaller than the grid size h is quite useless because all things bellow the size h cannot be resolved anymore. If the thickness  gets bigger than the distance between grid points the limit is incomplete as the numerical method is in principle capable to account for the volume integrals.

This consideration leads to the central question if it is not possible to take the two limits otherwise than separately. It seems to be possible to take them together but for practical reason another assumption is first made, i.e.

Z Ω ∂(ρϕ)α ∂t dV = Z Ω ˙ QαdV. (3.39)

3.7 The incomplete limit method 35

produced within the element, its only and immediate effect will be a rise of the temperature.

With this assumption and the assumption that the volume Ω does not depend on time t, equation (3.37) becomes

Z Ω ρα ∂ϕα ∂x vidV + Z ∂Ω (ρϕ)α(vα− vi) · ndS + Z ∂Ω jα· ndS = 0. (3.40)

The x–axis is here assumed to be pointing in the same direction as the interface normal n. It can be further assumed that the thickness  of the volume element is so small that the surface integral of tangential component of the flux and velocity, i.e. the part not having the direction of the interface normal n, can be neglected. This then yields to an equation similar to equation (3.38)

Z Ω

ρα ∂ϕα

∂x vidV + ρlϕl(vl− vi) + jl
n − ρsϕs(vs− vi) + js
n = 0. (3.41) The only remaining difficulty is to find a reasonable estimate for the first integral term. Maybe the clearest method is to equal the volume thickness  to the grid point distance h up to a certain factor γ,

 = γ · h. (3.42)

With this last proposition the general conservation equation (3.38) becomes γh

2 ρlϕ 0

l+ ρsϕ0s− ρlϕl− ρsϕs) + ρlϕl0(vl− vi) + jl
n − ρsϕ0s(vs− vi) + js
n ' 0. (3.43) with the understanding that ϕ0_land ϕ0sare the new interface values and ϕland ϕs are already known either from a previous time step or from the initial conditions. Note as well that the volume element is thought to be symmetric around the solid– liquid interface. When the limit of an infinitesimal small grid spacing h is taken, the first term of Equation (3.41) will become zero and one retains Equation (3.38) as should be.

With Equation (3.43) one can proceed as has been described above. For example, the interface condition for the concentration profile (2.25) can be rewritten in the form γ(k+Ex) 2 C 0 l  _χ=I −Cl  _χ=I  + ν(k − Ex)C0 l  _χ=I = FoDhl _∂C0 l ∂χ   _χ=I+  − FoDhs _∂C0 s ∂χ   _χ=I−  , (3.44)

From this Equation (3.44) it can be seen that this method leads to nearly the same interface condition as the other methods introduced in previous sections. Only the first term is new. Its value depends on the factor γ which can be chosen adequately.

Chapter 4

Numerical stability and

accuracy of the solidification

model

4.1

Introduction

As numerical solutions of partial differential equations are approximations of the exact analytical solutions, numerical errors are unavoidably. These errors indicate all deviations of the results from the exact values. There are three common sources for the deviations, i.e. intrinsic errors, inadequate choices of the model parameters values and bad relations between model parameters.

An intrinsic error is an error directly related to the kind of method in considera-tion. Usually no sensible choice of the model parameters allows to diminish the deviations. Sometimes, this error is small or occurs only for a certain range of parameter values. In extreme situations this error can make a method useless. This chapter is not concerned about this kind of errors. In Chapter 3, however, it was stated how an interface condition should not be included. Failure of the resulting method, when the advice is disregarded, is most likely due to an intrinsic error.

Quite another situation arises if the model parameters are chosen inadequately. In most numerical methods one or more parameters can be found that influence the accuracy of the approximation without changing the physical system. For example, in the FDM the results can often be improved by lowering the spacing of the mesh and increasing the duration of the time integration. Theoretically, the error connected with the model parameters can be lowered infinitely with