Open access

Simulation of Liquid Flow Permeability for Dendritic Structures during Solidification Process

Written By

S. M. H. Mirbagheri, H. Baiani, M. Barzegari and S. Firoozi

Submitted: 28 February 2011 Published: 05 July 2011

DOI: 10.5772/19647

From the Edited Volume

Computational Fluid Dynamics Technologies and Applications

Edited by Igor V. Minin and Oleg V. Minin

Chapter metrics overview

3,859 Chapter Downloads

View Full Metrics

1. Introduction

Alloy solidification is a two stage process, which starts by nucleation and ends by growth of solid phases. Subsequently, number, distribution and morphology (dendritic or non-dendritic) of the grains are formed during the solidification. Some critical defects such as micro/macro segregation, micro/macro porosities and micro/macro shrinkage take place in the solidification stage. The micro-defects are located in the interdendritic space, which are micro-channels that fluid flow through them in the last stage of the solidification. Herein, the region in the grain growth stage is introduced as a mushy zone (or porous media), where the solid phase is constantly progressing; and the ability of fluid to flow into the mushy zone is known as permeability of interdendritic liquid. Therefore, formation of the micro-defects depends on controlling of the permeability factor. In a great number of studies micro/macro solidification models have been simulated based on the permeability factor using Darcy’s law (Ganesan & Poirier, 1990; Nandapurkar et al., 1991; Poirier, 1987; Worster, 1991).

Interdendritic flow, in many CFD documents, is described using Darcy’s law, which relates the fluid flow rate to the pressure gradient, fluid viscosity, and permeability of the porous medium. To obtain an expression for the permeability as a function of the porosity of the porous medium, one generally considers flow through an idealized medium geometry, since it is impractical to solve the flow equations for the complex flow between the particles. Fig. 1 presents two viewpoints for investigation of the permeability in the porous media: metallurgical view and non-metallurgical (or common) view.

As shown in Fig. 1a two of the most commonly used geometries for analytical models are capillaries (Carman, 1937; Chen et al. 1995; Williams et al., 1974) and an array of spaced particles. A more realistic approach will be introduced that assumes geometry of a periodic or random array of cylinders. Since it is not possible to solve analytically for this type of flow over the full range of porosities, two limiting closed form solutions are used for lubrication and point-particle (dilute) models in low and high porosities, respectively. Analysis of permeability for Stokes flow through periodic arrays of cylinders were done by Sangani & Acrivos (1982), Sparrow & Loefler (1959) and, Larson & Higdon (1986).

The effect of fluid inertia on pressure drop required to drive the flow is a function of Reynolds number. Several authors computed the fluid flow through periodic arrays of cylinders as the function of Reynolds numbers for three ranges of low, moderate, and high Reynolds number (Cai & Berdichevsky, 1993; Edwards et al., 1990; Eidsath et al., 1983; Ghaddar, 1995; Koch & Ladd, 1997; McCartney, 1994; Nagelhout et al., 1995; Sangani & Mo, 1994; Sangani & Yao, 1988; Thom & Aplelt, 1961). Particularly noteworthy is the work of Koch & Ladd (1997) for modeling permeability and drag force based on Reynolds number using a Lattice-Boltzmann formulation. Also McCartney (1994) calculated the permeability in the range of low Reynolds numbers up to about 150 by Lattice Gas Cellular Automat (LGCA) model.

Figure 1.

Two categories of permeability models; a) Common permeability model; b) Metallurgical permeability model (solidification)

In the field of metallurgy, however, the models mentioned above had an important application in the mathematical modeling of flow through arrays of dendrites during the solidification of mushy alloys. Fig. 1b shows a schematic category for modeling of permeability during solidification process. There are many investigations for experimental measuring of permeability by variant methods during solidification process (Apelian et al., 1974; Duncan et al., 1999; Murakami et al., 1984; Murakami et al., 1983; Nielsen et al., 2001; Poirier & Ocansey, 1993; Streat & Weinberg, 1976). However, most of the experimental methods had two goals;

  1. an experimental goal is to measure the permeability of an alloy at a constant temperature and to correlate the permeability with microstructures quantities. In this case, it is necessary to obtain a quantitative stereological characterization of the solid-phase morphology. This has been achieved by quenching of the mushy sample during the permeability experiment and subsequent metallography and image analysis, which distinguish between the phases that were solid and liquid during the experiment.

  2. Another experimental goal is to measure the permeability and the temperature of an alloy during solidification/re-melting and to correlate the permeability with the solid fraction, which can be estimated from the temperature-time curves.

In this case, the composition of the interdendritic liquid varies throughout the experiment, which complicates the experimental design, the experiment method, and the data analysis (Nielsen et al., 2001).

The mathematical method includes two categories;

  1. flow through a network of equiaxed grains,

  2. flow through columnar dendrite networks with flows parallel and normal to the primary dendrites.

Piwonka & Flemings (1966), Apelian et al. (1974), Streat & Weinberg (1976), Liu et al. (1989) and Murakami et al. (1984) reported the permeability in equiaxed dendritic structures. In these investigations, however, the microstructural length scales were not reported or they could not be used to estimate permeability in a solidification model designed to calculate macro segregation. Ganesan & Poirier (1990 ) and Ocansey & Poirier (1993) measured and reported the permeability with inverse of the specific area of the solid selected as the length scale in equiaxed microstructures based on Kozeny-Carman model. There are at least two length scales associated with equiaxed dendritic solidification: the secondary dendrite arm spacing (or interdendritic flow) and the grain size (or extra-dendritic flow) (de Groh et al., 1993; Wang et al., 1995). Brown et al. (2002) developed a numerical model for the simulation of 3D flow through equiaxed dendrites of an Al13Cu3Si alloy and the determined the variation in permeability of structure as solidification progressed. The model involved the evolution of an equiaxed dendrite and the application of a CFD program to calculate permeability from Darcy’s law.

Streat & Weinberg (1976), Poirior (1987), Ganesan & Poirier (1990), Nandapurkar et al. (1991) Worster (1991) reported the permeability in columnar dendritic structures for macroscopic scale. However, due to the very complex microstructure of the dendrites, permeability determination still remains a challenge. Indeed, the dendritic columnar region is characterized, first, by a strong anisotropy, which requires knowledge of the different components of the permeability tensor (Murakami et al., 1984; Murakami et al., 1983; Poirier, 1987) and, second, by the non-uniform macroscopic properties such as the liquid volume fraction, which continuously varies from unity in the melt to zero in the solid region. Ganesan et al. (1992) showed that the permeability for a flow parallel to primary dendritic arms is hardly dependent on the secondary dendrite arm spacing. In the columnar configuration, data in different dendritic structures have been summarized by Poirier (1987) and quantitative relationships for parallel and normal permeability have been derived using a regression analysis. In the range of liquid fraction considered (>0.66) the regressions were in good agreement with the classical physical models, but the extrapolation beyond the upper limit failed. Numerical experiments for parallel and normal flow to the primary dendritic arms in columnar structures with a high liquid volume fraction (>0.6) have been found to be more successful (Ganesan et al., 1992). Mirbagheri (2008; Mirbagheri & Khajeh, 2008) measured and simulated interdendritic flow for mushy alloy based on permeability factor and modeled some micro structural factors on the mushy alloys permeability. Bahat et al. (1995) calculated from digitized images, the permeability for flow normal to primary dendrite arm using a Navier-Stokes finite element solver.

As mentioned above, permeability for flow through transverse sections of columnar-dendritic alloy have modeled as flows through array of circles, rhombi, cruciform and schematic dendrites. These models used surface area to volume ratio of the solid for normalizing permeability without predicting the effect of protuberance and radius of dendrite for a fixed fraction of solid as well the effect of angle of dendrite radius via liquid streamlines. In spite of these valuable researches in this field, determination of dendritic structural permeability due to the complex microstructure of the dendrites still remains a challenge in both fields of mathematical and experimental methods.

In this chapter, a numerical model has been introduced for determination of liquid flow permeability through dendritic solid phases during growth. The model includes three stages; first, numerical simulation of nucleation and growth of the equiaxed grains using a novel Cellular Automation Finite Difference (CAFD) method, and second, numerical simulation of micro fluid flow for interdendritic liquid alloys using Computational Fluid Dynamics (CFD) technique, and third, calculation of permeability based on Darcy’s law by the pressure and velocity results of CFD code. Finally determining and modeling the permeability variations versus the cooling rate and the solidification rate during growth. This model can be linked as a module into a commercial macro fluid flow code in order to predict the micro-defects such as the micro porosities, shrinkages or micro-segregations.


2. Computation models

In the present work, two separate computation models of the nucleation and grain growth, and interdendritic liquid flow have been developed and coupled for modeling of permeability in mushy alloys.

This is achieved by combining sub-models for each of these processes, i.e., computation of nucleation and grain growth by using the CAFD, and the micro fluid flow by using the CFD model for calculation of the interdendritic permeability. The governing equations are described in details in the next sections.

2.1. Solid phase generation code

The solid phase generation code, which has been developed in this investigation is based on the CA and KGT model (Atwood & Lee, 2000; Lee et al., 2001). In other words, the solidification model is non-constrained nucleation and growth. The model comprises of:

  1. Stochastic grain density based on local under-cooling, and

  2. grain growth based on local thermal and under-cooling (Kurz et al., 2001; Kurz et al., 1986).

However, in the present work, for a binary alloy system (Fig. 2) in the mushy zone, it is assumed that there is no constitutional under-cooling during growth phenomenon. Number, distribution, and morphology of nuclei as well as the growth rate are controlled by only thermal under-cooling, which is produced by Newtonian’s heat transfer in a two dimensional (2D) space. In this condition, at the beginning, all the liquid have the same under- cooling (i.e., the gradient of the under-cooling or temperature is equal to zero). As shown in Fig. 3, once the nucleation takes place, temperature around of the nuclei is raised, because of the liberation of latent heat. Therefore, the temperature gradient around of the grains is negative.

After the solidification, the microstructure of the solidified alloy will consist of fully equiaxed grains and no columnar grains. Therefore, in this investigation, initial and boundary conditions for the simulation of heat transfer and liquid flow are based on the non-constrained solidification.

Figure 2.

Phase diagram of an Al-Si binary alloy system. The regions between the liquidus and the solidus temperatures are the mushy zone

2.1.1. Governing equations

  1. Heat transfer equations

    In solidification process, there are two terms of heat transfer and latent heat (H f ). Interaction of these two terms affects the domain of thermal distribution. The heat transfer equation for the mushy zone may be written as:

    ρ C p T t = ( K T ) + ρ H f f s t E1
    f s = T l i q T T l i q T s o l E2

    where ρ, C p and K are the density, heat capacity and thermal conductivity, respectively; f s is solid fraction and T sol , T liq are solidus and liquidus temperature, respectively. ρH f (∂f s /∂t) is the heat source term, which is a function of temperature in the mushy zone, and is written as follows:

    ρ H f f s t = ρ H f ( f s T ) ( T t ) E3

    By differentiating Eq. 2, and substitution in Eq. 3:

    ρ H f f s t = ρ H f f s T T t = ρ H f ( 1 T l i q T s o l ) ( T t ) E4

    The fraction of solid in the mushy zone is estimated by Eq. 2. The release of latent heat between liquidus and solidus temperature is calculated by substituting Eq. 4 into the second term of Eq. 1. Therefore heat transfer equation is given by:

    ρ ( C P H f ( 1 T l i q T s o l ) ) T t = ( K T ) E5
    ρ [ C P e q ] T t = ( K T ) E6

    where C P e q can be considered as a quasi-specific heat given by:

    C P e q = [ C P H f ( 1 T l i q T s o l ) ] E7

    Figure 3.

    Temperature distribution around the grains

    Physical properties of the liquid and solid are assumed to be constant above T liq and below T sol, respectively. However, in the mushy zone, coefficients of heat conductivity and thermal capacity are presented as k mu = f L k L + f S k S , and C p m u = f L C P l + f S C P S (Mirbagheri & Silk, 2007). In a case that the alloy composition has no solidification range (i.e., ΔT 0 = T liqT sol = 0), such as eutectic composition, the right hand side in Eq. 4 approaches infinity and as a result, a virtual solidification range of 0.1-1 C is assumed.

  2. Nucleation equations

    When the temperature falls below the liquidus temperature, nucleation begins. In this condition the number of nuclei at each temperature and time are calculated from Eq. 7, as follows:

    N s = N t o t a l exp ( C m T l i q 2 T ( T l i q T ) 2 ) E8

    After determining the number of nuclei in each time step, by assigning a random distribution function, nuclei are distributed in the liquid domain. In the next step, the latent heat of solidification is calculated to adjust the temperature.

  3. Growth equations

The final stage of solidification process is the growth simulation. As mentioned before, once a negative temperature gradient is present in liquid adjacent to the grains, equiaxed grains grow as shown in Fig 3. The direction of primary arms of the equiaxed dendrites depend on crystal structure. Here, a B.C.C. crystal structure is assumed, in which each equiaxed dendrite has four perpendicular primary arms in 2D space, where they can grow in 48 crystalline directions (An et al., 2000). In order to simulate the morphology of the growth, a simple shape was used for grains based on Eq. 8 in polar coordinates.

r = R d + P d C o s ( 4 ( θ + θ 0 ) ) i f 0 P d R d 1 ; ( R * ) h R * = R d + P d E9

Fig. 4 shows a “cloverleaf” morphology for a dendrite section created based on Eq. 8, where P d is perturbation, R d radius of spherical nuclei prior to perturbation, and θ angle between primary arm direction and stream line.

Figure 4.

Shape of a nucleus at the beginning of the growth

A function (Eq. 9) needs to be defined for the dendrites radius growth rate (dr), which is added to the surface of existing grains in each time step of solidification stage.

r n e w = r o l d + d r ; d r = d f s 2 π r E10

Finally, results of nucleation and growth simulation at each df s, is used in the CFD code to calculate the permeability in domain.

2.1.2. CA-FD model

  1. CA-FD model

    The finite difference approximation of heat transfer equation is:

    T i , j n + 1 = T i , j n Δ t { ( U T X + V T Y ) + ( D Q X + D Q Y ) } E11

    where UTX, DTXL, DTXR and DQX are defined as followed:

    U T X = 0.5 { ( 1 + α ) U i 1 , j ( D T X L ) + ( 1 α ) U i , j ( D T X R ) } E12
    D T X L = T i , j T i 1 , j Δ x D T X R = T i + 1 , j T i , j Δ x D Q X = Q X R Q X L ρ C Ρ Δ x E13

    The terms of VTY and DQY are found in a similar way to UTX and DQX, respectively. If cells (i, j) and (i+1, j) are liquid:

    Q X R = k R T i + 1 , j n T i , j n Δ x Q X L = k l T i , j n T i 1 , j n Δ x E14

    K R and K L are thermal conductivity and found through the following equations:

    K R 1 = 0.5 ( 1 k i , j + 1 k i + 1 , j ) K L 1 = 0.5 ( 1 k i , j + 1 k i 1 , j ) E15

    If cell (i, j, k) contains a mixture of solid and liquid metal, the thermal conductivity of this cell is:

    k i , j = F S i , j k S i , j + F L i , j k L i , j , F S i , j + F L i , j = 1 , F L i , j = T i , j T s T l T s ; T s T T l E16

    In the freezing range the specific heat and liquid fraction of the mushy metal is found through the following equation:

    C p L S = Δ H f T l T s ; T s ( T ( T l E18

    where T s , T l are solidus and liquidus temperatures respectively. In finite difference form, C p is calculated as follows:

    C Ρ = C P L ; T i , j n ) T i , j n + 1 ) T E19
    C Ρ = C P S ; T i , j n + 1 ) T i , j n ) T S E20
    C Ρ = C Ρ L T i , j n T L T i , j n T i , j n + 1 + C Ρ L S T L T i , j n + 1 T i , j n T i , j n + 1 ; T i , j n ) T L ) T i , j n + 1 ) T S E21
    C Ρ = C Ρ L S T i , j n T S T i , j n T i , j n + 1 + C Ρ S T S T i , j n + 1 T i , j n T i , j n + 1 ; T L T i , j n ) T S ) T i , j n + 1 E22

    It should be noted that iteration is required here, because not only does T i,j depends on C p , but also C p depends on T i,j .

  2. Direction of crystalline growth

For a solidifying cell, as the solid fraction within the cell becomes greater than zero, the local temperature of the particles is obtained using the phase diagram and the under-cooling is calculated accordingly. In each solidifying cell, the change in solid fraction is primarily determined by KGT model (Kurz et al., 1986), which calculates the maximum growth rate based on a given under-cooling at near absolute stability limit. The solid fraction is further corrected by diffusion-controlled growth once the solid fraction reaches a critical value; it can grow into its neighboring liquid cells, providing the cells are under-cooled. Where the solid fraction of cells approaches unity, the cell is considered as fully solidified and the grain growth ends. A captured liquid cell by a growing neighboring cell is assigned the same grain orientation as its growing neighbor (Lee et al., 2001). Fig.5 shows the CA-FD algorithm of heat transfer during nucleation and growth

2.2. Inter-dendritic liquid flow code

2.2.1. Governing equations

  1. Fluid flow equations

    The Navier-Stokes and continuity equations are used to simulate flow of the interdendritic liquid through the network of dendritic solids. The Navier-Stokes equation for the incompressible liquid is given by the following equation (Bahat et al., 1995):

    ρ L D V D t = P + ρ L g + μ 2 V E23

    Also, the mass continuity equation for the incompressible liquid is as follows:

    ρ t + . ( ρ V ) = 0 E24

    It is often desirable to reproduce large scale physical experiments in scaled-down and more manageable laboratory settings. Information on the flow is contained in the parameters which characterize it, such as dynamic viscosity, velocity and density. If these parameters are combined in a suitable way to yield dimensionless quantities, then these enable one to make the desired statements relating the flows on the large and small scales. By introducing the following dimensionless variables (Griebel et al., 1998):

    x * = x L , t * = u t L , u * = u u , g * = L u 2 g , p * = p p ρ u 2 , a n d Re = ρ u L μ E25

    Navier-Stokes equation can be rewritten as the dimensionless form:

    D V D t = P + g + 1 Re 2 V E26

    Thus by solving Eqs. 16a, 16b and 16c, one can get the, velocity and pressure fields for an incompressible viscous fluid:

    u t + p x = 1 Re ( 2 u x 2 + 2 u y 2 ) ( u 2 ) x ( u v ) y + g x E27
    v t + p x = 1 Re ( 2 v x 2 + 2 v y 2 ) ( u v ) x ( v 2 ) y + g x E28
    u x + v y = 0 E29

    Figure 5.

    CA-FD algorithm of heat transfer during nucleation and growth

    Equation 15 can be written as two independent equations based on projection model as follows:

    ( δ V ) t + V . V = 1 Re 2 V ; δ V = V * V E30
    ( δ V ) t + P = 0 E31

    Equation 17a is independent of pressure and by using the Poisson’s equation (Eq. 18) the Navier-Stokes equation is solved implicitly.

    2 P = 1 Δ t . V * E32

    After solution of the Navier-Stokes equations, the pressure and velocity fields can be calculated to obtain the permeability of the mushy zone in each growth sequence.

  2. Permeability equations

Permeability is a measure of the ability of a porous material to transmit fluids. Darcy's law is a simple proportional relationship between the instantaneous discharge rate through a porous medium (V), the viscosity of the fluid (μ) and the pressure drop over a given distance. The law was formulated by Henry Darcy based on the results of experiments on the flow of water through beds of sand (Darcy, 1856).

V = K μ P E33

where, K is permeability. The problem with this law is that it is not valid at high Reynolds number. As the velocity increases, inertia effects appear, which complicates the calculation of permeability. To overcome this weakness, several corrections have been applied to the Darcy’s law by various researchers. In this investigation, the effective Darcy’s law equation was used for calculation of permeability (Eq. 19).

V = K μ f l P E34

2.2.2. Numerical solution of the fluid flow governing equations

The main purpose of the permeability simulation in the solidification process is calculating the velocity profile for inter-dendritic space during equiaxed grains growth. Therefore, governing equations are solved by FDM. The numerical solution method can be considered in four steps:

  1. Meshing of the system,

  2. Converting differential equations to finite difference approximation,

  3. Solution of the finite difference approximations of momentum in order to calculate velocity profile, and

  4. Calculation of permeability of the mushy zone by adding Darcy’s law in each growth sequences.

  1. Meshing

    The computational domain is divided into a number of cells with ΔX, ΔY dimensions, and all cell dimensions are equal for all calculations. Curved boundaries of system are approximated by stepped boundaries. For example, if the solid phase fraction occupies more than 0.9 volume of a cell, the cell is considered to be a complete solid cell. As shown in Fig. 6, domain is discretized as a staggered grid, i.e. the scalar and the vector variables are located at the centre and sides of the computational cell, respectively.

    Figure 6.

    Schematic of a staggered mesh and location of vector and scalar variables

    The diffusive terms of the momentum equation is discretized using central differencing. The discretization of the convective terms of momentum equation has been done using the donor-cell scheme (Mirbagheri et al., 2003)

  2. Finite difference approximations

    The finite difference approximation of momentum equations are:

    U i , j n + 1 = U i , j n + Δ t [ g x F U X F U Y + 1 Re ( V I S X ) ] + Δ t ( P i , j n + 1 P i + 1 , j n + 1 Δ x ) U i , j n + 1 = A x n + Δ t ( P i , j n + 1 P i + 1 , j n + 1 Δ x ) E35
    V i , j n + 1 = V i , j n + Δ t [ g Y F V X F V Y + 1 Re ( V I S Y ) ] + Δ t ( P i , j n + 1 P i , j + 1 n + 1 Δ y ) V i , j n + 1 = A y n + Δ t ( P i , j n + 1 P i , j + 1 n + 1 Δ y ) E36
    A x n = U i , j n + Δ t [ g x F U X F U Y + 1 Re ( V I S X ) ] A y n = V i , j n + Δ t [ g Y F V X F V Y + 1 Re ( V I S Y ) ] E37
    F U X = U i , j n 2 Δ x [ D U L + D U R + α S i g n ( U α ) ( D U L D U R ) ] ; D U L = U i , j n U i 1 , j n ; D U R = U i + 1 , j n U i , j n ; U α = ( D U R + D U L ) / 2 ; V I S X = ( U i + 1 , j n 2 U i , j n + U i 1 , j n Δ x 2 + U i , j + 1 n 2 U i , j n + U i , j 1 n Δ y 2 ) E38

    Other terms of the flux and the viscosity in x and y directions such as FUY, FVX, FVY, and VISY are obtained in a similar way as FUX and VISX terms, respectively. The finite difference approximation of Poisson’s equation (Eq. 23) is

    p i + 1 , j n + 1 2 p i , j n + 1 + p i 1 , j n + 1 ( δ x ) 2 + p i , j + 1 n + 1 2 p i , j n + 1 + p i , j 1 n + 1 ( δ y ) 2 = 1 δ t ( A x , i , j n A x , i 1 , j n δ x + A y , i , j n A y , i , j 1 n δ y ) E39
  3. Solving and computing procedures

    Various methods are available to adjust the pressure term in the Navier-Stokes equation (Hong, 2004; Versteeg & Malalasekra, 1995). Methods like SMAC (Amsden & Harlow, 1970) and SOLA (Hirt et al., 1975) use an explicit scheme to solve the pressure adjustment equation, which is based on divergence of pressure in each computational cell. The computational time are relatively long in these methods and are seldom used today due their low performance. Instead, other methods such as projection method are used to adjust the pressure term and compute new velocities that satisfy the continuity equation. Semi-implicit methods such as SIMPLE, SIMPLER and PISO are also frequently used (Versteeg & Malalasekra, 1995).

    To solve the governing equations of fluid flow in a mushy zone for binary alloys, a new CFD code has been developed based on fundamentals presented by (Griebel, 2011). Therefore, to correct the pressures calculated from the Navier-Stokes momentum equation, the Chorin’s projection method is used, and the resulting Poisson’s equation is solved using Successive Over-Relaxation (SOR) iterative method.

    Projection method uses an auxiliary velocity V * to obtain a Poisson’s equation for pressure. This equation can be solved using any solution algorithm such as the SOR or Gauss-Seidel methods. The momentum and continuity equations in time derivative form are:

    V n + 1 V n Δ t + V n . V n + P n + 1 = 1 Re 2 V n E40
    V n + 1 = 0 E41

    The superscripts (n) and (n+1) denote old and new time level, respectively. In projection method V n+1 domain is calculated at each new time step. Using the auxiliary velocity, the momentum equation can be split to two independent equations with Eq. 26 with no pressure term and Eq. 27 with pressure term.

    V * V n Δ t + V n . V n = 1 Re 2 V n E42
    V n + 1 V * Δ t + P n + 1 = 0 E43

    The divergence of Eq. 26 takes the form:

    . V n + 1 . V * Δ t + 2 P n + 1 = 0 E44

    Continuity equation (Eq. 25), requires that V n + 1 to be zero, thus

    2 P n + 1 = 1 Δ t . V * E45

    The overall solution process is using Eq. 26 to obtain V * , then to solving the Poisson’s equation (Eq. 29) to find the adjusted pressure values and finally solving Eq. 27 to obtain Vn+1. Fig. 7 describes the algorithm used in the present code.

  4. Calculation of permeability

After obtaining the pressure and velocity fields, which is called herein as the temporary permeability subroutine, and using Darcy’s law, coefficient of temporary permeability is calculated for domain in each solid fraction and their changes are saved until the end of solidification. This subroutine is showed at the lower part of the flowchart in Fig 7.

2.3. Validation of the present CFD code

The developed CFD code was validated by comparing the pressure gradient for an equiaxed grain with the Fluent code at the same initial and boundary condition. Fig. 8 shows four pressure fields adjacent to a single dendrite at solid fractions of 0.02, 0.08, 0.19, and 0.34, for and Al-6%wt.Si alloy. Fig. 9 shows comparison of the pressure gradient results between the present CFD code and the Fluent code based on data of Table 1. The CFD code predictions at high and low solid fractions are in good agreement with the corresponding Fluent code. Therefore, the present code could be valid for simulation of other domains.

Table 1.

Boundary conditions, thermo-physical properties, and initial condition for Fig 8 and Fig 10 to 13 domains

Figure 7.

Algorithm of the projection method for the present CFD code

Figure 8.

Simulation of the pressure fields for 4 solid fractions during Al-6%wt.Si alloy grain growth based on Table 1 data

Figure 9.

Comparison of the pressure gradient results of Fig. 8 domain between Fluent and present code based of Table 1 data


3. Temporary permeability and effective parameters

It is clear that the permeability is a function of pressure gradient based on Darcy’s law; therefore, the simulation of the pressure field in the mushy zone and the factors that affect it is our target. These factors such as, cooling rate, nucleation and growth rate, distribution of nucleus and grains, and morphology of grains, subsequently affects the permeability.

Since, in this investigation a Newtonian thermal condition is assumed; nucleation phenomena are entirely random. Therefore, the distribution of nuclei location could affect the behavior of the pressure field due to the drag force of nucleus on fluid flow. Figs. 10 to 12 show the effect of distribution of nuclei and grains on the pressure field for three constant solid fractions, 0.02, 0.10, 0.24, and 0.43 at a fixed cooling rate. Stream lines are substantially different in these figures, especial at high solid fractions (i.e. fs=0.43), which are affected from the number and size of the grains as well as the distribution of nuclei and grains.

Figure 10.

Simulation of the pressure fields adjacent to 4 types of nucleus distribution with 1% solid fraction and cooling rate of 100 J in each time step.(dQ/dt = 1000 J/s, fs=0.1%)

Figure 11.

Simulation of the pressure field adjacent to 4 types of nucleus distribution with 35% solid fraction and cooling rate of 100 J in each time step.(dQ /dt= 1000 J/s, fs=0.35%)

Fig 13 shows the effect of cooling rate on the pressure field for a fixed solid fraction (fs=0.43%). Results show that by increasing the rate of heat extraction from the domain (dQ), some parameter such as size, numbers, and distribution of nucleus and grains have changed and subsequently the pressure field has changed. For example at dQ = 70, grains are coarse and large, however at dQ = 300, grains are fine and small.

Fig 14 shows the temporary permeability that was calculated from the simulated pressure fields. In low solid fraction (fs< 0.06) and high solid fraction (fs>0.78), the behavior of permeability is as asymptotic function, because of f l factor in Eq. 19a. However between two solid fractions, behavior of the temporary permeability is almost liner. Markers, ■, ●,▲, and ◊, demonstrate the effect of different grains distributions on the permeability coefficient, which do not show a significant effect, especially at low and high solid fractions.

Figure 12.

Simulation of the pressure field adjacent to 4 types of nucleus distribution with 72% solid fraction and cooling rate of 100 J in each time step.(dQ/dt= 1000 J/s, fs=0.72%)


4. Conclusions

It is possible, especially in commercial codes, to obtain the thermal history in casting or solidification codes. In fact all numerical fluid flow and heat transfer software can save the cooling rate ( T t ) , temperature gradient ( G = Q k = T x ), and solidification rate ( v = f s t ) for each location of the meshed domain. Therefore, in the present code if the interdendritic permeability is defined as a function of cooling rate, the micro-structure will be related to the thermal history in a macro scale, and the formation of micro defects could be predicted. In other words, the present code is capable of predicting the micro-defects based on thermal history. However, the cooling rate is not an independent variable and should be determined by the temperature gradient and the solidification rate, as follows:

T t = ( T x ) ( x t ) =  G . v ,  for a single phase ,  and E46
T t = ( k . G ) / ( 1 + ΔHm . v ) ,  during phase change E47

Figure 13.

Simulation of the pressure field adjacent to 4 types of nucleus distribution with 0.45% constant solid fraction and dQ equals to 70, 130, 250, and 300 J in each time step

It means that increasing the cooling rate results in an increase in the number of nuclei and the change of dendrite morphology from faceted surface to dendritic surface (cloverleaf grain) and that in turn affects the interdendritic permeability. Therefore, the effect of temporary solid fraction (df s ) and heat extraction rate (dQ) on the temporary permeability is modeled as an asymptotic function and plotted in Fig. 15 and 16. Results shown in Fig. 16 show that after formation of 0.72% solid, grains are inter-connected and join together. Consequently, there is no fluid stream into the mushy zone, but only some micro fluid flow as a vortex of the entrapped liquid between the grains. It seems that the critical permeability for the formation of micro-porosity, in a mushy alloy such as Al-6%wt.Si, is less than 1E-10 to 1E-11 (m2). This is true, even at very high cooling rates of over 3000 J/s and different distributions of the equiaxed grains.

Figure 14.

Calculation of permeability based of simulation results of the present CFD and CA-FD codes during grain growth with 4 kinds of nucleus distribution (■●▲◊) versus solid fraction

In this investigation, the temporary permeability results, which was obtained based on the cooling and solidification rate can be modeled as a 3-dimensional surface. The resulting equation of this 3D-surface as given in Eq. 30 can be used in all casting simulation software.

L o g ( P e r m e a b i l i t y ) = 7.175 f s + 0.002 d Q + 5.185 E48

Finally, the significant findings of this investigation can be phrased as follows:

  1. In this investigation an algorithm was developed to calculate the permeability of mushy alloys. To simulate the permeability, CFD and CA-FD codes were coupled. The CA-FD code was used for nucleation and growth of the equiaxed dendrites and the CFD code was used for simulation of the pressure and velocity fields adjacent to the nuclei and grains.

  2. The temporary permeability was calculated based on the present codes (CFD and CA-FD)

    Figure 15.

    Changes in permeability results in the Al-6%Si alloy

    Figure 16.

    Modeling of temporary permeability as a 3-D surface during equiaxed grain growth based on simulation results of the present code

  3. Permeability simulation results showed that apart from low and high solid fractions, f S <0.06 and f S >0.72, respectively, the temporary permeability versus the solid fraction during solidification has a linier behavior.

  4. Results showed that the temporary permeability during nucleation and growth of the equiaxed grains depends on cooling and solidification rate.

  5. At high cooling and solidification rate fluctuations on temporary permeability was observed, which was presumably due to the perturbation of the grain surfaces or the evolution of semi-spherical grains to cloverleaf dendrites.

  6. Permeability results of present code can be utilized in commercial casting and solidification software.


  1. 1. Amsden A. A. Harlow F. H. 1970 The SMAC method: A Numerical Technique for Calculating Incompressible Fluid Flows. Los Alamos Scientific Laboratory Report LA-4370 New Mexico.
  2. 2. An S. U. Larionov V. Monastyrski V. Monastyrskaia E. Grafas I. Oh J. M. Lim O. D. Kim S. H. Lee J. H. Seo D. Y. 2000 The thermal analysis of the mushy zone and grain structure changes during directional solidification of superalloys. In: Superalloys 2000. Green, K. A., Pollock, T. M. & Kissinger, R. D., TMS.
  3. 3. Apelian D. Flemings M. C. Mehrabian R. 1974 Specific permeability of partially solidified dendritic networks of Al-Si alloys. Metallurgical and Materials Transactions B. 5 2533 2537 .
  4. 4. Atwood R. C. Lee P. D. 2000 A Combined cellular automaton and diffusion model for the prediction of porosity formation during solidification. Modelling of casting welding and advanced solidification processing IX, Aachen, Germany.
  5. 5. Bahat M. S. Poirier D. R. Heinrich J. C. 1995 Permeability for cross flow through columnar-dendritic alloys. Metallurgical and Materials Transactions B. 26 5): 1049 1056 .
  6. 6. Brown S. G. R. Spittle J. A. Walden-Bevan R. 2002 Numerical Determination of Liquid Flow Permeabilities for Equiaxed Dendritic Structures. Acta Materialia. 50 1559 1569 .
  7. 7. Cai Z. Berdichevsky A. L. 1993 Numerical simulation on the permeability variations of a fiber assembly. Polymer Composites. 14 6): 529 539 .
  8. 8. Carman P. C. 1937 Fluid flow through a granular bed. Transactions of the Institution of Chemical Engineers. 15 150 166 .
  9. 9. Chen Y. T. Davis H. T. Macosko C. W. 1995 Wetting of fiber mats for composites manufacturing: I. Visualization experiments. 41 10): 2261 2273 .
  10. 10. Darcy H. 1856 Les Fontaines Publiques De Ville de Dijon, Paris.
  11. 11. de Groh H. C. Weidman P. D. Zakhem R. Ahuja S. Beckermann C. 1993 Calculation of dendrite settling velocities using a porous envelope Metallurgical and Materials Transactions B. 24 5: 749 753 .
  12. 12. Duncan A. J. Han Q. Viswanathan S. 1999 Measurement of liquid permeability in the mushy zones of aluminum-copper alloys Metallurgical and Materials Transactions B. 30 4): 745 750 .
  13. 13. Edwards D. A. Shapiro M. Bar-Yoseph P. Shapira M. 1990 The influence of Reynolds number upon the apparent permeability of spatially periodic arrays of cylinders. Physics of Fluids A. 2 45 60 .
  14. 14. Eidsath A. R. G. C. Whitaker S. Herrmann L. R. 1983 Dispersion in pulsed systems- III : Comparison between theory and experiments for packed beds. Chemical Engineering Science. 38 1803 1816 .
  15. 15. Ganesan S. Chan C. L. Poirier D. R. 1992 Permeability for flow parallel to primary dendrite arms. Materials Science and Engineering A. 151 97 105 .
  16. 16. Ganesan S. Poirier D. R. 1990 Metallurgical and Materials Transactions. 21 B): 173 181 .
  17. 17. Ganesan S. Poirier D. R. 1990 Conservation of Mass and Momentum for the Flow of Interdendritic Liquid during Solidification. Metallurgical and Materials Transactions B. 21 173 181 .
  18. 18. Ghaddar C. K. 1995 On the permeability of unidirectional fibrous media: A parallel computational approach. Physics of Fluids A. 7 2563 2585 .
  19. 19. Griebel M. 2011 Retrieved 2011, from
  20. 20. Hirt C. W. Nichols B. D. Romero N. C. 1975 SOLA- A Numerical Solution Algorithm for Transient Fluid Flows. Los Alamos Scientific Laboratory Report LA-5852 New Mexico.
  21. 21. Hong C. P. 2004 Computer modelling of heat and fluid flow in materials processing, Taylor & Francis.
  22. 22. Koch D. L. Ladd A. J. C. 1997 Moderate Reynolds number flows through periodic and random arrays of aligned cylinders. J. Fluid Mech. 349 31 66 .
  23. 23. Kurz W. Bezencon C. Gäumann M. 2001 Columnar to equiaxed transition in solidification processing. Science Technology and Advance Materials. 2 1): 185 191 .
  24. 24. Kurz W. Giovanola B. Trivedi R. 1986 Theory of Microstructrual Development during Rapid Solidification. Acta Metallurgica. 34 5): 823
  25. 25. Larson R. E. Higdon J. J. L. 1986 Microscopic flow near the surface of two-dimensional porous media, Part 1. Axial flow. Journal of Fluid Mechanics. 166 449 472 .
  26. 26. Lee P. D. Atwood R. C. Dashwood R. J. Nagaumi H. 2001 Modeling of Porosity Formation in Direct Chill Cast Aluminium-Magnesium Alloys. Materials Science and Engineering. Vol.: pp. Submitted.
  27. 27. Lee P. D. Chirazi A. See D. 2001 Modeling Microporosity in Aluminum-Silicon Alloys: a Review. Journal of Light Metals. 1 1): 15
  28. 28. Liu C. Y. Murakami K. Okamoto T. 1989 Permeability of dendrite network of cubic alloys. Materials Science and technology. 5 11): 1148 1152 .
  29. 29. Mc Cartney J. F. 1994 Flow through arrays of cylinders: Lattice gas cellular automata simulations. Physics of Fluids. 6 2): 435 437 .
  30. 30. Mirbagheri S. M. H. 2008 Modeling of the equiaxed dendrite coarsening based on the interdendritic liquid permeability during alloy solidification. Metallurgical and Materials Transactions B-Process Metallurgy and Materials Processing Science. 39 3): 469 483 .
  31. 31. Mirbagheri S. M. H. Ashuri H. Varahram N. Davami P. 2003 Simulation of mould filling in lost foam casting process. International Journal of Cast Metals Research. 16 6): 554 565 .
  32. 32. Mirbagheri S. M. H. Khajeh E. 2008 Modelling and simulation of equiaxed dendritic structures permeability for Pb-Sn alloys. Materials Science and Technology. 24 12): 1444 1451 .
  33. 33. Mirbagheri S. M. H. Silk J. 2007 Simulation of Si concentration effect on the permeability for columnar dendrite structures during solidification of Al-Si alloy. Materials & Design. 28 1): 356 361 .
  34. 34. Murakami K. Shiraishi A. Okamoto T. 1983 Interdendritic fluid flow normal to primary dendrite-arms in cubic alloys. Acta Metallurgica. 31 1417 1424 .
  35. 35. Murakami K. Shiraishi A. Okamoto T. 1984 Fluid flow in interdendritic space in cubic alloys. Acta Metallurgica. 32 1423 1428 .
  36. 36. Nagelhout D. Bhat M. S. Heinricha J. C. Poirier D. R. 1995 Permeability for flow normal to a sparse array of fibres. Materials Science and Engineering A. 191 1-2): 203 208 .
  37. 37. Nandapurkar P. Poirier D. R. Heinrich J. C. 1991 Momentum equation for dendritic solidification. Numerical Heat Transfer. 19A 297 311 .
  38. 38. Nielsen O. Mo A. Applolaire B. Combeau H. 2001 Measurements and modeling of the microstructural morphology during equiaxed solidification of Al-Cu alloys. Metallurgical and Materials Transactions A. 32 8): 12049 12060 .
  39. 39. Piwonka T. S. Flemings M. C. 1966 Part VIII- Pore Formation in Solidification. Transactions of the Metallurgical Society of AIME. 236 1157 1165 .
  40. 40. Poirier D. R. 1987 Permeability for Flow of interdendritic liquid in columnar-dendritic alloys. Metallurgical and Materials Transactions B. 18 245 255 .
  41. 41. Poirier D. R. Ocansey P. 1993 Permeability for flow of liquid through equiaxial mushy zones. Materials Sciens and Engineering A. 171 231 240 .
  42. 42. Sangani A. S. Acrivos A. 1982 Slow flow past periodic arrays of cylinders with application to heat transfer. International Journal of Multiphase Flow. 8 193 206 .
  43. 43. Sangani A. S. Mo G. 1994 Inclusion of lubrication forces in dynamic simulations. Physics of Fluids. 6 1653 1662 .
  44. 44. Sangani A. S. Yao A. 1988 Transport processes in random arrays of cylinders. II. Viscous flow. Physics of Fluids. 31 9): 2435 2442 .
  45. 45. Sparrow E. M. Loefler A. L. 1959 Longitudinal laminar flow between cylinders arranged in regular array. AIChE Journal. 5 325 330 .
  46. 46. Streat N. Weinberg F. 1976 Interdendritic fluid flow in a lead-tin alloy. Metallurgical and Materials Transactions B. 7 3): 417 423 .
  47. 47. Thom A. Aplelt C. J. 1961 Field Computations in Engineering and Physics Van Nostrand, London.
  48. 48. Versteeg H. Malalasekra W. 1995 An Introduction to Computational Fluid Dynamics: The Finite Volume Method Approach, Longman Group Ltd., Harlow.
  49. 49. Wang C. Y. Ahuja S. Beckermana C. de Groh H. C. 1995 Multiparticle interfacial drag in equiaxed solidification. Metallurgical and Materials Transactions B. 26 1): 111 119 .
  50. 50. Williams J. G. Morris C. E. M. Ennis B. C. 1974 Liquid flow through aligned fiber beds. Polymer Engineering & Science. 14 6): 413 419 .
  51. 51. Worster M. G. 1991 Natural convection in a mushy layer. Journal of fluid mechanics. 224 325 359 .

Written By

S. M. H. Mirbagheri, H. Baiani, M. Barzegari and S. Firoozi

Submitted: 28 February 2011 Published: 05 July 2011