Simulation of Liquid Flow Permeability for Dendritic Structures during Solidification Process

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:

   $$\rho C_p\frac{\partial T}{\partial t}=\nabla \cdot \left(K\nabla T\right)+\rho H_f\frac{\partial f_s}{\partial t}$$ E1
   $$f_s=\frac{T_{liq}-T}{T_{liq}-T_{sol}}$$ 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:

   $$\rho H_f\frac{\partial f_s}{\partial t}=\rho H_f(\frac{\partial f_s}{\partial T})(\frac{\partial T}{\partial t})$$ E3

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

   $$\rho H_f\frac{\partial f_s}{\partial t}=\rho H_f\frac{\partial f_s}{\partial T}\frac{\partial T}{\partial t}=\rho H_f(\frac{-1}{T_{liq}-T_{sol}})(\frac{\partial T}{\partial 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:

   $$\rho \left(C_P-H_f(\frac{-1}{T_{liq}-T_{sol}})\right)\frac{\partial T}{\partial t}=\nabla \cdot \left(K\nabla T\right)$$ E5
   $$\rho \left[C_P^{eq}\right]\frac{\partial T}{\partial t}=\nabla \cdot \left(K\nabla T\right)$$ E6

   where $$C_P^{eq}$$ can be considered as a quasi-specific heat given by:

   $$C_P^{eq}=\left[C_P-H_f(\frac{-1}{T_{liq}-T_{sol}})\right]$$ 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^{mu}=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 ⁰ = T _liq – T _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_{total}\exp \left(\frac{C_m{T}_{liq}^2}{T{\left(T_{liq}-T\right)}^2}\right)$$ 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}Cos(4(\theta +{\theta }_0))if0\frac{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_{new}=r_{old}+dr;dr=\frac{d{f}_s}{2\pi {\displaystyle \sum 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-\Delta t\left\{\left(UTX+VTY\right)+\left(DQX+DQY\right)\right\}$$ E11

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

   $$UTX=0.5\left\{(1+\alpha )U_{i-1,j}(DTXL)+(1-\alpha )U_{i,j}(DTXR)\right\}$$ E12
   $$\begin{array}{l}DTXL=\frac{T_{i,j}-T_{i-1,j}}{\Delta x}\\ DTXR=\frac{T_{i+1,j}-T_{i,j}}{\Delta x}\\ DQX=\frac{QXR-QXL}{\rho C_{{\rm P}}\Delta x}\end{array}$$ 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:

   $$\begin{array}{l}QXR=-k_R\frac{T_{i+1,j}^n-T_{i,j}^n}{\Delta x}\\ QXL=-k_l\frac{T_{i,j}^n-T_{i-1,j}^n}{\Delta x}\end{array}$$ E14

   K _R and K _L are thermal conductivity and found through the following equations:

   $$\begin{array}{l}{K}_R^{-1}=0.5\left(\frac{1}{k_{i,j}}+\frac{1}{k_{i+1,j}}\right)\\ K_L^{-1}=0.5\left(\frac{1}{k_{i,j}}+\frac{1}{k_{i-1,j}}\right)\end{array}$$ E15

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

   $$\begin{array}{l}{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}}=\frac{T_{i,j}-T_s}{T_l-T_s};\\ T_s\le T\le T_l\end{array}$$ E16

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

   $$\begin{array}{l}{C}_p^{LS}=\frac{\Delta H_f}{T_l-T_s};\\ T_s\left(T\right(T_l\end{array}$$ E18

   where T _s , T _l are solidus and liquidus temperatures respectively. In finite difference form, C _p is calculated as follows:

   $$C_{{\rm P}}=C_P^L;T_{i,j}^n\left)T_{i,j}^{n+1}\right)T$$ E19
   $$C_{{\rm P}}=C_P^S;T_{i,j}^{n+1}\left)T_{i,j}^n\right)T_S$$ E20
   $$\begin{array}{l}{C}_{{\rm P}}=C_{{\rm P}}^L\frac{T_{i,j}^n-T_L}{T_{i,j}^n-T_{i,j}^{n+1}}+C_{{\rm P}}^{LS}\frac{T_L-T_{i,j}^{n+1}}{T_{i,j}^n-T_{i,j}^{n+1}};\\ T_{i,j}^n\left)T_L\right)T_{i,j}^{n+1})T_S\end{array}$$ E21
   $$\begin{array}{l}{C}_{{\rm P}}=C_{{\rm P}}^{LS}\frac{T_{i,j}^n-T_S}{T_{i,j}^n-T_{i,j}^{n+1}}+C_{{\rm P}}^S\frac{T_S-T_{i,j}^{n+1}}{T_{i,j}^n-T_{i,j}^{n+1}};\\ T_L{T}_{i,j}^n\left)T_S\right)T_{i,j}^{n+1}\end{array}$$ 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):

   $${\rho }_L\frac{D\overrightarrow{V}}{Dt}=-\nabla P+{\rho }_L\overrightarrow{g}+\mu {\nabla }^2\overrightarrow{V}$$ E23

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

   $$\frac{\partial \rho }{\partial t}+\overrightarrow{\nabla }.\left(\rho \overrightarrow{V}\right)=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):

   $$\begin{array}{l}{\overrightarrow{x}}^{*}=\frac{\overrightarrow{x}}{L},\\ t^{*}=\frac{u_{\infty }t}{L},\\ {\overrightarrow{u}}^{*}=\frac{\overrightarrow{u}}{u_{\infty }},\\ {\overrightarrow{g}}^{*}=\frac{L}{u_{\infty }^2}\overrightarrow{g},\\ p^{*}=\frac{p-p_{\infty }}{{\rho }_{\infty }{u}_{\infty }^2},and\\ \mathrm{Re}=\frac{{\rho }_{\infty }{u}_{\infty }L}{\mu }\end{array}$$ E25

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

   $$\frac{D\overrightarrow{V}}{Dt}=-\nabla P+\overrightarrow{g}+\frac{1}{\mathrm{Re}}{\nabla }^2\overrightarrow{V}$$ E26

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

   $$\frac{\partial u}{\partial t}+\frac{\partial p}{\partial x}=\frac{1}{\mathrm{Re}}\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)-\frac{\partial (u^2)}{\partial x}-\frac{\partial (uv)}{\partial y}+g_x$$ E27
   $$\frac{\partial v}{\partial t}+\frac{\partial p}{\partial x}=\frac{1}{\mathrm{Re}}\left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)-\frac{\partial (uv)}{\partial x}-\frac{\partial (v^2)}{\partial y}+g_x$$ E28
   $$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial 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:

   $$\frac{\partial (\delta \stackrel{\rightharpoonup }{V})}{\partial t}+\stackrel{\rightharpoonup }{V}.\nabla \stackrel{\rightharpoonup }{V}=\frac{1}{\mathrm{Re}}{\nabla }^2\stackrel{\rightharpoonup }{V};\delta \stackrel{\rightharpoonup }{V}=V^{*}-\stackrel{\rightharpoonup }{V}$$ E30
   $$\frac{\partial (\delta \stackrel{\rightharpoonup }{V})}{\partial t}+\nabla 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.

   $${\nabla }^{2}P=\frac{1}{\Delta t}\nabla .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=-\frac{K}{\mu }\nabla 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=-\frac{K}{\mu f_l}\nabla 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:

   $$\begin{array}{l}{U}_{i,j}^{n+1}=U_{i,j}^n+\Delta t\left[g_x-FUX-FUY+\frac{1}{\mathrm{Re}}(VISX)\right]+\Delta t\left(\frac{P_{i,j}^{n+1}-P_{i+1,j}^{n+1}}{\Delta x}\right)\\ U_{i,j}^{n+1}=A_x^n+\Delta t\left(\frac{P_{i,j}^{n+1}-P_{i+1,j}^{n+1}}{\Delta x}\right)\end{array}$$ E35
   $$\begin{array}{l}{V}_{i,j}^{n+1}=V_{i,j}^n+\Delta t\left[g_Y-FVX-FVY+\frac{1}{\mathrm{Re}}(VISY)\right]+\Delta t\left(\frac{P_{i,j}^{n+1}-P_{i,j+1}^{n+1}}{\Delta y}\right)\\ V_{i,j}^{n+1}=A_y^n+\Delta t\left(\frac{P_{i,j}^{n+1}-P_{i,j+1}^{n+1}}{\Delta y}\right)\end{array}$$ E36
   $$\begin{array}{l}{A}_x^n=U_{i,j}^n+\Delta t\left[g_x-FUX-FUY+\frac{1}{\mathrm{Re}}(VISX)\right]\\ A_y^n=V_{i,j}^n+\Delta t\left[g_Y-FVX-FVY+\frac{1}{\mathrm{Re}}(VISY)\right]\end{array}$$ E37
   $$\begin{array}{l}FUX=\frac{U_{i,j}^n}{2\Delta x}\left[DUL+DUR+\alpha Sign(U_{\alpha })(DUL-DUR)\right];\\ DUL=U_{i,j}^n-U_{i-1,j}^n;\\ DUR=U_{i+1,j}^n-U_{i,j}^n;\\ U_{\alpha }=(DUR+DUL)/2;\\ VISX=\left(\frac{U_{i+1,j}^n-2U_{i,j}^n+U_{i-1,j}^n}{\Delta x^2}+\frac{U_{i,j+1}^n-2U_{i,j}^n+U_{i,j-1}^n}{\Delta y^2}\right)\end{array}$$ 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

   $$\frac{p_{i+1,j}^{n+1}-2p_{i,j}^{n+1}+p_{i-1,j}^{n+1}}{{(\delta x)}^2}+\frac{p_{i,j+1}^{n+1}-2p_{i,j}^{n+1}+p_{i,j-1}^{n+1}}{{(\delta y)}^2}=\frac{1}{\delta t}\left(\frac{A_{x,i,j}^n-A_{x,i-1,j}^n}{\delta x}+\frac{A_{y,i,j}^n-A_{y,i,j-1}^n}{\delta y}\right)$$ 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:

   $$\frac{V^{n+1}-V^n}{\Delta t}+V^n.\nabla V^n+\nabla P^{n+1}=\frac{1}{\mathrm{Re}}{\nabla }^2{V}^n$$ E40
   $$\nabla V^{n+1}=0$$ E41

   The superscripts (n) and (n+1) denote old and new time level, respectively. In projection method V ⁿ⁺¹ 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.

   $$\frac{V^{*}-V^n}{\Delta t}+V^n.\nabla V^n=\frac{1}{\mathrm{Re}}{\nabla }^2{V}^n$$ E42
   $$\frac{V^{n+1}-V^{*}}{\Delta t}+\nabla P^{n+1}=0$$ E43

   The divergence of Eq. 26 takes the form:

   $$\frac{\nabla .V^{n+1}-\nabla .V^{*}}{\Delta t}+{\nabla }^2{P}^{n+1}=0$$ E44

   Continuity equation (Eq. 25), requires that $$\nabla V^{n+1}$$ to be zero, thus

   $${\nabla }^2{P}^{n+1}=\frac{1}{\Delta t}\nabla .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 Vⁿ⁺¹. 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 $$(\frac{\partial \text{T}}{\partial \text{t}})$$ , temperature gradient ( $$G=\frac{Q}{k}=\frac{\partial \text{T}}{\partial \text{x}}$$ ), and solidification rate ( $$v=\frac{\partial {\text{f}}_{\text{s}}}{\partial \text{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:

$$\frac{\partial \text{T}}{\partial \text{t}}=\left(\frac{\partial \text{T}}{\partial \text{x}}\right)\left(\frac{\partial \text{x}}{\partial \text{t}}\right)=\text{ G}.v,\text{ for a single phase},\text{ and}$$ E46

$$\frac{\partial \text{T}}{\partial \text{t}}=\left(\text{k}.\text{G})/(1+\text{ΔHm}.v\right),\text{ 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 (m²). 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.

$$-Log(Permeability)=7.175f_s+0.002dQ+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.