- Email: [email protected]
Prediction of the effective thermal conductivity of three-dimensional dendritic regions by the finite element method
Materials Science and Engineering A269 (1999) 90 – 97 www.elsevier.com/locate/msea
Prediction of the effective thermal conductivity of three-dimensional dendritic regions by the finite element method K. Ravindran, S.G.R. Brown, J.A. Spittle * Department of Materials Engineering, Uni6ersity of Wales Swansea, Swansea SA2 8PP, UK Received 5 October 1998; accepted 15 March 1999
Abstract The finite element method has been used to predict the variation in effective thermal conductivity, keff, in evolving three-dimensional dendritic mushy zones. The model demonstrates the influence of a dendrite-like geometry and volume fraction of solid on keff assuming that the conductivities of the solid and liquid phases remain constant and uniform as the structure evolves. The three-dimensional dendrite-like shapes have been generated by a cellular automaton-finite difference model. Numerical solutions are presented for both columnar and equiaxed dendritic zones. © 1999 Elsevier Science S.A. All rights reserved. Keywords: Thermal conductivity; Solidification; Dendrite growth; Numerical modelling
1. Introduction Numerical solution of the macroscopic transport of heat and fluid is becoming a common practice to analyse the various problems arising in metal casting. Significant advances have been made in the simulation of the free surface movement during mould filling, heat transfer during mould filling and the subsequent solidification and the development of thermal stresses. Efforts are also under way to enhance the predictive capability of these macroscopic models with microscopic models based on the understanding of nucleation, growth kinetics and other related phenomena. The ultimate goal is to develop tools which would be able to solve the entire casting process from the mould filling stage to the prediction of the microstructure and the mechanical properties at various sections of the cooled casting. Accurate numerical simulation of solidification transformation in a commercial multi-component alloy, leading to the prediction of microstructures and properties, requires accurate values of thermo-physical properties such as density, thermal conductivity, specific heat, * Corresponding author. Tel.: +44-1792-295589; fax: +44-1792295693. E-mail address: [email protected] (J.A. Spittle)
etc. Efforts initiated by the Department of Trade and Industry (DTI), UK, reveal some of the difficulties associated with measuring these properties in metals/alloys at high temperature. The measured property values for aluminium alloys showed a marked dependence on the thermal and mechanical histories of the samples and this could be due to the changes in the microstructure. It is rather difficult to isolate the latent heat effect in order to measure the thermal conductivity in the mushy stage accurately and hence efforts are being made to derive the thermal conductivity values from electrical conductivity measurements [1]. However, experimental measurement of all the phenomena accompanying solidification, as a function of time and varying fraction solid, within a continuously evolving dendrite structure is impossible. A numerical approach is, therefore, the only viable alternative way for predicting the thermal conductivity in such cases. Thermal conductivity of the mushy zone is a function of the morphology of the two phase region, the thermal conductivities of the solid and liquid (ks and kl), and the volume fraction of solid fs. Analytical expressions may be derived when some assumptions are made about the distribution of the solid and liquid in the mushy zone. The thermal resistance (which is inversely proportional to the thermal conductivity) of the solid and liquid
0921-5093/99/$ - see front matter © 1999 Elsevier Science S.A. All rights reserved. PII: S 0 9 2 1 - 5 0 9 3 ( 9 9 ) 0 0 1 6 0 - 4
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
regions may be lumped in such a way that the effective thermal resistance is either in ‘series’ or in ‘parallel’. The effective thermal conductivity for these two cases, kser and kpar, can be derived using the analogy of an equivalent electrical circuit:
kser =
 kskl à fskl + (1− fs)ks Ì
(1)
kpar =fsks + (1− fs)klà Škser and kpar are also known as the Wiener bounds. For any complex geometry, the effective thermal conductivity in a given direction should fall between these two limiting values. kpar gives the upper bound and kser the lower bound. In macroscopic alloy solidification analysis, due to the absence of other data, the thermal conductivity of the mushy zone is usually assumed to be a phase-average varying linearly with the volume fraction-solid, fs. This is also known as the ‘mixture model’ for effective thermal conductivity [2 – 6]: kmix = ks fs +kl(1− fs)
(2)
The value given above is the same as the equivalent thermal conductivity for a parallel circuit (kpar) and, as discussed earlier, this is an upper bound for keff and the actual value of keff could be anywhere between kpar and the lower bound, kser, depending on the morphology of the dendritic zone. In the case of columnar dendrites, keff depends on the direction as well, i.e. the effective thermal conductivity is anisotropic. Poirier and McBride [7] recommend the average of kpar and kser for keff. Since the conductivities in the solid and liquid phases of specific alloys are usually unknown, modellers frequently have to resort to using the solid and liquid conductivities of the pure solvent metal. Lee and Yang [8] have recently predicted the effective thermal conductivity of the columnar dendritic zone assuming that the dendrites are wedge-shaped. They predicted the thermal conductivity of the mushy zone along the growth direction (k xeff) as a function of fs using a steady state temperature field which has been averaged in the normal direction. This method is not
91
valid for a solidification analysis since the mushy zone usually encompasses many dendrites and the averaging has to be done for the whole dendrite and not just in one direction. Lee and Yang do not consider the effect of secondary arms and they suggested the use of kser as the thermal conductivity in the direction normal to growth (k yeff). Laouadi et al. [9] have added the effect of latent heat into a one-dimensional model for the effective thermal conductivity. None of these models considers the effect of the actual dendrite morphology on the effective thermal conductivity. In Section 2, we propose a general numerical approach for the prediction of thermal conductivity of the mushy zone when the dendrite morphology is known. In Section 3, a three-dimensional cellular automatonfinite difference model which is used for generating the dendrite shapes is briefly described. Finally numerical results obtained for generated dendritic structures are presented in Section 4.
2. Mathematical formulation In the present work, a periodic array of dendrites is assumed. The domain isolated for thermal-conductivity evaluation, a parallelepiped of dimensions L × W × B, encloses a single primary arm that has branched as shown in Fig. 1. The size of the computational domain is assumed to be much smaller than that of the domain where the effective property would be used for a macroscopic solidification analysis. A one-dimensional temperature field is imposed in a desired direction (say, x) and the steady state temperature field in the computational domain is obtained numerically by the finite element method. The unknown effective thermal conductivity k xeff of the isolated domain in the desired direction is evaluated by imposing the condition that the steady state heat flux is the same if the domain were filled with a medium having a constant thermal conductivity equal to k xeff. A similar formulation can be used for predicting the effective thermal conductivity in the y and z directions. The heat conduction equation and the boundary conditions prescribed for evaluating k xeff are given below: 9·k9T = 0
(3)
Boundary conditions:
 à T =T2 at x= L Ì (T = 0 at y=0, y=W, z =0 and z= B à (n Å T= T1 at x= 0
Fig. 1. Computational domain enclosing a branched primary arm of a columnar dendrite.
(4)
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
92
where n is the outward normal to the boundary of the computational diagram. The above equations may be nondimensionalised using the non-dimensional variables u, K, X, Y and Z as defined below: T−T2 k  K= ks à T1 −T2 Ì x y z X = Y= Z= à L L LÅ
u=
(5)
In the above, the length of the computational domain in the x direction (L) is used as the reference length and the thermal conductivity of the solid (ks) as the reference value for thermal conductivity. The non-dimensional form of the heat conduction equation and the boundary conditions are given below: 9·K9u =0
(6)
u=1 at X=0
 u=0 at X=1 Ã Ì (u à =0 at Y=0, Y=W/L, Z=0 and Z=B/L Å (n
(7)
With the nondimensionalisation defined above, the computed average non-dimensional wall flux is equal to the non-dimensional effective thermal conductivity (K xeff): k xeff L 2 = ks BW
(u q= − K = 1 at X= 0 (n
 à à u= 0 at X= 1 Ì Ã (u =0 at Y=0, Y=W/L, Z=0 and Z=B/L à (n Å
&& n (u (X
K
dY dZ
(8)
X=0
The above equation shows that the predicted effective thermal conductivity (k xeff) is independent of the temperature values (T1 and T2) prescribed on the boundary. It depends on the aspect ratio (L/W and L/B) of the domain isolated for averaging and not on the actual size of the domain (L, W and B). The effective thermal conductivity is the same for two dendrites having the same shape, but different sizes as long as the averaging domain is much smaller in size compared to the macroscopic domain where the averaged property would be used. The non-dimensional effective thermal conductivity (K xeff) does not depend directly on the thermal conductivity values of the solid or liquid (kl or ks) but on the thermal conductivity ratio (kl/ks). Therefore, numerical solutions can be obtained in terms of kl/ks and not on the actual values kl and ks which have to be provided by some other means such as from a database or an experiment. The effective thermal conductivity k xeff may also be obtained by prescribing a unit heat flux (non-dimen-
(9)
In this case, k xeff is evaluated by imposing the condition that for the same value of prescribed wall heat flux (q), the predicted average temperature at the X=0 wall should be the same if the domain were filled with a medium having a constant thermal conductivity equal to k xeff. The computed average wall temperature u( (X = 0) is related to k xeff by the following relation: 1 k 1 = xs = x K eff k eff u( (X = 0)
Boundary conditions:
K xeff =
sional) on the X=0 wall. The boundary conditions for this approach are given below:
(10)
The advantage of this alternate approach is that only the average of the predicted wall temperature is needed and not the temperature gradient which is not continuous in the nodes/element boundaries. Further the first derivatives of the variables are known to be accurate at the numerical (Gaussian) integration points and not on element boundaries. The heat conduction Eq. (6) is solved by the Galerkin finite element method. Eight node isoparametric brick elements have been used for discretizing the computational domain. The discrete form of the final equations may be written in matrix form as [10]: [A]{X}= %e [A e]{X}= {F}
(11)
where [A e] is the local stiffness matrix for cell e and the summation has to be performed over all the cells to form the global stiffness matrix, [A]. {X} is the array of unknown nodal temperature values and {F} is the nodal external heat flux vector which has nonzero values only for the nodes which fall on the boundary edges with a nonzero heat flux boundary condition. The stiffness matrix of a cell with a thermal conductivity k e may be obtained by numerical integration using the following relation: [A e]= k e
&&&
(9N)T9N dx dy dz
(12)
where 9N are the derivatives of the element interpolation (or shape) function. It may be noted from the above equation that the local stiffness matrix ([A e]) is the same for all cells but for a multiplying factor, the element thermal conductivity (k e). The cellular au-
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
tomata domains have a large number of cells and the storage requirements for the assembled global matrix ([A]) are quite high. As the discrete equations could not be solved by in-core solvers, an iterative solution approach based on the conjugate gradient method has been used for the solution of nodal unknowns [11]. No pre-conditioning is used as it involves an additional storage overhead of a large multi-dimensional array of the pre-conditioned element matrices. Also, there is no need to generate and maintain a three-dimensional matrix of cell nodal connectivity as the cells and their nodes can be easily addressed with a finite differencelike notation with three subscripts, i, j and k. Only a single element stiffness matrix (8 × 8) is stored and the solutions for a problem with 60×60 × 60 cells could be obtained within 12 min on a 200 MHz Pentium PC.
3. Cellular automaton-finite difference (CAFD) model for dendrite growth Modelling dendritic growth is complex as it requires a coupled analysis of the transport of heat and mass. Latent heat release during the phase change process and the need to consider other microscopic factors such as nucleation and growth make it a further challenge. In the CAFD model proposed by Brown and Spittle [12,13], the heat transfer phenomena is decoupled by assuming that the growth occurs under steady conditions (or at a constant cooling rate). For a given solidification problem, dendritic-like patterns can be simulated for different sections of a casting using the growth parameters obtained experimentally or numerically from a macroscopic modelling code. In the present work, a three-dimensional CAFD model, which is similar to the two-dimensional model [12,13], has been used to generate the evolving dendritic structures. A periodic array of dendrites is assumed and a single branched primary arm of a dendrite is isolated for modelling (as in Fig. 1). The computational domain is divided into a three-dimensional array of cubical cells. The cell size (edge length) is chosen such that the domain covers at least 10 secondary dendrite arm branches. Typical values of the cell size are 0.5–5 mm. All the cells are assumed to be liquid initially and a group of cells which are solid, i.e. the nucleus, is placed in an appropriate site at time t =0. In columnar dendritic growth, the left side (x = 0 boundary) is the wall on which the nucleus is placed. It is also assumed that no solute leaves or enters the x= L boundary. A periodic condition is imposed on the top and bottom boundaries (y =0 and y =W) and also on the back and front boundaries (z =0 and z =B). In the case of equiaxed dendritic growth, a periodic condition is imposed on all the boundaries and the nucleus is placed in the centre of the domain.
93
Initially all the cells are assigned a constant value of solute concentration. A temperature field which increases linearly along the growth direction (x) is imposed on the domain for columnar dendritic growth. The time step (Dt) is selected such that the temperature profile is advanced one cell during each step according to the prescribed cooling rate ((T/(t). In the case of equiaxed dendritic growth, temperature is assumed to be constant in the computational domain with the initial value slightly below the liquidus temperature according to the prescribed undercooling. The temperature of the domain is reduced gradually to the solidus temperature according to the prescribed cooling rate and the time step is selected such that there are at least 500 temperature steps during the analysis. Each temperature step includes five subcycles involving solidification, solute rejection and diffusion on a smaller time-step (Dt/5). A liquid cell is assumed to solidify in a given step if the local temperature is below the liquidus temperature (Tliq) corresponding to the composition of that cell if the growth rule is satisfied. The growth rule could be a complex condition based on the curvature of the solid–liquid boundary and a knowledge of related microscopic parameters for the alloy or a simple requirement that a minimum number of neighbouring cells should have already solidified. When all the cells which would solidify in a step are identified, they are marked as solid cells and the excess solute from the new solid sites is rejected into the neighbouring liquid cells. The diffusion solver employs an implicit point-by-point iterative method to obtain the new solute concentration field. Successive over relaxation (SOR) is used to accelerate the convergence of the diffusion step. The analysis is repeated until all the cells become solid or until the composition of the remaining liquid reaches the eutectic composition. The CAFD model is illustrated as a flow chart in Fig. 2. The main advantage of the present approach is that the continuous evolution of the dendrite can be modelled as a function of various growth parameters for a given alloy on a length scale which is sufficiently long to encompass several secondary arms. Formation of the secondary arms occurs naturally and the predicted values of the secondary arm spacing compare with published experimental results [12,13]. If the dendrites (or primary arms in the case of columnar dendrites) are uniformly spaced, the resulting shapes are symmetrical and this can be exploited to speed-up the numerical simulation. Only 1/4 of the dendrite has to be modelled in the case of columnar dendrites (B= W) and 1/8 of an equiaxed dendrite (L= B= W). A typical columnar dendrite simulated for a binary alloy system for a cooling rate of 187 K/s using a grid of 100 × 30× 30 cells having an edge length of 1 mm is given in Fig. 3 at different stages of solidification. The liquidus temperatures were obtained
94
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
Fig. 2. Flow chart of cellular automaton-finite difference model for dendrite growth.
from a database for the aluminium – copper system. The dendrite shapes shown in Fig. 3 are for an alloy with an average composition of 5% copper at 5, 10, 15 and 20% volume fraction solid. These results
were generated using a simple growth rule that at least five of the twenty-six neighbouring cells should be solid. An equiaxed dendrite simulated for the same alloy with a 60× 60× 60 grid is given in Fig. 4 for 5, 15, 15 and 20% solid. The edge length in this case is 2.5 mm which corresponds to a grain size of 300 mm. The growth rule required that at least six of the neighbouring cells should be solid. It should be noted that the Al–5wt.%Cu alloy above has been employed solely for the purpose of generating columnar and equiaxed dendrite-like geometries using the CAFD modelling method. Thereafter, in determining the effective thermal conductivity of the two phase solid+ liquid region as a function of fraction solid and conductivity ratio kl/ks (see next section), it is assumed that as the dendritic structure evolves the values of kl and ks remain constant and uniform. It is appreciated that the latter assumption can only reasonably apply to a pure metal. It is also recognised that pure metals cannot freeze in a columnar dendritic manner with a positive temperature gradient ahead of the advancing interface. However, the purpose of this paper is to demonstrate the application of a numerical approach to the determination of the effective conductivity of columnar and equiaxed type two-phase mushy regions. In the event that 3D computer models of the evolution of alloy dendrites improve, and that experimental knowledge of the variation of thermal conductivity with composition and temperature in specific alloys in both the solid and liquid states becomes available, the method outlined in this paper could be applied to estimate the values of effective thermal conductivities in actual alloys.
Fig. 3. Columnar dendrite generated by cellular automaton-finite difference model (120 ×30 ×30 grid).
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
95
Fig. 4. Equiaxed dendrite generated by cellular automaton-finite difference model (60 ×60 × 60 grid).
4. Numerical results A finite element code has been developed for solving the heat conduction Eq. (6) with either the temperature boundary condition (Eq. (7)) or the flux boundary condition (Eq. (9)). This code has been validated with the available exact/approximate solutions for simple geometry. Comparisons have also been made with the finite difference method for a complex dendrite shape 60× 60 ×60. The results summarised in Table 1 agree within 0.3% for both the methods with temperature/flux prescribed condition at X = 0. Due to the iterative solvers used in the codes, prescribing a temperature boundary condition results in less computer time than prescribing the flux. These results suggest that the finite element method is computationally more efficient for predicting the effective thermal conductivity of the CAFD generated dendrites. In the following paragraphs, numerical results obtained for both columnar and equiaxed dendrite shapes are presented.
As in the CAFD model, symmetry is exploited and only 1/4 of the domain around a dendrite is modelled for thermal conductivity estimation of columnar dendritic zones. Due to symmetry k zeff = k yeff and therefore only two principal thermal conductivity values need to be evaluated, k xeff and k yeff. The liquid thermal conductivity is usually less than that of the solid and therefore computations were carried out for values of Kl =kl/ks less than one only. The effective thermal conductivity Table 1 Predictions of FEM and FDM for temperature and flux boundary conditions Method
Boundary conditions at X =0
keff (60×60×60 CPU time (min) grid)
FEM
Temperature, Eq. (7) Flux Eq. (9) Temperature, Eq. (7) Flux Eq. (9)
0.5219
11.5
0.5218 0.5202
21 55
0.5202
170
FEM FDM FDM
96
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
Fig. 5. Effect of thermal conductivity ratio (k1/ks) on the predicted non-dimensional effective thermal conductivity values for columnar dendrites.
Fig. 6. Variation of effective thermal conductivity as a function of fraction solid for columnar dendrites (k1/ks = 0.5).
Fig. 7. Variation of effective thermal conductivity as a function of fraction solid for columnar dendrites (k1/ks = 0.01).
values are obtained for the dendrite-like structure shown in Fig. 3 for various values of fraction solid ( fs). In order to study the effect of the non-dimensional liquid thermal conductivity (Kl) on the Keff of complex
dendrite shapes, the following values of Kl were used: 0.5, 0.2, 0.1, 0.05 and 0.01. Numerical results for various dendrite shapes having volume fraction solid ( fs) values in the range 10–90% as a function of the liquid thermal conductivity are shown in Fig. 5. The curves for k xeff and k yeff are quite close (or overlap) for values of kl/ks \ 0.3 suggesting that the anisotropic nature of Keff may be ignored in this range. For low and high values of fs, the anisotropic nature may be ignored for all values of kl/ks. The lower the value of kl/ks, the higher is the difference between K xeff and K yeff. It should be noted that K yeff \ K xeff in almost all the cases. In order to compare the predicted effective thermal conductivity with the series and parallel thermal resistance models, the results are also plotted as a function of the volume fraction solid along with the curves for kser and kpar. Figs. 6 and 7 show the variation of normalized values of k xeff and k yeff with respect to fraction solid for kl/ks = 0.5 and 0.01, respectively. k xeff and k yeff are normalized with respect to (ks − kl) so that the normalized value would fall between 0 and 1 for all values of kl/ks. These figures show that the predicted effective thermal conductivity values are closer to kpar than kser. As pointed out earlier, k yeff is greater than k xeff and hence k yeff is closer to kpar. This is quite different to Ref. [8] where it is suggested to use k yeff = kser. These results do not mean that more heat is conducted in the y (or z) direction, but only suggest that the dendritic zones have a capacity to conduct more heat in the y and z direction than in the x direction. The actual amount of heat conducted in a direction depends also on the temperature gradient in that direction, which is obviously greater in the growth (x) direction for columnar dendrites. In the case of equiaxed dendritic regions, only 1/8 of the domain is used for predicting keff. The effective thermal conductivity is isotropic due to symmetry (k xeff = k yeff = k zeff). The predicted results for the dendrite in Fig. 4 are shown in Figs. 8 and 9 for kl/ks =0.5 and 0.01, respectively. The predictions of K xeff for the columnar dendrites are also included in these figures to enable comparisons between the columnar and equiaxed dendritic zones for a given fraction solid. For kl/ks =0.5, the values of k xeff are slightly higher for equiaxed dendrites (maximum about 0.6%) than the columnar dendrites. The values of k xeff of equiaxed dendrites are considerably higher than those of the columnar dendrites when the thermal conductivity of the liquid (kl) is 100 times lower than that of the solid, ks (kl/ks =0.01). The difference is quite high for low values of fraction solid (about 50% for fs = 0.5 and 4% for fs =0.8). Comparisons have also been made between the average of series and parallel thermal conductivities (as recommended in Ref. [7]) and the predicted keff. The magnitude of the maximum error between k xeff and ka = (kser + kpar)/2 is reported in Table 2 for various
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
97
function of kser and kpar for all values of fs and kl/ks. 5. Conclusions
Fig. 8. Variation of effective thermal conductivity as a function of fraction solid for columnar and equiaxed dendrites (k1/ks =0.5).
values of kl/ks. It may be noted that the errors could be quite high when the thermal conductivity of the liquid is considerably smaller than that of the solid. However, ka provides a reasonable estimate of keff for 1]kl/ks ] 0.2. ka gives a higher estimate of keff for small values of fs and kl/ks and therefore the error is negative in these cases. The maximum error does not occur at the same value of fs for the various values of kl/ks considered in this work. This clearly suggests that it is not possible to recommend a simple relation for keff which is only a
The finite element method has been used to predict the effective thermal conductivity of three-dimensional dendritic-like mushy zones as a function of fraction solid when the morphology of the mushy zone is known. The proposed approach is general and can be used to evaluate the effective thermal conductivity of dendrites generated by other numerical models as well. Prescribing the temperature boundary condition on all the walls yields efficient computations. The present work shows that the effective thermal conductivity values are closer to kpar than to kser. In the case of columnar dendrites, the anisotropic nature of thermal conductivity cannot be ignored when the liquid to solid thermal conductivity ratio (kl/ks) is less than 0.2. The results suggest that k yeff \ k xeff for the columnar and equiaxed dendrite shapes used and therefore is closer to kpar. The proposed approach could be used to generate different thermal conductivity curves for a variety of mushy zone configurations. Acknowledgements The authors are grateful to the Engineering and Physical Sciences Research Council (Grant No. GR/ K79147) for financially supporting this research project. References
Fig. 9. Variation of effective thermal conductivity as a function of fraction solid for columnar and equiaxed dendrites (k1/ks =0.01). Table 2 Maximum error between k xeff and the average of kser and kpar kl/ks
0.5 0.2 0.1 0.05 0.01
Error (%) =100×(k xeff−ka )/k xeff Columnar
Equiaxed
2.2 11 20 29 163
2.9 13 22 33 222
[1] P.N. Quested, K.C. Mills, R.F. Brooks, B. Monoghan, A.T. Dinsdale, A. Day, M.J. Richardson, R.J.L. Andon, R. Taylor, H. Szelagowski, in: M. Cross, J. Campbell (Eds.), Modeling of Casting, Welding and Advanced Solidification Processes V11, TMS, Warrendale, PA, 1995, p. 407. [2] W.D. Bennon, F.P. Incropera, Int. J. Heat Mass Transfer 30 (1987) 2161. [3] M. Rappaz, V. Voller, Met. Trans. 21A (1990) 749. [4] S.K. Sinha, T. Sundararajan, V.K. Garg, Int. J. Heat Mass Transfer 35 (1992) 2865. [5] P.W. Emms, A.C. Fowler, J. Fluid Mech. 262 (1994) 111. [6] S. Chang, D.M. Stefanescu, Acta Mater. 44 (1996) 2227. [7] D.R. Poirier, E. McBride, Mater. Sci. Eng. A224 (1997) 48. [8] S.L. Lee, J.H. Yang, Int. J. Heat Mass Transfer 41 (1998) 931. [9] Lauoadi, M. Lacroix, N. Galanis, Int. J. Num. Methods Heat Fluid Flow 8 (1998) 265. [10] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method— Volume 1: Basic Formulation and Linear Problems, 4th edn, McGraw-Hill, New York, 1989. [11] Jennings, J.J. McKeown, Matrix Computation, 2nd edn, Wiley, Chichester, 1992. [12] S.G.R. Brown, J.A. Spittle, in: M. Cross, J. Campbell (Eds.), Modeling of Casting, Welding and Advanced Solidification Processes V11, TMS, Warrendale, PA, 1995, p. 541. [13] S.G.R. Brown, J.A. Spittle, in: B.G. Thomas, C. Beckermann (Eds.), Modeling of Casting, Welding and Advanced Solidification Processes V111, TMS, Warrendale, PA, 1998, p. 179.
Prediction of the effective thermal conductivity of three-dimensional dendritic regions by the finite element method K. Ravindran, S.G.R. Brown, J.A. Spittle * Department of Materials Engineering, Uni6ersity of Wales Swansea, Swansea SA2 8PP, UK Received 5 October 1998; accepted 15 March 1999
Abstract The finite element method has been used to predict the variation in effective thermal conductivity, keff, in evolving three-dimensional dendritic mushy zones. The model demonstrates the influence of a dendrite-like geometry and volume fraction of solid on keff assuming that the conductivities of the solid and liquid phases remain constant and uniform as the structure evolves. The three-dimensional dendrite-like shapes have been generated by a cellular automaton-finite difference model. Numerical solutions are presented for both columnar and equiaxed dendritic zones. © 1999 Elsevier Science S.A. All rights reserved. Keywords: Thermal conductivity; Solidification; Dendrite growth; Numerical modelling
1. Introduction Numerical solution of the macroscopic transport of heat and fluid is becoming a common practice to analyse the various problems arising in metal casting. Significant advances have been made in the simulation of the free surface movement during mould filling, heat transfer during mould filling and the subsequent solidification and the development of thermal stresses. Efforts are also under way to enhance the predictive capability of these macroscopic models with microscopic models based on the understanding of nucleation, growth kinetics and other related phenomena. The ultimate goal is to develop tools which would be able to solve the entire casting process from the mould filling stage to the prediction of the microstructure and the mechanical properties at various sections of the cooled casting. Accurate numerical simulation of solidification transformation in a commercial multi-component alloy, leading to the prediction of microstructures and properties, requires accurate values of thermo-physical properties such as density, thermal conductivity, specific heat, * Corresponding author. Tel.: +44-1792-295589; fax: +44-1792295693. E-mail address: [email protected] (J.A. Spittle)
etc. Efforts initiated by the Department of Trade and Industry (DTI), UK, reveal some of the difficulties associated with measuring these properties in metals/alloys at high temperature. The measured property values for aluminium alloys showed a marked dependence on the thermal and mechanical histories of the samples and this could be due to the changes in the microstructure. It is rather difficult to isolate the latent heat effect in order to measure the thermal conductivity in the mushy stage accurately and hence efforts are being made to derive the thermal conductivity values from electrical conductivity measurements [1]. However, experimental measurement of all the phenomena accompanying solidification, as a function of time and varying fraction solid, within a continuously evolving dendrite structure is impossible. A numerical approach is, therefore, the only viable alternative way for predicting the thermal conductivity in such cases. Thermal conductivity of the mushy zone is a function of the morphology of the two phase region, the thermal conductivities of the solid and liquid (ks and kl), and the volume fraction of solid fs. Analytical expressions may be derived when some assumptions are made about the distribution of the solid and liquid in the mushy zone. The thermal resistance (which is inversely proportional to the thermal conductivity) of the solid and liquid
0921-5093/99/$ - see front matter © 1999 Elsevier Science S.A. All rights reserved. PII: S 0 9 2 1 - 5 0 9 3 ( 9 9 ) 0 0 1 6 0 - 4
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
regions may be lumped in such a way that the effective thermal resistance is either in ‘series’ or in ‘parallel’. The effective thermal conductivity for these two cases, kser and kpar, can be derived using the analogy of an equivalent electrical circuit:
kser =
 kskl à fskl + (1− fs)ks Ì
(1)
kpar =fsks + (1− fs)klà Škser and kpar are also known as the Wiener bounds. For any complex geometry, the effective thermal conductivity in a given direction should fall between these two limiting values. kpar gives the upper bound and kser the lower bound. In macroscopic alloy solidification analysis, due to the absence of other data, the thermal conductivity of the mushy zone is usually assumed to be a phase-average varying linearly with the volume fraction-solid, fs. This is also known as the ‘mixture model’ for effective thermal conductivity [2 – 6]: kmix = ks fs +kl(1− fs)
(2)
The value given above is the same as the equivalent thermal conductivity for a parallel circuit (kpar) and, as discussed earlier, this is an upper bound for keff and the actual value of keff could be anywhere between kpar and the lower bound, kser, depending on the morphology of the dendritic zone. In the case of columnar dendrites, keff depends on the direction as well, i.e. the effective thermal conductivity is anisotropic. Poirier and McBride [7] recommend the average of kpar and kser for keff. Since the conductivities in the solid and liquid phases of specific alloys are usually unknown, modellers frequently have to resort to using the solid and liquid conductivities of the pure solvent metal. Lee and Yang [8] have recently predicted the effective thermal conductivity of the columnar dendritic zone assuming that the dendrites are wedge-shaped. They predicted the thermal conductivity of the mushy zone along the growth direction (k xeff) as a function of fs using a steady state temperature field which has been averaged in the normal direction. This method is not
91
valid for a solidification analysis since the mushy zone usually encompasses many dendrites and the averaging has to be done for the whole dendrite and not just in one direction. Lee and Yang do not consider the effect of secondary arms and they suggested the use of kser as the thermal conductivity in the direction normal to growth (k yeff). Laouadi et al. [9] have added the effect of latent heat into a one-dimensional model for the effective thermal conductivity. None of these models considers the effect of the actual dendrite morphology on the effective thermal conductivity. In Section 2, we propose a general numerical approach for the prediction of thermal conductivity of the mushy zone when the dendrite morphology is known. In Section 3, a three-dimensional cellular automatonfinite difference model which is used for generating the dendrite shapes is briefly described. Finally numerical results obtained for generated dendritic structures are presented in Section 4.
2. Mathematical formulation In the present work, a periodic array of dendrites is assumed. The domain isolated for thermal-conductivity evaluation, a parallelepiped of dimensions L × W × B, encloses a single primary arm that has branched as shown in Fig. 1. The size of the computational domain is assumed to be much smaller than that of the domain where the effective property would be used for a macroscopic solidification analysis. A one-dimensional temperature field is imposed in a desired direction (say, x) and the steady state temperature field in the computational domain is obtained numerically by the finite element method. The unknown effective thermal conductivity k xeff of the isolated domain in the desired direction is evaluated by imposing the condition that the steady state heat flux is the same if the domain were filled with a medium having a constant thermal conductivity equal to k xeff. A similar formulation can be used for predicting the effective thermal conductivity in the y and z directions. The heat conduction equation and the boundary conditions prescribed for evaluating k xeff are given below: 9·k9T = 0
(3)
Boundary conditions:
 à T =T2 at x= L Ì (T = 0 at y=0, y=W, z =0 and z= B à (n Å T= T1 at x= 0
Fig. 1. Computational domain enclosing a branched primary arm of a columnar dendrite.
(4)
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
92
where n is the outward normal to the boundary of the computational diagram. The above equations may be nondimensionalised using the non-dimensional variables u, K, X, Y and Z as defined below: T−T2 k  K= ks à T1 −T2 Ì x y z X = Y= Z= à L L LÅ
u=
(5)
In the above, the length of the computational domain in the x direction (L) is used as the reference length and the thermal conductivity of the solid (ks) as the reference value for thermal conductivity. The non-dimensional form of the heat conduction equation and the boundary conditions are given below: 9·K9u =0
(6)
u=1 at X=0
 u=0 at X=1 Ã Ì (u à =0 at Y=0, Y=W/L, Z=0 and Z=B/L Å (n
(7)
With the nondimensionalisation defined above, the computed average non-dimensional wall flux is equal to the non-dimensional effective thermal conductivity (K xeff): k xeff L 2 = ks BW
(u q= − K = 1 at X= 0 (n
 à à u= 0 at X= 1 Ì Ã (u =0 at Y=0, Y=W/L, Z=0 and Z=B/L à (n Å
&& n (u (X
K
dY dZ
(8)
X=0
The above equation shows that the predicted effective thermal conductivity (k xeff) is independent of the temperature values (T1 and T2) prescribed on the boundary. It depends on the aspect ratio (L/W and L/B) of the domain isolated for averaging and not on the actual size of the domain (L, W and B). The effective thermal conductivity is the same for two dendrites having the same shape, but different sizes as long as the averaging domain is much smaller in size compared to the macroscopic domain where the averaged property would be used. The non-dimensional effective thermal conductivity (K xeff) does not depend directly on the thermal conductivity values of the solid or liquid (kl or ks) but on the thermal conductivity ratio (kl/ks). Therefore, numerical solutions can be obtained in terms of kl/ks and not on the actual values kl and ks which have to be provided by some other means such as from a database or an experiment. The effective thermal conductivity k xeff may also be obtained by prescribing a unit heat flux (non-dimen-
(9)
In this case, k xeff is evaluated by imposing the condition that for the same value of prescribed wall heat flux (q), the predicted average temperature at the X=0 wall should be the same if the domain were filled with a medium having a constant thermal conductivity equal to k xeff. The computed average wall temperature u( (X = 0) is related to k xeff by the following relation: 1 k 1 = xs = x K eff k eff u( (X = 0)
Boundary conditions:
K xeff =
sional) on the X=0 wall. The boundary conditions for this approach are given below:
(10)
The advantage of this alternate approach is that only the average of the predicted wall temperature is needed and not the temperature gradient which is not continuous in the nodes/element boundaries. Further the first derivatives of the variables are known to be accurate at the numerical (Gaussian) integration points and not on element boundaries. The heat conduction Eq. (6) is solved by the Galerkin finite element method. Eight node isoparametric brick elements have been used for discretizing the computational domain. The discrete form of the final equations may be written in matrix form as [10]: [A]{X}= %e [A e]{X}= {F}
(11)
where [A e] is the local stiffness matrix for cell e and the summation has to be performed over all the cells to form the global stiffness matrix, [A]. {X} is the array of unknown nodal temperature values and {F} is the nodal external heat flux vector which has nonzero values only for the nodes which fall on the boundary edges with a nonzero heat flux boundary condition. The stiffness matrix of a cell with a thermal conductivity k e may be obtained by numerical integration using the following relation: [A e]= k e
&&&
(9N)T9N dx dy dz
(12)
where 9N are the derivatives of the element interpolation (or shape) function. It may be noted from the above equation that the local stiffness matrix ([A e]) is the same for all cells but for a multiplying factor, the element thermal conductivity (k e). The cellular au-
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
tomata domains have a large number of cells and the storage requirements for the assembled global matrix ([A]) are quite high. As the discrete equations could not be solved by in-core solvers, an iterative solution approach based on the conjugate gradient method has been used for the solution of nodal unknowns [11]. No pre-conditioning is used as it involves an additional storage overhead of a large multi-dimensional array of the pre-conditioned element matrices. Also, there is no need to generate and maintain a three-dimensional matrix of cell nodal connectivity as the cells and their nodes can be easily addressed with a finite differencelike notation with three subscripts, i, j and k. Only a single element stiffness matrix (8 × 8) is stored and the solutions for a problem with 60×60 × 60 cells could be obtained within 12 min on a 200 MHz Pentium PC.
3. Cellular automaton-finite difference (CAFD) model for dendrite growth Modelling dendritic growth is complex as it requires a coupled analysis of the transport of heat and mass. Latent heat release during the phase change process and the need to consider other microscopic factors such as nucleation and growth make it a further challenge. In the CAFD model proposed by Brown and Spittle [12,13], the heat transfer phenomena is decoupled by assuming that the growth occurs under steady conditions (or at a constant cooling rate). For a given solidification problem, dendritic-like patterns can be simulated for different sections of a casting using the growth parameters obtained experimentally or numerically from a macroscopic modelling code. In the present work, a three-dimensional CAFD model, which is similar to the two-dimensional model [12,13], has been used to generate the evolving dendritic structures. A periodic array of dendrites is assumed and a single branched primary arm of a dendrite is isolated for modelling (as in Fig. 1). The computational domain is divided into a three-dimensional array of cubical cells. The cell size (edge length) is chosen such that the domain covers at least 10 secondary dendrite arm branches. Typical values of the cell size are 0.5–5 mm. All the cells are assumed to be liquid initially and a group of cells which are solid, i.e. the nucleus, is placed in an appropriate site at time t =0. In columnar dendritic growth, the left side (x = 0 boundary) is the wall on which the nucleus is placed. It is also assumed that no solute leaves or enters the x= L boundary. A periodic condition is imposed on the top and bottom boundaries (y =0 and y =W) and also on the back and front boundaries (z =0 and z =B). In the case of equiaxed dendritic growth, a periodic condition is imposed on all the boundaries and the nucleus is placed in the centre of the domain.
93
Initially all the cells are assigned a constant value of solute concentration. A temperature field which increases linearly along the growth direction (x) is imposed on the domain for columnar dendritic growth. The time step (Dt) is selected such that the temperature profile is advanced one cell during each step according to the prescribed cooling rate ((T/(t). In the case of equiaxed dendritic growth, temperature is assumed to be constant in the computational domain with the initial value slightly below the liquidus temperature according to the prescribed undercooling. The temperature of the domain is reduced gradually to the solidus temperature according to the prescribed cooling rate and the time step is selected such that there are at least 500 temperature steps during the analysis. Each temperature step includes five subcycles involving solidification, solute rejection and diffusion on a smaller time-step (Dt/5). A liquid cell is assumed to solidify in a given step if the local temperature is below the liquidus temperature (Tliq) corresponding to the composition of that cell if the growth rule is satisfied. The growth rule could be a complex condition based on the curvature of the solid–liquid boundary and a knowledge of related microscopic parameters for the alloy or a simple requirement that a minimum number of neighbouring cells should have already solidified. When all the cells which would solidify in a step are identified, they are marked as solid cells and the excess solute from the new solid sites is rejected into the neighbouring liquid cells. The diffusion solver employs an implicit point-by-point iterative method to obtain the new solute concentration field. Successive over relaxation (SOR) is used to accelerate the convergence of the diffusion step. The analysis is repeated until all the cells become solid or until the composition of the remaining liquid reaches the eutectic composition. The CAFD model is illustrated as a flow chart in Fig. 2. The main advantage of the present approach is that the continuous evolution of the dendrite can be modelled as a function of various growth parameters for a given alloy on a length scale which is sufficiently long to encompass several secondary arms. Formation of the secondary arms occurs naturally and the predicted values of the secondary arm spacing compare with published experimental results [12,13]. If the dendrites (or primary arms in the case of columnar dendrites) are uniformly spaced, the resulting shapes are symmetrical and this can be exploited to speed-up the numerical simulation. Only 1/4 of the dendrite has to be modelled in the case of columnar dendrites (B= W) and 1/8 of an equiaxed dendrite (L= B= W). A typical columnar dendrite simulated for a binary alloy system for a cooling rate of 187 K/s using a grid of 100 × 30× 30 cells having an edge length of 1 mm is given in Fig. 3 at different stages of solidification. The liquidus temperatures were obtained
94
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
Fig. 2. Flow chart of cellular automaton-finite difference model for dendrite growth.
from a database for the aluminium – copper system. The dendrite shapes shown in Fig. 3 are for an alloy with an average composition of 5% copper at 5, 10, 15 and 20% volume fraction solid. These results
were generated using a simple growth rule that at least five of the twenty-six neighbouring cells should be solid. An equiaxed dendrite simulated for the same alloy with a 60× 60× 60 grid is given in Fig. 4 for 5, 15, 15 and 20% solid. The edge length in this case is 2.5 mm which corresponds to a grain size of 300 mm. The growth rule required that at least six of the neighbouring cells should be solid. It should be noted that the Al–5wt.%Cu alloy above has been employed solely for the purpose of generating columnar and equiaxed dendrite-like geometries using the CAFD modelling method. Thereafter, in determining the effective thermal conductivity of the two phase solid+ liquid region as a function of fraction solid and conductivity ratio kl/ks (see next section), it is assumed that as the dendritic structure evolves the values of kl and ks remain constant and uniform. It is appreciated that the latter assumption can only reasonably apply to a pure metal. It is also recognised that pure metals cannot freeze in a columnar dendritic manner with a positive temperature gradient ahead of the advancing interface. However, the purpose of this paper is to demonstrate the application of a numerical approach to the determination of the effective conductivity of columnar and equiaxed type two-phase mushy regions. In the event that 3D computer models of the evolution of alloy dendrites improve, and that experimental knowledge of the variation of thermal conductivity with composition and temperature in specific alloys in both the solid and liquid states becomes available, the method outlined in this paper could be applied to estimate the values of effective thermal conductivities in actual alloys.
Fig. 3. Columnar dendrite generated by cellular automaton-finite difference model (120 ×30 ×30 grid).
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
95
Fig. 4. Equiaxed dendrite generated by cellular automaton-finite difference model (60 ×60 × 60 grid).
4. Numerical results A finite element code has been developed for solving the heat conduction Eq. (6) with either the temperature boundary condition (Eq. (7)) or the flux boundary condition (Eq. (9)). This code has been validated with the available exact/approximate solutions for simple geometry. Comparisons have also been made with the finite difference method for a complex dendrite shape 60× 60 ×60. The results summarised in Table 1 agree within 0.3% for both the methods with temperature/flux prescribed condition at X = 0. Due to the iterative solvers used in the codes, prescribing a temperature boundary condition results in less computer time than prescribing the flux. These results suggest that the finite element method is computationally more efficient for predicting the effective thermal conductivity of the CAFD generated dendrites. In the following paragraphs, numerical results obtained for both columnar and equiaxed dendrite shapes are presented.
As in the CAFD model, symmetry is exploited and only 1/4 of the domain around a dendrite is modelled for thermal conductivity estimation of columnar dendritic zones. Due to symmetry k zeff = k yeff and therefore only two principal thermal conductivity values need to be evaluated, k xeff and k yeff. The liquid thermal conductivity is usually less than that of the solid and therefore computations were carried out for values of Kl =kl/ks less than one only. The effective thermal conductivity Table 1 Predictions of FEM and FDM for temperature and flux boundary conditions Method
Boundary conditions at X =0
keff (60×60×60 CPU time (min) grid)
FEM
Temperature, Eq. (7) Flux Eq. (9) Temperature, Eq. (7) Flux Eq. (9)
0.5219
11.5
0.5218 0.5202
21 55
0.5202
170
FEM FDM FDM
96
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
Fig. 5. Effect of thermal conductivity ratio (k1/ks) on the predicted non-dimensional effective thermal conductivity values for columnar dendrites.
Fig. 6. Variation of effective thermal conductivity as a function of fraction solid for columnar dendrites (k1/ks = 0.5).
Fig. 7. Variation of effective thermal conductivity as a function of fraction solid for columnar dendrites (k1/ks = 0.01).
values are obtained for the dendrite-like structure shown in Fig. 3 for various values of fraction solid ( fs). In order to study the effect of the non-dimensional liquid thermal conductivity (Kl) on the Keff of complex
dendrite shapes, the following values of Kl were used: 0.5, 0.2, 0.1, 0.05 and 0.01. Numerical results for various dendrite shapes having volume fraction solid ( fs) values in the range 10–90% as a function of the liquid thermal conductivity are shown in Fig. 5. The curves for k xeff and k yeff are quite close (or overlap) for values of kl/ks \ 0.3 suggesting that the anisotropic nature of Keff may be ignored in this range. For low and high values of fs, the anisotropic nature may be ignored for all values of kl/ks. The lower the value of kl/ks, the higher is the difference between K xeff and K yeff. It should be noted that K yeff \ K xeff in almost all the cases. In order to compare the predicted effective thermal conductivity with the series and parallel thermal resistance models, the results are also plotted as a function of the volume fraction solid along with the curves for kser and kpar. Figs. 6 and 7 show the variation of normalized values of k xeff and k yeff with respect to fraction solid for kl/ks = 0.5 and 0.01, respectively. k xeff and k yeff are normalized with respect to (ks − kl) so that the normalized value would fall between 0 and 1 for all values of kl/ks. These figures show that the predicted effective thermal conductivity values are closer to kpar than kser. As pointed out earlier, k yeff is greater than k xeff and hence k yeff is closer to kpar. This is quite different to Ref. [8] where it is suggested to use k yeff = kser. These results do not mean that more heat is conducted in the y (or z) direction, but only suggest that the dendritic zones have a capacity to conduct more heat in the y and z direction than in the x direction. The actual amount of heat conducted in a direction depends also on the temperature gradient in that direction, which is obviously greater in the growth (x) direction for columnar dendrites. In the case of equiaxed dendritic regions, only 1/8 of the domain is used for predicting keff. The effective thermal conductivity is isotropic due to symmetry (k xeff = k yeff = k zeff). The predicted results for the dendrite in Fig. 4 are shown in Figs. 8 and 9 for kl/ks =0.5 and 0.01, respectively. The predictions of K xeff for the columnar dendrites are also included in these figures to enable comparisons between the columnar and equiaxed dendritic zones for a given fraction solid. For kl/ks =0.5, the values of k xeff are slightly higher for equiaxed dendrites (maximum about 0.6%) than the columnar dendrites. The values of k xeff of equiaxed dendrites are considerably higher than those of the columnar dendrites when the thermal conductivity of the liquid (kl) is 100 times lower than that of the solid, ks (kl/ks =0.01). The difference is quite high for low values of fraction solid (about 50% for fs = 0.5 and 4% for fs =0.8). Comparisons have also been made between the average of series and parallel thermal conductivities (as recommended in Ref. [7]) and the predicted keff. The magnitude of the maximum error between k xeff and ka = (kser + kpar)/2 is reported in Table 2 for various
K. Ra6indran et al. / Materials Science and Engineering A269 (1999) 90–97
97
function of kser and kpar for all values of fs and kl/ks. 5. Conclusions
Fig. 8. Variation of effective thermal conductivity as a function of fraction solid for columnar and equiaxed dendrites (k1/ks =0.5).
values of kl/ks. It may be noted that the errors could be quite high when the thermal conductivity of the liquid is considerably smaller than that of the solid. However, ka provides a reasonable estimate of keff for 1]kl/ks ] 0.2. ka gives a higher estimate of keff for small values of fs and kl/ks and therefore the error is negative in these cases. The maximum error does not occur at the same value of fs for the various values of kl/ks considered in this work. This clearly suggests that it is not possible to recommend a simple relation for keff which is only a
The finite element method has been used to predict the effective thermal conductivity of three-dimensional dendritic-like mushy zones as a function of fraction solid when the morphology of the mushy zone is known. The proposed approach is general and can be used to evaluate the effective thermal conductivity of dendrites generated by other numerical models as well. Prescribing the temperature boundary condition on all the walls yields efficient computations. The present work shows that the effective thermal conductivity values are closer to kpar than to kser. In the case of columnar dendrites, the anisotropic nature of thermal conductivity cannot be ignored when the liquid to solid thermal conductivity ratio (kl/ks) is less than 0.2. The results suggest that k yeff \ k xeff for the columnar and equiaxed dendrite shapes used and therefore is closer to kpar. The proposed approach could be used to generate different thermal conductivity curves for a variety of mushy zone configurations. Acknowledgements The authors are grateful to the Engineering and Physical Sciences Research Council (Grant No. GR/ K79147) for financially supporting this research project. References
Fig. 9. Variation of effective thermal conductivity as a function of fraction solid for columnar and equiaxed dendrites (k1/ks =0.01). Table 2 Maximum error between k xeff and the average of kser and kpar kl/ks
0.5 0.2 0.1 0.05 0.01
Error (%) =100×(k xeff−ka )/k xeff Columnar
Equiaxed
2.2 11 20 29 163
2.9 13 22 33 222
[1] P.N. Quested, K.C. Mills, R.F. Brooks, B. Monoghan, A.T. Dinsdale, A. Day, M.J. Richardson, R.J.L. Andon, R. Taylor, H. Szelagowski, in: M. Cross, J. Campbell (Eds.), Modeling of Casting, Welding and Advanced Solidification Processes V11, TMS, Warrendale, PA, 1995, p. 407. [2] W.D. Bennon, F.P. Incropera, Int. J. Heat Mass Transfer 30 (1987) 2161. [3] M. Rappaz, V. Voller, Met. Trans. 21A (1990) 749. [4] S.K. Sinha, T. Sundararajan, V.K. Garg, Int. J. Heat Mass Transfer 35 (1992) 2865. [5] P.W. Emms, A.C. Fowler, J. Fluid Mech. 262 (1994) 111. [6] S. Chang, D.M. Stefanescu, Acta Mater. 44 (1996) 2227. [7] D.R. Poirier, E. McBride, Mater. Sci. Eng. A224 (1997) 48. [8] S.L. Lee, J.H. Yang, Int. J. Heat Mass Transfer 41 (1998) 931. [9] Lauoadi, M. Lacroix, N. Galanis, Int. J. Num. Methods Heat Fluid Flow 8 (1998) 265. [10] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method— Volume 1: Basic Formulation and Linear Problems, 4th edn, McGraw-Hill, New York, 1989. [11] Jennings, J.J. McKeown, Matrix Computation, 2nd edn, Wiley, Chichester, 1992. [12] S.G.R. Brown, J.A. Spittle, in: M. Cross, J. Campbell (Eds.), Modeling of Casting, Welding and Advanced Solidification Processes V11, TMS, Warrendale, PA, 1995, p. 541. [13] S.G.R. Brown, J.A. Spittle, in: B.G. Thomas, C. Beckermann (Eds.), Modeling of Casting, Welding and Advanced Solidification Processes V111, TMS, Warrendale, PA, 1998, p. 179.