A mathematical model for surface segregation in aluminum direct chill casting

A mathematical model for surface segregation in aluminum direct chill casting

Thevik, Havard J

A two-dimensional mathematical model for the development of macrosegregation at and close to the ingot surface during direct chill (DC) casting of aluminum rolling sheet ingots is presented. The model accounts for macrosegregation caused by exudation of interdendritic melt and macrosegregation associated with solidification shrinkage. Equations for the conservation of energy, solute, momentum, and mass during the stationary phase of the process are solved numerically by a finiteelement method. The solution domain corresponds to a vertical cross section at the middle of the longest side of the slab. The main simplifications in the modeling concept are to assume that the solid in the mushy zone moves with the casting speed, and that the alloy is binary and solidifies according to the lever rule. The thickness and solute concentration of the surface layer and the macrosegregation close to the surface are calculated, and modeling results are compared to measurements on full-scale castings.


THE surface quality of direct-chilled (DC) cast aluminium rolling sheet ingots is often reduced by a segregated layer of exudations (Figure 1). For example, on AA5182 ingots with a nominal magnesium concentration of 4.7 pct cast with the spout-and-pin regulating technology, the layer can be more than 1-mm thick and have an average magnesium concentration of up to about 10 pct.[‘] The removal of the surface layer before further processing of the ingot entails high costs.

The exudated layer is caused by interdendritic melt flow through a partly solidified (mushy) and remelting shell close to the mold. This remelting and the resulting metal flow are due to an air gap between the shell and the mold caused by the global solidification contraction of the ingot[26] (Figure 2). The metallostatic head then forces interdendritic liquid through the mushy shell toward the ingot surface. In addition to exudation, the solidification shrinkage also leads to transport of interdendritic liquid.[7,8,9] While the latter flow phenomenon leads to positive macrosegregation (i.e., the solute concentration becomes larger than the nominal alloying concentration) close to the ingot surface, exudation results in negative segregation in the zone through which the interdendritic melt flow has taken place.[5’0 In the AA5182 example referred to previously, it should be noted that the negative macrosegregation in the solute-depleted zone caused by the exudation is much more severe than the positive contribution associated with the solidification shrinkage. On the other hand, when the exudation is very small, the solidification shrinkage becomes the leading mechanism behind the macrosegregation formation close to the surface.[11]

This study is directed toward the mathematical modeling of the macrosegregation formation associated with exudation and solidification shrinkage. This segregation will be referred to as surface segregation. Former macrosegregation studies focusing on either exudation or solidification shrinkage can be found in References 5 and 10, and 8 and 9, respectively. Haug et al.[11] carried out a one-dimensional study in which both macrosegregation mechanisms were addressed, and a case study relevant for aluminum DC casting was incorporated.

The purpose of the present article is to present an extension of the modeling concept in Reference 11 to a twodimensional situation relevant to the DC casting process. Section II describes briefly the physical basis of the model, the governing equations, and their numerical solution. In Section III, model predictions are compared to surface segregation measurements on full-scale castings, and weak points in the model are discussed.


A half-vertical cross section perpendicular to the wide side of the rolling slab, starting at the ingot center, constitutes the rectangular two-dimensional solution domain (Figure 2). The Cartesian coordinates are denoted by x and y, respectively. Only the stationary part of the casting process is considered. This means that solidified metal leaves the solution domain at y = 0 with the casting speed V^sub c^. It should be noted that the total volume flux through the boundary y = h is somewhat larger than the flux through the boundary y = 0, because of the solidification shrinkage and the exudation.

A. Governing Equations

The governing equations are based upon the general volume-averaged conservation equations discussed in Reference 12. The following set of assumptions was discussed by Haug et al.[11] and is invoked to simplify the general equations.

(1) The solid phase is rigid and moves downward with the casting velocity.

(2) The momentum transfer between the liquid and solid phases in the solidification interval is modeled by a Darcian term. The associated mushy zone permeability is isotropic and related to the liquid fraction by the Kozeny-Carman relation.

(3) There is no pore formation.

(4) There is no macroscopic diffusion of solute.

(5) The dispersion fluxes are set to zero.

(6) The thermophysical properties are constant within each phase.

(7) The alloy is binary.

(8) There is local thermodynamical equilibrium in the mushy zone (lever rule), and the liquidus and solidus lines in the thermodynamic equilibrium-phase diagram are linear.

Here, t, p, h, lambda, T, c, mu, p, and g, denote the time, density, enthalpy per unit mass, heat conductivity, temperature, solute concentration, dynamic viscosity in the liquid, pressure, and acceleration due to gravity, respectively. The subscripts I and s indicate quantities associated with the liquid and solid phases, respectively, and the superscript tr stands for transpose of a tensor. Although only the stationary solution is of interest here, the condition (partial d)/(partial d)t = 0 is not introduced in the equations. This is due to the numerical solution method discussed in Section II-.

B. Boundary Conditions

The height of the solution domain is chosen sufficiently large such that a zero diffusive energy flux can be imposed along the boundary y = 0. Along this boundary, the ingot is completely solidified and, thus, Nu = 0. At y = h, liquid melt enters at a given casting temperature (T^sub c^). The pressure at y = h is equal to the atmospheric pressure, and the volume-averaged viscous stress tensor tau = (mu)[inverted delta u + (inverted deltau)^sup w^] is zero. Along the symmetry line x = 0, the horizontal component of the temperature gradient, the normal velocity, and the tangential viscous stress all equal zero.

As indicated in Figure 2, the boundary along the vertical surface, x = b, is subdivided into five regions: a secondary cooling region with “low” heat transfer (y element of [0, ysub 1^]); a secondary cooling region with “high” heat transfer, reflecting boiling in the cooling water film just below to the water impingement point (y element [y^sub 1^ y^sub 2^]); the air gap region (y element [y^sub 2^, y^sub 3^]); the primary cooling region (y element [y^sub 3^, y^sub 4^]); and a region reflecting a free surface meniscus or a hot-top (y element [y^sub 4^, h]). The latter region avoids solidification at the inflow boundary. The heat-transfer coefficients in these regions are denoted by alpha^sub 1^ alpha^sub 2^, alpha^sub 3^, alpha^sub 4^, and alpha^sub 5^, respectively. For simplicity, the outside reference temperature is, in the case studies in Section III, given the same value of T^sub 0^ in all regions. Concerning the mechanical boundary conditions, there are two regions: the air gap and secondary cooling regions (y element [0, y^sub 3^]), in which the pressure is equal to the atmospheric pressure and the viscous stress is zero, and the contact region y element [y^sub 3^, h], where v = 0.

Notice that the model does not include any mechanism for determining where the ingot contracts from the mold. This means that the position of the air gap (if any) must be defined a priori.

C. Numerical Solution

The numerical analysis and solution is based upon a dimensionless version of the governing equations. The new equations have the same structure as their dimensional counterparts.

The spatial discretization is based upon a weighted residual method using the finite-element grid with rectangular elements (Figure 3). While the Galerkin formulation is employed for the energy, mass, and momentum equations, weight functions for the solute equation are chosen according to the streamline upwind Petrov-Galerkin formulation [13] in order to avoid spatial oscillations in the solution. Piecewise bilinear polynomials are used as nodal basis functions for the primary unknowns pc, ph, and v.

In evaluating a, the values for T and ph from the previous iteration are employed.

The discretization described previously yields a set of linear equations for the energy, solute, and combined mass and momentum equations. The system of linear equations originating from the energy equation has a sparse and nonsymmetric coefficient matrix and is solved by a Bi-cgstab method.[16] The linear system of equations obtained from the solute equation (Eq. [4]) and the combined momentum and mass equations (Eq. [19]) are solved by Gauss elimination.

The stationary solution is obtained by an implicit time integration scheme. At each time step, a loop in which the governing equations are solved sequentially is performed in order to obtain the overall solution.

(1) The value of ph is calculated from Eq. [3].

(2) The value of pc is calculated from Eq. [4].

(3) The values of T, g^sub 1^,c^sub 1^ h^sub 1^ and p are updated by using Eqs. [7] through [15].

(4) The value of v is calculated from Eq. [19].

Values from the most recent iteration or time step, if the previous loop is performed only once per time step, are used for h^sub l^, h, c^sub i^, c, p, and v in the last term on the lefthand side, as well as in the last term on the right-hand side, of Eqs. [3] and [4] (notice that ph and pc are regarded as primary unknowns in the convective terms on the left-hand side of these equations).

Initially, only the energy equation is solved (step 1 in the iteration loop above) in order to obtain a reasonable temperature field. The start value for the temperature is then 757 K, and c and u are set to c^sub 0^ and zero, respectively. After about 600 seconds, the temperature is nearly stationary, and all the four steps in the iteration loop are invoked. In the modeling cases discussed in Section III, the time step was 0.5 seconds initially and 0.3 seconds when all equations were solved. The iteration loop was performed once for each time step. The calculation was terminated when stationary conditions were obtained in the solution domain (after approximately 1000 seconds). The global conservation of mass and solute were checked, and it was found that |(J^sub in^ – J^sub out^)/J^sub in^|

The implementation of the numerical procedure is based upon the object-oriented C+ + class library Diffpack.[17,18] The typical central processing unit time required the case studies presented in Section III is about 12 hours on a HEWLETT-PACKARD Work Station 9000/770 (based on a grid with 7130 nodes, Figure 3). Approximately 60 pct of the computation time can be associated with the direct method (Gaussian elimination) required to solve the linear system of equations for the velocity.

D. Verification of the Numerical Solution

The numerical implementation has been compared to the one-dimensional model of Haug et al.,[11] in which a finitedifference, control-volume based discretization was employed. As a validation case, an Al-4.8 pct Mg alloy with an initial temperature of 950 K was cooled unidirectionally. The material data used for this alloy are given in Table I. At the cold end of the alloy, the heat-transfer coefficient is 300 W/m^sup 2^K and the ambient temperature is 303 K, while the hot end is isolated. The solution domain is taken to be 0.1 m. The final inverse segregation profile close to the chill, for the present model, and that of Reference 11 are displayed in Figure 4(a). The solute concentration profiles are nearly identical, and do also correspond well to the analytically calculated surface value.[19]

The present numerical implementation has also been tested in a case with exudation. The case is similar to the one-dimensional study previously mentioned, the only differences being that alpha = 1000 W/m^sup 2^K for 0 = 27 s (alpha is then 300 W/m^sup 2^K). The pressure difference over the solution domain, which is the driving force for the exudation, corresponds to a metallostatic height of 5 cm. As shown in Figure 4(b), the correspondence between the the present numerical code and that of Haug et al.[11] is very good. The exudated layer is not included in Figure 4(b), but both the thickness and solute concentration of the layer are (nearly) identical for the two implementations. Notice that Haug et al. modeled the momentum balance by Darcy’s law, while the full momentum equation is employed in the present model. Therefore, the good correspondence between the two models also reflects the fact that inertial terms and viscous stresses can be neglected in the momentum balance when predicting exudation for slowly solidifying, unidirectional model problems (cf also the scaling analysis of Krane and Incropera).[20]


In order to illustrate the capabilities of the present model, predictions will be compared to experimental data from the casting trials of rolling sheet ingots reported in Reference.[21] The alloy used in these casting trials was a nearly binary Al-3 pct Mg alloy, and the solute concentration was measured from the surface to the center of the ingot. The solute concentration in the surface region for “ingot C” is shown in Figure 5. The first 1 to 2 mm within the as-cast ingot surface have a magnesium concentration substantially higher than the nominal value. Then, there is a region of approximately 3 to 5 mm in which the solute concentration is less than 3 pct Mg. Beyond this solute-depleted region, the concentration levels out at the nominal concentration. During the stationary period of the casting trials, the temperatures were logged close to the surface within the ingot. These measurements are shown in Figure 6.

In the modeling of the casting experiment associated with the solute concentration in Figure 5, the material and process data summarized in Table I are used. The grid employed in this calculation is shown in Figure 3. It contains 7130 nodes and is refined in the mold region in order to adequately resolve the fluid-flow field within the thin mushy zone.

In determining the heat-transfer coefficients listed in Table I, the temperature measurements close to the surface have been used. This choice of heat-transfer coefficients gives a reasonable correspondence between measured and calculated temperatures (Figure 6).

Using the data of Table I, the mathematical model predicts the solute concentration distribution referred to as “case 1”. For this case, the calculated melt flow close to the mold is shown in Figure 7. The measured (mean) thickness and solute concentration of the exudated layer are compared to the predicted values (cf Eqs. [16] and [17]) in Table II. Furthermore, the solute concentrations within the ingot are displayed in Figure 8. In this figure, x = 0 is the position of the initial casting surface, i.e., the vertical, right boundary in the model. In order to facilitate a comparison between the measured and calculated values, this position corresponds to the transition between the depleted zone and the exudated layer in the measurements, i.e., x = 0 in Figure 8 corresponds to x == 1.9 mm in Figure 5.

Qualitatively, the model is able to reproduce the observed solute concentration profiles. Both a highly segregated surface layer and a solute-depleted zone within the ingot are predicted. Quantitatively, however, it is seen from Table II that the model calculates a somewhat thinner exudated layer than measured, while the average concentration is overestimated. Furthermore, the calculated solute-depleted zone extends over a larger region than the measured solute depletion (Figure 8). The quantitative discrepancy is believed to be mainly due to the quite simple treatment of the boundary conditions in the mold region. Although a contact zone with a constant (high) heat-transfer coefficient and an air-gap zone in which the heat-transfer coefficient is also constant (and low) seem to reflect the exudation mechanism qualitatively, the real heat-transfer conditions are probably more complicated. For the melt flow, the transition between no exudation and exudation is assumed to be at the same position as the change in thermal boundary conditions. However, since the ingot surface is rather rough, it might be such that exudation occurs while there are still some regions of the surface in good (thermal) contact with the mold.

The model predictions are, furthermore, very sensitive to the length of the primary cooling zone. This can be demonstrated by reducing and increasing the contact length by, e.g., 1 mm to 4 and 6 mm, respectively. The corresponding modeling results are labeled “case 2” when the contact length is 4 mm (v^sub 3^ = 846 mm) and “case 3” when the the contact length is 6 mm (y^sub 3^ = 844 mm). All other data are as given in Table I. It is seen from Table II and Figure 8 that the increase in the primary cooling length results in less macrosegregation, while decreasing the contact length leads to more macrosegregation. The reason for this behavior is that the longer contact length allows more heat to be extracted from the alloy, and the temperatures will be considerably lower in the surface region, especially at the start of the air-gap region. This means that the liquid fraction close to the surface in the air-gap region also will decrease when the contact length increases, and the exudation consequently becomes less. It is apparent that the calculated macrosegregation depends strongly on the length of the primary cooling region, and this illustrates the need for a criterion or a more sophisticated model which can predict the contact length.

As discussed in Reference 22, there is also considerable uncertainty related to the permeability model. Since the “friction forces” associated with the interdendritic melt flow are inversely proportional to the permeability, the extent of exudation depends strongly on the permeability. Moreover, the solidification is approximated by a lever-rule description, which also might contribute to the differences in simulated and measured temperatures and solute concentrations.[23]

It should be noted that free convection and the possible turbulence associated with this phenomenon are not accounted for in the present model (cf assumption no. 6 in Section Il-A). Including this phenomena could affect the macrosegregation in two ways. (1) The thermal field and the position of the liquid-mush transition might change significantly, resulting in different interdendritic flow patterns. 2) Solutal and thermal gradients might induce transport of liquid within the mushy zone. Concerning the first point, it is reasonable to believe that the temperature distribution close to the surface in the mold-contact and air-gap region is mainly governed by the thermal boundary conditions. Concerning the second point, it is believed that the driving force for free convection is not able to induce a significant amount of liquid flow within the mushy zone. These assumptions are supported by findings in the work of Reddy and Beckermann.[24] They performed calculations on a DC casting of an Al-Cu round ingot and found that macrosegregation due to free convection is almost neglectable compared to shrinkage-induced macrosegregation. In view of this, it is expected that the free convection would have a very small influence on the surface segregation predictions in the present case studies.


A mathematical model for the development of surface segregation due to exudation and solidification shrinkage in a DC casting is presented. The model is used to predict the surface segregation in a DC casting of an Al-3 pct Mg sheet ingot. Calculated values are in good qualitative agreement with measured values. Both a highly segregated surface layer and a solute-depleted zone within the ingot are predicted.


This research is part of the Brite-EuRam Project No. BE1112 under Contract No. BRPR-CT95-0112, which includes the following partners: Alusuisse-Lonza Services AG, Switzerland; Calcom SA, Switzerland; Ecole Polytechnique Federale de Lausanne, Switzerland; Elkem Aluminium ANS, Norway; Hoogovens Corporate Services BV, The Netherlands; Hydro Aluminum ASA, Norway; Institut National Polytechnique de Lorraine, France; Pechiney Recherche, France; Technische Universiteit Delft, The Netherlands; and Vereinigte Aluminiumwerke AG, Germany. SINTEF participates in the EMPACT Project as a major subcontractor to Hydro Aluminium ASA and Elkem Aluminium ANS. The authors thank the European Commission and the project partners for financial support.

*HEWLETT-PACKARD is a trademark of Hewlett-Packard Company, Colorado Springs, Co.


1. I.H. Hove, B. Andersson, and D. Voss: Proc. 3rd Conf on Aluminium Alloys-Their Physical and Mechanical Properties, L. Arnberg, O. Lohne, E. Nes, and N. Ryum, eds., Trondheim, Norway, June 1992, vol. 2, pp. 264-69.

2. E.F. Emley: Int. Met. Rev., 1976, pp. 75-115.

3. K. Buxmann: Metallurgy, 1977, vol. 31 (2), pp. 163-10.

4. L. Ohm and S. Engler: Metallurgy, 1989, vol. 43 (6), pp. 520-24. 1989.

5. A. Mo: Int. J Heat Mass Transfer, 1993, vol. 36 (18), pp. 4335-40.

6. B.R. Henriksen and E.K. Jensen: Light Metals, TMS-AIME, Warrendale, PA, 1993, pp. 969-77.

7. M.C. Flemings and G.E. Nereo: Trans. TMS-AIME, 1967, vol. 239, pp. 1449-61.

8. S. Minakawa, I.V. Samarasekera, and F. Weinberg: Metall. Trans. B, 1985, vol. 16B, pp. 595-604.

9. J.H. Chen and H.L. Tsai: Int. J Heat Mass Transfer, 1993, vol. 36 (12), pp. 3069-75.

10. A. Mo, T.E. Johnsen, B.R. Henriksen, E.K. Jensen, and O.R. Myhr: in Light Metals, U. Mannweiler, ed. TMS-AIME, Warrendale, PA, 1994, pp. 889-96.

11. E. Haug, A. Mo, and H.J. Thevik. Int. J Heat Mass Transfer, 1995, vol. 38 (9), pp. 1553-63.

12. J. Ni and C. Beckemla: Metal. Trans. B, 1991, vol. 22B, pp. 349-61. 13. N. Brookes and J.R. Hughes: Computer Methods Appl. Mech. Eng., 1982, vol. 32, pp. 199-259.

14. C.A.J. Fletcher: Computational Techniques for Fluid Dynamics, 2nd ed., Springer-Verlag, New York, NY, 1991, vol. II. IS. Genevieve Amiez and Pierre-Alain Gremaud: Numer. Math., 1991,

vol. 59, pp. 71-89.

16. H.A. van der Vorst: Siam J Sci. Comput., 1992, vol. 1, pp. 631-44. 17. A.M. Bruaset and H.P. Langtangen: Numerical Methods and Software Tools in Industrial Mathematics, M. Daehlen and A. Tveito, eds. Birkh*user, 1997, pp. 61-90.

18. The Diffpack home page. http://www.nobjects.com/prodserv/diffpack/. 19. E. Haug and H.J. Thevik: Modeling of Casting, Welding and Advanced Solidification Processes VIII, B. Thomas and C. Beermann, eds., Warrendale, PA,1998, pp. 321-28. 20. M.J.M. Krane and F.P. Incropera: Int. J Heat Mass Transfer, 1996, vol. 39, pp. 3567-79.

21. H.J. Thevik, A. Mo, A. Hkonsen, S. Benum, E.K. Jensen, and B.R. Henriksen: in Modeling of Casting, Welding and Advanced Solidification Processes VIII, B. Thomas and C. Beckermann, eds., TMS, Warrendale, PA, 1998, pp. 273-80.

22. A. Mo, T. Rusten, H.J. Thevik, B.R. Henriksen, and E.K. Jensen: in Light Metals, Reidar Huglen, ed., TMS-AIME, Warrendale, PA, 1997, pp. 667-74.

23. H.J. Thevik and A. Mo: Int. J. Heat Mass Transfer, 1997, vol. 40

(9), pp. 2055-65.

24. A.V. Reddy and C. Beckermann: Metall. Mater. Trans. B, 1997, vol. 28B, pp. 479-89.

HAVARD J. THEVIK, formerly Research Scientist with SINTEF, Oslo, Norway, is now with Det Norske Veritas, N-1322, Hovik, Norway. ASBJORN MO, Research Director, and TORGEIR RUSTEN, Senior Scientist, are with SINTEF, N-0314 Oslo, Norway.

Manuscript submitted April 29, 1998.

Copyright Minerals, Metals & Materials Society and ASM International Feb 1999

Provided by ProQuest Information and Learning Company. All rights Reserved