**1. Introduction**

20 Numerical Simulations / Book 1

512 Computational Simulations and Applications

[34] V. P. Solovjov, B. W. Webb, "SLW Modeling of Radiative Transfer in Multicomponent Gas Mixtures," *Journal of Quantitative Spectroscopy & Radiative Transfer*, vol. 65, pp. 655–672,

[35] C. P. Thurgood, A. Pollard, A. B. Becker, "The TN Quadrature Set for the Discrete

[36] L. S. Rothman, R. R. Gamache, R. H. Tipping, C. P. Rinsland, M. A. H. Smith, D. C. Benner, V. M. Devi, J. M. Flaud, C. Camy-Peyret, A. Perrin, A. Goldman, S. T. Massie, L. R. Brown, "The HITRAN Molecular Database: Editions of 1991 and 1992," *Journal of*

[37] J. M. Hartmann, R. Levi Di Leon, J. Taine, "Line-by-line and Narrow-band Statistical Model Calculations for H2O," *Journal of Quantitative Spectroscopy & Radiative Transfer*, vol.

Ordinates Method," *Journal of Heat Transfer*, vol. 117, pp. 1068–1070, 1995.

*Quantitative Spectroscopy & Radiative Transfer*, vol. 48, pp. 469–507, 1992.

2000.

32, pp. 119–127, 1984.

This chapter is devoted to the problems connected with numerical modeling of moving boundary ones. In particular the solidification and cooling processes proceeding in the system casting-mould are considered. The subject matter of solidification process modeling is very extensive and only the selected problems from this scope will be here discussed. From the mathematical point of view the thermal processes proceeding in the domain considered (both in macro and micro/macro scale) are described by a system of partial differential equations (energy equations) supplemented by the geometrical, physical, boundary and initial conditions. The typical solidification model bases on the Fourier-Kirchhoff type equations, but one can formulate the more complex (coupled with the basic model) ones concerning the heat convection in a molten metal sub-domain, the changes of local chemical constitution of solidifying alloy (segregation process) etc. Here we limit oneself only to the tasks connected with the predominant heat conduction process, and (as the examples of various problems) the macro models of alloys solidification and the micro/macro models of pure metals crystallization will be presented, at the same time the direct and inverse problems will be discussed.

In order to construct a numerical model and an adequate computer program simulating the course of problem considered, one must accept a certain mathematical description of the process (the governing equations). The next step is the transformation of this mathematical model into a form called the re-solving system constructed on the basis of a selected numerical method.

After transformation of the algorithm developed into a computer program and supplementing it with suitable pre- and post-processing procedures (input data loading, graphic presentation of results, print-outs etc.), and carrying out of computations, one obtains the results including information concerning the transient temperature field, kinetics of solidification process, the temporary shapes of sub-domains etc. They may have a form of numerical print-outs, e.g., giving the temperature field at distinguished set of points or the volumetric fraction of solid state at the neighbourhood of these points.

Numerical modeling of the heat and mass transfer in solidifying metal is a typical interdisciplinary problem and requires particular knowledge in the field of foundry practice, mathematics (a course in mathematical analysis offered in technical schools is quite sufficient here), thermodynamics (in particular heat transfer), numerical methods and, in the final stage, also programming and operation of computer equipment.

Numerical Modeling of Solidification Process 515

domain. In the case of alloys solidification the segregation effects can be also taken into

In this sub-chapter the model called the one domain approach will be discussed (Mochnacki & Suchy, 1995). Let us denote the temperatures corresponding to the beginning and the end of solidification process as *TL* and *TS* . If one considers the solidification of pure metals or eutectic alloys then the border temperatures can be introduced in an artificial way substituting the solidification point *T*\* by a certain temperature interval [*T*\* Δ*T*, *T*\* + Δ*T*]. Numerical experiments show that the assumed values of Δ*T* (in reasonable limits) are not very essential, and the results of numerical

One assumes the knowledge of temperature-dependent function *fS* for interval [*TS*, *TL* ] and

(,) d (,) d *S S f x t f Txt t Tt*

( ,) ( ) ( ) ( ,) *Txt C T T Txt*

where *C*(*T* ) = *c*(*T* ) *L* df*S*/d*T* is called 'a substitute thermal capacity'. This parameter can be

*S P SL*

where *cL*, *cP*, *cS* are the volumetric specific heats of molten metal, mushy zone and solid state sub-domains, one can use the equation (9) as the model of thermal processes proceeding in the whole, conventionally homogeneous, casting domain. It is the reason that the approach

The function *fS* must fulfill the conditions *fS*(*TL* ) = 0 and *fS*(*TS* ) = 1, additionally for

( )d ( )

corresponds to the change of volumetric physical enthalpy for [*TS*,*TL*] gives the additional

*PL S*

*CT T c T T L* (10)

*<sup>f</sup> CT c L T T T T*

*L L*

*c T T*

 

*S S*

*c T T*

*t*

<sup>d</sup> ( ) <sup>d</sup>

*L*

*T*

*S*

*T*

Introducing the following definition of substitute thermal capacity

(7)

(8)

(9)

account, etc.

then

**2.2 Macro models of solidification** 

simulations of solidification process are similar.

Introducing this formula to equation (1) one obtains

defined in the different ways.

presented is called 'a one domain method'.

information assuring the proper choice of *fS* .

As an example the following function can be considered

*T* < *TS*: *fS* = 0, for *T* > *TL*: *fS* = 0.

The integral

### **2. Heat transfer in the system casting – mould**

#### **2.1 Description of the process**

At first, we will consider the heat transfer processes in a system casting-mould-environment The solidification and cooling of casting material (alloy or pure metal) can be described by the following energy equation

$$\mathcal{L}(T)\frac{\partial \mathcal{T}(\mathbf{x},t)}{\partial t} = \nabla \left[\mathcal{L}(T)\nabla T(\mathbf{x},t)\right] + L\frac{\partial f\_S(\mathbf{x},t)}{\partial t} \tag{1}$$

where *c*(*T*) is a volumetric specific heat of casting material, λ(*T*) is a thermal conductivity, *L* is a volumetric latent heat, *T* = *T*(*x, t*), *fS* = *fS*(*x, t*) denote the temperature and the local volumetric fraction of solid state, *x* denotes the spatial co-ordinates, *t* is a time. One can see, that only heat conduction in a casting volume is considered. The energy equation can constitute a base both in a case of macro-scale modeling and micro/macro one (Stefanescu, 1999; Mochnacki & Majchrzak, 2007a; Mochnacki & Suchy, 1995). The differences appear on a stage of solidification rate *fS*/*t* computations.

The equation determining the course of thermal processes in a mould sub-domain is of the form

$$\mathcal{L}\_m(T) \frac{\partial T\_m(\mathbf{x}, t)}{\partial t} = \nabla \left[ \lambda\_m(T) \,\nabla T\_m(\mathbf{x}, t) \right] \tag{2}$$

where the index *m* identifies the mould sub-domain, the non-homogeneous mould can be also considered.

On the external surface of mould the continuity condition (3rd type of boundary condition)

$$-\lambda\_m \frac{\partial T\_m(\mathbf{x}, t)}{\partial \mathbf{n}} = \mathbf{a} \left[ T\_m(\mathbf{x}, t) - T\_a \right] \tag{3}$$

is, as a rule, accepted. Here α is a heat transfer coefficient, *Ta* is an ambient temperature, /*n* denotes a normal derivative.

On the contact surface between casting and mould the continuity condition is given

$$-\lambda \frac{\partial T(\mathbf{x},t)}{\partial \mathbf{n}} = \frac{T(\mathbf{x},t) - T\_m(\mathbf{x},t)}{R(\mathbf{x},t)} = -\lambda\_m \frac{\partial T\_m(\mathbf{x},t)}{\partial \mathbf{n}}\tag{4}$$

where *R* is a thermal resistance. For *R* = 0 (a such assumption can be done in the case of sand mix mould) the last equation takes a form

$$\begin{cases} -\lambda \frac{\partial \mathcal{T}(\mathbf{x}, t)}{\partial n} = -\lambda\_m \frac{\partial \mathcal{T}\_m(\mathbf{x}, t)}{\partial n} \\ T(\mathbf{x}, t) = T\_m(\mathbf{x}, t) \end{cases} \tag{5}$$

The initial temperature distribution for *t* = 0 is also known

$$\mathbf{x}, t = \mathbf{0}: \quad T(\mathbf{x}, \mathbf{0}) = T\_0(\mathbf{x}) \; , \quad T\_m(\mathbf{x}, \mathbf{0}) = T\_{m0}(\mathbf{x}) \tag{6}$$

The mathematical model presented above can be more complicated. One can consider the convectional component of heat transfer which appear in the molten metal sub-

At first, we will consider the heat transfer processes in a system casting-mould-environment The solidification and cooling of casting material (alloy or pure metal) can be described by

> (,) (,) ( ) () (,) *Txt Sf xt c T T Txt L t t*

where *c*(*T*) is a volumetric specific heat of casting material, λ(*T*) is a thermal conductivity, *L* is a volumetric latent heat, *T* = *T*(*x, t*), *fS* = *fS*(*x, t*) denote the temperature and the local volumetric fraction of solid state, *x* denotes the spatial co-ordinates, *t* is a time. One can see, that only heat conduction in a casting volume is considered. The energy equation can constitute a base both in a case of macro-scale modeling and micro/macro one (Stefanescu, 1999; Mochnacki & Majchrzak, 2007a; Mochnacki & Suchy, 1995). The differences appear on

The equation determining the course of thermal processes in a mould sub-domain is of the

(,) ( ) <sup>λ</sup> () (,) *<sup>m</sup> <sup>m</sup> m m T xt c T T T xt*

where the index *m* identifies the mould sub-domain, the non-homogeneous mould can be

On the external surface of mould the continuity condition (3rd type of boundary condition)

 (,) λ α (,) *<sup>m</sup> m ma T xt T xt T*

is, as a rule, accepted. Here α is a heat transfer coefficient, *Ta* is an ambient temperature,

(,) (,) (,) (,) λ λ (,)

where *R* is a thermal resistance. For *R* = 0 (a such assumption can be done in the case of sand

(,) (,) λ λ

The mathematical model presented above can be more complicated. One can consider the convectional component of heat transfer which appear in the molten metal sub-

*Txt T xt n n*

*<sup>m</sup> <sup>m</sup>*

*Txt Txt T xt T xt n Rx t n*

*m m <sup>m</sup>*

(4)

0 0 *t Tx T x T x T x* 0 : ( , 0) ( ) , ( , 0) ( ) *m m* (6)

*t*

*n*

On the contact surface between casting and mould the continuity condition is given

(,) (,)

 

*Txt T xt*

The initial temperature distribution for *t* = 0 is also known

*m*

(1)

(2)

(3)

(5)

**2. Heat transfer in the system casting – mould** 

a stage of solidification rate *fS*/*t* computations.

**2.1 Description of the process** 

the following energy equation

form

also considered.

/*n* denotes a normal derivative.

mix mould) the last equation takes a form

domain. In the case of alloys solidification the segregation effects can be also taken into account, etc.

#### **2.2 Macro models of solidification**

In this sub-chapter the model called the one domain approach will be discussed (Mochnacki & Suchy, 1995). Let us denote the temperatures corresponding to the beginning and the end of solidification process as *TL* and *TS* . If one considers the solidification of pure metals or eutectic alloys then the border temperatures can be introduced in an artificial way substituting the solidification point *T*\* by a certain temperature interval [*T*\* Δ*T*, *T*\* + Δ*T*]. Numerical experiments show that the assumed values of Δ*T* (in reasonable limits) are not very essential, and the results of numerical simulations of solidification process are similar.

One assumes the knowledge of temperature-dependent function *fS* for interval [*TS*, *TL* ] and then

$$\frac{\partial f\_S(\mathbf{x},t)}{\partial t} = \frac{\mathbf{d} \, f\_S}{\mathbf{d}T} \frac{\partial T(\mathbf{x},t)}{\partial t} \tag{7}$$

Introducing this formula to equation (1) one obtains

$$\mathcal{C}(T)\frac{\partial \mathcal{T}(\mathbf{x},t)}{\partial t} = \nabla \left[\mathcal{A}(T)\nabla T(\mathbf{x},t)\right] \tag{8}$$

where *C*(*T* ) = *c*(*T* ) *L* df*S*/d*T* is called 'a substitute thermal capacity'. This parameter can be defined in the different ways.

Introducing the following definition of substitute thermal capacity

$$\mathbf{C}(T) = \begin{cases} c\_L & T > T\_L \\ c\_P - L \frac{\mathbf{d} \cdot f\_S}{\mathbf{d}T} & T\_S \le T \le T\_L \\ c\_S & T < T\_S \end{cases} \tag{9}$$

where *cL*, *cP*, *cS* are the volumetric specific heats of molten metal, mushy zone and solid state sub-domains, one can use the equation (9) as the model of thermal processes proceeding in the whole, conventionally homogeneous, casting domain. It is the reason that the approach presented is called 'a one domain method'.

The function *fS* must fulfill the conditions *fS*(*TL* ) = 0 and *fS*(*TS* ) = 1, additionally for *T* < *TS*: *fS* = 0, for *T* > *TL*: *fS* = 0. The integral

$$\int\_{T\_S}^{T\_\perp} \mathbf{C} \{ T \} \, \mathrm{d}T = c\_P (T\_L - T\_S) + L \tag{10}$$

corresponds to the change of volumetric physical enthalpy for [*TS*,*TL*] gives the additional information assuring the proper choice of *fS* .

As an example the following function can be considered

Numerical Modeling of Solidification Process 517

*r*

 ( ,) ( ) ( ) ( ,), ( ) ( ) *Hxt <sup>T</sup> aT H xt aT t C T*

In Figure 2 the course of enthalpy function for carbon steel (0.44%C) under the assumption that the substitute thermal capacity *C* (*T*) is approximated by piece-wise constant function

The transformation of typical boundary and initial conditions supplementing the equation

The mathematical model of solidification basing on the enthalpy function is often used at the stage of numerical modeling. In particular, the approach discussed seems to be very effective in a case of more complex courses of phase changes (e.g. solidification of cast iron, the process proceeds partially at the constant temperature and partially at the interval of temperature and then the introduction of enthalpy function and the application of numerical procedure called the Alternating Phase Truncation Method (Mochnacki & Suchy, 1995; Majchrzak & Mochnacki, 1995; Majchrzak & Mendakiewicz, 1993) leads to the simple

The other form of energy equation determining the thermal processes in domain of

*r*

 

(19)

*T*

*T*

solidifying alloy can be obtained using the Kirchhoff transformation, this means

*U T*

(18) can be found, among others, in (Mochnacki & Suchy, 1995).

( ) ( )d

 

(17)

(18)

*T HT C*

where the lower limit of integration corresponds to the optional reference level.

*T*

**2.3 Introduction of enthalpy function and Kirchhoff transformation**  Let us introduce the physical enthalpy of alloy defined in the following way

The energy equation (8) written using the function (17) takes a form

( ) ( )d

(equation (16)) is shown.

Fig. 2. Enthalpy function.

and effective numerical algorithm).

$$\mathcal{L}\_S(T) = \left(\frac{T\_L - T}{T\_L - T\_S}\right)^n \tag{11}$$

and then

$$\frac{\mathbf{d}\,f\_S(T)}{\mathbf{d}\,T} = -\frac{n}{T\_L - T\_S} \left(\frac{T\_L - T}{T\_L - T\_S}\right)^{n-1} \tag{12}$$

Finally

$$\text{CC}(T) = c\_P + \frac{L}{T\_L - T\_S} \text{ m} \left( \frac{T\_L - T}{T\_L - T\_S} \right)^{n-1} \tag{13}$$

The quotient *L*/(*TL TS* ) = *csp* is called 'a spectral latent heat'. Introducing this parameter one has

$$\mathcal{C}(T) = c\_P + c\_{sp} \quad n \left( \frac{T\_L - T}{T\_L - T\_S} \right)^{n-1} \tag{14}$$

It is easy to check that

$$\int\_{T\_S}^{T\_L} \left[ c\_P + c\_{sp} \cdot n \left( \frac{T\_L - T}{T\_L - T\_S} \right)^{n-1} \right] dT = c\_P (T\_L - T\_S) + L \tag{15}$$

Above formulas are very often used for exponent *n* = 1 and then

$$\text{CC}\{T\} = \mathcal{c}\_P + \frac{L}{T\_L - T\_S} = \mathcal{c}\_P + \mathcal{c}\_{sp} \quad \text{ } T \in \left[T\_S, T\_L\right] \tag{16}$$

So, the linear course of *fS* leads to the constant value of mushy zone thermal capacity. Assuming additionally the constant values of *cS* and *cL* one obtains the staircase function as an approximation of *C* (*T* ). In Figure 1 this type of function is shown. If the changes of border temperatures are taken into account (macrosegregation effect) then the families of functions *C*(*T* ) should be considered (Szopa, 1999).

Fig. 1. Substitute thermal capacity

*L*

*T T f T T T* 

*S L*

*L TT CT c n*

*P sp*

*L*

*T T*

*T T*

*L S*

*T T*

*T T CT c c n*

*f T n TT T TTTT* 

*L S*

*LSLS*

*LS LS*

*L*

1

*c cn T cT T L*

*P sp PL S*

( ) , [,] *<sup>P</sup> P sp S L*

So, the linear course of *fS* leads to the constant value of mushy zone thermal capacity. Assuming additionally the constant values of *cS* and *cL* one obtains the staircase function as an approximation of *C* (*T* ). In Figure 1 this type of function is shown. If the changes of border temperatures are taken into account (macrosegregation effect) then the families of

*<sup>L</sup> CT c c c T TT*

*L S*

*T T*

*TT TT* 

The quotient *L*/(*TL TS* ) = *csp* is called 'a spectral latent heat'. Introducing this parameter

*L*

*n*

1

1

*n*

1

*n*

d( )

(15)

(16)

*n*

(11)

(12)

(13)

(14)

( )

*S*

d () d

( )

( )

*T n*

*T L S*

Above formulas are very often used for exponent *n* = 1 and then

*L*

*S*

functions *C*(*T* ) should be considered (Szopa, 1999).

Fig. 1. Substitute thermal capacity

*P*

and then

Finally

one has

It is easy to check that

#### **2.3 Introduction of enthalpy function and Kirchhoff transformation**

Let us introduce the physical enthalpy of alloy defined in the following way

$$H(T) = \bigcap\_{T\_r}^{T} \mathcal{C}(\mu) \text{ d}\mu \tag{17}$$

where the lower limit of integration corresponds to the optional reference level. The energy equation (8) written using the function (17) takes a form

$$\frac{\partial \mathcal{H}(\mathbf{x},t)}{\partial t} = \nabla \left[ a(T) \nabla H(\mathbf{x},t) \right], \qquad a(T) = \frac{\mathcal{A}(T)}{\mathcal{C}(T)} \tag{18}$$

In Figure 2 the course of enthalpy function for carbon steel (0.44%C) under the assumption that the substitute thermal capacity *C* (*T*) is approximated by piece-wise constant function (equation (16)) is shown.

Fig. 2. Enthalpy function.

The transformation of typical boundary and initial conditions supplementing the equation (18) can be found, among others, in (Mochnacki & Suchy, 1995).

The mathematical model of solidification basing on the enthalpy function is often used at the stage of numerical modeling. In particular, the approach discussed seems to be very effective in a case of more complex courses of phase changes (e.g. solidification of cast iron, the process proceeds partially at the constant temperature and partially at the interval of temperature and then the introduction of enthalpy function and the application of numerical procedure called the Alternating Phase Truncation Method (Mochnacki & Suchy, 1995; Majchrzak & Mochnacki, 1995; Majchrzak & Mendakiewicz, 1993) leads to the simple and effective numerical algorithm).

The other form of energy equation determining the thermal processes in domain of solidifying alloy can be obtained using the Kirchhoff transformation, this means

$$\mathcal{U}(T) = \bigcap\_{T\_r}^{T} \mathcal{A}(\mu) \text{ d}\mu \tag{19}$$

Numerical Modeling of Solidification Process 519

where *u* = *R*/*t* is a crystallization rate (*R* is a grain radius), *N* is nuclei density, *v* is a coefficient which equals 1 in case of spherical grains or *v* < 1 (e.g. dendritic growth, as in

> ( , ) ( ,) *Sf x t xt*

( , ) ( ,) *Sf x t xt t t*

The model above discussed determines the geometrical volume (volume fraction) and it is

In order to take into account the limitations of growth in the final stages of the process one

( , ) ( ,) 1 (,) *<sup>S</sup>*

*f x t xt <sup>f</sup> x t*

d () <sup>d</sup> 1 () *S S f f* 

*t t*

*S*

(23)

(24)

(25)

(26)

Figure 3 (Majchrzak & Piasecka, 1997)).

Fig. 3. Dendritic growth.

Then

In the case of so-called linear model we have

and if *fS* = 1 then the crystallization process stops. The derivative of *fS* with respect to time equals

the correct assumption on the first stages of crystallization.

introduces the following modification of equation (24)

It is easy to check that the introduction of function (19) Leeds to the following form of equation (8)

$$\phi\left[T(\mathsf{U})\right]\frac{\partial \mathsf{U}\mathcal{U}(\mathsf{x},t)}{\partial t} = \nabla^2 \mathsf{U}(\mathsf{x},t) \; \; \; \; \qquad \phi\left(\mathsf{U}\right) = \frac{\mathsf{C}\left[T(\mathsf{U})\right]}{\lambda\left[T(\mathsf{U})\right]}\tag{20}$$

One can see, that the right hand side of energy equation (20) becomes a linear one and it can be essential at the stage of numerical modeling. The details concerning this approach and also the transformed form of typical boundary and initial conditions can be found in (Mochnacki & Suchy, 1995; Szopa, 1999).

#### **3. Micro/macro models of solidification**

A generalized form of solidification model belonging to the second generation ones (Stefanescu, 1999) will be here discussed. The equation describing the thermal processes proceeding in domain of solidifying metal contains the term (source function) controlling the course of latent heat evolution. A capacity of internal heat sources is proportional to crystallization rate, more precisely, to time derivative of function *fS*(*x*, *t*) corresponding to the local and temporary volumetric fraction of solid state at the point considered (Majchrzak et al., 2006; Mochnacki & Szopa, 2007; Fraś et al., 1993; Lupa et al. 2004) see: equation (1). If the micro/macro approach is taken into account, then the changes of *fS* result from the crystallization laws in a micro scale (nucleation and nuclei growth). In this way the macro model basing on the Kirchhoff-Fourier-type equation is coupled with the model describing the processes proceeding on the level of single grains.

In literature one can find two basic models determining the mutual connections between *fS* and temporary volume of grains. In particular, the following function is introduced

$$
\rho(\mathbf{x},t) = \mathcal{N}(\mathbf{x},t) \, V(\mathbf{x},t) \tag{21}
$$

where *N* is a nuclei density [nuclei/m3], *V* is a single grain volume. If one assumes *fS*(*x*, *t*) = ω(*x*, *t*) then the linear model is considered (Majchrzak et al., 2006), while if *fS*(*x*, *t*) = 1 − exp[−ω(*x*, *t*)] then the exponential one is taken into account (Fraś et al., 1993; Lupa et al. 2004). In this place one can see that for the small values of ω, both models lead to the same results because exp(−ω) 1−ω and 1 − exp(−ω) ω. During the final stages of solidification the models discussed give the different results. The assumption of exponential changes of ω allows, in a certain way, to take into account the mutual geometrical interactions between the grains.

In this sub-chapter the generalization of the models previously presented will be shown. The differential equation from which linear and exponential models result, will be modified and in this way the new possibilities of *fS* definition will appear.

As was mentioned, the micro/macro models of solidification discussed below require the introduction of function being the product of nuclei density and single grain volume. In a practical realization the exponent ω (equation (21)) is determined by the following formula (Mochnacki & Szopa, 2007; Fraś et al., 1993)

$$\rho(\mathbf{x},t) = \frac{4}{3}\pi \,\mathrm{v} \,\mathrm{N}(\mathbf{x},t) \left[\int\_0^t \mu(\tau) \,\mathrm{d}\tau\right]^3 \tag{22}$$

It is easy to check that the introduction of function (19) Leeds to the following form of

<sup>2</sup> ( ,) [ ( )] [ ( )] ( , ), ( ) [ ( )] *Uxt C TU T U Uxt U*

One can see, that the right hand side of energy equation (20) becomes a linear one and it can be essential at the stage of numerical modeling. The details concerning this approach and also the transformed form of typical boundary and initial conditions can be found in

A generalized form of solidification model belonging to the second generation ones (Stefanescu, 1999) will be here discussed. The equation describing the thermal processes proceeding in domain of solidifying metal contains the term (source function) controlling the course of latent heat evolution. A capacity of internal heat sources is proportional to crystallization rate, more precisely, to time derivative of function *fS*(*x*, *t*) corresponding to the local and temporary volumetric fraction of solid state at the point considered (Majchrzak et al., 2006; Mochnacki & Szopa, 2007; Fraś et al., 1993; Lupa et al. 2004) see: equation (1). If the micro/macro approach is taken into account, then the changes of *fS* result from the crystallization laws in a micro scale (nucleation and nuclei growth). In this way the macro model basing on the Kirchhoff-Fourier-type equation is coupled with the model describing

In literature one can find two basic models determining the mutual connections between *fS*

where *N* is a nuclei density [nuclei/m3], *V* is a single grain volume. If one assumes *fS*(*x*, *t*) = ω(*x*, *t*) then the linear model is considered (Majchrzak et al., 2006), while if *fS*(*x*, *t*) = 1 − exp[−ω(*x*, *t*)] then the exponential one is taken into account (Fraś et al., 1993; Lupa et al. 2004). In this place one can see that for the small values of ω, both models lead to the same results because exp(−ω) 1−ω and 1 − exp(−ω) ω. During the final stages of solidification the models discussed give the different results. The assumption of exponential changes of ω allows, in a certain way, to take into account the mutual geometrical

In this sub-chapter the generalization of the models previously presented will be shown. The differential equation from which linear and exponential models result, will be modified

As was mentioned, the micro/macro models of solidification discussed below require the introduction of function being the product of nuclei density and single grain volume. In a practical realization the exponent ω (equation (21)) is determined by the following

<sup>4</sup> ( , ) ( , ) ( )d <sup>3</sup>

 *x t vN x t u* 

0

 

*t*

and temporary volume of grains. In particular, the following function is introduced

and in this way the new possibilities of *fS* definition will appear.

formula (Mochnacki & Szopa, 2007; Fraś et al., 1993)

*t T U*

 

( ,) ( ,) ( ,) *xt Nxt Vxt* (21)

3

(22)

(20)

equation (8)

(Mochnacki & Suchy, 1995; Szopa, 1999).

**3. Micro/macro models of solidification** 

the processes proceeding on the level of single grains.

interactions between the grains.

where *u* = *R*/*t* is a crystallization rate (*R* is a grain radius), *N* is nuclei density, *v* is a coefficient which equals 1 in case of spherical grains or *v* < 1 (e.g. dendritic growth, as in Figure 3 (Majchrzak & Piasecka, 1997)).

Fig. 3. Dendritic growth.

In the case of so-called linear model we have

$$f\_S(\mathbf{x}, t) = \alpha(\mathbf{x}, t) \tag{23}$$

and if *fS* = 1 then the crystallization process stops. The derivative of *fS* with respect to time equals

$$\frac{\partial \, f\_S(\mathbf{x}, t)}{\partial t} = \frac{\partial \, oo(\mathbf{x}, t)}{\partial t} \tag{24}$$

The model above discussed determines the geometrical volume (volume fraction) and it is the correct assumption on the first stages of crystallization.

In order to take into account the limitations of growth in the final stages of the process one introduces the following modification of equation (24)

$$\frac{\partial f\_S(\mathbf{x},t)}{\partial t} = \frac{\partial \, \alpha(\mathbf{x},t)}{\partial t} [1 - f\_S(\mathbf{x},t)] \tag{25}$$

Then

$$\frac{\text{d}\,f\_S(o)}{1 - f\_S(o)} = \text{d}\, o\,\tag{26}$$

Numerical Modeling of Solidification Process 521

second power of undercooling below the solidification point *T* \* (Mochnacki & Szopa, 2007;

<sup>2</sup> 2 \* *Nxt Txt T Txt* ( ,) ( ,) ( ,)

where η is the nucleation coefficient. The nucleation stops when Δ*T*(*x*, *t* + Δ*t* ) < Δ*T*(*x*, *t*),

d ( ,) ( ,) <sup>d</sup> *Rxt <sup>m</sup> Txt*

 

where μ is the growth coefficient, *m* [1, 2]. In literature one can also find the following

d ( ,) ( ,) ( ,) ( ,) <sup>d</sup> *Rxt uxt Txt Txt*

Because the problem considered is non-steady one, so one introduces the time grid defined

Let us consider the control volume Δ*Vi* from casting domain. During a certain interval of time the temperature at central point of Δ*Vi* decreases below the solidification point and the

 

1 2

2 3

 

011 1 <sup>1</sup> 0 ... ... , *s s S ss t t t t t t tt t* (37)

*t*

*t*

Now, the certain aspects of function ω modelling will be discussed.

while the metal domain is divided into *m* control volumes .

We find a number of the first 'portion' of nuclei *Ni*

 

(34)

(35)

(1) (using equation(34)) and the final

(36)

Fraś et al., 1993; Lupa et al. 2004)

Fig. 4. Undercooling below point *T*\*

where μ1, μ2 are the growth coefficients.

crystallization process starts.

radius of grains (formula (35)). The first value of ω is equal to

equation

as follows

The solid phase growth (equiaxial grains) is determined by

additionally for *T*(*x*, *t*) > *T\**: *N*(*x*, *t*) = 0 (Figure 4).

next

$$f\_S(a) = 1 + C \exp(-a) \tag{27}$$

Because for ω = 0: *fS* = 0 therefore *C* = 1 and finally

$$f\_S(o) = 1 - \exp(-o) \tag{28}$$

The last equation corresponds to the well known exponential model (the Kolmogoroff formula).

Here, the following modification of equation (24) can be proposed (Mochnacki & Szopa, 2007; Mochnacki & Szopa, 2010)

$$\frac{\partial f\_S(\mathbf{x},t)}{\partial t} = \frac{\partial \, o(\mathbf{x},t)}{\partial t} [1 - f\_S(\mathbf{x},t)]^n \tag{29}$$

where *n*  0. From the physical point of view the same conditions as in previous case are fulfilled and the component 1 *fS* changes from 0 to 1. So, we have

$$\frac{\mathbf{d}\,f\_S(o)}{\left[1-f\_S(o)\right]^\eta} = \mathbf{d}\,o\,o\tag{30}$$

It is easy to check that the solution fulfilling the condition ω = 0: *fS* = 0 is of the form

$$f\_S(o) = 1 - \left[ (n-1)o + 1 \right]^{\frac{1}{1-n}} \tag{31}$$

One can see that the last power-type formula constitutes the generalization of linear and exponential models. For *n* = 0 one obtains the linear one, while for *n* = 1 one has

$$\lim\_{n \to 1} 1 - [(n-1)\,\alpha + 1]^{\frac{1}{1-n}} = 1 - \exp(-\alpha) \tag{32}$$

The others values of *n* can be also introduced. For example

$$\begin{aligned} n=2: \quad f\_S &= 1 - \frac{1}{\alpha+1} \\ n=3: \quad f\_S &= 1 - \frac{1}{\sqrt{2\alpha+1}} \\ n=\frac{1}{2}: \quad f\_S &= 1 - \left(1 - \frac{\alpha}{2}\right)^2 \end{aligned} \tag{33}$$

It should be pointed out that if *n*  1 then for ω : *fS* 1. This property is not fulfilled for *n* < 1 and it must be taken into account on the stage of numerical simulation. For example if *n* = 1/2 then solidification process takes place only for ω 2, the values ω > 2 are not physically correct.

Now, the way of function ω modeling will be discussed.

A driving force of nucleation and nuclei growth is an undercooling below solidification point *T* \*. We assume that a local and temporary number of nuclei is proportional to the

( ) 1 exp( ) *Sf C*

( ) 1 exp( ) *Sf*

The last equation corresponds to the well known exponential model (the Kolmogoroff

Here, the following modification of equation (24) can be proposed (Mochnacki & Szopa,

( , ) ( ,) 1 (,) *<sup>n</sup> <sup>S</sup>*

*f x t xt <sup>f</sup> x t*

where *n*  0. From the physical point of view the same conditions as in previous case are

*n*

1 1

ω 1

2ω 1

 

 d () <sup>d</sup> 1 () *S*

It is easy to check that the solution fulfilling the condition ω = 0: *fS* = 0 is of the form

*n*

*n f*

*n f*

*n f*

1

The others values of *n* can be also introduced. For example

*n*

Now, the way of function ω modeling will be discussed.

exponential models. For *n* = 0 one obtains the linear one, while for *n* = 1 one has

*f f* 

*S*

( ) 1 ( 1) 1 <sup>1</sup> *<sup>n</sup> Sf n*

One can see that the last power-type formula constitutes the generalization of linear and

lim 1 ( 1) 1 1 exp( ) *<sup>n</sup>*

<sup>1</sup> 2: 1

*S*

*S*

*S*

<sup>1</sup> 3: 1

<sup>1</sup> <sup>ω</sup> : 11 2 2

It should be pointed out that if *n*  1 then for ω : *fS* 1. This property is not fulfilled for *n* < 1 and it must be taken into account on the stage of numerical simulation. For example if *n* = 1/2 then solidification process takes place only for ω 2, the values ω > 2 are not

A driving force of nucleation and nuclei growth is an undercooling below solidification point *T* \*. We assume that a local and temporary number of nuclei is proportional to the

*S*

1

2

(31)

(32)

(29)

(27)

(28)

(30)

(33)

*t t*

Because for ω = 0: *fS* = 0 therefore *C* = 1 and finally

fulfilled and the component 1 *fS* changes from 0 to 1.

next

formula).

So, we have

physically correct.

2007; Mochnacki & Szopa, 2010)

second power of undercooling below the solidification point *T* \* (Mochnacki & Szopa, 2007; Fraś et al., 1993; Lupa et al. 2004)

$$N(\mathbf{x}, t) = \eta \,\,\Delta T(\mathbf{x}, t) \,\,^2 = \eta \left[\,\,T\,^\* - T(\mathbf{x}, t)\right]^2\tag{34}$$

where η is the nucleation coefficient. The nucleation stops when Δ*T*(*x*, *t* + Δ*t* ) < Δ*T*(*x*, *t*), additionally for *T*(*x*, *t*) > *T\**: *N*(*x*, *t*) = 0 (Figure 4).

Fig. 4. Undercooling below point *T*\*

The solid phase growth (equiaxial grains) is determined by

$$\frac{\text{d}\,R(\mathbf{x},t)}{\text{d}\,t} = \mu \,\Delta T(\mathbf{x},t)^m\tag{35}$$

where μ is the growth coefficient, *m* [1, 2]. In literature one can also find the following equation

$$
\Delta \mathbf{u}(\mathbf{x},t) = \frac{\mathbf{d} \, \mathbf{R}(\mathbf{x},t)}{\mathbf{d}t} = \mu\_1 \, \Delta \mathbf{T}(\mathbf{x},t)^2 + \mu\_2 \, \Delta \mathbf{T}(\mathbf{x},t)^3 \tag{36}
$$

where μ1, μ2 are the growth coefficients.

Now, the certain aspects of function ω modelling will be discussed.

Because the problem considered is non-steady one, so one introduces the time grid defined as follows

$$0 = t^0 < t^1 < t^1 < \dots < t^{s-1} < t^s < \dots < t^S \quad , \quad \Delta t = t^s - t^{s-1} \tag{37}$$

while the metal domain is divided into *m* control volumes .

Let us consider the control volume Δ*Vi* from casting domain. During a certain interval of time the temperature at central point of Δ*Vi* decreases below the solidification point and the crystallization process starts.

We find a number of the first 'portion' of nuclei *Ni* (1) (using equation(34)) and the final radius of grains (formula (35)).

The first value of ω is equal to

Numerical Modeling of Solidification Process 523

0 0 *e V e*

where Δ*H*0 is a change of control volume enthalpy during the time interval Δ*t*, *Qe* is a heat conducted at the time Δ*t* from the adjoining nodes to the node *x*0, *qV* is a mean capacity of internal heat sources resulting from the evolution of latent heat in the control volume Δ*V*0. If one assumes that the heat fluxes flowing to the element Δ*V*0 are proportional to the

> 0 0

<sup>0</sup> 2 2 *<sup>f</sup> e e <sup>e</sup> f f*

boundary volumes (Mochnacki & Suchy, 1995). In a case of no-flux boundary condition (see:

The change of enthalpy of the control volume Δ*V*0 during the time Δ*t* equals (Domański

where *C*0 *f* is the volumetric specific heat (substitute thermal capacity), *f*, *f*+1 denotes two

 <sup>1</sup> 0 0 00 0

*h h <sup>R</sup>* 

*e*

*<sup>e</sup>* are the thermal conductivities in the control volumes Δ*V*0 and Δ*Ve* at the

*<sup>f</sup> <sup>R</sup> <sup>e</sup>* (in numerical realization e.g. <sup>10</sup>

*ff f H CT T V* (44)

. The other definition of thermal resistance should be introduced for the

*<sup>f</sup> <sup>R</sup> <sup>e</sup>* is the thermal resistance between points *x*0 and *xe* (Majchrzak & Mendakiewicz, 1993), Δ*Ae* surface limiting the domain Δ*V*0 in the direction *e*. If we denote by *he* the distance

*f f <sup>e</sup> e e <sup>f</sup> e T T Q At <sup>R</sup>*

0

surface limiting the domain Δ*V*0 in direction *e* is a part of 'adiabatic' boundary.

et al., 2009a; Domański et al., 2009b; Mochnacki & Ciesielski, 2007)

Let us write the balance equation using the explicit scheme

*H QV <sup>q</sup> <sup>t</sup>* (41)

, then we shall obtain a solving system of the

(42)

(43)

<sup>0</sup> <sup>10</sup> *<sup>f</sup> <sup>R</sup> <sup>e</sup>* ) if the

The energy balance for the control volume Δ*V*0 can be written in the form

Fig. 5. Control volume *V*<sup>0</sup>

type 'explicit scheme'. So

between the nodes *x*0, *xe* then

 and *<sup>f</sup>* 

example presented) one assumes 0

where <sup>0</sup>

where 0

*f* 

successive time levels.

moment *t* = *t f*

temperature differences at the moment *t* = *t f*

$$
\Delta \alpha\_i^{(1)} = \frac{4}{3} \pi \,\mathrm{v} \,\mathrm{N}\_i^{(1)} \,\Delta \,\mathrm{R}\_i^{(1)} \tag{38}
$$

This value is introduced into the macro model for transition corresponding to the next time interval.

In the second stage of crystallization process modeling we find the quantity *Ni* using the equation (34) and next we can estimate the size of the second generation, this means *Ni* (2) = *Ni Ni* (1). We can also find the new increment of the grains radius Δ*Ri* (2). It should be pointed out that the current radii of the first generation are equal to Δ*Ri* (1) + Δ*Ri* (2), while for the second generation: Δ*Ri* (2).

The new value of function ω is determined by the formula

$$\rho \alpha\_{i}^{(2)} = \frac{4}{3} \pi \upsilon \left[ N\_{i}^{(1)} \left( \Delta R\_{i}^{(1)} + \Delta R\_{i}^{(2)} \right)^{3} + N\_{i}^{(2)} \Delta R\_{i}^{(2)3} \right] \tag{39}$$

The next steps of crystallization process modeling result from the generalization of considerations above presented. It should be pointed out that after the transition by the maximum undercooling, the number of nuclei at control volume Δ*Vi* is constant and this fact must be taken into account in adequate numerical procedure. The knowledge of two successive values of ω*<sup>i</sup>* (*<sup>k</sup>*) allows to estimate the source term in equation (1) using the differential quotient, in particular

$$\operatorname{L} \frac{\partial f\_S(\mathbf{x}, t)}{\partial t} = \operatorname{L} \frac{\operatorname{d} f\_S}{\operatorname{d} \, \alpha \nu} \frac{\partial \operatorname{o}(\mathbf{x}, t)}{\partial t} \approx \operatorname{L} [(n - 1)\alpha + 1]^{\frac{n}{1 - n}} \frac{\Delta \, \alpha \nu}{\Delta t} \tag{40}$$

### **4. Examples of numerical methods application**

In this sub-chapter only the selected problems connected with numerical modelling of solidification process will be discussed. In particular a certain version of control volume method and next the algorithms basing on the boundary element method will be here presented.

So, the first part of considerations is devoted to the application of control volume method (CVM) for numerical simulation of solidification process both in the macro and micro/macro scale. The control volume method (Mochnacki & Suchy, 1995; Domański et al., 2009a; Domański et al., 2009b; Mochnacki & Ciesielski, 2007) constitutes an effective tool for numerical computations of the heat transfer processes. The domain analyzed is divided into *N* volumes. The CVM algorithm allows to find the transient temperature field at the set of nodes corresponding to the central points of the control volumes. The nodal temperatures can be found on the basis of energy balances for the successive volumes. In order to assure the correctness and exactness of the algorithm proposed the control volumes are generated in the shape of the Voronoi polygons (see: Fig. 5).

Let us consider the control volume Δ*V*0 with the central node *x*0 . It is assumed here that the thermal capacities and capacities of the internal heat sources are concentrated in the nodes representing elements, while thermal resistances are concentrated in the sectors joining the nodes

(1) 4 (1) (1) 3 *i ii*

This value is introduced into the macro model for transition corresponding to the next time

In the second stage of crystallization process modeling we find the quantity *Ni* using the equation (34) and next we can estimate the size of the second generation, this means

(1). We can also find the new increment of the grains radius Δ*Ri*

 <sup>3</sup> (2) <sup>4</sup> (1) (1) (2) (2) (2)3 3 *i i i i ii*

The next steps of crystallization process modeling result from the generalization of considerations above presented. It should be pointed out that after the transition by the maximum undercooling, the number of nuclei at control volume Δ*Vi* is constant and this fact must be taken into account in adequate numerical procedure. The knowledge of two

<sup>1</sup> ( ,) d ( ,) [( 1) 1] <sup>d</sup>

In this sub-chapter only the selected problems connected with numerical modelling of solidification process will be discussed. In particular a certain version of control volume method and next the algorithms basing on the boundary element method will be here

So, the first part of considerations is devoted to the application of control volume method (CVM) for numerical simulation of solidification process both in the macro and micro/macro scale. The control volume method (Mochnacki & Suchy, 1995; Domański et al., 2009a; Domański et al., 2009b; Mochnacki & Ciesielski, 2007) constitutes an effective tool for numerical computations of the heat transfer processes. The domain analyzed is divided into *N* volumes. The CVM algorithm allows to find the transient temperature field at the set of nodes corresponding to the central points of the control volumes. The nodal temperatures can be found on the basis of energy balances for the successive volumes. In order to assure the correctness and exactness of the algorithm proposed the control volumes are generated

Let us consider the control volume Δ*V*0 with the central node *x*0 . It is assumed here that the thermal capacities and capacities of the internal heat sources are concentrated in the nodes representing elements, while thermal resistances are concentrated in the sectors joining the

 

*tt t*

*S S <sup>n</sup> f xt f xt L L L n*

*vN R R N R*

(*<sup>k</sup>*) allows to estimate the source term in equation (1) using the

(40)

*n*

*vN R* (38)

(2). It should be

(39)

(2), while for

(1) + Δ*Ri*

 

pointed out that the current radii of the first generation are equal to Δ*Ri*

(2). The new value of function ω is determined by the formula

**4. Examples of numerical methods application** 

in the shape of the Voronoi polygons (see: Fig. 5).

 

interval.

(2) = *Ni Ni*

the second generation: Δ*Ri*

successive values of ω*<sup>i</sup>*

presented.

nodes

differential quotient, in particular

*Ni*

Fig. 5. Control volume *V*<sup>0</sup>

The energy balance for the control volume Δ*V*0 can be written in the form

$$
\Delta H\_0 = \sum\_{\varepsilon} Q\_{\varepsilon} + \left| \Delta V\_0 \right| q\_V \Delta t \tag{41}
$$

where Δ*H*0 is a change of control volume enthalpy during the time interval Δ*t*, *Qe* is a heat conducted at the time Δ*t* from the adjoining nodes to the node *x*0, *qV* is a mean capacity of internal heat sources resulting from the evolution of latent heat in the control volume Δ*V*0. If one assumes that the heat fluxes flowing to the element Δ*V*0 are proportional to the temperature differences at the moment *t* = *t f* , then we shall obtain a solving system of the type 'explicit scheme'. So

$$Q\_e = \frac{T\_e^f - T\_0^f}{R\_{0e}^f} \Delta A\_e \Delta t \tag{42}$$

where <sup>0</sup> *<sup>f</sup> <sup>R</sup> <sup>e</sup>* is the thermal resistance between points *x*0 and *xe* (Majchrzak & Mendakiewicz, 1993), Δ*Ae* surface limiting the domain Δ*V*0 in the direction *e*. If we denote by *he* the distance between the nodes *x*0, *xe* then

$$R\_{0e}^f = \frac{h\_e}{2\lambda\_0^f} + \frac{h\_e}{2\lambda\_e^f} \tag{43}$$

where 0 *f* and *<sup>f</sup> <sup>e</sup>* are the thermal conductivities in the control volumes Δ*V*0 and Δ*Ve* at the moment *t* = *t f* . The other definition of thermal resistance should be introduced for the boundary volumes (Mochnacki & Suchy, 1995). In a case of no-flux boundary condition (see: example presented) one assumes 0 *<sup>f</sup> <sup>R</sup> <sup>e</sup>* (in numerical realization e.g. <sup>10</sup> <sup>0</sup> <sup>10</sup> *<sup>f</sup> <sup>R</sup> <sup>e</sup>* ) if the surface limiting the domain Δ*V*0 in direction *e* is a part of 'adiabatic' boundary.

The change of enthalpy of the control volume Δ*V*0 during the time Δ*t* equals (Domański et al., 2009a; Domański et al., 2009b; Mochnacki & Ciesielski, 2007)

$$
\Delta H\_0 = \mathbb{C}\_0^f \left( T\_0^{f+1} - T\_0^f \right) \left| \Delta V\_0 \right| \tag{44}
$$

where *C*0 *f* is the volumetric specific heat (substitute thermal capacity), *f*, *f*+1 denotes two successive time levels.

Let us write the balance equation using the explicit scheme

Numerical Modeling of Solidification Process 525

As an example the results of computations concerning the solidification of aluminium casting (2D problem) are presented. The macro control volume shown in Figure 8 (Mochnacki & Ciesielski, 2007) has been divided into micro control volumes. The heat fluxes (the Neumann boundary conditions) given on the periphery of macro control volume result from the solution of macroscopic problem (here the finite difference method has been applied). The micro model has been solved using the CVM under the assumption that the nuclei density is a constant value (this assumption is often used), the nuclei growth results from equation (35) for *m* = 2, while the changes of solid state fraction are determined by the Kolmogoroff formula. The result (cooling curve) shown in Figure 9 concerns the central point of macro control volume being in the direct thermal contact with sand mould subdomain. The micro/macro approach allows one to observe the undercooling below the solidification point and next, the growth of local temperature. This effect cannot be observed

when the numerical solution bases on the macro approach.

Fig. 6. Example of structure and its discretization.

$$\left|C\_{0}^{f}\left(T\_{0}^{f+1} - T\_{0}^{f}\right)\right| \Delta V\_{0}\big| = \sum\_{\varepsilon} \frac{T\_{\varepsilon}^{f} - T\_{0}^{f}}{R\_{0\_{\varepsilon}}^{f}} \Delta A\_{\varepsilon} \Delta t + \left|\Delta V\_{0}\right| q\_{V} \Delta t \tag{45}$$

or

$$T\_0^{f+1} = \sum\_{\varepsilon=0}^{n} \mathcal{W}\_{\varepsilon} T\_{\varepsilon}^{f} + \frac{q\_V \Delta t}{\mathcal{C}\_0^{f}} \tag{46}$$

where

$$\mathcal{W}\_{\varepsilon} = \frac{\Delta t \,\Delta A\_{\varepsilon}}{\mathcal{C}\_{0}^{\prime} \mathcal{R}\_{0\epsilon}^{\prime} \left| \Delta V\_{0} \right|^{\prime}} \quad \varepsilon = 1, \ldots, n \tag{47}$$

$$\mathcal{W}\_0 = 1 - \sum\_{\iota=1}^n \mathcal{W}\_{\iota} \tag{48}$$

In order to assure the stability of the above explicit scheme the coefficient *W*0 must be positive.

It was mentioned above that the control volumes have been constructed in the form of Voronoi polygons (Thiessen cells). A single polygon is defined by the lines that bisect the lines between the central point and its surrounding points. The bisecting lines and the connection lines are perpendicular to each other and this property is very essential when the control volume method is used.

Many algorithms to construct the Voronoi polygons can be found in literature. One popular method is based on the Delaunay triangulation (Domański et al., 2009b; Mochnacki & Ciesielski, 2007) and this method has been used at the stage of 'in house' computer program construction. The aim of computations was the analysis of thermal processes in domain of solidifying cast composite made from Al-Si alloy (metal matrix) and silicon (fibres). In particular the problem of 'spontaneous' solidification of matrix due to presence of fibres has been analyzed.

In Figure 6 an example of the structure of cast composite with 40% fibres is shown. The geometrical parameters of a matrix-fibres system are chosen on the basis of the optical micrographs presented in (Domański et al., 2009b). It can be seen that the fibres diameters are different, and also their mutual positions are rather incidental. The only unquestionable information concerning the geometry of the system results from the volumetric fraction of the fibres in the domain analyzed. So, the numerical procedure realizing the mesh generation, fibres localization, values of fibres radii bases on the application of random numbers generation. The sub-domain presented in Figure 6 corresponds to the central part of composite (side of square equals 100 m) and on the external surface the no-flux conditions have been assumed. It turned out that the spontaneous solidification takes place for 55% fraction of fibres. The kinetics of solidification in this case is shown in Figure 7. The input data concerning the parameters of cast composite components can be found in (Domański et al., 2009b).

A similar algorithm basing on the control volume method can be used also in a case of macro/micro modelling of solidification. One can see that the shapes of Thiessen cells are similar to the primary structure of casting and it was a reason of investigations connected with the application of presented version of CVM in a scope of macro/micro problems.

*e e T T CT T V A t Vq t <sup>R</sup>*

*<sup>n</sup> f f <sup>V</sup>*

*q t T WT*

*t A W en CR V*

0

0 0

1

*e W W*

In order to assure the stability of the above explicit scheme the coefficient *W*0 must be

It was mentioned above that the control volumes have been constructed in the form of Voronoi polygons (Thiessen cells). A single polygon is defined by the lines that bisect the lines between the central point and its surrounding points. The bisecting lines and the connection lines are perpendicular to each other and this property is very essential when the

Many algorithms to construct the Voronoi polygons can be found in literature. One popular method is based on the Delaunay triangulation (Domański et al., 2009b; Mochnacki & Ciesielski, 2007) and this method has been used at the stage of 'in house' computer program construction. The aim of computations was the analysis of thermal processes in domain of solidifying cast composite made from Al-Si alloy (metal matrix) and silicon (fibres). In particular the problem of 'spontaneous' solidification of matrix due to presence of fibres has

In Figure 6 an example of the structure of cast composite with 40% fibres is shown. The geometrical parameters of a matrix-fibres system are chosen on the basis of the optical micrographs presented in (Domański et al., 2009b). It can be seen that the fibres diameters are different, and also their mutual positions are rather incidental. The only unquestionable information concerning the geometry of the system results from the volumetric fraction of the fibres in the domain analyzed. So, the numerical procedure realizing the mesh generation, fibres localization, values of fibres radii bases on the application of random numbers generation. The sub-domain presented in Figure 6 corresponds to the central part of composite (side of square equals 100 m) and on the external surface the no-flux conditions have been assumed. It turned out that the spontaneous solidification takes place for 55% fraction of fibres. The kinetics of solidification in this case is shown in Figure 7. The input data concerning the parameters of cast composite components can be found in

A similar algorithm basing on the control volume method can be used also in a case of macro/micro modelling of solidification. One can see that the shapes of Thiessen cells are similar to the primary structure of casting and it was a reason of investigations connected with the application of presented version of CVM in a scope of macro/micro problems.

*e e f*

, 1,..., *<sup>e</sup>*

*e*

*C*

*f f*

*f e V*

(45)

(46)

(47)

(48)

00 0 0 0

<sup>1</sup> <sup>0</sup>

1 0

*e f f e*

*e*

00 0

1 *n*

0

*ff f e*

or

where

positive.

been analyzed.

(Domański et al., 2009b).

control volume method is used.

As an example the results of computations concerning the solidification of aluminium casting (2D problem) are presented. The macro control volume shown in Figure 8 (Mochnacki & Ciesielski, 2007) has been divided into micro control volumes. The heat fluxes (the Neumann boundary conditions) given on the periphery of macro control volume result from the solution of macroscopic problem (here the finite difference method has been applied). The micro model has been solved using the CVM under the assumption that the nuclei density is a constant value (this assumption is often used), the nuclei growth results from equation (35) for *m* = 2, while the changes of solid state fraction are determined by the Kolmogoroff formula. The result (cooling curve) shown in Figure 9 concerns the central point of macro control volume being in the direct thermal contact with sand mould subdomain. The micro/macro approach allows one to observe the undercooling below the solidification point and next, the growth of local temperature. This effect cannot be observed when the numerical solution bases on the macro approach.

Fig. 6. Example of structure and its discretization.

Numerical Modeling of Solidification Process 527

Now, the problems connected with the boundary element method (Brebbia et al., 1984; Brebbia & Dominguez, 1992; Majchrzak, 2001) application in the scope of solidification process modeling will be discussed. So, the boundary element method constitutes the effective tool for numerical simulation of solidification processes (Majchrzak & Mochnacki, 1995; Majchrzak & Mochnacki, 1996; Majchrzak & Szopa, 2001; Mochnacki & Majchrzak, 2007b). The method assures the exact approximation of real shapes of the boundaries and also the very good exactness of boundary conditions approximation (Brebbia et al., 1984; Brebbia & Dominguez, 1992). These features of the BEM are very essential in the case of solidification processes modeling because the big gradients of temperature near the boundary require the adequate exact method of temperature field computations. The main disadvantage of the method is the fact that it can be used only for linear tasks (the constant values of metal and mould thermophysical parameters) because the fundamental solutions are only available for linear equations. Nevertheless, several techniques have been developed for the numerical solution of nonlinear problems. In the case of solidification problem the basic BEM algorithm for linear Fourier's equation should be coupled with procedures correcting the temporary solutions for successive levels of time. These procedures are called the temperature field correction method (Mochnacki & Suchy, 1995; Mochnacki & Szopa, 2008), the alternating phase truncation method (Mochnacki & Suchy, 1995; Majchrzak & Mochnacki, 1995; Majchrzak & Mendakiewicz, 1993) and the artificial

Fig. 9. Temperature history at the central point of macro control volume.

heat source method (Mochnacki & Majchrzak, 2010; Mochnacki & Szopa, 1998).

quite effective version of the BEM) will be presented. Let us consider the following parabolic equation

The basic idea of the BEM can be realized in different ways. In the case of parabolic equations the I and II scheme of the BEM, the dual reciprocity BEM, the BEM using discretization in time etc. are used. The detailed information concerning the BEM application in the scope of heat transfer problems can be found in the numerous books and papers published by Ewa Majchrzak from the Silesian University of Technology (Poland). Below, the short discussion of the BEM using discretization in time taken from (Majchrzak & Szopa, 2001; Mochnacki & Szopa, 1998; Mochnacki & Szopa, 2002) (it is not very popular but

> <sup>2</sup> ( ,) ( ,) ( ,) *T x t Q x t a T x t tc*

(49)

Fig. 7. Kinetics of composite solidification.

Fig. 8. The inner structure of macro control volume.

Fig. 7. Kinetics of composite solidification.

Fig. 8. The inner structure of macro control volume.

Fig. 9. Temperature history at the central point of macro control volume.

Now, the problems connected with the boundary element method (Brebbia et al., 1984; Brebbia & Dominguez, 1992; Majchrzak, 2001) application in the scope of solidification process modeling will be discussed. So, the boundary element method constitutes the effective tool for numerical simulation of solidification processes (Majchrzak & Mochnacki, 1995; Majchrzak & Mochnacki, 1996; Majchrzak & Szopa, 2001; Mochnacki & Majchrzak, 2007b). The method assures the exact approximation of real shapes of the boundaries and also the very good exactness of boundary conditions approximation (Brebbia et al., 1984; Brebbia & Dominguez, 1992). These features of the BEM are very essential in the case of solidification processes modeling because the big gradients of temperature near the boundary require the adequate exact method of temperature field computations. The main disadvantage of the method is the fact that it can be used only for linear tasks (the constant values of metal and mould thermophysical parameters) because the fundamental solutions are only available for linear equations. Nevertheless, several techniques have been developed for the numerical solution of nonlinear problems. In the case of solidification problem the basic BEM algorithm for linear Fourier's equation should be coupled with procedures correcting the temporary solutions for successive levels of time. These procedures are called the temperature field correction method (Mochnacki & Suchy, 1995; Mochnacki & Szopa, 2008), the alternating phase truncation method (Mochnacki & Suchy, 1995; Majchrzak & Mochnacki, 1995; Majchrzak & Mendakiewicz, 1993) and the artificial heat source method (Mochnacki & Majchrzak, 2010; Mochnacki & Szopa, 1998).

The basic idea of the BEM can be realized in different ways. In the case of parabolic equations the I and II scheme of the BEM, the dual reciprocity BEM, the BEM using discretization in time etc. are used. The detailed information concerning the BEM application in the scope of heat transfer problems can be found in the numerous books and papers published by Ewa Majchrzak from the Silesian University of Technology (Poland). Below, the short discussion of the BEM using discretization in time taken from (Majchrzak & Szopa, 2001; Mochnacki & Szopa, 1998; Mochnacki & Szopa, 2002) (it is not very popular but quite effective version of the BEM) will be presented.

Let us consider the following parabolic equation

$$\frac{\partial T\left(\mathbf{x},t\right)}{\partial t} = a\left(\nabla^2 T\left(\mathbf{x},t\right) + \frac{Q\left(\mathbf{x},t\right)}{c}\right) \tag{49}$$

Numerical Modeling of Solidification Process 529

Applying the second Green formula to the first component of criterion (53) and taking into account the property (55) of fundamental solution one obtains the following boundary

\*

*T t T x q x t* 

( ξ , ) ( , ) d <sup>λ</sup>

*q xT xt*

1 ( , ) ( ξ , ) d

*T xt T x*

*f*

( , ) ( ξ , ) d <sup>λ</sup>

For example, in a case of 1D problem (plate of thickness *L*) the last equation takes a form

\*

*f f*

 

<sup>1</sup> (, ) (,)(, )

0

0

11 12 21 22

> *f f*

 

*g g qt g g qLt h h Tt p z h h TLt pL zL*

 

*L*

*a t*

*L*

*T t T xqxt*

*Q xt T x*

*f*

<sup>1</sup> ( <sup>ξ</sup> ,) ( <sup>ξ</sup>, ) ( , ) d <sup>λ</sup>

*f f*

1 \*

*f*

0

 

(59)

(60)

one obtains the system of equations which can be written in the

*L*

(57)

(58)

(61)

\*

0

\* 1

\*

(0, ) (, ) (0, ) (0) (0) () () (, )

*f f*

 

*L f*

<sup>1</sup> (,)(, ) () ()

*q xTxt p z*

<sup>1</sup> ( ) ( , ) ( , )d

*<sup>f</sup> p T x Txt x*

<sup>1</sup> ( ) ( , ) ( , )d

*<sup>f</sup> z Qxt T x x*

 

\*

1

\*

11 12 21 22

1

*a t*

can be calculated in analytical way.

integral equation (ξ Γ)

where

and

For ξ 0+ and for ξ *L*

matrix form

where

\* \* *q x n T x* ( ξ , ) λ ( ξ , ) (56)

where *a* is a thermal diffusivity, *c* is a volumetric specific heat, *Q* is a capacity of internal heat sources. The equation (49) is supplemented by the assumed boundary and initial conditions. At first the time grid

$$0 = t^0 < t^1 < \dots < t^{f-2} < t^{f-1} < t^f < \dots < t^F < \infty \tag{50}$$

with constant step Δ*t* = *t f t f*  1 is introduced.

We substitute the energy equation (10) for time *t*  [*t f* 1, *t f*] by the equation in which the time derivative is approximated by the differential quotient, in particular

$$\frac{T\left(\mathbf{x},t^{f}\right) - T\left(\mathbf{x},t^{f-1}\right)}{\Delta t} = a\,\,\,\nabla^{2}T\left(\mathbf{x},t^{f}\right) + \frac{Q\left(\mathbf{x},t^{f}\right)}{c} \tag{51}$$

The last equation can be written in the form

$$\begin{aligned} \nabla^2 T \left( \mathbf{x}, t^f \right) - \frac{1}{a \Delta t} T \left( \mathbf{x}, t^f \right) + \\ \frac{1}{a \Delta t} T \left( \mathbf{x}, t^{f-1} \right) + \frac{Q \left( \mathbf{x}, t^f \right)}{\lambda} = 0 \end{aligned} \tag{52}$$

Using the weighted residual method criterion (WRM) (Stefanescu, 1999; Szopa, 1999) one obtains

$$\begin{aligned} &\int\_{\Omega} \left[ \nabla^2 T(\mathbf{x}, t^f) - \frac{1}{a \, \Delta t} T(\mathbf{x}, t^f) + \\ &\frac{1}{a \, \Delta t} T\left(\mathbf{x}, t^{f-1}\right) + \frac{Q(\mathbf{x}, t^f)}{\lambda} \right] T^\*(\boldsymbol{\xi}, \mathbf{x}) \, \mathrm{d}\Omega = 0 \end{aligned} \tag{53}$$

where *T\**(ξ, *x*) is the fundamental solution and for the objects oriented in rectangular coordinate system it is a function of the form

$$T^\*(\boldsymbol{\xi}, \boldsymbol{x}) = \begin{cases} \frac{\sqrt{a \, \Delta t}}{2} \exp\left(-\frac{r}{\sqrt{a \, \Delta t}}\right), & \text{for } \text{D problem} \\\\ \frac{1}{2\pi} \, \mathop{\,\mathrm{K}\_0}\left(\begin{array}{c} r \\ \hline \sqrt{a \, \Delta t} \end{array}\right), & \text{for } \text{2D problem} \\\\ \frac{1}{4\pi r} \, \exp\left(-\frac{r}{\sqrt{a \, \Delta t}}\right), & \text{for 3D problem} \end{cases} \tag{54}$$

where K0() is the modified Bessel function of zero order, *r* is the distance between the points *x* and ξ.

The fundamental solution fulfills the following equation

$$-\nabla^2 T^\*(\xi, \mathbf{x}) - \frac{1}{a\Delta t} T^\*(\xi, \mathbf{x}) = -\mathcal{S}\left(\xi, \mathbf{x}\right) \tag{55}$$

Additionally, on a basis of formula (54) the heat flux resulting from the fundamental solution

$$q^\*(\xi, \mathbf{x}\_\prime) = -\lambda \, n \cdot \nabla T^\*(\xi, \mathbf{x}\_\prime) \tag{56}$$

can be calculated in analytical way.

Applying the second Green formula to the first component of criterion (53) and taking into account the property (55) of fundamental solution one obtains the following boundary integral equation (ξ Γ)

$$\begin{aligned} \left(T\left(\xi,t^{f}\right)\right) &+ \frac{1}{\lambda} \int\_{\Gamma} T^\*\left(\xi,x\right) q\left(x,t^{f}\right) \mathrm{d}\,\Gamma = \\ &\frac{1}{\lambda} \int\_{\Gamma} q^\*\left(\xi,x\right) T\left(x,t^{f}\right) \mathrm{d}\,\Gamma + \\ &\frac{1}{a\,\Delta t} \int\_{\Omega} T\left(x,t^{f-1}\right) T^\*\left(\xi,x\right) \mathrm{d}\,\Omega + \\ &\frac{1}{\lambda} \int\_{\Omega} Q\left(x,t^{f}\right) T^\*\left(\xi,x\right) \mathrm{d}\,\Omega \end{aligned} \tag{57}$$

For example, in a case of 1D problem (plate of thickness *L*) the last equation takes a form

$$\begin{aligned} \left[T(\boldsymbol{\xi},t^{f}) + \left[\frac{1}{\lambda}T^{\*}(\boldsymbol{\xi},\mathbf{x})q(\mathbf{x},t^{f})\right]^{L}\right]\_{0}^{L} &= \\ \left[\frac{1}{\lambda}q^{\*}(\boldsymbol{\xi},\mathbf{x})T(\mathbf{x},t^{f})\right]\_{0}^{L} + p(\boldsymbol{\xi}) + z(\boldsymbol{\xi}) \end{aligned} \tag{58}$$

where

528 Computational Simulations and Applications

where *a* is a thermal diffusivity, *c* is a volumetric specific heat, *Q* is a capacity of internal heat sources. The equation (49) is supplemented by the assumed boundary and initial

We substitute the energy equation (10) for time *t*  [*t f* 1, *t f*] by the equation in which the

<sup>2</sup> (, ) (, ) (, ) (, ) *f f <sup>f</sup> T x t T x t <sup>f</sup> Q x t a T x t* 

<sup>1</sup> (, ) (, )

*f f*

*<sup>f</sup> <sup>f</sup>*

<sup>1</sup> (, ) (, ) <sup>0</sup>

Using the weighted residual method criterion (WRM) (Stefanescu, 1999; Szopa, 1999) one

1 \*

where *T\**(ξ, *x*) is the fundamental solution and for the objects oriented in rectangular co-

exp , 2

exp , 4

*r a t*

where K0() is the modified Bessel function of zero order, *r* is the distance between the

<sup>2</sup> \* \* <sup>1</sup> *T x T x x* ( <sup>ξ</sup> ,) (<sup>ξ</sup> , ) <sup>δ</sup> ( <sup>ξ</sup> , ) *a t*

Additionally, on a basis of formula (54) the heat flux resulting from the fundamental

*Q x t T x t T x* 

<sup>1</sup> (, ) (, ) (<sup>ξ</sup> , )d 0 <sup>λ</sup>

 

*a t*

*a t r*

K for 2 D problem

λ

 *tc* 

1

1

<sup>1</sup> (, ) (, )

*T x t T x t a t*

*<sup>f</sup> <sup>f</sup>*

*a t r*

<sup>1</sup> ( , ) , 2

1

*<sup>r</sup> T x* 

*f f*

 *T x t T x t a t Q x t T x t*

 

time derivative is approximated by the differential quotient, in particular

2

*a t*

2

\* <sup>0</sup>

*a t*

ordinate system it is a function of the form

The fundamental solution fulfills the following equation

0 1 2 1 0 ... ... *fff F t t t t t t*  (50)

(51)

for 1 D problem

for 3 D problem

(55)

(52)

(53)

(54)

conditions. At first the time grid

with constant step Δ*t* = *t f t f*  1 is introduced.

The last equation can be written in the form

obtains

points *x* and ξ.

solution

$$p(\xi) = \frac{1}{a\Delta t} \Big| \int\_0^L T^\*(\xi, \mathbf{x}) \, T(\mathbf{x}, t^{f-1}) \, \mathbf{dx} \tag{59}$$

and

$$z(\xi) = \frac{1}{\lambda} \bigcup\_{0}^{L} \mathcal{Q}(x, t^{f}) \; T^{\*}(\xi, x) \; \text{d}x \tag{60}$$

For ξ 0+ and for ξ *L* one obtains the system of equations which can be written in the matrix form

$$
\begin{bmatrix}
\mathcal{g}\_{11} & \mathcal{g}\_{12} \\
\mathcal{g}\_{21} & \mathcal{g}\_{22}
\end{bmatrix}
\begin{bmatrix}
q(0,t^f) \\
q(L,t^f)
\end{bmatrix} = \\
$$

$$
\begin{bmatrix}
h\_{11} & h\_{12} \\
h\_{21} & h\_{22}
\end{bmatrix}
\begin{bmatrix}
T(0,t^f) \\
T(L,t^f)
\end{bmatrix} + 
\begin{bmatrix}
p(0) \\
p(L)
\end{bmatrix} + 
\begin{bmatrix}
z(0) \\
z(L)
\end{bmatrix}
\tag{61}
$$

where

Numerical Modeling of Solidification Process 531

The second example concerns the 3D problem solution (Szopa, 1999). Figure 12 shows the

The boundary surface was divided into constant boundary elements (squares), whereas the interior was divided into constant internal cells (cubes). The boundary condition of the form *x*  Γ: *T*(*x*, *t*) = 1465°C was assumed (this value results from well known Schwarz's solution and was treated as a constant one). The initial temperature *T*(*x*, 0) = 1550°C whereas the thermophysical parameters correspond to the cast steel 0.35%C. The solidification process was taken into account by the correction of temperature field for successive levels of time (Temperature Field Correction Method). In Figure 13 the temporary temperature field for

shape of considered steel casting domain (dimensions are marked in mm).

Fig. 11. Source function evolution.

*t* = 180s and two selected vertical sections is shown.

Fig. 12. Temperature field in domain considered.

$$g\_{11} = -g\_{22} = -\frac{\sqrt{\Delta t}}{2\sqrt{\lambda c}} \qquad g\_{12} = -g\_{21} = \frac{\sqrt{\Delta t}}{2\sqrt{\lambda c}} \exp\left(-\frac{L}{\sqrt{a\Delta t}}\right) \tag{62}$$

while

$$hh\_{11} = h\_{22} = -\frac{1}{2} \qquad h\_{12} = h\_{21} = \frac{1}{2} \exp\left(-\frac{L}{\sqrt{a\Delta t}}\right) \tag{63}$$

Equation (61) allows to determine the 'missing' boundary values (temperatures and heat fluxes for *x* = 0 and *x* = *L*), while in the second step of computations the temporary temperatures in the set of internal points ξ (0, *L*) can be found on the basis of equation (58).

As the example of 1D solution concerning the solidification process (micro/macro approach) a problem of heat transfer in domain of aluminium plate (*L* = 0.02 [m]) is discussed. On the outer surface of the plate the boundary temperature *TB* = 655°C, while for the axis of symmetry *qB* = 0 are assumed. The growth coefficient is equal to μ = 310−<sup>6</sup> [m/sK2], at the same time the nucleation coefficient: ψ = 1010 [m−3 K−2], exponent *n* = 1.5 (equation (31)). Pouring temperature: *T*0 = 700°C, solidification point: *Tc* = 660°C. The thermophysical parameters of Al are assumed to be constant. The solidification model corresponds to description previously presented (see: formulas (38), (39), (40)). In domain of plate 20 internal points has been distinguished and in Figure 10 the cooling curves at points 10, 14 and 18 are presented, while in Figure 11 the evolution of source function at the same points are shown.

Fig. 10. Cooling curves.

1 1 2 2 *<sup>L</sup> h h h h exp a t* 

Equation (61) allows to determine the 'missing' boundary values (temperatures and heat fluxes for *x* = 0 and *x* = *L*), while in the second step of computations the temporary temperatures in

As the example of 1D solution concerning the solidification process (micro/macro approach) a problem of heat transfer in domain of aluminium plate (*L* = 0.02 [m]) is discussed. On the outer surface of the plate the boundary temperature *TB* = 655°C, while for the axis of symmetry *qB* = 0 are assumed. The growth coefficient is equal to μ = 310−<sup>6</sup> [m/sK2], at the same time the nucleation coefficient: ψ = 1010 [m−3 K−2], exponent *n* = 1.5 (equation (31)). Pouring temperature: *T*0 = 700°C, solidification point: *Tc* = 660°C. The thermophysical parameters of Al are assumed to be constant. The solidification model corresponds to description previously presented (see: formulas (38), (39), (40)). In domain of plate 20 internal points has been distinguished and in Figure 10 the cooling curves at points 10, 14 and 18 are presented, while in Figure 11 the evolution of source function at the same

*t tL*

 *c c at*

(62)

(63)

11 22 12 21 2 2

11 22 12 21

the set of internal points ξ (0, *L*) can be found on the basis of equation (58).

while

points are shown.

Fig. 10. Cooling curves.

*g g g g exp*

Fig. 11. Source function evolution.

The second example concerns the 3D problem solution (Szopa, 1999). Figure 12 shows the shape of considered steel casting domain (dimensions are marked in mm).

The boundary surface was divided into constant boundary elements (squares), whereas the interior was divided into constant internal cells (cubes). The boundary condition of the form *x*  Γ: *T*(*x*, *t*) = 1465°C was assumed (this value results from well known Schwarz's solution and was treated as a constant one). The initial temperature *T*(*x*, 0) = 1550°C whereas the thermophysical parameters correspond to the cast steel 0.35%C. The solidification process was taken into account by the correction of temperature field for successive levels of time (Temperature Field Correction Method). In Figure 13 the temporary temperature field for *t* = 180s and two selected vertical sections is shown.

Fig. 12. Temperature field in domain considered.

Numerical Modeling of Solidification Process 533

problem solution with respect to parameter *zk* is defined as a partial derivative <sup>1</sup> *U xt Txtz z z <sup>k</sup>* ( , ) ( , , ,... ) / *n k* (Kleiber, 1997; Dems & Rousselet, 1999). As it was mentioned below the methods of sensitivity analysis allows one 'to rebuilt' the solution of problem considered to the solution concerning the other value of parameter *zk*. So, the set of parameters introduced to the computer program as the input data we denote by *z*10, *z*<sup>2</sup>

*zn*0. Additionally we are interested in the others solutions corresponding to changed input data (e.g. new pouring temperature). The methods of sensitivity analysis give the possibilities to transform the basic solution on the solution for disturbed initial parameters.

00 0 0

*k n kk k*

*k kn*

0 (the second order sensitivity) result from the differentiation of equations

00 0 0 0 0

where *T* is a searched function (e.g. temperature), *x* is a spatial coordinate, *t* is a time, *Uk*

and conditions describing the physical process considered with respect to parameter *zk* (direct approach) We can notice that formula (64) concerns to the disturbance of distinguished parameter, but one can analyze also the disturbance of the certain set of

Sensitivity analysis with respect to boundary and initial conditions (e.g. heat transfer coefficient α, ambient temperature *Ta* , initial temperatures *T*0, *Tm*0) or mould parameters (*cm*, λ*m*) will be discussed, at the same time the direct approach will be applied. The problem is more complicated in comparison with the others ones, because the mathematical model of solidification is strongly non-linear (Szopa, 2005; Szopa & Wojciechowska, 2003; Szopa et al.,

At first, the energy equation (8) should be differentiated with respect to parameter *zk*, this

*k k*

*t*

*z tz*

(,) ( ) () (,)

 (66)

*t*

*Txt C T T Txt*

() d() () d() , *k kk k* d d *CT CT T T T T z Tz z Tz* 

(,) ( ) [ ( ) ( , )]

*Uxt C T T Uxt*

(,) [ '( ) ( , ) ( , )] '( ) ( , )

*Txt TUxt Txt C TUxt*

( , , , ,..., ,..., ) ... <sup>2</sup>

( , , , ,..., ,..., )

*Tx t z z z z z*

*<sup>z</sup> Tx t z z z z U z V*

0

*k*

1 2

respect to *zk* at the same point. The equations determining the values of *Ui*

The transformation results from the Taylor formula, namely

1 2

denotes the derivative *T*/*zk* at the start point 0, whereas *Vk*

sensitivity) and *Vi*

means

Because

parameters *zk*, simultaneously.

2004; Mochnacki & Szopa, 2009).

**5.1 Sensitivity analysis of macro models** 

therefore (using additionally Schwarz theorem)

0, ...,

0

(64)

(65)

(67)

0 (the first order

0 is the second derivative with

Fig. 13. Temperature field in domain considered.
