The University of Sheffield: Computer Modelling of Solidification of Pure Metals and Alloys
The University of Sheffield: Computer Modelling of Solidification of Pure Metals and Alloys
The University of Sheffield: Computer Modelling of Solidification of Pure Metals and Alloys
By
June, 1995
Computer Modelling of Solidification of Pure Metals and Alloys
Summary
Two numerical models have been developed to describe the volumetric changes
during solidification in pure metals and alloys and to predict shrinkage defects in the
castings of general three-dimensional configuration. The first model is based on the
full system of the Continuity, Navier-Stokes and Enthalpy Equations. Volumetric
changes are described by introducing a source term in the Continuity Equation
which is a function of the rate of local phase transformation. The model is capable
of simulating both volumetric shrinkage and expansion.
The second simplified shrinkage model involves the solution of only the Enthalpy
Equation. Simplifying assumptions that the feeding flow is governed only by
gravity and solidification rate and that phase transformation proceeds only from
liquid to solid allowed the fluid flow equations to be excluded from consideration.
Several modifications to FLOW-3D have been made to improve the accuracy and
efficiency of the metal/mould heat transfer and solidification algorithms.
M. R. Barkhudarov
Contents
1 Introduction 6
2
2.8 Modelling the Latent Heat Release 51
3.1 Introduction .. 66
3
5.2.1 Source Term Calculation . . . . . . . . 97
4
7 Discussion 128
7.1 Physical Models. · 128
Bibliography 152
Appendices 164
5
Chapter 1
Introduction
Casting technologies differ greatly in casting and die materials used, filling methods,
cooling rates and final product properties. Ingots and billets are castings of simple
shape, usually suitable for mechanical working, e.g. forging, rolling, extrusion and
drawing. In investment casting (or lost wax, or precision casting) a pattern of wax
or other fusible material is die cast and a plaster of Paris mould is then made round
it. When the plaster of Paris hardens the wax is burnt away, leaving a clean cavity
for the metal casting with precise dimensions and good surfaces. These and many
others processes share the same physical phenomena, that is: fluid flow of the liquid
metal, heat flow from the hot metal to the cold mould and finally liquid to solid
phase transformation of the metal in the mould [7].
Each of these processes and their interactions depend not only on the thermophysical
properties of the materials, but also on the starting and boundary conditions, i.e.
the initial metal and mould temperatures, pouring method and mould geometry.
All these factors together are called casting design.
A castings design is aimed at producing the part with the desired properties, at the
lowest cost, using available materials. One of the objectives is to produce uniform,
equi-axed, fairly fine crystals in large castings, for this structure gives metal with
the best strength and other mechanical properties [6,7,8].
Apart from grain structure, several other factors determine the quality of a cast
metal. In all casting processes care must be taken to avoid splashing the metal on the
sides of the mould, otherwise splashed solid films form and then become coated with
6
oxide which prevents them from bonding with the rest of the casting (cold shuts).
Oxide films on the liquid surface may also become folded into the metal unless the
pouring is done with care. The oxide film on aluminium is particularly troublesome,
partly for this reason and partly because, being a continuous and mechanically
strong film, it prevents the metal from filling sharp corners and narrow channels in
the mould [2]. These effects can be predicted given an appropriate numerical model
for fluid flow and free surface.
Other main factors are segregation and porosity [2,6,7,8]. Porosity in a casting is
largely defined by gas evolution and shrinkage in the solidifying metal. Predict-
ing macro-porosity due to shrinkage is the aim of the mathematical and numerical
modelling presented in this work.
Until recently casting industry relied mostly on experience, intuition and a kind of
esoteric knowledge passed on from generation to generation. But over the last 10 to
15 years there has been an increase in interaction between foundries and computer
software developing research groups. Due to commercial pressure from other metal
forming techniques and the necessity of having high quality products, the foundry
designer is looking for up-to-date analytical tools to improve the casting techniques.
The need for a more sophisticated approach to the design of the manufacturing pro-
cess is also dictated by the increased range of new alloys and products manufactured
7
using casting technology for which no previous experience is available.
It is widely accepted that the use of a computer model which simulates the physi-
cal phenomena of a casting process could save time and resources during the design
stage. It helps to reduce the risk of misruns, defect formation etc. Due to the breath-
taking achievements of the computer industry and equally wonderful developments
in the computational fluid dynamics there are a number of commercially available
software packages capable of reliable, robust and sufficiently fast simulations of a
remarkable range of fluid flows: viscous and inviscid, compressible and incompress-
ible fluid flows with transient free surfaces and shock waves; non-Newtonian fluids,
stress analysis, heat transfer are also included into some numerical models. More
and more computer codes are produced specifically to describe solidification and
they are successfully used on a regular basis.
However, by no means the task of modelling the casting process fully, e.g. mI-
crostructure and defect formation, can be considered completed. Accurate simula-
tion of filling stage is still considered a difficult task and in many cases skipped,
assuming an instantaneous filling. Accuracy of cooling and solidification solution
requires an adequate description of the mould geometry, metal/mould interface and
a knowledge of the material properties. Finally, the phase transformation and defect
formation often involves microscale phenomena and pose a problem of coupling the
corresponding micromodels to the standard macromodels used in computer simula-
tions.
The existing computer tools differ in capabilities, accuracy, applicability, and none
of them models the process in full. There is still an increasing demand for a uni-
versal and reliable analytical tool which would both increase the productivity of the
industry and help to create an insight into physical phenomena occurring during
the solidification process. Here is an incomplete list of what a foundryman would
expect from a computer model to do according to Piwonka [1]:
• establish the temperature distribution in the mould and the liquid metal at
the start of solidification;
8
• predict the grain size and shape;
• predict the distribution of gas content within the casting;
- mould filling with viscous fluid including free surface and solidification effects;
Attention is being specifically paid to volumetric shrinkage since it is one of the main
causes of macro-defects in castings and there is still no adequate computer model
which would reliably predict these defects. Fluid flow is inlcuded into consideration
to investigate its influence on the defect formation.
The proposed mathematical model will be able to predict the formation of casting
defects due to shrinkage of the solidifying material involving the formation of macro-
cavities. The work is based on a proprietary general purpose software FLO W-3Dl.
The code is capable of modelling the first four items listed immediately above though
some improvements in the heat transfer and solidification models have been made
as a result of the present work.
Numerical modelling of mould filling and solidification can be divided into four
major parts:
• free surface tracking which concerns mainly the filling stage when there are
multiple and essentially transient free surfaces;
lThe code is a property of FLOW SCIENCE Inc., Los Alamos, U.S.A.
9
• heat conduction and heat transfer during and after filling; here the major
computational effort goes into resolving the metal/mould interface, given that
almost all the heat from the metal is lost to the mould. Besides, if the metal
and mould thermal properties differ greatly, then steep temperature gradients
are expected near the interface which also need to be accurately accounted for;
• solidification during and after filling; here attention is paid to the way the
latent heat is accounted for.
Fluid flow modelling is far more difficult than heat transfer calculation since fluid
motion requires a solution of three momentum equations (in three dimensions) and
a continuity equations which are all coupled, even though in terms of the real time
filling takes only a small fraction of the whole casting process. Nevetherless it is
widely accepted that filling affects the subsequent solidification and possible defect
formation because [2]
• heat transfer during mould filling affects the temperature distribution in the
casting after the filling. This is especially true for thin wall castings where a
great deal of the superheat may be removed during filling [2]. Without going
through the mould filling process, it is impossible to ascertain correct initial
conditions for subsequent solidification analysis;
• fluid flow modelling of the liquid melt plays a critical role in optimising the
design of a casting system which consists of a mould, sprue, runner and ingates,
in order to eliminate mould erosion and porosity problems. It is important
to know how the filling proceeds, i. e. whether the metal splashes or any
premature freezing occurs. Cold shuts and misruns can only be avoided by
knowing the detailed flow field and heat transfer characteristics of the liquid
melt during the mould filling process. Large castings, though, in which filling
takes only a small fraction of the whole time, may be less sensitive to the
initial fluid flow.
• fluid flow and heat transfer have a profound effects on grain size, macro- and
micro-porosity and segregation of alloying elements of a casting. Fluid flow
occurs during solidification after the mould is full due to residual circulation,
thermal convection and feeding.
• foreign inclusions, oxide films, introduced during filling stage, may also be a
cause of a poor quality casting.
10
The shrinkage formation model developed in this work is based on full Navier-Stokes
and continuity equations coupled with a solidification model. Pressure and velocities
are dependent on the volumetric changes occurring during the phase transformation.
The model would allow one to simulate the effects of residual circulation and/or ther-
mal convection on the solidification and to predict velocities in the feeding flow. The
use of the full fluid flow equations makes the model most general, e.g. volumetric
expansion due to remelting can also be easily handled. The numerical implementa-
tion of the model is based on the VOF method which uses the fluid fraction function
to represent free surfaces and voids.
A brief overview of the numerical methods for solving Navier-Stokes equations with
free moving boundaries, heat transfer and phase transformation is given in Chapter
2.
The outline of the general numerical method, including equation discretisation, free
surface representation and the solution procedure, is given in Chapter 3. Those are
the basis of FLOW-3D and form the framework in which the numerical modelling
is carried out. The latter includes modifications of the metal/mould heat transfer
algorithm which improved the accuracy of the calculations and enhancements to the
solidification model, presented in Chapter 4.
Chapter 5 contains mathematical and numerical equations for two shrinkage models:
a full model (Ml) which is based on Navier-Stokes equations, and a simplified model
(M2) which employs the energy conservation equation and does not describe fluid
flow during solidification.
Simulations results are presented in Chapter 6 and are compared with the results of
experimental castings available from literature. A detailed discussion of the results
and numerical problems encountered during calculations is given in Chapter 7.
11
Chapter 2
We shall be dealing with the flow, heat transfer and phase change of gases and liquids
(melts) commonly termed as fluids. The Eulerian description of the evolution of
fluid at each point in time and space of the physical space will be used rather than
the Lagrangian method where the evolution of individual physical fluid volumes
is described. The main difference between the two methods is that an Eulerian
coordinate system is fixed in the physical space while a Lagrangian coordinate system
is attached to the fluid particles and moves and deforms with the fluid [3].
The Eulerian hydrodynamical equations describing a combined fluid and heat flow
can be written in the following general differential forrnt, e.g. [3,4]:
Continuity Equation
-ap
at = - V'. (pv) (2.1)
IThat is assuming that all functions in these equations are continuous and have continuous
first derivatives. In a more general case of a flow with discontinuities the integral form of the
conservation equations will have to be used [3].
12
Momentum Equation - Newton's Second Law
(2.2)
at
conservation of energy
aE - (v . V)E
1
+ - (P . V) . v + -
A d qe
(2.3)
at P dt
where
The Continuity Equation, Eq. (2.1), states that a change of the density in a fluid
element is equal to the total flux of the fluid advected (or convected) into the element
and it is also a statement of the conservation of mass.
Momentum Equation, Eq. (2.2), states that the fluid momentum in a fluid element
changes due to the advection of the momentum into the element, action of external
sur face forces P, and of external body forces F.
p = p(p,T)
13
where
eij = ! (av j + aVi)
2 aXi aXj
is the shear rate tensor; p the normal stress, or pressure, and f the shear stress tensor.
In general, p and f may depend on some other physical parameters, e.g. temperature
and phase concentrations, but we are not considering these complications here.
The dependence of f on e defines the viscous model of the fluid. If the values of
eij are small then Tij can be expanded by Taylor series about eij = 0 and retaining
only linear terms. In this case the shear stress is a linear function of the shear rate
and coefficients Bijkl are dependent on temperature and other physical parameters.
All gases and simple liquids, molten metals and slags obey the linear viscosity model
which are termed Newtonian fluids.
Assuming further that the fluid is isotropic (that is, its properties at a point are
independent of the geometrical direction), components of tensor 13 are reduced to
only two independent parameters and the final form of the stress tensor P is
(2.5)
where>. and J.L are constant and positive coefficients, J.L is the fluid dynamic viscosity
coefficient and TJ = >. + ~J.L is the second viscosity coefficient. Momentum Eq. (2.2)
for constant >. and J.L can be written as
with v = J.L/ P being the fluid kinematic viscosity. Eq. (2.6) is called the Navier-
Stokes equation.
Eqs. (2.1) and (2.2) may be simplified further by assuming the incompressibility of
the fluid, i. e.
p = Po = const (2.7)
Then \1 . v = 0 in Eqs. (2.5) and (2.6).
The density of the fluid does change due to thermal contraction/expansion (linear
contraction) which results in thermal, or natural, convection in the fluid in the
presence of body forces. The thermally induced flow can be an important factor,
14
for example, in redistributing heat, and can be taken into account by the use of the
Boussinesq approximation:
p = Po + 8p, 8p ~ 1
po
Then P = Po in all the terms of Eq. (2.6) except for the pressure gradient term,
where variable density is retained. A simple linear dependence of 8p on temperature
can be used to complete the model:
8P = - {3 Po (T - To)
where To is the reference temperature and
(3 = J:.
V
(av)
aT p
The Energy Conservation Eq. (2.3) states that a change of the internal energy E
in a fluid element is the result of the advection of E into and out of the element,
work due to the internal surface forces, -!CP. \7) ·V, and a flux of the external heat
dqe /dt. Using the Continuity Equation and Eq. (2.4) we can write
1 dl/ p 1
-Pdt + pT;je;j
A
The second term on the righthand side of Eq. (2.8) constitutes the heating of the
fluid due to the viscous friction.
We will assume that viscous heating is negligibly small compared with the advection
and the external heat flux terms in Eq. (2.3). This is true for all fluids for the range
of flows considered in this work. Viscous friction may become important in high
pressure die-casting flows [5].
External heat flux can be written in terms of a vector quantity using simple geo-
metrical considerations:
dqe _ Id·
dt -
-- wq
p (2.9)
According to Fourier's law, which is derived similarly to the linear viscous fluid
model [3], the heat flux is expressed as a linear function of the the temperature
gradient in the fluid:
q = -k\7T (2.10)
where k is the fluid thermal conductivity; it is positive and it may be a function of
temperature and other physical parameters.
15
Eq. (2.3) can be rewritten now as
aE - (v. 'V)E - p (d1/p) + ~'V. (k'VT) (2.11)
at dt p
Strictly speaking a system can reach the equilibrium state only if in does not interact
with other systems and is given sufficient time to allow the system to be in the
equilibrium at each intermediate point of the process. In reality the interaction
processes, such as metal-mould heat flow, occur at finite speeds. In that case we can
define a local equilibrium so that each infinitesimally small volume of the system
is at equilibrium with its own temperature and pressure, though the system as a
whole is not. This assumption is valid if heat fluxes across these volumes are not
too large. According to Fourier's law, Eq. (2.10), the fluxes are small if temperature
gradients are sma1l 2 •
Large velocities during filling stages can lead to the development of turbulence in
the metal flow. Reynolds number in such casting processes as gravity and pressure
die castings can reach values of 10 5 implying that turbulence levels may be very
high. Turbulence in the metal affects velocity distribution in the flow, filling time,
free surface motion and heat transfer [7].
16
au") -)
au'·
f.=v
( -ax·
au~' +_J
J ax' ax· , I
where Vo is the molecular value of the kinematic viscosity and VT the turbulent
viscosity
k2
VT = Cv - f.
In the k - f. model typically C v = 0.09, as derived from experimental studies, while
in the RNG k - f. model Cv is a function of local shear rates. In the latter case the
model is more accurate and suitable for metal flow conditions in complicated runner
and mould systems of typical industrial castings.
dH = TdS + dp (2.12)
p
17
or
(2.14)
where dHdm is the heat given to or taken from the surroundings of the element, called
latent heat. Eq. (2.14) shows that though phase transformation occurs isothermally,
it is accompanied by heat flow leading to a change in the enthalpy. The latent heat
of solidification L per unit mass is
L = dHdm (2.15)
dm
The latent heat of melting can be thought of as the energy required to pull the
atoms apart to the more openly packed structure of the liquid.
It follows from Eqs. (2.11) and (2.12) that the fluid specific heat at a constant
pressure, Cp , is needed to evaluate the enthalpy:
H = {T CpdT (2.16)
iTo
while for the internal energy the specific heat at a constant volume, Cv , is required:
E = {T CvdT
}yo
Obviously, for liquids and solids, it is much easier to measure Cp than Cv , and the
advantage of using the enthalpy instead of the internal energy is apparent. Taking
into account the latent heat we have:
(2.17)
The fraction of solid function fs has to be defined so that Eq. (2.17) can be solved for
temperature. If solidification is assumed to be in equilibrium, then for pure metals it
will proceed at a constant temperature Tm and fraction of solid is uniquely defined
by the enthalpy. For alloys the equilibrium occurs over a range of temperatures.
The lever rule for binary alloys then relates the initial alloy composition, Co, and the
current values of compositions of the liquid, CII and solid, Ca, to is [6,7]:
Co - Cl
fa = - -
Cs - c, (2.18)
If the phase diagram of the alloy is known then compositions Cl and c, can be found
from the diagram given a value of the current temperature.
18
There are several factors, however, that may complicate this procedure.
When a new phase is being formed in the system, this phase often appears first as
small nuclei in the old phase, which then grow by the addition of more material
from the old phase. The surface energy dominates the nucleation process because
of the small size of the nuclei and its high surface curvature. Experimental studies
show that in a pure metal free from inclusions and away from the walls a substantial
undercooling is required to promote solidification - approximately 0.2Tm, e.g. 295°K
for pure iron [6].
Once an undercooled droplet is nucleated, freezing is then very rapid. This differs
from the casting of metals under practical conditions in a foundry where nucleation
occurs at temperatures within lOOK of the ideal freezing point due to the presence
of foreign inclusions and mould walls, which serve as nucleation centres. The rate
of growth of the solid on the nucleus in this case is much smaller, as it is controlled
by the rate of removal of the latent heat released at the solid-liquid interface.
Recently a number of 'macro-micro' models have been proposed to predict the mi-
crostructural evolution of a casting [11]-[16]. Comparisons of the predicted grain
structures, using a rather sophisticated probabilistic model, with the experimental
results for dendritic and eutectic alloys, given by Gandin et al [17] are very impres-
sive. However, the main difficulty in applying this model in a general case is the
requirement of a very fine spatial resolution: the maximum computational cell size
for 2-D and 3-D simulations was 10 J.Lm which is unacceptable for a 3-D industrial
scale fluid flow and solidification modeling. Besides, modelling of the solidification
kinetics requires additional input parameters, e.g. the relationship between the un-
dercooling!:!:..T and the number of grains nAT in different parts of the casting. These
19
are either defined from some assumptions and/or obtained empirically. The latter
makes kinetics models case-dependent and complicates its validation. Besides, it is
difficult to take into account the influence of the metal flow on the grain evolution,
the importance of which has been shown by Ortega [18].
In alloy systems solidification speed may be comparable to the rate of solute diffusion
and advection. In this case the distribution of the solute in both solid and liquid
phase will be non-uniform and this will affect the local temperature of solidification.
To model these phenomena the solute transport equation will need to be solved for
each phase [19]
-ac
at + ('\7 . v) c = '\7. (D'\7 c) (2.20)
where c is the solute concentration and D the diffusion coefficient. Given the values
of the composition and temperature the value of the solid fraction can be found
from the alloy phase diagram [7].
divv = 0 (2.21 )
Bv 1
- v· '\7v - - '\7p + v~v + -1 F (2.22)
at po po
aH
at = - (v . '\7)H +
1
po'\7· (k'\7T) (2.23)
Elliptic Eq. (2.21), parabolic Eqs. (2.22), (2.23) and Eq. (2.17) are a system of
non-linear, coupled second order partial differential equations written in terms of
the fluid velocity, pressure and temperature as the unknown variables. Only one
boundary condition for the velocity and for the temperature at each fixed boundary
will be sufficient to define the solution of the these equations; the conditions can
be either of Dirichlet type (a specified velocity/temperature value) or of Neumann
type (a specified value of the velocity gradient/heat flux normal to the boundary)
or a linear combination of both [20,21]. Velocity, pressure and temperature are then
defined uniquely away from the boundaries given the appropriate initial conditions.
Boundary conditions specific to the viscous fluid model are kinematic no-slip (Dirich-
let) conditions at the walls
v=O (2.24)
20
In the presence of a free surface, which is defined as a boundary between two fluids
with very different densities (for example, for water and air Pw/ pa ~ 1000), the
lighter fluid (gas) can be substituted by a constant pressure boundary condition
at the free surface. The exclusion of the gas phase from consideration poses an
additional problem of setting the boundary condition at the generally transient free
surface.
- pressure in the fluid PIs results from the sum of the gas pressure, surface tension
force and viscous stresses and can be written as
(2.26)
where
Vn = V· nj
Zero flux of the tangential momentum is assumed through the free surface. For
tangential viscous stress this gives
TT = 0 (2.27)
where the subscript 'T indicates the projection onto the plane tangent to the free
surface. This assumption is well justified for casting problems due to generally large
Reynolds numbers in the liquid metal flow. An example of a case when tangent
stresses are important is the interaction between ocean surface waters and atmo-
spheric winds.
Heat flow boundary conditions at the walls and at the free surface can be specified
in terms of the heat flux qbc per unit area. The same expression can be used in both
cases, that is:
(2.28)
21
where
h is the heat transfer coefficient which characterises the thermal resistance of the
interface and
The value of h varies with temperature and geometrical properties of the interface,
such as the wall surface finish and the size of an air gap between the metal and the
mould. At the free surface h is defined by the convection conditions in the gas, such
as turbulence intensity, whether it is natural or forced and the free surface orienta-
tion. Radiation heat transfer coefficient is strongly temperature dependent [22] and
is influenced by the presence of other radiating surfaces. Ludley and Szekely showed
that for a molten steel held in a ladle the heat loss from the unprotected surface
was of the order of 20 to 30% of that lost to the ladle walls by conduction [23].
There are few cases for which analytical expressions for h are available [22] and
normally experimentally defined values are used [24].
Numerical methods have to be used to solve Eqs. (2.21)-{2.23) since these are non-
linear, transient, coupled second-order equations. There are four widely adopted
basic approaches to obtain a numerical solution of these equations: the finite dif-
ference, finite element, the finite volume and the boundary element method [25].
22
The finite difference method (FDM) is based on the properties of Taylor expansions
and on the straight-forward application of the definition of derivatives. It is the
oldest of the methods applied to obtain numerical solutions of differential equations,
and the first application is considered to have been developed by Euler in 1768.
The idea of finite difference methods is quite simple, since it corresponds to an
estimation of a derivative by the ratio of two differences according to the definition
of the derivative.
For a function u(x) the derivative at a point x can be approximated, for example,
by forward differencing:
au =u(x+box)-u(x) + O(box)
ax box
indicating that the truncation error O(x) goes to zero as the first power of box.
Most finite difference models used in solidification simulation have employed regu-
lar spaced rectangular meshes, which are not particularly suitable for representing
complex shaped castings, although curved surfaces can be approximated by using a
fine node spacing and step approach [28J (see also FAVOR method in Section 3.3.2).
For curvilinear meshes the discretization of the equations can be performed after a
transformation from the physical space (x, y, z) to a Cartesian, computational space
(e, 7], () [29,30]. This method of body fitted coordinates combines the flexibility of
the finite element domain discretisation and the simplicity of the finite difference
equation approximations [31,32].
This is the technique in which the integral formulation of the conservation laws are
discretized directly in the physical space. The finite volume method (FVM) takes
advantage of an arbitrary mesh, where large number of options are open for the
definition of the control volumes around which the conservation laws are expressed.
23
The method uses the conservation equations in their general integral form. For
a scalar quantity U, with volume sources Q, and a flux vector F, which contains
only convective contributions, the integral equation for a discrete volume 0 with
boundary S is given by
(2.29)
The essential significance of this formulation lies in the presence of the surface
integral and the fact that the time variation of U inside the volume only depends on
the surface values of the fluxes. The advantage, especially in the absence of source
terms, is that the fluxes are calculated only on two-dimensional surfaces. Eq. (2.29)
is replaced by the discrete form
(2.30)
index J indicates that the quantity is averaged over the volume OJ.
The finite element method (FEM) originated from the field of structural analysis
as a result of research carried out mainly between 1940 and 1960. The concept of
'elements' can be traced back to the techniques used in stress calculations, whereby a
structure was subdivided into small substructures of various shapes and re-assembled
after each 'element' had been analysed. The development of this technique and its
formal elaboration led to the introduction of what is now called the finite element
method by Turner et al [33] in a 1956 paper dealing with the properties of a triangular
element in plane stress problems.
24
After its successful applications to structural mechanics it soon appeared that the
method could also be used to solve continuous field problems [34] including solidifi-
cation [35].
The main advantage of the finite element method is the ease with which complex
geometries can be handled and the wall boundary conditions can be set more ac-
curately. The flexibility regarding node size and shape means that curved surfaces
can be accurately meshed.
Apart from setting the complex task of generating the best mesh, finite element
method has other disadvantages. Generally, it requires more computer time for
calculations than the FDM. Besides, any physical problem written in finite element
formulation loses its clarity, the physical nature of the initial equations becomes
obscured by the mathematical symbolism. Despite this and supported by successful
developments in computer technology, the finite element method is widely used
at present, both in commercial and scientific programming [38]-[40]. It has been
successfully applied to model casting process in 2-D [35,41] and 3-D [42,43] and
commercial 3-D FEM packages, such as ProG AST [44,45], have also been marketed.
An alternative to the more familiar finite element, finite difference and control vol-
ume approaches is the boundary element method (BEM). With this method, the
differential equations describing the flow are recast into integral representations so
that the solution away from the boundaries is written in terms of the boundary
values. The solution of the flow problem is then based on the numerical quadrature
of the integrals. It has been shown that this method can be successfully applied to
phase transformation and free surface problems with moving boundaries, as well as
to general viscous flows [25,39,46,47]. Detailed description of the boundary element
method can be found in these references.
25
2.3 Discretisation of the Equations of Motion
The finite volume approach in structured grids will be used in this work to solve the
fluid motion and energy Eqs. (2.21)-(2.23). The choice is dictated by the use of a
commercial CFD code FLOW-3D which employs this method (see Chapter 3).
26
2.3.1 Advection and Diffusion Terms Approximations
and
The evaluation of
OU 2 ovu n uk - uL (VU)T - (VU)B
n:,x,i,i = (ox + Oy) Ii,] = L\x + ~y
Numerical stability analysis shows that the MAC method is unconditionally unstable
for inviscid flows and requires a non-zero viscosity to maintain stability [27]. More-
over, central differencing for the advection and viscous terms introduces a limitation
to the cell Reynolds number which for a one-dimensional flow is
uL\x
Re c = - - , (2.33)
11
The control volume, upwind differencing approach to solving Eqs. (2.17), (2.21)-
(2.23) is widely used [5,49,50,51]. The upwind differencing scheme was first used by
27
Gentry et al [52] and Hirt et al developed it for incompressible flows and named the
method SOLA (SOLution Algorithm, see Section 2.3.2) [53].
Fig. 2.1c shows a computational cell and its control volume (CV) for calculating
the x-direction momentum fluxes. Assuming for simplicity a uniform mesh spacing
~x, at position (i,j) we have:
and
(2.36)
are the estimates by averaging of the velocities at the control volume boundaries.
Viscous terms are discretised by central differencing, in the same way as in MAC
method.
The upwind differencing for the advection terms removes the cell Reynolds number
restriction imposed by Eq. (2.33) for viscous flows and is conditionally stable for
inviscid fluids [27]. The essence of one-sided differencing scheme given by Eqs.
(2.34)-(2.36) is to "advect" the value of a variable (velocity) taken upstream from
the considered node. Sometimes the method is also called the donor cell method.
If for example UR in Eq. (2.34) is negative then Ui+1,j, instead of ui,i" is used as the
quantity advected through the right boundary of the control volume.
• it possesses the transportive property which means that the effect of a pertur-
bation is only advected in the direction of the velocity. All methods which use
centred-space derivatives for the advection terms do not possess this property
and filters have to be designed to avoid spurious solutions.
• compared with the central differencing method, the upwind method is not
stability limited by cell Reynolds number. Vanka [54] used the advantages
of both schemes by applying the upwind method only in the cells which had
Re c ~ 2 and the central differencing for all other cells, thus obtaining a
"hybrid" method.
28
Roache [27] mentions a surprising agreement of the upwind "donor cell" method
with some second-order calculations of the driven cavity problem. In fact the term
"accuracy" of a numerical method is rather arbitrarily defined. In most cases it is
understood in the formal Taylor series expansion sense. It is possible to represent
spatial derivatives "more accurately" in this sense but not to retain for example the
conservative property.
and
1
Vy = 2v~y(1 - cy )
This disadvantage of the first order upwind differencing makes it inapplicable for
the advection of variables with sharp discontinuities in the flow region, such as
free surfaces. Section 2.5 describes two numerical methods specially developed to
preserve the discontinuities.
It has been shown that the "donor cell" upwind method demonstrate good com-
parisons with physical experiments at high Reynolds numbers [27]. In calculation
of driven cavity flow, Torrance et al [55] have shown that the upwind differencing
method applied to the conservative equations is considerably more accurate than
second-order differencing applied to the non-conservative equations.
If advection and/or viscous terms on the left-hand side of Eqs. (2.31) and (2.32)
are estimated at time t n +b i.e. implicitly, then these equations become coupled and
have to be solved at all nodes simultaneously to obtain the first-guess velocities u(j+l
and V(j+l. This is achieved by inner iterations on the velocities. The inner iterations
usually take much less iterations to converge then the outer, i.e. pressure iterations.
Roache [27] suggests that ideally the advection terms should be approximated ex-
plicitly and the viscous terms implicitly. This would reflect the physical nature of
these processes, that is the finite velocity of advection and the infinite velocity of
diffusion 3 •
3The infinite diffusion velocity means that if there is a local perturbation in a uniform distri-
bution of a variable, then it will instantaneously affect all areas due to the diffusion.
29
In the present work all terms on the right-hand sides of Eqs. (2.31), (2.32), except
for the pressure gradients, are approximated explicitly, that is at time tn. In this
case the solution matrix for uo+I and vo+ 1 is diagonal and the first-guess velocities
are easily calculated at each node. The explicitness introduces some limitations on
the time-step size ilt which are discussed in Section 2.10.
As has already been mentioned, the upwind differencing method has first order
formal accuracy. It is well known that higher order methods, both in space and time,
pose a problem of resolving the boundary conditions without losing accuracy [27].
Often only first order equations can be used at the boundaries. The first order
truncation errors propagate from the boundaries into the computational domain,
depriving the solution of its higher-order accuracy. One of the specific features of
casting is that the ratio of the internal mould walls area, at which the boundary
conditions are set for both the heat and fluid flows, to the volume of the open-to-
flow domain is large. This fact reduces the advantages of employing higher-order
accurate discretisation methods for the heat and fluid flow equations and justifies
the use of the first order numerical method.
Among all other numerical methods, the Leith's approximation scheme for the ad-
vection terms is worth of mention since it is fully explicit but is unconditionally
stable [27]. It is also called the semi - Lagrangian method [56,57]. The details and
an application of the method are given in Appendix B.
Due to the incompressibility of the fluid pressure in the Navier-Stokes equations and
velocity components in the Continuity Equation have to be approximated implicitily.
This leads to the following set of algebraic equations (in two dimensions)
1 n+l n+I
An Pi+I,j - Pi,i
u· . - - ---'--"'----'''-- (2.37)
1,1 PO ilx
(2.38)
where uf,j and v~i include terms approximated at time level n. Pressure gradient
terms in Eqs. (2.37) and (2.38) are approximated using the cell-centred pressures
as they lie exactly on the control volume boundaries (see Fig. 2.1c). If we choose
a computational cell (i,j) of a staggered grid (see Fig. 2.1b) as the control volume
30
for the Continuity Equation, then its discretised form at time t n +! is
L.. -
I,) -
_1_(u n+1 _
~x i,j
n+l)
U i - 1,j
+ _1_( n+l
~y Vi,j
_ n+l) - 0
Vi,j_l- (2.39)
The iteration procedure to bring the velocity divergence to zero in every cell and
update pressures is a major part of the computational effort in obtaining numerical
solution. Pressure iterations may take as much as 90% of the total CPU-time [5J
and much attention is paid to its optimisation.
1+1 I ~t
Ui_l,j = Ui_l,j ~flpI+l (2.43)
POUX
I+! I ilt
v··
I,)
= v·I,). +- - ~PI+l
Poily (2.44)
ilt
~flPI+l (2.45)
Pouy
For each cell encountered, the divergence Li,j is computed using the most current
velocity values available. The iterations continue until I L 1< f. This technique is
called Gauss-Seidel (GS) iteration method. The technique which uses only previous
iteration values to compute the new values through the whole iteration cycle is called
Jacobi method.
Lipinski et al [5] argue that one should use Jacobi type pressure iterations to assure
symmetric results. When LL; is calculated in a control volume using GS method
and the mesh is swept in x-direction from i = 1 to i = imax, then Ui-l,j has already
been iterated in i - 1,j cell while ui,i has not. This introduces an asymmetry in
symmetric flows. GS method is symmetric at the level of convergence, but when the
results on the way to convergence are allowed to be non-symmetric problems can
occur due to truncation errors eventually leading to a non-symmetric final solution.
31
The effect of this asymmetry could be reduced by alternating the direction of the
mesh sweeping at different time steps.
However, no such problems have been encountered using GS method in the present
work. This method has also the advantage of smaller computer memory required
for storage, though Jacobi method would be easier to vectorise and/or parallelise.
The two iteration methods described above represent the point-by-point (PP) iter-
ation procedure4 • The main drawback of these iteration schemes is that the lowest
wavelength error components damp most slowly. Thus, regardless of the initial er-
ror distribution these components will dominate for large iteration number 1 [27].
Frankel [58] showed that, for the most resistant error components, asymptotically 1
GS iterations are worth 21 Jacobi iterations. In other words the latter converges at
half the speed of GS method [59].
The are a few variations of the above procedures aiming at achieving faster conver-
gence.
In the line - by - line (LL) method a line of cells is chosen. It is assumed that
pressures along the neighbouring lines are known from their "latest" values. The
correction values for p's along the chosen line are calculated by directly solving the
tridiagonal simultaneous linear equations using available efficient solvers.
In the block - by - block (BB) iteration method, developed by Wang [49], four cells,
2 x 2, are chosen as a group. Assuming again that p's in the surrounding cells are
known from their "latest" values, the equations for the four pressures in the chosen
cells are solved simultaneously by Cramer's rule or direct matrix inversion. In cells
which cannot be organized into groups a point-by-point method may be used.
o< w < 2
The optimum value of w which gives the fastest convergence lies in the interval
between 1.0 and 2.0 and depends on the coefficient matrix of the linear system. For
the incompressible fluid problems a value of w between 1.7 and 1.8 is often optimum.
4The Jacobi method is also known as Richardson's method or method of simultaneous dis-
placements, and the Gauss-Seidel algorithm is also known as Liebman's method or successive
replacements method.
32
It can be shown that the convergence acceleration due to over-relaxation is similar
to choosing a maximum time step size in explicit methods for solving transient
equations [27].
The PP, LL and BB methods were used to solve a simple two-dimensional problem in
a square domain covered by N x N square mesh cells [60]. On one of the boundaries
of the domain pressure was fixed at a constant value while the pressure inside the
domain was given an arbitrary initial value. Then one of the three iteration methods
were used to obtain the solution, which is obviously a uniform pressure equal to the
boundary value.
• The BB method takes the shortest computing time to obtain the solution; and
the LL method needs the longest computing time among the three methods.
• The number of iterations required are almost the same for LL and BB methods.
The iteration number needed for the PP method was always larger, and more
so for smaller values of w. The reason for this is that the information diffuses
slowly relative to the other techniques because the pressure value of a node is
influenced only by its neighbours.
• Different values of w were tested. For each iteration method there is an op-
timum value of w with which minimum number of iterations is needed. This
value is about 1.6 to 1.7.
In a flow problem which involves transient free surfaces and moving obstacles, use of
the BB and LL methods becomes difficult because the coefficient arrays at different
locations are not the same. Additional computing time is required to determine
these coefficients. Therefore, the PP method appears to be the most efficient for
these problems.
A widely-used variation of the line-by-line method for solving the linear algebraic
equations is alternating-direction-implicit (ADI) scheme [61]. This method con-
verts the line-by-line routine into a direct solver of two-dimensional linear problems
through splitting the time step in halves (for 2-D, or in three equal parts for 3-D)
and judicious alternation of the direction in which the triagonal matrix algorithm
is employed. If the ADI method is used in the y-direction, the iteration consists
33
of a sweep through all i cells solving for the pressures and velocities in the j-index
direction.
In SOR methods, the number of iterations required for convergence, lmate' increases
with N, the total cell number. For ADI methods applied to square regions, lmate is
almost independent of N, so that for large enough N, ADI methods are preferable.
In the numerical experiments of Birkoff et al [63] ADI methods were nearly four times
faster than the optimum SOR method in a 40 x 40 grid. But for non-rectangular
regions the ADI methods are not certainly known to be faster, and the SOR methods
are easier to program.
Another comparison of the PP, LL and ADI methods was carried out by Chuan et
al [64]. They considered the Stephan problem.
The finite volume, fully implicit method was used to obtain the system of algebraic
equations that represent transient 3-D heat conduction with phase change, modelled
by the enthalpy method. The material used was water, freezing at Tm = DOC, though
an artificial "mushy zone" was used for the reasons discussed in Section 2.8.
The computational domain was a cube with three adjacent cold walls and the other
three being adiabatic. The resulting equations were solved by fully implicit al-
gorithms in one, two and three dimensions. The iteration involves the unknown
temperatures and the position of the solidification front at the new time step.
• For this type of non-linear problem the LL and PP solvers perform equally
well in terms of CPU efficiency and accuracy. Only for the three-dimensional
case does the ADI method substantially outperform the other two procedures.
• The relation between CPU time and time step size is non-linear. The use of
too large a time step affects accuracy while providing little CPU savings.
Chuan et al also concluded that the enthalpy method in combination with the finite
volume method produces an effective model to simulate freezing.
34
The simplicity and flexibility of the point-by-point SOR method in combination
with the finite volume upwind differencing scheme produced an efficient solver for
incompressible transient fluid flows called SOLA [53]. Its extension to include free
surfaces, SOLA-VOF [51] (see Sections 2.4 and 2.5), is now widely used in mould
filling simulations [49], [65]-[67].
Vanka [54] applied the block - implicit multigrid solution algorithm for incompress-
ible flows. The multigrid technique is based on the idea that each frequency range
of error must be iterated on an optimised grid. For a given grid the error compo-
nent wavelengths comparable with the mesh size are smoothed out most efficiently.
The multigrid technique cycles between coarse and finer grids until all the frequency
components are appropriately smoothed. For a predefined set of grids, the solution
is initiated on the coarsest grid M. GS iterations are performed on grid M until
it is fully converged. The converged solution is then interpolated to the next, finer
grid and is iterated on until the required convergence level. When the finest level
is solved to the desired accuracy, then the whole solution cycle is terminated. In
this formulation, when iterations start on a given grid the long wavelength error
components have been already smoothed on the coarser grids.
The multigrid technique has been shown to be one of the most efficient solvers
for linear and nonlinear equations [25,68] though the procedure of choosing the
grid hierarchy and the interpolation of the variables between them requires some
delicate programming [69]. It is unclear, though, whether this method will speed
up the smoothing of the error components with wavelengths larger than the scale of
the computational domain, e.g. uniform fluid pressure adjustments.
35
by line segments. This method removes the limitation of single-value surfaces,
but the extension to three dimensions is extremely difficult.
3. Marker Particles Method [48]. This is also called Marker And Cell method
(MAC). In this method marker particles are spread over all fluid-occupied
regions, with each particle specified to move with the fluid velocity at its
location. The MAC offers the advantage of eliminating all logic problems
associated with intersecting surfaces and can readily be extended to three
dimensions. However, it requires a significant increase in computer storage
and running time.
A number of variations of the MAC method exists, e.g. Simplified MAC
method (SMAC) and Leading Marker Method (LMM). These variations aim
at saving the computer storage and cutting the computer running time. Lin
and Hwang [71], for example, store and compute only particles in the surface
cells and those interior cells that are next to the surface. In this formulation the
method cannot describe the opening of new voids in the fluid. An algorithm
has to be designed for introducing new particles if the number of surface cells
and their neighbours increases sharply during filling. Both are more likely to
occur in a three-dimensional flow.
4. Volume of Fluid Method [51]. The Volume Of Fluid (VOF) method uses a
continuous function to describe fluid configuration. It evolved from the original
MAC method. The F(x, y, z, t) function has the unity value at any point
occupied by fluid and zero otherwise; it preserves its discontinuous nature to
track a free surface boundary. When averaged over the cells of a computing
mesh, the average values of F is equal to the fractional volume of the cell
occupied by fluid. Cells with F value between zero and unity contain a free
surface. The VOF method requires much less computer storage space and
running time, and has none of the logic problems of handling intersecting
surfaces in three dimensions.
The flexibility of the VOF method has enabled its use in combination with a
variety of numerical algorithms [5,49,65,67,72,74,75,76]' including those based
on the finite element approach [35,43].
5. Variable Density or Two Phase Method This method models free surface by
introducing a density discontinuity: regions free of fluid are assumed to be
filled with air [50]. In this case the fluid equations are solved for both fluid
and air. The advantage of this method over the VOF method is that the
boundary conditions at the free surface do not have to be set explicitly. It
36
allows problems with two or more fluids of different densities in contact with
each other to be modelled.
However, for most of the casting problems air can be modelled as a region
with constant pressure without the need to resolve the velocity field in it.
This assumption can be made because of the large density differences between
metal and air, and velocities of air are small compared with the speed of sound
in the air.
Therefore solving the full density equations for metal and air may be a waste
of computer resources. Besides, a higher order numerical scheme is required
for the solution of this equation to preserve sharp interfaces between the fluids
(see Section 2.5).
One way to accelerate the computation is to consider a fictitious air, called
"numerical air" whose properties are such that the linear system issued from
the discretised problem is solved rapidly. Careful attention has to be paid to
interfacial conditions between the fictitious air and the fluid. The basic idea
of the method is to solve a flow model in the air which gives very shallow
pressure gradients. Consequently, the constant pressure condition at the free
surface will be almost satisfied. This procedure is used in a commercial CFD
package SIMULOR developed to model mould filling [77,78]. A shortcoming
of the method is a relative imprecision of the exact location of the free surface
if the mesh is coarse.
6. Lagrangian Method. This method appeared only recently and is still being
developed (not to be confused with the Lagrangian method of tracking dis-
continuities developed by Hirt et al [79]). It is based on the representation
of the flowing or deforming medium as a collection of mutually interacting
parcels of material [80]. Individual parcels move as solid bodies subjected
to the boundary conditions. The main problem in this method is to define
the forces acting between the parcels as well as energy exchange and inter-
action with walls. These definitions depend on the number (or size) of the
parcels used and at present are chosen posteriori making the resulting motion
look more realistic. However, if the problem of arbitrariness is overcome, the
method offers flexibility and control in introducing properties of the materials.
The VOF method is employed in FLOW-3D and is used in the present work for
tracking free surfaces. The transport equation solved for F is
of + ('\1. v) F
at = 0 (2.46)
37
Velocities are obtained from the momentum and continuity equations prior to ad-
vecting the fraction of fluid function.
Once the fluid configuration is defined, pressure boundary conditions, Pi,i' have to be
set at the centres of the cells containing the free surface. In the original SOLA-VOF
code this pressure is obtained by linear interpolation (or extrapolation, depending
on whether the surface cell centre is covered with fluid on not) between the surface
pressure ps and the pressure Pi,i-l in the full fluid cell neighbour in the direction
normal to the free surface as shown in Fig. 2.2a:
where TJ = D..y / d and d is the distance from the centre of cell (i, j - 1) to the free
surface. As pressure Pi,i-l takes part in the iteration procedure, Pi,i also changes
from iteration to iteration according to Eq. (2.47). The following relation must be
satisfied to ensure convergence [81]
1] < 2 (2.48)
which is always true in uniform meshes [60]. But it may not be the case in a non-
uniformly spaced mesh. For 1] > 2, if Pi,i-l changes by D..Pi,i-l during an iteration,
then Pi,i will change by a greater quantity (1 - TJ )D..Pi,i-l' During the next iteration
Pi,i-l will try to compensate for the change in the boundary pressure causing a swing
of Pi,i in the other direction but with an even larger amplitude, and no convergence
will be reached.
To avoid such problems (and to speed up the calculation) Anzai and Niyama [82]
used Ps as the cell-centred value of the pressure in the surface cell, i.e. put Pi,i = Ps'
The results of their calculations for a die mould filling show little difference from
the original method. This is not surprising given that the velocity scale was 10
m/ s and the shape of the free surface is defined by the bulk liquid flow rather than
the surface gravity waves. The Froude number NFr = u 2 /gd for the mesh size
d = 5 mm is around 50. This method, however, is not acceptable for low Froude
numbers characteristic for shrinkage feeding flow, for example, (NFr ~ 0.005 for
velocities of the order of 1 mmj s) because the shape of the free surface in this case
is mainly defined by the differences in pressures in the surface cells.
38
where h is the distance between the free surface and the cell centre. It is impor-
tant that h is negative when the cell centre is outside the fluid. This procedure
ensures continuous and smooth change of pressure in a cell as the free surface passes
through it. Eq. (2.49) is a good approximation in a low Froude number flow and
works equally well in fast flows with the velocity magnitude of the order of 1m/ s [84J.
In this formulation pressures in the surface cells are fixed during the iteration pro-
cedure.
A shortcoming of using the body forces to define the surface cell pressures is clear
in a situation when body forces are absent but a non-zero pressure gradient exists
in the fluid due to a difference in pressures applied at the fluid boundaries [85J. In
that case, according to Eq. (2.49), the surface pressure is assigned to the cell centre
without any adjustment, ignoring the net pressure gradient in the fluid.
By setting surface cell pressures using Eq. (2.49) surface gravity waves are generated,
and care must be taken to avoid possible numerical instability arising from their
propagation (see Section 2.10).
In Section 2.3.1 the upwind (donor-cell) differencing methods was described to calcu-
late the advection the fluid momentum. Unfortunately, it is only first order accurate
and is highly diffusive and thus unsuitable for the free surface advection. A direct
application of this method would result in smearing of the interface, leading to an
non-physical result.
Usmani et al [35J used a continuous, rather than a step-like, fluid fraction function
to avoid the numerical difficulties in advecting the step function. In their approach
the free surface is associated with the value of F( t, x, y) = 0.0, instead of the range
of values between 0.0 and 1.0 as used in the VOF method. F(t, x, y) < 0.0 indicates
the fluid region and F(t, x, y) > 0.0 indicates the voids. Furthermore, at every time
step function F is smoothed in a layer on both sides of the free surface to obtain
a linear profile of F in the direction normal to the free surface. This serves two
purposes: first, a linear distribution of F reduces the numerical diffusion as it is
proportional to the second spatial derivative of F; second, the smoothing helps to
avoid developing undesirable steep gradients of the F-function.
However, this approach is only suitable for tracking an already existing free surface.
Opening new void is beyond its capabilities, in other words the full cells are always
39
assumed to stay full unless the free surface passes though them. The original VO F
method allows one not only to track existing free surfaces, but also to model the
opening of new voids and free surfaces. The latter is crucial in solidification shrinkage
modelling.
In the YOF method the donor - acceptor flux approximation is used to advect free
surface. The essential idea is to use information about F downstream as well as
upstream of a flux boundary to establish an approximate interface shape and then
to use this shape in computing the flux. In each computational cell containing free
surface, the free surface is represented by a plane parallel to one of the coordinate
planes which depends on the distribution of F in the neighbouring cells.
When solving scalar advection Eq. (2.46) in a staggered grid, control volumes are
chosen to coincide with the cells so that the convective fluxes can be easily calculated
at the cell faces.
Consider the amount of F to be fluxed in the x-direction through a cell face during
a time step ~t. The flux of volume crossing this face per unit cross-sectional area is
£ = u~t, where u is the normal velocity at the cell face. The sign of u determines
the donor and acceptor cells, i. e. cells loosing and gaining volume, respectively. The
amount of F fluxed across the cell face is defined as 8F times the cell face area
and
Fe = max{(FDM - FAD)' I £ I- (FDM - FD~XD, O.O}
Here single subscripts denote the acceptor (A) and donor (D) cells. The double
subscript, (AD) refers to either A or D, depending on the orientation of the interface
relative to the direction of flow as explained below. FDM is the maximum of FD and
the F value in the cell upstream of the donor cell.
The min feature in Eq. (2.50) prevents the fluxing of more F from the donor cell
than it has to give, while the max feature accounts for an additional F flux, Fe, if
the amount of void (1.0 - F) to be fluxed exceeds the void amount available. Fig.
2.3 provides an explanation of Eq. (2.50). The donor and acceptor cells are defined
in Fig. 2.3a for fluxing across a vertical cell face. When AD = D, the flux is a
donor cell value,
8F = FD ·£
in which the value in the donor cell is used to define that portion of the cell face
area exposed to fluid (Fig. 2.3b). Satisfying Courant stability criterion (Section
2.10) guarantees that it is not possible to empty the donor cell in this case.
40
When AD = A, the value of F in the acceptor cell is used to define that portion
of the cell face area which F is flowing through. In Fig. 2.3c, all the fluid in the
donor cell is fluxed because everything lying between the dashed line and the cell
boundary moves into the acceptor cell. This is an example of the min test in Eq.
(2.50). In Fig. 2.3d more fluid than the amount FA' £, must be fluxed, so this is an
example of the max test.
Whether the acceptor or donor cell is used depends on the mean surface orienta-
tion. The acceptor cell is used when the surface is advected in a direction normal
to itself, otherwise, the donor cell value is used. The reason for testing on surface
orientation is that an incorrect steepening of surface waves will occur if the acceptor
cell is always used to compute the fluxes. In effect, the acceptor method is nu-
merically unstable because it introduces a negative diffusion of F. Instabilities do
not grow to unbounded values, however, because of the min and max tests used
in the flux definition. On the other hand, when the surface is advected normal to
itself, a steepening that maintains the step-function character of F is exactly what
is required.
Once the flux has been computed by the above method, it is multiplied by the
flux boundary area to get the amount of the fluid to be subtracted from the donor
cell and added to the acceptor cell. In this way the fluid volume defined by F is
conserved. When the advection process is repeated for all cell boundaries in the
mesh, the resulting F values correspond to the time-advanced values satisfying Eq.
(2.46) and still define sharp interfaces.
Youngs [86] developed a 2-D VOF advection model in which the reconstructed free
surface is allowed to have a non-zero slope in a cell, rather than forcing it to be
parallel to a coordinate plane. The Youngs' method is presumably more accurate,
but excessive calculations may be required to define the orientation of the free surface
in every cell when the method is applied in three dimensions.
Referring to Fig. 2.3e, the flux through the right-hand face of the cell can be
computed as [50]
for positive u e • The van Leer approach determines the face values in terms of local
41
t ;".. ,.'-
,
c-
gradients of F:
where
and
if be > 0
if be < 0
To ensure monotonicity, if be X 8w < 0 then (aF / ax}p = 0 and it reverts to upwind
conditions.
The van Leer method has been proven particularly successful in preserving sharp
discontinuities and was applied to advecting free surface in mould filling simula-
tions [50,88].
The use of the VOF method with either the donor-acceptor or the van Leer advection
schemes limits the time step size so that the free surface cannot be moved by more
than one cell in one cycle. This restriction can be inconvenient, especially if it
is significantly more severe than other time step limitations, e.g. viscous or heat
transfer time step stability limits.
Swaminathan and Voller [72,73] suggested an "enthalpy type" model for advecting
free surfaces. The analogy is used with the enthalpy method of describing phase
transformations, according to which a computational cell undergoing the phase
change will remain fixed at the melting/freezing temperature until the associated la-
tent heat has been extracted/supplied (Section 2.8). In this situation a cell changing
phase is unable to transport heat to cells which have not yet changed phase.
In the 'enthalpy type' free surface advection model a new variable, G is introduced
such that
o ~ F < 1 when G = 0
42
o< G ~ 1 when F = 1
where F is the ordinary VOF function. The transport equation is then used in the
form
of
- + (V· v)G = 0 (2.51 )
at
In the context of the analogy with the phase change problem, F takes the role of
the enthalpy and G of the temperature. Eq. (2.51) is equivalent to an isothermal
phase change problem in which the specific heat is zero, i.e. the 'enthalpy' consists
only of the latent heat.
Eq. (2.51) is discretised in a fully implicit, first order form and solved in each
control volume simultaneously. The simultaneous solution procedure is significantly
simplified due to the relationship between F and G.
1. no numerical diffusion and smearing of the free surface is present due to the
introduction of variable G;
2. to take full advantage of the large time step size advection terms in the mo-
mentum equations must also be approximated in an implicit way. This usually
would require an interation method to solve the discretised equations thus in-
troducing additional calculations in the solution procedure.
3. if the time step is so large that the free surface is advanced by several cells,
a problem arises of setting velocities and pressures in the newly filled cells
43
which are required as the initial values for the pressure correction iteration
procedure. As the solution for the Laplace equation for pressure, derived from
the continuity equation (2.21), is defined by the boundary conditions [27],
the first-guess pressure values are unimportant for the converged solution,
if correct boundary conditions are given (though the number of iterations
necessary to achieve the convergence is affected by the initial value). However,
it is unclear how to calculate the velocities.
4. the condition for a computational cell to be full before passing the fluid to an
empty neighbour is a good approximation if the free surface is advected in the
direction normal to the free surface. In a case when the fluid flow is parallel
to the free surface this approximation will lead to an incorrect position of the
free surface.
At present the VOF approach, in conjunction with the donor-acceptor or van Leer
advection method, is the most powerful and general tool to describe transient free
surfaces.
Eq. (2.23) is to be solved numerically for the fluid enthalpy. A numerical analogue
of it for a control volume V, chosen to coincide with a cell, is
(2.52)
Here a first order forward differencing is used for the temporal derivative. f::.H:d is
the advective flux of the fluid enthalpy across the cell faces, estimated at time level n,
i. e. explicitly. It is closely related to the advection of the fraction of fluid function,
F (Section 2.5). Once the amount 8F of the VOF function that is advected through
a cell face per unit area has been calculated, the enthalpy advection is estimated by
8H = Hv8F (2.53)
where HD is the enthalpy of the upstream (donor) cell per unit fluid volume. The
total flux of the enthalpy into the cell is then equal to the sum of the fluxes through
each cell face times the corresponding face areas:
f::.H:d = L Ak t5H;:
face.
In that way the fluid enthalpy is conserved as the fluid moves through the mesh.
44
The discretised Fourier conduction terms in Eq. (2.52) are represented by QT which
is the sum of fluxes at control volume faces. In two-dimensions, a uniform grid and
constant heat conduction coefficient k, the central differencing for the Laplacian
term in Eq. (2.23) at cell (i,j) gives
QT = k (Ti+l,j - 2Ti,j + Ti-1,j + Ti,j+l - 2Ti,j + Ti,j-l) (2.54)
~X2 ~y2
Eq. (2.54) gives a second order accurate approximation for thermal diffusion. The
advantage of this approximation is that it only involves the immediate cell neigh-
bours, in other words, it uses the same expression for all internal fluid or mould cells.
For the boundary cells, either the temperature or the heat flux must be specified as
the boundary condition.
A fully explicit formulation for the conductive heat fluxes, i.c. with all temperatures
in Eq. (2.54) evaluated at time level tn, is employed and a time step stability limit
is introduced. For a one-dimensional flow it is (see also Section 2.10)
~X2
~t< -
20:
A fully implicit approximation of the conduction term would imply all temperature
on right-hand side of Eq. (2.54) to be taken at the new time level tn+!' This will
remove the time step restriction but the system of the discretised equations becomes
coupled and iterations are required at each time step to obtain the solution [52]. Note
that at each iteration Eq. (2.17) will have to be solved for the temperature. The
latter can be complicated by a non-linear dependence of the fraction of solid function
fs on the temperature (Section 2.8) and the whole solution procedure may prove to
be ineffective. To avoid this complication, the energy equation can be written solely
in terms of temperature and the phase transformation taken into account using the
"effective specific heat" [5] (see also Section 2.8).
A semi-implicit method was suggested by DuFort and Frankel [89]. The essence
of the approach is to use only the central node temperature in expression (2.54),
Tj,j, at the new time level t n+!; all the other temperatures are taken at time tn.
In this case there is no stability limit for the time step and the solution matrix for
the temperatures is purely diagonal. The cost of employing this method can be
demonstrated on a one-dimensional heat conduction problem without phase change
and fluid flow. Eq. (2.52) then can be written for a node i as
T'!'+!
. = T'!'. + ~ t a Ti'tl - 2Tr+!
~X2
+ Tt'-l (2.55)
45
or
where
!:::.t
!:::.t'
and
2
· uAt' =!:::.x
11m --
~t-+cx> 2a
Thus Eq. (2.55) is equivalent to a fully explicit expression, with an erroneously
small time step. It is erroneous because the user may think that the result of n
time steps applies to the time n!:::.t, when it actually should apply to the time n6it'.
The errors result from the absence of energy conservation property in this method
since temperatures at different time levels are used to evaluate the heat flux in the
right-hand side of Eq. (2.55).
A fully explicit first order numerical approximation to the thermal diffusion term
(2.54) will be used in the present work due to its simplicity with regard to the
solution matrix inversion and to the development of the shrinkage model.
The importance of the interfacial heat transfer strongly influences the choice of the
numerical technique to solve Eqs. (2.21 )-(2.23). Often the solution accuracy for
fluid flow is sacrificed in favour of the heat transfer.
The use of unstructured meshes in the finite element approach allows one to solve
successfully both problems of the geometry and heat transfer modelling. The latter
46
is achieved by placing the mesh nodes on the interface (boundary fitted mesh) and
using temperatures, the metal and the mould, at each of those nodes so that the
temperature difference in Eq. (2.56) can be accurately estimated [90]. Besides, a
formation of an air gap can be described if the nodes are fixed to the solid material
as it moves off the mould wall [50]. The effect of the gap size is then taken into
account by adjusting the heat transfer coefficient h. An experimental procedure
for measuring the gap size and the heat transfer coefficient has been given by Hou
and Pehlke [24]. Evans et al presented an experimental techniques of measuring the
interfacial heat flux [91].
The use of structured mesh in the finite volume approach introduces the following
problem. Generally the metal/mould interface is represented by a stepwise surface
because each computational cell can be either completely open to flow or completely
blocked by the mould, and the approximate interface passes along the cell faces. A
finer mesh is not a solution because (a) there is always a limit to how fine a mesh
one can use and (b) refining the mesh does not make a better approximation for
the interface area. Figs. 2.4a, b show a finite difference approximation to a straight
line by two meshes, a coarse and a fine. Although the finer mesh representation
seems to be closer to the straight line the length of the discretised line in each
case is the same and depends only on the line inclination. As a consequence, finite
difference cell nodes do not lie on the interface, and the cell-centred temperatures
cannot represent the interfacial temperatures.
Another way of describing solid walls is called Fractional Area and VOlume Ratio
(FAVOR) technique [88,123]. In this method a cell volume and faces are allowed to
be partially blocked (Fig. 2.4c) Four additional arrays are required then for each cell
(three in two dimensions): the fractional open cell volume and fractional open areas
for the right, back and top cell faces as shown in Figure 2.4d for a two-dimensional
cell. In this case an interface containing cell has two temperatures: metal and mould.
FAVOR method allows much smaller numbers of cells to be used to represent curved
surfaces and gives better wall boundary condition approximations not only for the
heat transfer, but for the viscous stresses as well. The price paid for the improved
accuracy are the additional storage requirement and possibly a more severe time
step stability limit (Section 2.10 and Chapter 3)
It is shown in Section 4.1 how, using the FAVOR technique, temperature nodes
can be efficiently introduced at the interface on both sides without any significant
increase in the computation. The actual interfacial area in each cell can be estimated
at the preprocessing stage and used as dA in Eq. (2.56) [92].
47
Structured finite difference and finite volume meshes are fixed in space and cannot be
used to model efficiently the movement of the solid material and the gap formation
due to thermal stresses. An attempt to model the gap formation was made by
Huang et al [94]. Their model is based on the assumption that the magnitude of an
interfacial gap size primarily depends on the thermal contraction of the casting as
it solidifies. For a point P at the metal/mould interface, its movement relative to
the geometrical centre of the casting and the movement of the casting as a whole
due to gravity is estimated. As a result the point can move away from the interface
creating a gap of size .6./. The heat transfer coefficient is adjusted according to the
formula
h = kg(T)
B+Cf).[
where kg(T) is the gas thermal conductivity at temperature T, and Band Care
fitted from experimental data. This model is simple and efficient and has been
implemented in two-dimensions.
There have been numerous attempts to simplify and speed up the solution algorithm
for the energy equation. Dantzig and collaborators [95,96] developed a boundary cur-
vature method based on the fact that, generally, information on mould temperature
is not required. Hence, the mould can be excluded completely from the calculation
given that the correct heat fluxes at the metal/mould interface are known. The
latter is calculated by correlating the actual interface shape with simple geometries
for which an analytical solution exists and applying this solution to calculate the
heat fluxes and temperatures in the mould without discretising the mould domain
and solving the energy equation in it.
The boundary curvature method is especially attractive for modelling sand castings,
because large differences in metal and sand thermal diffusivities require substantial
computational effort to accurately calculate heat flux in the sand. A much finer
mesh would have to be used in the sand than in the casting itself.
However, the boundary curvature method is limited by its lack of generality in de-
scri bing complex geometries and by lack of physical transparency in the formulation
of the method. The method has been implemented in SPIDER, a finite element
code [97].
An analytical solution can be used to describe the penetration of heat into the
mould for a short time after the contact between the metal and mould wall, when
the thermal gradients in the sand are very steep and cannot be resolved by the
finite volume mesh [93]. Stoehr and Wang [66] used an analytical expression for the
48
interfacial heat flux per unit area when modelling thin wall sand casting:
q =
/kPC
-;t (T metal - Tsand )
in which k, P and C refer to the relevant properties of the sand and time is measured
from the arrival of metal to that point. This is an attractive alternative to mesh
refinement at the interface. However, care must be taken in using the analytical
functions since physical conditions at the interface do not generally correspond to
those assumed by the analytical solutuion. Besides, a special algorithm has to be
derived to make a smooth transition from the analytical solution to the numerical
one at later times of simulation.
Two-domain methods have been also developed in which the casting and mould are
meshed independently and the two solutions of Eq. (2.52) are found in each of the
domains separately. The interface boundary conditions are then used to match the
metal and mould temperatures. This technique allows one to take advantage of the
fact that the energy equation for the mould does not contain advection and phase
transformation terms:
aT = _1 V( kVT) (2.57)
at pC
Different time step sizes can be used in the two domains and linear temporal inter-
polation is then employed to bring the solutions together [98J.
A similar approach is used by Chen and Tsai [99]. At a given time step Eq. (2.52) is
solved together with the interfacial condition, Eq. (2.56), to obtain only the casting
temperature distribution. The mould temperature is taken from the previous time
step. Then in the next time step Eq. (2.57) is solved with the known temperature
from the previous time step solution for metal. At each time step fully implicit
formulation is used for the corresponding temperatures, and the solution is found
iteratively. The method allowed optimum time steps and iteration parameters to be
used for each material.
The general shortcoming of the two-domain method with separate time stepping is
the absence of the conservation property: the net energy change in the mould may
not be equal to the loss of the energy in the metal. It can be shown that for the
Chen and Tsai method this results in a similar effect to that of the DuFort-Frankel
method described earlier in this Section. Neglecting advection, phase transformation
and assuming for simplicity a lump temperature model for both metal and mould,
the Chen-Tsai model becomes
T n+l Tn-l
1 - 1 = -Ql (Tr+I - T;) (2.58)
2~t
49
r.2n +
2
-
r.2n = 0:' (T n+1 _ r.n+2)
2 I:lt 2 1 2
(2.59)
The same time step size was chosen for domains 1 (metal) and 2 (mould) with
0:'1 = hi PI C1 and 0:'2 = hi P2 C 2 . The metal and mould temperatures are calculated
at times tn+l and tn+2 respectively. The loss of the conservation property arises
from two factors. Firstly, the heat fluxes out of the metal and into the mould are
evaluated at different times. Secondly, the two temperatures used to evaluate these
fluxes in both equations are taken from different times. If Tt+1 on the right-hand
side of Eq. (2.58) is substituted by an estimate for the metal temperature at time
(2.60)
then the heat flux and time tn will be estimated in a more conserving manner. Eqs.
(2.58), (2.59) can be rewritten after some algebra as:
(2.61 )
(2.62)
where
~(T.n+2
2 2
+ T.n)
2
and
~t
I:lt' = - - - -
1 + O:'l~t'
Similarly to the DuFort-Frankel method, the actual time steps, ~t' and I:lt", at which
the solutions for the metal and mould advance are smaller than ~t. Furthermore,
the times, to which those solutions actually correspond, are not the same for metal
and mould and the difference grows as the calculation proceeds because ~t' =j:. ~t"
due to the differences in thermal diffusivities. The error in the Chen-Tsai method
(as well as the DuFort-Frankel) is a truncation error, i.e. it is proportional to a
power of the time step size.
A two-domain method is also used in the present work, but energy equations for
the metal and mould are coupled at the same time level and solved simultaneously.
The latter is simplified by the use of a fully explicit formulation of the numerical
analogue equations.
50
2.8 Modelling the Latent Heat Release
There are several widely used numerical methods to take into account the latent
heat in the energy equation.
Enthalpy Method. The enthalpy method of accounting for the latent heat release
uses the Energy Equation in the form of Eq. (2.23). Once a numerical solution for
H at a certain time level is found, using one of the methods described in Sections
2.6 and 2.7, the current released latent heat and metal temperature are found from
Eq. (2.17).
The advantage of the enthalpy method (Eq. (2.23)) is that it is not possible for a
node to drop through the freezing range without the latent heat being accounted
for, since it is incorporated at every stage of the computation. In other words the
enthalpy method is energy conserving.
An artificial freezing range is usually introduced to prevent the cooling curves from
exhibiting false plateaus, However, this will cause a diffusion of the solid/liquid
interface. In general, the smaller the artificial mushy zone, the better the solution
accuracy [99].
Another enthalpy-based method, called enthalpy diffusion model, has been suggested
by Mundin and Fortes [104]. Here the energy equation is written in the form
8H
8i=V.(aVH) (2.63)
where a varies from k/ pC to zero (in the moving front). Eq. (2.67) gives an energy
conserving numerical solution. The temperature oscillations, typical in the enthalpy
method, are absent but a stretching of the phase changing region is introduced since
enthalpy is continuous across the liquid/solid interface.
51
for alloys with large freezing range T/ - Ts [105]. As was mentioned in Section 2.1.2
it is not among the objectives of this work to include solidification kinetics. We will
only note that the enthalpy method assumes the enthalpy of a cell to be a unique
function of the cell temperature. This would introduce a difficulty in modelling un-
dercooling, characteristic to the kinetics models, as it implies at least three enthalpy
values for one temperature when the cell is, first, undercooled, then the temperature
recovers due to solidification and, finally, decreases as the solidification progresses.
Several authors used experimentally measured evolution of the metal enthalpy versus
temperature to incorporate the latent heat release [100,106]. If a metal sample, V,
is cooled at a sufficiently slow rate to assume it to be at a uniform temperature T,
then the change of the total metal enthropy is
dH =V
& p
(edT _L dis )
& & =~(T.
A 00
-T)
where h, A, and Too are the metal/mould heat transfer coefficient, the interface area
and the ambient temperature respectively. Since
dis = dis dT
dt dT dt
the total latent heat released to a certain time can be expressed as a function of the
temperature and cooling rate
dis ( T - Too) hA
L dt = e 1 + e dT / dt ' e= -
pVc
and all that has to be done is to accurately measure T( t).
The experimental method is perhaps the most accurate and flexible way to describe
the latent heat release, especially for multi-component alloys. Most solidification
modelling packages have an input option to include a table of enthalpy-versus-
temperature values.
Modified Specific Heat Method. This method involves artificially raising the specific
heat within the freezing range to account for the release of latent heat. The simplest
method employs a constant higher value of specific heat within the solidification
range [66,93]:
L
eef f = e + T/ - Ts ' (2.64)
where T/ and Ts are liquidus and solidus temperatures, respectively. This is only
suitable for narrow freezing range alloys. A more accurate expression is given by [100]
eeff = e - L. dT
di.
52
where dfs/dT can be defined by using an appropriate solidification model.
If the liquidus and solidus lines are assumed to be straight, together with a complete
mixing of the solute in both liquid and solid phases, then the lever rule applies [7]
(see also Eq. (2.18) in Section 2.1.2):
(2.65)
where J( is the partition coefficient and Tm is the melting point ofthe solvent metal.
If the conditions of the solidification are such that there is no diffusion of the solute
in solid and complete diffusion in the liquid phase, then the fraction of solid is
determined by Scheil equation:
fs = 1, (2.66)
A problem with the modified specific heat method is that it is possible for a node to
'jump' over the freezing range in a single time step, thereby missing out the latent
heat evolution. This is illustrated in Fig. 2.5a, where the effective heat capacity is
calculated using Eq. (2.63). A post-iterative check is necessary to ensure that the
latent heat has been accounted for. If this not the case then the nodal temperature
must be readjusted accordingly [101].
It may be difficult to justify the application of these models in the form of Eqs.
(2.64) and (2.65) in numerical modelling. For example, ScheH equation, Eq. (2.65),
assumes an infinite diffusivity of the solute in the liquid, no diffusion in the solid
phase of the material and that the total amount of the solute in the casting, Me, is
constant
Me = const (2.67)
It implies that the liquid phase in the casting always has a uniform concentration,
and the fraction of solid function in Eq. (2.65) is the net mass solid to liquid ratio
in the casting.
In the control volume approach Eq. (2.65) is applied to every control volume, with
fs being the fraction of solid in the control volume. This effectively means that Eq.
(2.66), instead of being reasonably applied to the whole casting volume, is used in
53
each control volume, thus, becoming much stronger and, generally, inaccurate and
making the numerical solution mesh-dependent. In reality the amount of solute in
each control volume may change due to fluid advection and diffusion.
Similar problems can be seen in the application of the lever rule. Despite that, Eqs.
(2.64) and (2.65) and their variations have produced useful predictions.
Temperature Recovery Method. In the case of pure metals the latent heat is divided
by the specific heat to give a virtual temperature, ~Tv [102]. This may be regarded
as the temperature change over which an amount of specific heat equal to the latent
heat would be evolved. Using one of the above solidification models the variation of
enthalpy within the solidification range can be obtained. The modified specific heat
capacity at a particular temperature can then be determined from the slope of the
enthalpy versus temperature curve.
The liberation of latent heat is accounted for by holding the temperature of a node at
the freezing point until a number of excess degrees equal to the virtual temperature
have been accumulated after which the temperature is allowed to fall.
In alloy solidification the latent heat can be accounted for in a stepped manner
according to the percentage of solid formed at different temperatures in the freezing
range [103]. The disadvantage of the temperature recovery method is that it requires
that additional post-iterative calculations be made.
In the present work the enthalpy method will be employed to model solidification,
using Eqs. (2.17) and (2.23).
Shrinkage induced porosity results when there is an inadequate supply of liquid metal
to counter the volumetric shrinkage on solidification exhibited by most metals.
Gas induced porosity may result from the evolution of gas during solidification by
either one of the two mechanisms:
54
• by the formation of a compound gas during solidification (e.g. CO in steel).
A further possibility is that air bubbles may be introduced into the liquid metal
at the filling stage. If these bubbles are unable to escape then gas porosity will
result [2J.
Usually in castings the gas content is kept as low as possible by using either inert
gas degassing techniques, vacuum treatment or by the addition of an element which
form a solid compound with the dissolved gas. Exceptions to this are steels in which
CO is allowed to form to compensate for solidification shrinkage. A similar technique
is sometimes used in die casting when dispersed porosity is considered less harmful
than a concentrated shrinkage cavity.
If during the solidification of a pure metal part of the casting freezes over, leaving a
substantial pocket of liquid metal isolated from the feeding top, then gross internal,
or secondary porosity results. In the case of alloys solidifying with a dendritic front,
it is not necessary for the metal to freeze over to 100% solid to stop the feeding flow.
This is because when the fraction of solid in the mushy zone is high, the remaining
interdentritic liquid regions do not make a continuous passage connecting the feeder
and the rest of the liquid casting.
Microporosity generally occurs randomly between the dendrite arms of alloys so-
lidifying in a pasty manner. Uram et al [107] et alobserved that the amount of
microporosity between dendrites in a columnar region is considerably less than in
the equiaxed region. The reason for this is that the interdendritic liquid regions in
the columnar region are orientated along the crystal axis and they are more likely
to connect the shrinking material in the mushy zone with the bulk liquid metal.
55
In the formation of the internal porosity the contributions of dissolved gas and the
pressure drop across the mushy zone due to incompressibility of metal are additive.
Assuming no barrier to pore nucleation the condition for the formation of a pore of
radius r is [7]
20-
P :S Pg - -
r
where p is the local metal pressure and Pg is the equilibrium gas pressure.
At high gas contents pores can form early during solidification and remain spherical
in shape. Gas bubbles forming at later stages will become trapped in the mushy
region, their shape distorted to the shape of the interdendritic space. The total void
volume is the sum of shrinkage and the evolved gas volume.
The present research is primarily concerned with shrinkage porosity, though a uni-
formly distributed gas porosity could be taken into account by decreasing accord-
ingly the solid phase density (Chapter 7). The exclusion of detailed gas micro-
porosity modelling is mainly due to the complexity of the physical phenomena.
Experimental results indicate that pore sizes in alloys varies markedly with cooling
rate, gas content and, to a lesser extent, with grain refining [108].
56
porosity is calculated using the Solid Fraction Gradient Method, which assumes
that pores are generated by the solidification shrinkage of a large quantity of molten
metal remaining at the time of the fs.cr isosurface disappearance.
- molten metal flows downwards due to the gravity, z. e. effects of pressure and
viscosity are neglected;
- the flow speed is much greater than that of the solidification front; hence feeding
can be assumed to occur instantaneously;
Whilst this method is not able to predict the formation of dispersed porosity in
regions of low temperature gradients it can predict the actual shape and position
of gross shrinkage cavities. This is not the case with methods which ignore any
consideration of fluid flow and simply attempt to locate thermal centres within the
casting.
- no heat flux due to the metal flow is calculated; this can be a combination of the
residual circulation from the filling stage, the feeding flow and thermal convection.
The latter may be increasingly important for low conductivity metals, such as steels.
- the use of the Solid Fraction Gradient Method requires a high accuracy numerical
heat conduction/transfer model, and may be sensitive to the method of calculating
latent heat.
- it is not clear from the description of the method [111] whether the computational
elements are allowed to be partially filled during the calculation and how is it taken
into account in terms of the heat flow.
A similar method was used by Nagasaka et al [112]. Here the authors use the
solid fraction gradient, 'V fa, as the key parameter that controls the driving force for
feeding. The movement of liquid was considered only in the direction of the gradient,
57
i.e. normal to the solidification front. In addition a critical value of V' f, is applied
as a criterion for the feeding. At each time step all possible feeding paths (consisting
of cells with fs < fs,er) to a solidifying cell are determined and the solid fraction
gradient is calculated along each of them. The cell then can only be fed if at least
in one feeding path V' fs does not fall below the critical value. The total shrinkage
volume generates a macro-cavity in the elements with the minimum pressure head
thus taking into account gravity.
This method also requires a high degree of numerical accuracy in evaluating the
fraction of solid gradients. Furthermore, the determination of all possible feeding
paths for each solidifying cell may be inefficient in 3-D problems.
The agreement with the experiments obtained by Imafuku and Chijiiwa and Na-
gasaka et ai, is good, although the accuracy is limited by the size of the individual
elements. The choice of the critical solid fraction value is rather arbitrary. A value
of /6,cr = 0.9 has been suggested by Davies [113] and Spittle and Brown [114].
The method of Niyama et al [115] is based on the experimental observation that cen-
treline porosity results if the temperature gradient, G, along the centreline is below
a certain value. It was confirmed experimentally that this technique can success-
fully predict centreline shrinkage in 100 mm diameter steel casting. A temperature
gradient of 0.2 DC jmm was determined as the critical temperature gradient required
to avoid centreline shrinkage, but its value was shown to vary with the section
thickness.
The same authors have found that if the criterion G j v'Ii is used, where R is the
cooling rate, then the same critical value of 0.8 DCl/2s1/2mm-l can be applied to all
castings of similar compositions regardless of the section size. Although this criterion
is essentially empirical, the authors showed that it has a theoretical basis [116].
However, it is unclear how to apply this criterion to general casting configurations
since Gjv'R is a dimensional parameter.
58
This method was applied to predict the centreline porosity in a cylindrical steel
casting, for which VTS cr =3 mm/ s was chosen.
The authors claimed that their method is more general than the G /.JR method
because the value of VTS cr does not depend on the general cooling rate conditions
and casting geometry. However, VTS cr depends on the freezing range on a particular
alloy and ideally it should be determined experimentally. Furthermore, the VTS
method cannot be used to predict macroshrinkage.
An advantage of criterion methods is their simplicity and ease of use. They can be
easily added to an appropriate computer code.
The main disadvantage of the criterion methods is that they require experimental
verification of one or more of the critical parameters, making these models geom-
etry and materials dependent. Criterion methods also do not model the dynamic
interaction between developing macro- and micro-porosity regions with the rest of
the casting. An internal cavity would affect the fluid and heat flow in the casting.
The liquid fraction II is found from a finite volume solution of the macroscopic
energy equation. The velocity is found from Darcy's law
1
v= --(\1p - pg) (2.69)
]{ pil
where
1
]{
is the permiability coefficient and d is the dendrite cell size. Gas pressure in the
porosity region is calculated using equation
20-
p = Pg -
r
(2.70)
59
The diameter of the porosity forming first is assumed to be the same as the root
diameter of the dendrite cells.
The conservation equation for gas content, Q, in this case hydrogen in an aluminium
alloy, is
(2.71 )
where Qo is the initial gas content and Qs, QI are the gas contents in the solid and
liquid phases, respectively. The last term in Eq. (2.71) describes the amount of the
gas in pores.
When the volume element is in the mushy zone, the amount of porosity is calculated
from metal pressure p using Eqs. (2.68) and (2.69). The gas pressure is then
calculated from Eq. (2.70) and the new amount of porosity from Eq. (2.71).
In this model feeding flow velocities are estimated at each time step, though the
advection terms in the momentum and energy equations are neglected. The flow is
assumed to be dominated by the friction and gravity forces. However, this model
is more suitable for predictions of micro-porosity due to gas evolution than for
macroshrinkage modelling.
Stoehr and Wang [66] included a volumetric source term in the continuity equation
in their coupled 2-D heat transfer and fluid flow model. If the solid/liquid mixture
density in each control volume is estimated as the volume average5
(2.73)
then the Continuity Equation, Eq. (2.21), can be rewritten, for 0 < Is < 1, as
·
dwv = PI - ps
P
(aatls + 't'7f)s
V· v (2.74)
The source term in Eq. (2.74), evaluated from the solution of the heat transfer
equation, would generate the flow necessary to accommodate the volumetric change
and can describe both shrinkage and expansion (e.g. due to remelting) of the metal.
5 f. in Eq. (2.73) is the volume fraction of solid. If the mass fraction of solid was used, as
everywhere else in the text, then Eq. (2.73) would transform to
p = P.P' (2.72)
f,Pl+(l-f.)p,
60
However, the authors do not present any numerical technique to find the solution
of the coupled continuity and momentum equations, nor do they claim to have used
the model to simulate the formation of shrinkage cavities.
Evans et al [120] developed a shrinkage model based on heat flow only for a 2D
axisymmetric cylindrical geometry. They employed a value of the critical solid
fraction Is.er = 0.8 which was higher than that used by Imafuku and Chijiiwa.
Micro-porosity is estimated after the minimum solid fraction in the casting reaches
the critical value. This value of Is,er is used to estimate feeding in the vertical
direction. In the horizontal direction it is assumed that feeding occurs until the
metal is fully solid, i.e. effectively ISler = 1.0 for the horizontal flow. Cells are not
allowed to be partially filled. Instead, the liquid level is reduced by a full cell layer
every time when sufficient volumetric shrinkage has been accumulated. This requires
a fine mesh to reduce the resulting truncation error.
BEX BEY)
f7 xy = /'i, ( By + Bx
where Ex, Ey are the displacements, /'i, = E/(1 - J.L;), 11(1 is the Poissons ratio, E the
Youngs modulus and /1 the coefficient of linear thermal expansion.
The enthalpy equation is solved in the whole domain while the stress equations are
only solved in the solid phase. The boundary conditions for the latter are such that
the solid metal surface, including the liquid/solid interface, moves along the direction
of the interface normal and towards the solid phase due to thermal contraction of
the latter.
The authors use control volume approach with unstructured mesh in which the
shape of the individual control volumes is initially rectangular, but then it changes
to follow deformation of the material. The solution proceeds as follows:
2. estimate 6.T in each control volume and solve the stress equations;
61
3. using the newly obtained displacements estimate the size of the air gap and
the movement and deformation of each control volume;
4. calculate porosity within the liquid material by estimating the stretching of
the liquid control volumes due to the movement of the solid phase. The stress
boundary conditions at the liquid/solid interface ensure that the solid phase
pulls the liquid phase apart, forcing it to create voids. The void regions are
evenly distributed in the liquid phase at each time step, therefore, gravity is
not taken into account.
For the metal/mould heat transfer calculation coincident nodes are used as described
in Section 2.7.
The assumption inherent in step 4 above may be unrealistic. Cavities inside a casting
appear because of the volumetric changes during the phase transformation occurring
at the liquid/solid interface, and not because of the linear contraction of the solid
phase upon cooling. Fig. 2.7 shows the predicted cavity in a 2-D square casting. In
fact, given only thermal solid contraction, there should not be an internal cavity.
The cavity in Fig. 2.7 results from the way the boundary conditions are set at the
liquid/solid interface. The actual movement of the interface due to the thermal
contraction of the solid phase depends on its position relative to the geometrical
centre of the casting and the mould walls but it may not always move towards the
solid phase along the normal to the boundary.
There are several restrictions on time-step size must be observed to avoid numerical
instabili ties .
• for the explicit upwind differencing scheme, described in Section 2.3.1, the
time step size has to satisfy the Courant-Friedrichs-Lewy (CFL) condition
~t < ,
mm''.J,'k
(~Xi ~Yj ~Zk)
-- -- --
u"',J,k " v'·
(2.75)
',J,k W··
',J,k
where u, v and w are velocity components. Physically this condition means
that fluid must not flow across more than one computational cell per time
step.
62
• if diffusion terms are approximated explicitly, then the time step is limited
further:
1
f:lt < (2.76)
OOk.
2.max ',],Ok[a',],
o
(..,.l.,..+..,.l.,..+ 1 )]
~xt ~y; ~z~
where a = J-l/ P for the viscous terms in the momentum equation (Section 2.3.1)
and a = k/ pC for the heat diffusion terms in the energy equation (Section
2.6)6.
The restriction physically means that no quantity should diffuse more than
one cell in one time step .
• free surfaces introduce another type of stability condition associated with the
propagation of surface waves and the way the free surface cell pressure is
determined in the calculation (Section 2.4). If a body force G is applied to
the fluid in a direction normal to the free surface, there may be surface waves
with speeds of order of y'G h, where h is the depth of fluid or length of the
wave. The value of h is defined by the pressure interpolation procedure in the
surface cell, given by Eq. (2.49), and is equal to the cell size in the direction
normal to the free surface.
The actual condition is that surface waves should not propagate more than
one cell in one time step. For example, if z is the normal direction to the free
surface, then
(2.77)
Similar limits must be imposed in the x and y directions for each cell containing
a free surface.
The range of the stable time step sizes is given by the condition, Eqs. (2.75}-(2.77),
with the smallest right-hand side.
These time step size stability limits result from using the von Neumann stability
analysis of the linearised numerical analogue equations [27]. In this method the
influence of boundaries is excluded and the solution behaviour is investigated by
analysing the evolution of a single Fourier component of the numerical solution.
6In fact,
a',i,1: =
is a more correct expression. As dH/dT = C + Ldf./dT, the time step restriction becomes less
restrictive for partially solid cells and, therefore, the latent heat has a stabilising role in an explicit
numerical algorithm.
63
The mathematical foundations for the analysis of convergence and stability of nu-
merical schemes are well-developed only for linear systems. The results from the
linear theory are used as guidelines to nonlinear problems, the justification depend-
ing on numerical experiments. In practice, the stability limits given by right-hand
sides of Eqs. (2.75)-(2.77) are decreased further by multiplying them by so called
safety factors which usually do not exceed 0.5.
The main interest of the present work is the development of a hydrodynamic model
to simulate coupled fluid flow and shrinkage during solidification. Shrinkage defects
are among the main causes of poor quality castings and a computer model capable
of accurate predictions of theses defects.
Most of the existing numerical shrinkage models ignore fluid flow, including residual
circulation, thermal convection and shrinkage induced flow. Those models which do
include fluid flow have not been used to produce a comprehensive analysis of the
influence of the fluid flow on the shrinkage defect formation.
The shrinkage model based on full fluid flow equations and which is described in
Chapter 5 includes numerical models for:
- fluid flow;
- free surface;
The work is based on a commercially available numerical model of fluid and heat
flow with free surface motion described in Chapter 3. The metal/mould interfacial
heat transfer and solidification models are enhanced to improve the accuracy and
include physical models as described in Chapter 4.
An additional shrinkage model will be developed which is based only on heat flow
equations thus ingnoring the fluid flow aspects of the solidification and shrinkage
process. This model is used for comparisons of its simulation results with the results
of the full shrinkage model. It will be shown that in some situations the results of
the two models are close inducating that fluid flow may be not important in these
64
cases.
65
V ..
I,J
~
,,.
p ..
u.1- l,J' I,J u ..
(a) ...
~ • ~,
4'
I,J
~
, ..
V ••
I,j-
1
V ••
I,j
u.1- l,J' p .. u ..
I,J I,J
(b)
•
vI,j-
.. 1
(c)
Uj-l,j UL
.. I
I
.:
Ip.J+ l',j
Uj+lJ
Vi,j-l vB vi+l,j-l
Figure 2.1. (a) the staggered mesh arrangement for the primitive
variables; (b) the control volume for the continuity
equation discretisation7 (c) the control volume for
the x-component of the momentum equation discretisa-
tion in the upwind differencing method.
(8)
. .
.... :::::--~
(b)
Figure 2.2. setting the pressure in surface cell (i, j): (a) using
linear interpolation between Pi,j and Psi (b) assuming
hydrostatic distribution between the free surface and
the cell centre due to the body forces Gni h is negative
if the cell centre is outside the fluid.
u·Dt
1-
(a) donor _o+;--'"~acceptor
(b)
eN
,..,
n
Ue
(c) eW (
w • p
--(.~
e
e E
.S
(a) (b)
VF· .
J,J
AT-I,J.
(c) (d)
Figure 2.4. Line and surface disctretisation techniques: (a) and (b) show
the conventional finite difference technique transforming a
line (dotted line) into a stepwise line (heavy solid line);
in both cases the discretised line length is the same. (c) The
discretisation of a circle using the conventional technique
(light shaded area) and the FAVOR method (dashed line). (d)
Variables introduced by the FAVOR method; in this case
VFi.j=O.35, ARi-l,j=l.O, ARi,j=O.O, ATi,j-l=O.2 and ATi.j=O.5.
c
- •
(a)
------------;- - - - -.-
-
T
1000
(b)
900
800
10 20 30
t, seconds
Figure 2.S. (a) The modified specific heat method (equation (72».
If the cell temperature 'jumps' from Tn-l to Tn in one
time step, then the latent heat is not accounted for.
(b) The 'heat wave' effect due to the use of the enthalpy
method is shown on a simulated cell cooling curve. The
distance between peaks corresponds to the solidification
front passing through a cell. The solidification tempera-
ture is 933 0 K.
r
••
•
Interdendllhc feeding
04 ~ <r-----04 ~
Mass (elastic) (plastic)
Icedln g
.. Soltd feeding
..
(macro) (micro)
: .99
l .gs
191
1-'6
1'$
I
L _ __
This chapter contains a description of the parts of the fluid flow model employed in
the proprietary CFD code FLOW-3D that are relevant to the present work [83].
3.1 Introduction
Extensive tests have been performed on the isothermal fluid flow model. The tests
were supported by experimental simulations using water. It was found that FLOW-
3D was capable of modelling successfully fast moving flows typical for many casting
configurations [84,122].
The set of the hydrodynamical equations constituting the mathematical model for
casting simulation consists of continuity Eq. (2.21) for an incompressible fluid,
Navier-Stokes Eq. (2.22) for linear viscous fluid and the energy equation written in
66
terms of enthalpy to model heat flow and phase transformation, Eqs. (2.17) and
(2.23).
• undercooling and nucleation effects are also neglected assuming that solidifi-
cation is solely controlled by the removal of the latent heat. In other words,
dG = 0 (Eq. (2.19)) during phase transformation and liquid at the solidifi-
cation front is always at the ideal melting point. This assumption simplifies
greatly the mathematical modelling of solidification and is acceptable for sim-
ulations of macro-shrinkage effects in both pure metals and alloys.
• Eq. (2.20) for the solute evolution in a melt is not used in the present work.
Instead, simple models like lever rule or Scheil equation (Eqs. (2.64), (2.65))
will be used to determine fraction of solid function. This assumption is par-
tially justified since the choice of latent heat release model has little effect on
the metal temperature and freezing time for alloys with narrow freezing range,
e.g. low carbon steels [105].
Eqs. (2.17), (2.21)-(2.23) are solved with boundary conditions given by Eqs. (2.24)
and (2.25). In many cases of fluid flow during casting surface tension and viscous
forces in Eq. (2.26) can be neglected. This is due to high Weber numbers during
the filling stage (the ratio of the dynamic pressure to the surface tension force)
pRv 2
We = --~1(J
where R is the characteristic free surface curvature radius, and practically zero free
surface curvature at later stages of solidification, i. e. high Bond numbers
Bo = gpR ~1
(J
which is the ratio of the gravitational force to the surface tension force l .
1 In aluninium based alloys (T may be large because of the formation of an oxide film on the
metal surface [2].
67
Viscous stresses at the free surface are small compared to the fluid inertia, z. e.
Reynolds number is large
Lvp
Re = --» 1,
J1
and to the gravitational force, i.e. Galileo number is large
gp2L
Ga = --»
J12
1
After neglecting the viscous and surface tension forces Eq. (2.26) for the fluid
pressure at the free surface reduces to
p = pg
A more rigorous statement than Eq. (2.27) is used to express the tangential viscous
stresses at the free surface:
aVn = 0
az
aV r = 0
an
where a/az is a gradient along an arbitrary direction 1in the plane tangent to the free
surface. This assumption substantially simplifies the setting of the viscous boundary
conditions in the numerical model.
Further simplification of the mathematical model is made by assuming that the free
surface is an adiabatic boundary, i.e. neglecting radiation and convective heat losses
to the air,
qJs = 0
This may be a valid assumption for the following reasons:
a) the wall heat transfer coefficient is usually much larger than that of the free
surface 2 ;
b) the total area of the free surface is usually much smaller than that of the metal
mould interface;
c) often free surfaces only exist during the filling stage; if there is a free surface after
filling but the air is confined by the mould walls then the air is heated up quickly,
reducing the subsequent heat transfer by radiation and convection.
2Radiation may be important for high temperature metals and alloys, such as steels. Mould
erosion has even been observed due to the radiation.
68
Eq. (2.28) is used for the heat flux per unit area at the metal/mould interface
q = h (Tmetal - Tmould)
The following computational methods are used to solve numerically Eqs. (3.1)-(3.3)
in three dimensions:
2. all temporal derivatives are approximated by the first order accurate backward
differencing; all spatial derivatives, except for the pressure gradients, are esti-
mated at the 'old' time level, i.e. explicitly, introducing time step size limits
to maintain stability of the numerical solution.
3. momentum advection terms are approximated using the upwind (donor cell)
differencing method which, being of first order formal accuracy3, possesses
transportive and conservative properties though introduces numerical diffu-
sIon;
69
6. free surface is represented by the VOF method, which proved to be one of the
most accurate technique to describe transient fluid interfaces; donor-acceptor
advection scheme was chosen to solve the conservation equation for the fluid
fraction function, Eq. (2.46), so that the sharpness of the free surface is
preserved; all other scalar quantities, such as enthalpy, that are ~iscontinuous
across free surface are advected in the same way.
7. the enthalpy solidification model is completed either by lever rule, Eq. (2.64),
or Scheil equation, Eq. (2.65). For pure materials a solidifying cell is kept at
the melting point until it is fully solid;
Fluid velocities and pressures are located at staggered mesh locations as shown in
Fig. 3.2:
• u-velocities and fractional areas Ax at the centres of cell faces normal to the
x direction;
• v-velocities and fractional area Ay at the centres of cell faces normal to the Y
direction;
• w- velocities and fractional areas Az at the centres of cell faces normal to the
Z direction;
70
volume, surface fluxes, surface stresses and body forces can be computed in terms
of surrounding variable values.
Most terms in the equations are evaluated using the current time-level values of the
local variables. This produces a simple and efficient explicit computational scheme,
though it requires a restrictive time step size to maintain computatTonally stable
and accurate solution.
The procedure for advancing a solution through one time increment, l:l.t, consists of
three steps:
3. Finally, when there is a free surface, it must be updated using Eq. (3.5) and
the donor-acceptor advection method, described in Section 2.5, to give new
fluid configuration. Energy values must be updated to reflect advective and
diffusive processes.
Repetition of these steps will advance a solution through any desired time interval.
At each step, of course, appropriate boundary conditions must be imposed at all
mesh, obstacle and free-boundary surfaces.
Curved obstacles, wall boundaries, or other geometric features are embedded in the
mesh by defining the fractional face areas and fractional volumes of the cells that are
open to flow. The fractional areas and volumes are incorporated into flow equations
as shown below in Eqs. (3.1)-(3.4).
71
Flow equations can now be written in the following form using Cartesian coordinates:
(3.1)
ou
ot
1 [OU
+ VF u Ax Ox
ou ou
+ v Ay oy + w Az oz =
1 1 op
Po ox + Gx - Dx
+ /I
0 [Ax ou
{ 2 ox ox 1+ oy0 [All ( oxov + Ou)
oy 1+ OZ
0 [Az ( oz ox ) 1}
ou + ow (3.2)
(3.3)
(3.4)
where (G x, Gy , Gz ) denotes body forces and D = (Dx, Dy, Dz ) is the drag force
employed to model flow in the porous media and is written in a general form
(3.6)
and, finally, the equation for the volume of fluid fraction function F
of
VF""at"
0 0 0
+ ox(uAxF) + 8y(vA y F) + oz(wAzF) = 0 (3.7)
72
For a one-fluid flow F represents the open volume fraction occupied by the fluid.
For a cell with dimensions dx, dy, dz and total volume
Vo = dxdydz
Vfluid = F· Vapen = F· VF · dx dy dz
Thus, fluid of constant density exists where F = 1, and void regions correspond to
locations where F = O. "Voids" are regions without fluid mass that have a uniform
pressure assigned to them. Physically, they represent regions filled with a vapor of
gas whose density is insignificant with respect to the fluid density.
When Eqs. (3.1 )-(3. 7) are applied to a partially blocked rectangular control volume
cell, then VF , Ax, Ay and Az represent the fractional open-to-flow volume and face
areas of the cell, respectively. New terms have to be introduced on the right-hand
side of momentum Eqs. (3.2)-(3.4), W8 x , W8 y , W8 z , associated with the wall shear
stress. If these terms are omitted then there is no wall shear stress because the
remaining terms contain the fractional flow areas (Ax, A y , Az) which vanish at walls.
The energy equation for the solid obstacle (mould) contain the complements of VF
and area fractions:
(3.8)
The use of the FAVOR method in numerical modelling influences the time step size
stability limit described in Section 2.10. For example, the CFL criterion (Eq. (2.75))
transfers into
(3.9)
It can be seen from Eq. (3.14) that cells with small fractional volume and large
face areas open to the flow can introduce very small time step size limits. Normally,
however, the small VF/A ratio is balanced by small velocities at such faces.
73
Fractional volume and areas in each cell are estimated at the preprocessing stage
and then used during the calculation. Within each cell containing metal/mould in-
terface the mould surface is assumed to be flat. FAVOR method allows the program
to calculate viscous shear stresses and metal/mould heat transfer more accurately
than if a stepwise approach to geometry description was used, since the fluid flow
direction at the walls, the wall location and area in each cell are e;timated with
better accuracy.
A generic form for the finite-difference approximation of momentum Eq. (3.2), for
example, is
n+l n+l
n.
Uj 1 k
+ ~t n+l . ( _ Pi+l,j,k - Pi,j,k
A
" pL..l.X·1+ 1
2
The advective, viscous and body force terms have an obvious meaning, e.g. FU X
means the advective flux of U in the x direction; V I SX is the x-component viscous
force; Gx includes gravitational and other body forces and W8 x is the viscous wall
force in the x direction.
The advective and viscous terms are all evaluated using old-time level (n) values
for velocities (wall shear stresses are implicitly evaluated as described below in this
Section). Old-time pressure values pn are used to get a first guess for the new
velocities.
Specific approximations chosen for the various acceleration terms in Eq. (3.2) are
relatively unimportant and have been described in Section 2.3. In FLOW-3D a
modified donor-cell approximation has been developed that retains its accuracy in a
variable mesh and reduces to a conservative difference expression when the mesh is
uniform [51,123]. This method approximates advective fluxes in the nonconservative
form U· Vu. Then the general form of this approximation is
FUX
+ (3.11)
where
UR = 0.5· (Ui+1,j,k A FR,i+1,j,k + Ui,j,kAFR,i,j,k)
UL = 0.5· (Ui,j,kAFR,i,j,k + Ui-l,j,kAFR,i-l,j,k)
74
and
All other acceleration" terms in the momentum equations are approximated by stan-
dard central differences. The difference techniques for viscous diffusion..processes are
fairly straightforward. First order differences are applied to velocity components to
obtain local shear rates, which are then multiplied by averaged cell-faced area frac-
tions. The results are then differenced to approximate the net viscous stress. All
quantities are evaluated explicitly in these calculations.
The approach for the wall stress in the w-velocity equation, for example, is as follows.
Wall shears influencing w can arise from wall areas located on x and y cell faces
surrounding w. For anyone of these faces, if the fractional area A is less than unity,
the remaining area fraction (1 - A) is considered to be a wall on which a stress is
generated. On an x face to the right of w, for instance, the force due to wall shear,
wS
_~
z -
( aw) ~ _211(1-Az~x2
ax 11 ax
Ax)(w - wo)
Similar stresses are evaluated at each of the four surrounding cell walls and their
sum is taken as the total stress. Wall stresses are included in an implicit way to
avoid possible numerical instabilities arising in cells with large wall areas and small
flow volumes. Although this makes the momentum equation for w implicit, it is
trivial to solve since it is linear in wand is not coupled to equations for other cells.
n+1 _ -n + .u.Atn+ 1 . _
n+l n+l
Pi+1,j,k - Pi,j,k Kn n+l
)
(3.12)
'+2')' U"',),"k
U " "k - u"',3,"k - "1" k
',3, ( p~x" 1
'+2
n +1
v i,j,k -n
Vi,)",k
+ .u.Atn+l . (_ p~nl,k - pi,t.l _ .n"
T/n
J 1 k
vn+1)
"" k (3.13)
P~yj+t "+2' ',3,
75
(3.14)
where, e.g.
[{:" l ' = 2 J(!J-I,],'k . J(!J- 1 'k
1+ ,],
I+ 2 ,],k J(!J-,
I,],k
+ J(!J-I+l,],k,
and u, v and w include all explicitly evaluated terms and the wall shear stresses.
For incompressible flows the continuity Eq. (3.1) can be interpreted as an elliptic
condition on the cell pressures and velocities. Velocities computed from Eqs. (3.12)-
(3.14) must satisfy the following discretised approximation of continuity Eq. (3.1)
n+l A ,. - A FR,I-l,],k
" n+ l A n+ l A FB,i,j-l,k
Ui,j,k FR,I,],k
n+!
Ui-l,j,k + vi,j,k FB,i,j,k -
v i,j-l,k
~Xi ~Yj
+ w~tl A ,. -
I,],k FT,I,],k
W!ltl
I,],k-l
A FT,I,],k-l
,. = 0
(3.15)
~Zk
Pi,j,k
(3.16)
8L"k
1,1,
8Pi,j,k
+
(see also Eqs. (2.39)-(2.45)). Eq. (3.16) is simply a Newton type of relaxation
process that will produce the value of p needed to make L = O. Then for u-velocities,
for example,
~p. 'k ~tn+1
un+! ----+ u~tl + I,J, (3.18)
i,j,k I,],k Po ~X I'+12,J,'k 1 + ~tn+! K':'+l 'k
I 2,J,
76
.6.p'I,),'k
Ui-l,j,k ~ Ui-l,j,k - (3.19)
Po .6.X'_l 'k
I 2 ,),
1 + .6.tn +1 I<~ 1 ,),k'
.- 2
In each cell the velocity values used in evaluating L and in Eqs. (3.18), (3.19) are
most current values available during the iteration process. To start the iteration
process, the new estimated velocities from Eqs. (3.2)-{3.14) are used with the pres-
sures remaining from the previous time step. Velocities located at zero area faces
are not adjusted.
In cells containing a free surface, that is, a cell containing fluid, but with one or
more empty neighbours, a different procedure is used (Section 2.4). The boundary
condition wanted in these cells is that the pressure is a specified value, p" at the
surface. The surface pressure is set equal to the neighbouring void region pressure,
Pg ,
(3.20)
where n is the index of the adjacent void region. This pressure is then imposed on
the solution by extrapolating it to the pressure, Pi,j,k, at the centre of the surface cell
assuming a hydrostatic distribution within the cell (Eq. (2.49)). The hydrostatic
variation depends on the net body force in the direction normal to the surface. This
surface pressure is not changed during the pressure iteration, but treated as a fixed
boundary value 4 .
A complete iteration consists of adjusting pressures and velocities in all full cells
according to Eqs. (3.16), (3.18) and (3.19). Convergence of the iteration is achieved
when all cells have Li,i,k values whose magnitudes are below some small number,
3 1
t . VF,i,i,k' Typically, t is of order 10- 8- , although it can vary with the specific
problem being solved.
The donor-acceptor advection method (Section 2.5) with the newly obtained veloc-
ities is employed to find the new time level fluid configuration by solving Eq. (3.7).
Other scalar quantities are advected using the same method as described in Section
2.6, i.e. if 8F is the quantity of fluid advected in x-direction from cell (i, j, k) to
4The fact that surface cell pressures do not take part in the iteration process effectively means
that surface pressure is evaluated explicitly. This may lead to a numerical instability if a void
pressure varies with the void volume as p V'Y = const, where 'Y is a constant parameter [85]
77
cell (i + 1,j, k) over the time tlt per unit cross-sectional area, then the amount of
enthalpy advected through the same cell face is
The new F values occasionally have values slightly less than zero or slightly greater
than unity. Therefore, after the scalar advection calculations have been completed,
a pass is made through the mesh to reset values of F less than zero back to zero
and values of F greater than one back to one. Accumulated changes in fluid volume
introduced by these adjustments during a calculation are recorded and printed out.
The total fluid volume is also printed. Volume errors after hundreds of cycles are
usually a fraction of a percent of the total fluid volume.
For application of free surface boundary conditions and advection of the F function,
it is necessary to assign an approximate normal direction to the surface. The free
surface in each surface cell is approximated by a plane normal to one of the coordi-
nate directions so that the normal points in the direction of an empty neighbour cell.
If the surface cell has more than one empty neighbour, then the direction selected
is the one having the largest F value in the opposite neighbour cell.
Once the inward normal direction has been determined, the surface location is de-
fined by a flat surface normal to this coordinate direction that extends to the proper
height.
A two-domain method is also used in the present work, but energy Eqs. for the
metal and mould are coupled at the same time level and solved simultaneously. The
latter is simplified by the use of a fully explicit formulation of the numerical analogue
equations. For a control volume P (which in this case conincides with a mesh cell)
the heat flux bq per unit area at each face is calculated as
C _ k Tp - Tad]
oq - a
~x
where Tadj is the temperature of the adjacent control volume across the face, ka is
the harmonic average of thermal conductivities of the two cells
(3.22)
78
and L\x is the distance between the two cell centres. The use of the harmonic
average, rather than the arithmetic one, is generally a more accurate procedure,
more obvious so in a case when the values of kp and kadj differ substantially. In the
present work thermal conductivity varies between: solid and liquid phases in metal,
and different materials in the mould.
If a cell face includes a metal/mould interface than the heat flux is estimated by
using Eq. (2.56). The same flux hq (times the face area) is used for both cells
sharing the face, thus, ensuring the conservation property. The total heat flux into
the cell is equal to the sum of all the conductive and interfacial fluxes multiplied by
corresponding areas.
Thermal conduction and heat transfer terms appear in both the fluid and obstacle
energy Eqs. (3.6) and (3.8). These terms are treated in analogus ways. The nu-
merical analogue of Eq. (3.6), with forward time differencing and first order spatial
approximations, is
Fn+1 n+1 F n+1 Hn
Vi i,j,k H i,j,k - i,j,k i,j,k _ xn
F,i,j,k po L\t n+1 - i,j,k
where X represents advection terms, T and TW are cell-centred fluid and obstacle
temperatures respectively, h heat transfer coefficient, W A an interfacial area, A a cell
face area, k an averaged heat transfer coefficient (Eq. (3.23)) and L\ an appropriate
spatial increment.
The same procedure is employed when solving the solid energy Eq. (3.8) wall tem-
peratures.
79
2-D cylinder: Shrinkage Modelling. Probl em header
$xput
remark='computational parameters' ,
icolor=l, autot=2, epsi=O.OOOl, twfin=llOO.O,
delt=l, itb=l, pltdt=lOO, prtdt=lOOOO.O, Num erical and
ifenrg=2, ihtc=2, ipdis=l,
cyl=l.O, ishr=l, physi cal parameters
remark='physical properties', and properties,
gz=-980.0, pcav=-lOOOOO.O,
rhof=2.54, rhofs=2.71, - speci fication of the
remark='thermal properties', nume rical model
tll=933, tsl=933, clhtl=3.97+09, required for the
cvl=0.896+07, cvsl=1.17e+07,
simul ation.
thcl=1.Oe+07, thcsl=2.2Se+07,
remark='boundary conditions',
wl=l, wr=2, wf=l, wbk=l, wb=2, wt=2,
tbcd=293.0, tbct(1,2)=293.0, tbct(l,S)=293.0, tbct(l,6)=293 .0,
Send
$mesh
nxcelt=12,
px(l)=O.O, nxcell(1)=7, Mesh specifications
px(2)=4.9, sizex(2)=O.7,
px(3)=12.0,
nycelt=l,
py(l)=O.O,
py(2)=l.O,
nzcelt=2S,
pz(l)=6.0,
pz(2)=13.6, sizez(2)=O.7, nzcell(2)=20,
pz(3)=30.0,
Send
Sobs
nobs=2, Obsta cle (mould)
iofo(l,l)=l,
cc(l)=-l.O, rah(1)=4.9, zh(1)=13.6, geom etryand
iofo(1,2)=2, physi cal properties.
cc(2)=-1.O, ral(2)=4.9,
kobs(1)=60000, kobs(2)=60000,
rcobs(1)=1.7e+07, rcobs(2)=1.7e+07,
hobsl(l)=l.Se+06, hobs1(2)=l.Se+06,
Send
$f! Initial fluid
flht=28.6,
oonfi guration
Send
$bf Baffle specifications
Send
$temp Initial temperature
tempi=1023, distri bution
Send
$motn Non-i nertial forces
Send
$grafic
ncplts=3, Graph ical requests
yc1(1)=O.S, yc2(1)=O.5, kontyp(1)=6,
ycl(2)=O.5, yc2(2)=O.5, kontyp(2)=5,
ycl(3)=O.S, yc2(3)=O.S, kontyp(3)=S, ictyp(3)=1,
nvplts=l,
yvl(1)=0.5, yv2(1)=O.S,
xloc(l)=O.OOl, yloc(l)=O.OOl, zloc(1)=21.8,
Send
$parts Marke r particales
Send C!n"l' ifications
VF,r, F,p,H,T,K
Care must be taken to resolve areas where the flow variables have steep gradients.
Discontinuity surfaces and internal boundaries represent such areas. VOF method,
for example, allows to resolve free surface boundary inside a computational mesh
cell and set the boundary conditions in the cell according to the amount of fluid in
it and the surface orientation. Another example is the use of the FAVOR method to
describe the obstacle surfaces inside the mesh cells. The FAVOR method requires
much less cells to represent a sphere, for example, than if a stepwise approach was
employed.
Further use of the FAVOR method has been made to improve the accuracy of the
metal/mould heat transfer calculations, taking into account that cell-centred nodes
do not, in general, lie on the interface, Assume that temperature profile between
the metal/mould interface and the nearby cell node is linear on both sides of the
interface. In this case if dXl is the distance from a fluid interfacial cell node to the
interface and dX2 is the distance from a corresponding mould interfacial cell node
to the interface, then the effective heat transfer coefficient for those nodes is [92]
h _ 1
eJ J - !kl.
kl
+ !h + ~
k2
where kl and k2 are the metal and mould thermal conductivities respectively and
80
1/ h is the interfacial thermal resistance. Then the heat flux across the interface can
be calculated using the metal and mould cell temperatures Tl and T2 , located at the
nodes, as
q = heJdTl - T2 )
This procedure gives a higher-order accurate approximation to the interfacial bound-
ary conditions than if h was used to estimated the flux. As will be shown below in
this Section, this approach is equivalent to introducing two additional, zero-volume
cells on each side of the interface. The idea has been first introduced and successfully
applied in conjunction with finite element method [90,5].
Consider two neighbouring cells, one of which is completely open to the flow (VF = 1)
and the other completely blocked (VF = 0), so that their common face represent the
metal/mould interface. According to the control volume approach, Tl is the average
fluid temperature in cell 1 and T2 is the average mould temperature in cell 2 (Fig.
4.1a). Both temperatures are located at the cell centres. The x-axis is chosen to be
normal to the interface positioned at x = o.
In the basic method the flux per unit area between the two cells is calculated as
(4.1)
The actual interfacial boundary condition, given by Eq. (2.56), involves metal and
mould surface temperatures, T61 and Ta2 respectively
(4.2)
with 1/ h representing the interface thermal resistance. For the case shown in Fig.
4.1a, Tal is located at the right face centre of cellI and Ta2 is located at the left face
centre of cell 2. To simplify the analysis, suppose that
(4.3)
To estimate the accuracy of Eq. (4.1), we will expand T2 in Taylor series about
x = 0:
1 aT 2
1aT 2 3
T2 = Ta2 + 2 ax Ix=O Dox + -4-a
x Ix=O
Dox + O(Dox
2 ) (4.4)
where Dox is the size of cell 2. Substituting T2 from Eq. (4.4) into Eq. (4.1) and
taking into account Eqs. (4.2) and (4.3) yields
1 aT
ax Ix=O Dox + O(Dox )
2 )
h· ( Tsl - T82 - --
2
1 aT
Qo - h - -
2 ax Ix=O Dox + O(Dox 2 )
1 ~x)
Qo ( 1 + "2h T + O(Dox2) (4.5)
81
where k is the mould conductivity and
Qo = - k -
aT (4.6)
ax Ix=o
It is clear from Eq. (4.5) that boundary condition approximation given by Eq. (4.1)
is of the first order formal accuracy with respect to the cell size ~x. Fo.! low thermal
conductivity mould materials, such as sands, the first order term on the right-hand
side of Eq. (4.5) may appear to be larger than the zeroth order term 1 . The heat
flux Q1 is overestimated compared to Qo because the zeroth and first order terms
have the same sign.
The basic, first order heat transfer method represents metal/mould interface by a
row of cells with a non-zero thickness. Therefore, the interface has a thermal capacity
proportional to the layer thickness and the mould specific heat. Since interfacial cell
sizes and mould geometry vary across the computational domain, the value of Ql
becomes strongly mesh dependent.
The interfacial temperature T2' can be expressed as a function of Tal and T2 using
Eq. (4.9)
(4.11)
IFor silica sand the typical values are h = 1.5 X 103 Watt m- 2 K-l, k = 0.6 Watt m- 1 1(-1
and ~x = 5 mm, and the first order term on the right-hand side of Eq. (4.5) is over six times
larger than the zeroth order term.
82
and Eq. (4.10) transforms into
uxp
A C -dT2 (4.12)
dt
where
h1~
heff = 2 uX
k
(
4.13 )
h + 1A
2ux
ISthe effective heat transfer coefficient for the cell-centred temperature T 2 • The
expression for he!! takes into account the resistance of the interface and of the
mould material lying between the interface and the location of T 2 .
he!! . (Tal - T 2 )
2
1 aT 1a T 2 3 )
he!!' ( Tsl - Ts2 - -- ~x - --2 ~x - O(~x )
2 ax Ix=o 4 ax Ix=o
Qo - -1 h ef !· ~ Ta 2
A 2
ux + O( ux
A 3)
= Qo + O(~x 2 ) (4.14 )
4 ux Ix=o
Eq. (4.14) shows that Q2 is a second order accurate approximation to Qo. Besides,
the second order term is proportional to ~:r which is normally much smaller than
the first order derivative ~ in Eq. (4.5).
Furthermore, for the error function solution for a one-dimensional heat flow into a
semi-infinite media, ~:r Ix=O = 0 and the accuracy of Q2 increases to a third order.
In the example above this analytical solution exists if Tal = const and h --+ 00 or
The modified method is more accurate because it resolves the temperature profile
in the interfacial cell by a linear function, instead of the uniform temperature used
in the basic method. The main assumptions in the modified method are:
• the interface has no heat capacity, i.e. the heat fluxes on both sides of the
interface are equal;
• the temperature profile between the interface and the centre of the cell is linear
as a consequence of the first order approximation of the heat conduction terms.
83
The analysis of this Section suggests that the gain in accuracy of the modified heat
transfer method over the basic one is larger for larger values of t!:1xhj k.
h· [To - T(O,t)] = - k -
aT (4.15)
ax Ix=O
where hand k are constant parameters and To is a constant ambient temperature.
Initially the sand is at a uniform temperature T j • The analytical solution, T(x, t) is
given by [22]
For the numerical simulation the length of the sand is 200 mm, which is a good
approximation of a semi-infinity during the first hour of simulation for the mould
properties listed in Table 1.
First, a coarse uniform mesh is used consisting of N = 10 cells with !:1x = 20 mm,
as shown in Fig. 4.2a. The mesh is chosen in such a way that the left face of
the first cell in sand coincides with the metal/mould interface. The temperature at
the centre of this cell, that is at x = 10 mm, is recorded during the first hour of
simulation and the results of the first order and modified methods are compared to
the analytical solution in Fig. 4.3a. Fig. 4.3b shows the predicted interfacial heat
fluxes, Ql and Q2, normalised by the analytical solution and Fig. 4.3c shows the
total predicted energy fluxed through the interface, that is fci Ql dt and fci Q2 dt, also
normalised by the analytical solution.
It is clear from Fig. 4.3 that the modified heat transfer algorithm describes more
accurately both temperature and energy evolution in the mould. The total energy
flux, predicted by the modified method, is within 20% of the analytical solution
soon after 500 s, while the basic procedure gives an error of more than 30% even at
t=3600 s.
The second computer simulation is carried out with a non-uniform mesh with N=lO
cells and the interfacial cell size !:1x = 2 mm, the rest of the cell sizes gradually
increasing towards the right domain boundary as shown in Fig. 4.2h. Temperatures
at two positions are recorded: at x = 10 mm, as before, and at x = 1 mm which is
at the centre of the interfacial cell.
84
The same time step size, !::1t = 5 s, is used throughout the simulations, which is
sufficiently small to ensure numerical stability for the two meshes.
Fig. 4.4a shows the comparisons of the temperatures histories given by the two
algorithms and the analytical solution and Figs. 4.4b, c show the comparisons of the
interfacial heat flux and total fluxed energy respectively.
The use of the non-uniform mesh improves the accuracy dramatically. As expected,
the modified method shows less mesh dependence. In a general case, this would
make the heat extraction along the metal/mould interface more uniform.
In both Figs. 4.3c and 4.4c lines representing the total fluxed energy predicted by
each of the two methods lie on opposite sides of the exact solution: the first order
algorithm overestimates the flux and the modified one underestimates it.
If the assumption of Eq. (4.3) is dropped then the same treatment is applied to cell
1 in Fig. 4.1b to introduce the metal surface temperature Tt'. The result is (see also
Eqs. (2.61), (2.62))
( 4.17)
where
! h.Ja....!L
h ~Xl ~X2
eff = Ih...h-
2 ~Xl
+4 l...h...a..
4 ~Xl ~X2
+ Ih..a..
2 ~X2
(
4.18
)
and indexes 1 and 2 refer to cells 1 and 2 respectively and 1/ hef f is the effective
thermal resistance between locations of temperatures Tl and T2 , consisting of the
interface resistance 1/ h and resistances of materials on both sides of the interface,
!::1xd2kl and !::1x2/ 2k2.
Consider now a situation where the metal mould interface is inside a cell, i.e. the
interfacial cell is partially blocked. Let VF,l be the fractional open volume in it
and To and Tl the cell-centred metal and mould temperatures respectively. Suppose
that the interface is normal to the x-axis and cell 2 is the next fully blocked cell
neighbour along the axis, with temperature T2 (Fig. 4.1c).
Assuming first that To represents the metal surface temperature, the energy balance
between the two cells is
( 4.19)
(4.20)
85
In the basic heat transfer algorithm
is the distance between centres of cells 1 and 2 irrespective of the value of VF,I'
Therefore, only the thermal resistance between these two points, r = Sx/ k, is taken
into account in Eqs. (4.19), (4.20). The resistance between the interface and the
centre of cell 2, estimated by assuming linear temperature distribution between the
nodes, is
(4.21 )
Compared to r a , the basic method overestimates the resistance if VF,1 > 0.5, and
underestimates if VF,l < 0.5.
o - (4.22)
(4.23)
(4.24)
Eq. (4.22) constitutes the heat flux balance at the interface and
A _ (1 - VF,d ~Xl
uXI' - 2
It can be seen from Eqs. (4.22)-(4.24) that the total thermal resistance between the
interface and cell 2 centre is equal to ra'
d~ k
(1 - VF,I)pC-
dt
= hef! . (To - Td + -;:--.
uxv
(T2 - TI ) (4.25)
86
with
Applying a similar procedure on the metal side of the interface, the heat flux between
Xl' and Xl", the centre of the open volume VF,l' Vi of cell 1 with te~perature To,
is expressed as
(4.26)
where
(4.27)
kme , kmo being the metal and mould thermal conductivities respectively, and
Numerical solutions for the interfacial heat flux and total energy, fluxed through
the interface, at t = 30 s from the beginning of the simulation are shown in Figs.
4.5a and 4.5b, respectively. Each numerical result is normalised by the corresponding
analytical solution. Compared to the first order method, the modified method results
are more accurate and less dependent on the value of VF in the interfacial cell, though
it is obvious that the used mesh resolution is not sufficient to give a satisfactory
accuracy for all values of VF • It is interesting that their appear to be an optimum
value of VF for each method, which may be dependent on the material and interface
properties.
In a general case, when the interface in a cell is arbitrarily orientated, Xl' and Xl"
are associated with the geometrical centres of the blocked and open volumes of the
interfacial cell. 6.Xl' and 6.XI" are approximated by
87
where
( 4.30)
~Xn is the size of the cell along the coordinate axis in the direction nearest to the
direction of the interface normal and Al is the area of the interface in the cell.
The integral on the right-hand side of the expression for the enthalpy given by Eq.
(3.4) can be written as
If constant values of Cs and C/ are assumed then the first and third integrals on
the right-hand side of Eq. (4.31) are easily evaluated. The second integral contains
specific heat of the solid/liquid mixture and is variable across the freezing range if
C s =I- C/
( 4.32)
and
rT/ rl dT
iT. CdT = - Jo [JsC s + (1 - Is)Cd dIs dIs (4.33)
The effect of the contribution to the enthalpy given by Eq. (4.34) or (4.35) is
that during solidification both latent heat and specific heat have to be removed
88
as the temperature decreases between T/ and T s , unlike solidification at a constant
temperature. This generally leads to larger solidification times for alloys than for
pure materials, given the same latent heat. The effect is more apparent for larger
specific heat and freezing range.
The simplest way to model alloy solidification is to assuming linear variation of solid
fraction with temperature in the freezing range
(4.36)
This model can be enhanced by adding optional lever rule and Scheil Eqs. (2.64),
(2.65) (Section 2.8). It can be seen from Eq. (4.35) that, within the freezing range,
enthalpy becomes a nonlinear function of the temperature, even for the linear latent
heat release,
where fs(T) is defined by one of the solidification models. For the Scheil's model
Eq. (4.37) transforms into
1
H(T) GaTe + 2 [(fB,e + f,l(T))C + (2 -
B fSle - fs(T))Gd(T - Te)
+ (l-fB(T))L (4.38)
Fig. 4.6 shows the variation of the fraction of solid and enthalpy for an AI-4.5%Cu
alloy (Table 3) for each of the solidification models, given by Eqs. (2.64), (2.65),
(4.36), (4.37) and (4.38).
In the modified solidification algorithm Eq. (4.37) is solved by the Newton iteration
method in which at iteration n + 1
89
where I(Tn) is the difference between the left- and right-hand sides of Eq. (4.37).
The starting temperature is taken from the previous time step. Iterations stop when
T/-Ts
1Tn+! - Tn 1< f, f = 10000
It usually takes one to five iterations to satisfy the convergence criter.ia since H(T)
is a smooth function and dH/dT > 0 for all temperatures. Generally, the number
of iterations depends on the cooling rate and the time step size.
When T = Te (for the Scheil model), Eq. (4.38) becomes a linear function of Is
The solid fraction function, Is, is not explicitly calculated in the original code nor is
III used as an output variable for graphical and numerical analysis at the end of the
simulation. However, spatial and temporal evolution of the fraction of solid function
could provide important information about the progress of solidification.
A modification of the code has been made to include solid fraction as an additional
variable. It is either a function of temperature for solidification models defined by
Eqs. (2.64), (2.65) or (4.36) or is an independent variable for pure metal and eutectic
solidification processes. Solid fraction, stored in an additional array, serves as an
input for other parts of the program, such as drag function calculation (Section 4.3)
and shrinkage simulation (Chapter 5).
A modification has also been made to include solid fraction as an output quantity.
The modification is based on the existing routine for storing and displaying numerical
data, so that all display options for other variables are applicable to the solid fraction.
Those include:
• history plots for the evolution of the solid fraction with time at specified loca-
tions in the casting.
• 2-D spatial colour and contour line plots. This can be combined with velocity
vectors and marker particle plots.
90
• both 2-D and 3-D spatial plots can be stored with predefined time increments
for subsequent animation.
Additional dashed single contour line plot of fs = 0.5 is provided to visualise the
approximate position of the solidification front, especially useful for nC;lJ'row freezing
range alloys. This is done similarly to the existing FLOW-3D method of visualising
the position of the free surface.
Fig. 4.7 shows an example of a 2-D colour temperature plot with flow velocities and
solidification front. More examples of the use of the solid fraction can be seen in
other sections.
When metal solidifies, its position in space no longer changes, i.e. the velocity of the
solid phase is zero (or equal to a constant non-zero value in the case of continuous
casting process, for example). One of the ways to account numerically for the
corresponding change of fluid momentum is to increase the solid/liquid mixture
viscosity coefficient /lej I in a solidifying cell so that for a completely solid cell /leI I ~
00 [124]
fs ) -2.5 l.,p
/-leI f = /-l ( 1 - - (4.39)
fa ,1'
where fa,1' is the solid fraction at close packing of the solid (equals to 0.62 for solid
spherical particles), similar to the critical solid fraction fs,er at which any flow seizes
according to the shrinkage models described in Section 2.9.
As a result of the "effective viscosity" model, given by Eq. (4.39), the solid ma-
terial will adhere to the mould, including moving walls. The latter occurs during
continuous, twin-roll and melt-spinning casting processes.
Viscous terms in the discretised momentum equations will have to be included im-
plicitly since the explicit formulation will cause severe time step size limitations as
J.LeJ J increases. The implicitness will make momentum Eqs. (3.7)-(3.9) coupled with
each other, requiring an efficient solution for the first guess velocities.
In the present model the cell momentum change due to solidification is accounted
for by using the drag force concept. The solid phase material is assumed to be at
rest or moving at a uniform speed.
91
Here we will only consider the drag due to solidification. The main role of the drag
force is to reduce the fluid velocity to zero when it solidifies. Assuming that solid
material is at rest with respect to the computational mesh, the solidification/melting
process is approximated by using a drag coefficient, J«(T), that is a function of tem-
perature or fraction of solid. The drag should be effectively infinite when material is
in the solid phase. At intermediate states, consisting of a mush, the drig should also
assume an intermediate value. The behaviour of the drag coefficient for 0 < fs < 1
is very much the choice of the user since it is difficult to validate any particular drag
model. The only restriction is that J( must be non-negative.
For a Darcy-type solidification drag, the flow in the mushy zone is thought of as a
flow of the liquid fraction through a fixed matrix structure created by the dendrites.
Coefficient J( is then dependent on the volumetric solid fraction [125] (Eq. (2.70))
J( ftv (4.40)
rv (1 _ fsv)3
Variations of Eq. (4.40) include introducing a critical solid fraction, so that the
drag is equal to zero for fs < fs,cr, before the solid phase forms a rigid structure.
Such regions are often called "slurry zones" , where solid material floats freely in the
metal [2].
Huang et al suggested two expressions for the drag coefficient: one for a columnar
dendritic region, I<C, and the other, I<e, for an equiaxed [117]. The drag is dependent
on the specific surface area available for flow in the dendritic mushy region
(4.42)
Another equation for the drag force coefficient is presented in this Section. The
derivation is based on the simplifying assumption that there is no frictional inter-
action between solid and liquid phases. This assumption is more appropriate in the
situations when the liquid/solid interface can be represented with a surface rather
than a mushy zone.
92
Consider a cell full of fluid (Fig. 4.8a). At time t, ms and m/ are masses of the solid
and liquid phases in the cell respectively
(4.43)
and
-- m6
fs m
Since the solidified fluid velocity is zero, the total momentum of the cell is equal to
the liquid phase momentum
(4.44)
Suppose that there are no body and pressure forces acting on the cell and that the
liquid phase can flow freely, i.e. neglect viscous friction between the two phases. In
that case during solidification, assuming also that volumetric shrinkage effects are
negligible,
u/ = const, m = const (4.45)
The cell momentum change over the time dt due to the loss of mass by the liquid
phase to the solid phase, using Eqs. (4.44), (4.45), is
Far the solid/liquid mixture velocity in the cell, which appears in the discretised
equations,
(4.47)
dM = mdu m (4.48)
Now zqualising the right-hand sides of Eqs. (4.46) and (4.48), divided by dt, and
taking into account Eq. (4.47) yields
um df,
(4.49)
1- f, dt
the solution of which is
um = Uo (1 - f,) (4.50)
where Uo is the liquid velocity for fa = O. It is easily seen that U m = 0 when fs = 1.
93
Comparing the right-hand side of Eq. (4.49) with the definition of the drag force,
Eq. (3.10), gives an expression for the drag force coefficient
J( = 1 dfs
(4.51 )
1 - fs dt
In a more general case, when m and are not constant because of £he net inflow
u/
or outflow of fluid in the cell and the action of other forces, the drag coefficient
becomes
J( = _1 dms
(4.52)
ml dt
The right-hand side of Eq. (4.52) needs to be defined for fs 1 as both the
numerator and denominator vanish at that point. By putting
f{ = 00 for fs =1
the condition of an infinite drag in the solid phase is satisfied.
The drag is zero in a partially solidified cell if the solidification front in it does not
move. If part of the cell melts, i.e. dms < 0 then f{ = 0 as the cell momentum
does not change because the newly created liquid has zero velocity immediately after
melting. Thus, condition f{ ~ 0 is also satisfied.
A I-D simulation of pure aluminium solidification has been carried out to test the
drag model. The geometrical setup is shown in Fig. 4.8b. The metal is flowing
at a constant velocity Uo over a horizontal mould plate in the absence of pressure,
body and viscous forces. The results are shown in Fig. 4.9. As the metal in a cell
starts to solidify, the drag force is applied and the cell horizontal velocity changes
in accordance with Eq. (4.50), as shown in Figs. 4.9b - d for three consecutively
solidifying cells.
The drag force given by Eq. (4.52) does not represent any physical force acting on
the fluid. It is just a consequence of using a solid/liquid mixture velocity, defined
by Eq. (4.47), rather than modelling the flow of each phase separately.
Eq. (4.52) is applied here to both mushy and plane front solidification. However,
frictional drag in the mushy zone may dominate the flow and the developed drag
model will not be sufficiently accurate.
94
(a)
Dx
1-- ~I
::::;::::::::=:::::::=::::::;=;:;:;:;:;: ::::;
(b)
~ldX I~ Dx-dx ~ 1
(c)
IDX1" 1
~ ~
~
1~
.... DX2'
~
I
DX1'
(a)
200mm
Figure 4.2. uniform mesh (a), with 6x=20 nun, and non-uniform mesh (b),
with the smallest cell size 2 nun, used to simulate the one-
dimensional heat flow used as a test for the modified heat
transfer algorithm. White dots show where the t emperatures
werecompared.
coarse Me sn
TerrpeuIU'e aI x-l a rrrn
1050
: : : : : :
.....~!':"!~.t=~.~~.::'.~::"r.=.=·=:::·r.·:=::::+·:==-=·r.·· .···· .· ...... t................... .
esc
~
750
of
(0)
I 660
650
450
Co arse Mesn
Inet1acia1 Heal FlDc
:::::::::::::::::IL:::::::::::::t:::::::::::::.:::c::::::::::::::1::::::::::::::::t:::::::::::::::::1:::::::::::::::::t:::::::::::::::::
I 3.5
i: iii iii
1Jne. $
Co arse Mesn
TOIai lner1acia1 Ene<!1f FlDc
3.6
I
w
2.5
I
(c)
1.5
w • ._ 1 .
1
~ ••••• , " . , 0 •• • • "1····...... ...... ' .
Cl5
·················l;/::·:~······(··············I···· · ············1················-(···· ············1 ....·........····t·················
1
IJne. $
(n)
tine, S
Fine Mesh
Irter1acia1 Heal FlDC
1.1'T"---"',CY',----,------.,.-----..,.-----,----,-----..,------.
"
"""
:'
!~ . . . . ...........:". ....... ....:.................
105 :: : ......... .. ..... 1..- ................:..................:.................
_1 • • • •• • •• • •• •• • ••
(b)
OJ
I
I{'···~;··Li
}\,
.. · .............. H--....::..~.....,,.-:~==
.
:1
___+--___+-~~~~~~......"""d
····_··-t········r········r········
.;:;: ;1.... :·
., i
° ::
:
:I
:
i:
:
::~
1 l.';. i.
tine, S
Fine Mesh
Taallrter1acia1 Energ( FlDC
1.25;..,.------,..----...,-------,-----....,----:-----:----.,.-----.
1.2
i 1.1
1.1
w
tine. S
1 •• • • •• •••••••••• .• •• ••••• •••• • •••••••••••••••••••••• •• •••••••••• •••••••••• • •••• • • •• •••• •• ••••• ••••• •• • • ••••• •• • •••• ••••••••••• ••••• • ••• • •••••••• • • •
(a)
(]
o.
(11)
0.' 0.' OA
qlElI1 voUne,
0.1
VF
0.. 0.' G.' ...
O.
:--.-----------..
........ 1.4
0.6 -
........• .. -.-:~.
1.3
0.7
01
- 1.21 .,~
0.6
C/l ~
0.6
0.4
1.1
0.3
0.2
,+----.----y----:1::----:r---...--...:::;.-r---i-ll.e
~o ~o
temperature, K
temperature, oK
0.0
330.0 mm
1 2 3 4 5
(a)
me ta l
(b)
lOOmm
(a)
1.0 T'----~--------l-
(b)
60.0
time, s
1.0 -r--------..____- -__ " , . - - - - - - 1 -
(c)
60.0
time, s
1.0 1
(d)
60.0
time, s
Figure 4.9. 1-0 heat and fluid flow simulation result for the flow
shown in figure 4.8,(b): (a) velocities and solidification
front (dashed line) at t=23.0 s; (b)-(d) 1 - fraction of
solid, 2 - 1 K is the drag coefficient, 3 - hori-
l+K~t
zontal velocity for cells A, B, C, respectively.
Chapter 5
A solidification shrinkage/expansion model based on full fluid and heat flow equa-
tions is presented in this chapter.
Consider a computational cell (or control volume) of open volume Vo and boundary
E containing a mixture of liquid and solid of volumes VI and Va, respectively. The
mixture volume, V, and mass, m, are
m = pV (5.2)
respectively, where F is the cell fluid fraction. The mixture density, p, is
P=
PIV/ + P.. Va (5.3)
V
V=v,+V/ (5.4)
Over a time increment dt, the mixture volume in the cell changes due to solidifica-
tion, dV·, and a flux of the liquid phase through cell faces, dV('. The total change
in the mixture volume, dV, is
where
S = (1 _ps) PI
dY,
dt
~
Va
(5.7)
constitutes the change of volume per unit time per unit cell volume due to phase
transformation.
The corresponding change of the mixture density, from Eqs. (5.1) and (5.3), is
dP -- ( ps - - V,dV/ = dp*
pI ) V/dY, V2 + dp' (5.8)
where
(5.9)
is the change of density due to phase transformation only. Similarly, the fraction of
fluid function is
dF = dV = dF* + dF' (5.10)
Va
and
dF* = (1 _Ps) PI
dY,
Va
Sdt (5.11 )
The integral form of the continuity equation, Eq. (2.1), for the cell open volume Va
IS
PI' f
lr.
Vn du = - ~
at lvf P dr (5.12)
The liquid phase density is used in the surface integral because it is assumed that
only pure liquid phase can flow. At those parts ofthe cell boundary that are blocked
by the solidified material, velocity is zero and no contribution is made to the integral.
The volumetric integral on the right-hand side of Eq. (5.11) is calculated over fluid
volume V, since Vo - V is void.
Replacing the density on the right-hand side of Eq. (5.11) by the average cell-centred
value, given by Eq. (5.3), yields
96
Eq. (3.6) for the fluid fraction function must be rewritten to take into consideration
Eqs. (5.9) and (5.10)
at + ~.
aF Va Jr.f vnFda = S (5.14)
We will neglect the variation of the mixture density in the momentum~quation, Eq.
(2.22), so that only the drag force in it accounts for solidification.
where qr. denotes heat conduction and transfer fluxes through the cell boundary. HI
is the enthalpy of the liquid phase even if two phases are present along the boundary
since only the liquid phase flows. After replacing P and H by the cell-centred average
quantities, Eq. (5.14) can be rewritten as
Eqs. (5.6), (5.7), (5.12), (5.13), and (5.15), together with the momentum equa-
tions, represent the mathematical model of the solidification process with volumetric
changes.
dVa is defined by solving the fluid and heat flow equations. Therefore the calculation
of its source term and the solution of these equations are coupled. An implicit
representation of the source S, defined by Eq. (5.6), would require a simultaneous
solution of the fluid flow and enthalpy equations, making the numerical algorithm
complicated and inefficient.
At the beginning of n + 1 time cycle the value of dVa is estimated in each cell using
values of fluid fraction and solid fraction functions known at times t n - 1 and tn
(5.17)
97
where
PI' fs
fs'V = --....:.....;.--=--=-----
pds + Ps(1 - fs)
is the cell volumetric fraction of solid.
(5.1S)
The use of ~tn+l, rather than ~tn, ensures that the subtracted volume at cycle n + 1
exactly corresponds to the volumetric changes given by Eq. (5.16).
The corresponding change of the cell fluid fraction is given by Eq. (5.10). After
all the advection contributions to the cell fluid content are calculated (as described
in Section 2.5), giving an intermediate value Pi~/:l, the final adjustment is made to
include the volumetric change
(5.19)
The left-hand side of Eq. (5.15) can be represented in the following form
~( PF H) = aF H F H ap (5.20)
at Pat + at
Combining Eqs. (5.15) and (5.19) yields
PI .
aF H
--at" PI
= P Vo
(- PIfJ HI Vn du + qE )
1 ap
--PIFH (5.21)
E P at
The enthalpy equation is used in this form in the numerical calculations. It effec-
tively implies that P1F H, instead of pF H, is used as the cell total heat content.
This simplifies the process of specifying initial and boundary conditions though an
additional source term appears on the right-hand side.
(5.22)
(5.23)
9S
If Rn < 0 then stability analysis shows that I = n + 1 to ensure unconditional
numerical stability. Eq. (5.22) is easily resolved for un+! since the source term is
linear. The explicit formulation, i.e. I = n, is stable only if ~tn+l < l/Rn.
If Rn > 0 then both explicit and implicit methods give exponentially growing solu-
tions since the analytical solution is an exponent. The explicit methoa, however, is
preferable because it always gives zero overshoot, while the implicit method gives
zero overshoot only if At n+ l < -1/ Rn .
The drag coefficient at n + 1 time level is calculated from Eq. (4.55) in which
ml = P1F n (1 - f~)Vo
Combining expressions for ml and m5 with Eq. (5.17), the drag coefficient is given
by
(5.24)
while 8Liljlk/8pilj,k, in Eq. (3.17) remains unchanged since the source term sn is
formulated explicitly and does not depend on pn+1.
For a solidifying cell S < 0 if p. > PI, leading to a negative value of the cell velocity
divergence, div v. The latter means that there is a net flux of liquid metal into the
99
cell. Similarly, if S > 0, the velocity divergence is positive and the remaining liquid
is pushed out of the cell.
Special care must be taken if a cell (or a cluster of adjacent cells) containing both
phases are surrounded by solidified cells. In this case velocities at cell faces are equal
to zero and if S t= 0, then there is no solution for the cell pressure ifconditions of
incompressibility and continuity are not dropped.
If S > 0 then the remaining liquid must flow out of the cell. In practice pressure
in such confined region will increase causing solid phase deformations and, perhaps,
even partial remelting. Since these effects are not described by the present model,
such situations are resolved by setting to zero the net volume increase in the confined
region, i.e. S = 0, so that numerical solution does exist. This is equivalent to
assummg that the liquid becomes compressible to accommodate the increase in
volume.
The solution is more complicated ifthe source term, S, in a confined cell is negative.
In this case
LOOk
1,3, = -S-°k
'", > 0
and, according to Eq. (3.22), pressure will decrease indefinitely.
In practice this will cause dissolved gasses to evolve into bubbles and to open internal
voids which, after complete solidification, remain as shrinkage defects. A critical
pressure, Per, at which bubbles start to appear is employed in cavitation models [3].
The same approach is employed in the present numerical model. The value of the
critical pressure is unimportant so long as it is much smaller than the pressures in
the casting to avoid ambiguity on where a cavity should be created. This is because
cell pressure serves as the only flag to open a cavity. Since the actual value of Per
is arbitrary and internal cavity volumes at every time step are defined only by the
amount of the volumetric shrinkage, the gas evolution process is not described by
the model.
When
Pi,i,k ::; Per (5.27)
then the cell changes from an 'incompressible and continuous' cell, for which con-
100
vergence criterion is
I L""kl<t·ViF.""k
1,1, ,t,J,
Pi,i,k = Per
and the value of Li,i,k is defined by iterations. The condition of fluid continuity in
the cell is dropped and the cell then serves as a feeding basin for other solidifying
cells.
This algorithm as based on the cavitation model in the original version of FLOW-3D.
In the latter, input parameter PCAV defines the critical pressure.
In the FLOw-3D cavitation model, cells in which P < PC AV are defined at each
time step before the iteration begins, so that during iteration their pressure is driven
to PC AV and a net outflow of fluid from these cells is allowed to open a cavitation
void. The 'cavitation' status of the cells is not changed during iteration. The rate
of the cavity growth is a free parameter and is defined by the user.
In the present shrinkage model the status of a cell can be changed by the algorithm
during iterations. Suppose after I iterations pressure in a cell is greater than Per:
1
Pi,i,k > Per
and
L~I,J,"k > 0
If
1+1 _ 1
Pi,J,k - Pi,i,k
+ ~Pi,J,k
1+1
< Per
then the following adjustment is made
~/+1
Pi,i,k -
_Per -
1
Pi,j,k (5.28)
so that
1+1
Pi,i,k = Per
101
Velocity adjustments for the cell (i,j, k) (Eqs. (3.24) and (3.25) and corresponding
expressions for the other components) are made using 6.p given in Eq. (5.27). For
the next iteration cell (i,j,k) is marked as an 'internal cavity' cell.
In the subsequent iterations, [' > 1+ 1, the magnitude of LL,k is not required to be
below the value of the convergence criterion f if pL,k does not increast!. As soon as
the value of pL,k exceeds Per, the 'incompressible and continuous' status is returned
to the cell, even if a cavity has already been created in it.
• Situations are avoided where iterations are carried out on equations which do
not have a solution. The latter would happen if pressures and velocities are
iterated in a confined region and the continuity condition is imposed on every
cell in the region 1 .
Allowing cells with P ::; Per to have a non-zero velocity divergence ensures the
existence of a solution .
• More than one cell can appear to have pressure below Per. This allows the
openning of internal voids in several cells simultaneously, i.e. in one time step.
In the shrinkage model the critical pressure is substantially lower than casting (at-
mospheric) pressure. If the confined liquid region consists of several cells then a
uniform pressure decrease is required until the lowest cell pressure in the region
reaches the value of Per' This uniform pressure adjustment represents an error com-
ponent of an infinite wavelength which, as mentioned in Section 2.3.2, is the slowest
to converge.
102
solution for velocities in the region since only pressure gradient is present in the
incompressible fluid equations.
If the net source term S in a confined region is positive than the situation is resolved
as described in Section 5.2.2. If S is negative than it will take only one iteration for
the minimum pressure in the region to fall below the value of Per and the correspond-
ing cell (and possibly some of its neighbours) will be the location of the opening of
a new internal void.
As the solidification proceeds internal cavities grow. Eventually some cells may turn
into surface cells or become completely empty. Pressure in the empty cells is set to
be equal to PeT and treated similarly as all other void regions, and pressures in the
surface cells are set in the usual way, that is using .Eq. (3.26).
In many situations feeding criterion, such as critical solid fraction fs,er, remains a
useful albeit arbitrary simplification of shrinkage modelling.
Ideally, the model should adjust fluid pressures and velocities according to the drag
force by solving the flow equations. For example, if liquid has to flow through a high
drag area to feed a solidifying region then the pressure in the region should decrease
as the drag force increases. Eventually the pressure reaches the critical value and a
cavity starts to open in the solidifying region, though some feeding may still occur
despite the high drag.
To avoid solving flow equations in the presence of a high drag force, a critical value
of solid fraction is used as a cut-off parameter for flow.
For the drag force coefficient, given by Eq. (5.23), the value of the term K' = l+kAt
in Eqs. (3.25)-(3.27) lies usually above 0.9 and only for /s -+ 1 does this term plays
a substantial role (depending, of course, on the size of the time step as well).
Therefore, for pure materials, for which Eq. (5.23) has been derived, the drag force
103
poses little difficulty in achieving convergence. In accordance with this, it is assumed
that for pure metal a cell can feed and be fed if fs < 1. Flow through a cell seizes,
i.e. K' = 0, when the cell becomes fully solid. In other words, for pure metals
fll,er = 1.
For alloy systems flow in a mushy zone experiences much larger resistance. Instead
of using a high value of the drag coefficient in such locations, a critical solid fraction
fll,er < 1 is used so that flow into a cell ceases when fll ~ f.,er' Eq. (5.23) is still
used for solid fraction below fs,er' This treatment of alloys is an idealisation of the
flow in mushy zones, when for low fs values the melt represents a mixture of liquid
and freely floating solid crystals. As solidification progresses crystals starts forming
a rigid structure and the flow resistance increases steeply over a small range of fs
values.
• in alloys flow resistance can be anisotropic. Feeding is easier along the colum-
nar crystals than across, and a single value of fs,er is not sufficient to resolve
such situation. In that sense the critical solid fraction gradient [112], yo fs,
(Section 2.9) or anisotropic drag coefficient [117] (Eqs. (4.41), (4.42)) offer
some flexibility;
• the value of fs,er varies for different materials, casting geometries and solidifi-
cation conditions. The right choice of fs,er is almost an art.
The present shrinkage model has been developed to predict macro-defects in castings,
since micro-porosity simulation requires modelling of gas evolution, interdendritic
flow, nucleation and dendrite growth. The value of fll,er divides the total shrinkage
volume between micro- and macro-porosity. Here the critical fraction of solid method
is employed to predict macro-cavities, unlike most other shrinkage models. Micro-
porosity will inevitably be present as a consequence of using fll,er < 1 for alloy
systems, but its spatial distribution is not claimed to be accurate and will not be
compared with experimental data 2 •
The value of the critical solid fraction is included in the list of input parameters
specified in PREPIN.INP file as FRSLCR. However, for pure materials it is always
equal to one.
21f the average micro-porosity in the solidified casting is known a priori, e.g. from the liquid
metal gas content, then the solid phase density P. can be reduced correspondingly to account for
that. Modelling of micro-porosity can then be avoided by setting f.,er = 1.0.
104
5.3 Simplified Shrinkage Model
The shrinkage model described above is capable of modelling liquid metal volumetric
shrinkage/expansion due to phase transformation, including the induced flows and
cavity formation, occurring in both phase change directions. The model-is applicable
to most general 3-D mould and casting configurations.
The cost of that is the necessity to solve fluid flow equations not only during the
filling stage, but also during solidification making calculations CPU intensive. Fur-
thermore, the time step size is smaller than in a no-flow situation. Usually the most
severe restriction comes from the free surface gravity waves stability criterion, Eq.
(2.77) of Section 2.10 3. Since the solidification stage is normally much longer than
the filling stage, a significant increase in the CPU time is obvious.
A simplification of the full shrinkage model (hereafter referred to as model M1) can
be found by considering the following [110,126]:
• Often solidification occurs in such conditions that the liquid metal free surface
is horizontal and flat all the time.
With these assumptions, another shrinkage model (M2), which does not solve the
fluid flow equations, has been developed. In brief, at each time step
1. Isolated liquid regions are identified in the casting. This time an isolated
region is defined as one bounded by solidified cells or mould, and free surface.
2. The total volumetric change due to phase change over D.t is estimated in each
of the regions and removed from the top of it. The 'top' of a region is defined
by the direction of the gravity.
3The CFL criterion (Eq. (2.75)) is not a problem since normally flow velocities are small during
solidification.
105
5.3.1 Numerical Implementation of the Simplified Model
It is essential that gravity is normal to one of the coordinate planes and the free
surface is aligned to that plane.
At time step n + 1 the volumetric change in each cell dViJ,k, occurred during the
previous time step £ltn, is estimated using Eq. (5.16).
Eq. (5.13) is not solved in every cell directly. For an isolated liquid region R/ of
vol ume Vln the total change of volume is
dvt =
and
'"'
~
. . k) . F!'·
(1 - fn",v,'", I,),
k
(i,j,k·)eR,
where f:'v,i,j,k is the volumetric solid fraction in cell (i, j, k) at time tn.
pnf'kI. = fn.tI,'.J,
I"~,
.. k.
and
dVi n = dVin - Vi".
k* = k* - 1
and the procedure is repeated from step 2;
4. For VI". ~ dVin. Since gravity is normal to k* coordinate plane the liquid cells
in this plane should have the same value of the fluid fraction, Fn+l. However,
the fluid mixture level in the cells can only be reduced at the most by their
liquid phase content. Therefore,
106
(a) an estimate to Fn+! is made assuming that the mixture level is reduced
in every cell by the same amount
because the source term in Eq. (5.20) is balanced by the surplus of enthalpy delivered
by the net fluid flow into the cell. As the flow term is excluded from consideration
in the M2 model, so must be the source term.
The drag force coefficient calculation is also simplified because it is only necessary
to know when the drag becomes infinite to identify isolated liquid regions
• CPU time is saved as momentum and continuity equations are not solved and
no iteration is needed to obtain velocities and pressures. Enthalpy equation is
also simplified. The time step size, controlled by heat conduction/transfer, is
usually considerably larger that when the full system of equations is solved.
107
The limitations of the M2 model are:
• Free surface must be plane and normal to the gravity; its orientation must not
change during solidification.
• The mesh must be chosen in such a way that the free surfaceIS parallel to
one of the coordinate planes. In other words, body forces direction have to be
constant and parallel to one of the mesh directions.
The simplified shrinkage model is similar to the one developed by J .Evans [120]
(Section 2.9) with two major differences: first, the M2 model is applicable to general
3-D casting configurations (with the limitations stated in the preceding paragraph),
while Evans' model was designed to model only centre line porosity in cylindrical
ingots. Secondly, cells here are allowed to be partially filled, while in Evans' model
they can only be either full or empty: a cell is not emptied until enough shrinkage
volume has accumulated to empty it completely. Furthermore the same value of the
critical fraction of solid is used in all directions in the M2 model.
The choice between M1 and M2 models is made by the user through the input
parameter ISHR: ISHR=l initiates model M1 and a value of ISHR above 1 initiates
model M2.
108
liquid
feeding flow
semi-solid zone
(high drag zone)
r-L
I _ oJ - V,
,=(1+1
no
Fk*=fs,Vk*
m .!!.!!!!!!!!!!!
II I 1IIJ11 I! 11 III
1 dV=dV-V *
:iHmHHmm k
litiiiiiiiliiiiii yes
!ii!iliW!iiWi k*=k*-1
~=F
!:!:!!!::::!:: :::
111 11 1: 1111111111
ii iiiiiiiiiiiiil
lJ!!H:ii::::::l!:
!!!! !!!!!!!!!!!!
+
mllml;nl!l!!
(F' - first guess fo the
yes
for F",I,J,'k*=FI,J,
' 'k
p
f,s,V, I,J,
" k*'F'I,J,'k*>F' cell Fi,j,k*=Js, v, i,j,k*'F";,j,k*
i,j,k*
I VFk*= VFk*-F"i,j,k*
noJor all
k* cells dV=dV-(F\j,k*-Fi,j,k*)
B
Vt=Vt-Vi,j,k
This Section describes briefly experimental casting simulations which are used for
validation of numerical predictions given in the later sections. The experiments were
carried out independently of the present authour.
Sand moulds. Castings were moulded in 100 AFS silica sand bonded with 4% sodium
silicate [127]. The sand was hand rammed over the wooden patterns and CO 2
hardened at 80°C. The mould surface was dried prior to casting by exposing it to
high intensity lamp light overnight.
Aluminium castings. The aluminium used for castings was commercially pure. The
copper was weighed and put into the crucible after the aluminium had been melted
to make the AI-4.5%Cu alloy. The metals were heated in a 3 MKHz induction
furnace.
A porous plug degassing method using argon was employed for both pure aluminium
and the alloy. In order to improve degassing, a flux was used on the melt surface to
avoid air entrapment during degassing.
109
6.1a shows the overall geometry and dimensions of the casting. The metal
was poured directly into the mould at 90°f( superheat. The filling stage took
a fraction of a second so that initial conditions in the metal and mould are
close to an instantaneous filling situation. Initially the mould was at room
temperature of 293 0 f(. Immediately after the filling had been <:,ompleted the
open top of the metal was covered with an insulating board preheated to 770 0 f(
to minimise heat losses from the metal free surface.
Thermocouples were placed both in the mould and metal, as shown in Fig.
6.1b, and temperature readings were recorded for half an hour after the pour-
ing. According to the thermocouple located at the centre of the casting (th/ c
No. 1 in Fig. 6.31), the total solidification time, t s , is 754.0 s. Figure 6.3
shows the sectioned casting with a shrinkage cavity at the top.
2. Chill mould cylindrical castings. These castings were produced by J.Evans [126].
The mould geometry and dimensions are shown in Fig. 6.4 The top cylindri-
cal mould was made of grey cast iron. This was fixed on top of a pure iron
cylindrical mould of the larger diameter, positioned on a sand base, using a
18 mm thick steel plate.
One casting was made in pure iron (E2) and a second in 0.6%C, 1.1 %Mn,
0.2%Si steel (E3). Since any gas evolution would have affected the shape of
the shrinkage cavities, the castings were fully killed. This required the addition
of 0.15%Si and 0.05%AI to the pure iron and 0.05%AI to the alloy steel.
Both castings were top poured, the pure iron at 1843 0 f( and the alloy at
1873 0 f( corresponding to 34 0 f( and 64 0 f( degrees of superheat, respectively.
The top of the casting was then insulated using Kaowool board. Solidification
times are not known exactly since no thermocouple data is available.
Figs. 6.5a, b show the sectioned castings. The pure metal casting pas two
major centre-line shrinkage cavities: one at the top and the other, internal
cavity, at the location where the mould diameter changes. The alloy casting
shows a similar pattern but the two cavities are connected by a narrow passage
and the lower cavity extends farther down the casting axis than in the case of
the pure iron casting.
110
Fig. 6.7 shows the centre plane section of the solidified casting with a well
defined cavity in the feeder.
4. T-shaped sand casting. These simple shaped casting experiments allowed the
influence of changes in gating and feeding on the size and location of shrinkage
cavities to be assessed [129]. In particular, the position and size-of cavities at
the T-junction were observed. Casting dimensions and overview are shown in
Figs. 6.S and 6.9, respectively. The top of the feeder is open to atmosphere.
The runner system had been designed to ensure control of the filling conditions
of the casting. A reservoir at the top of the sprue was first filled with the melt,
keeping the sprue entrance closed with a stopper (Fig. 6.Sb). The stopper was
removed as soon as the melt temperature, measured by a single thermocouple,
reached a specified value.
The filling time was varied by changing the diameter of the choke, dch , at
the bottom of the sprue. Three values were used: dch = 8 mm, 15 mm, and
22 mm.
Feeding conditions at the T-junction were varied by changing the vertical to
horizontal section thickness ratio, R a / b = alb (Fig. 6.Sa), viz. R a / b = 0.6, 1.0,
and 1.67.
Two values of pouring temperature, T = 963 0 K and T = 1013° K, were used
corresponding to 30 0 K and SOD K of superheat, respectively.
(a) Pure aluminium (E5-E16). Figs. 6.lD and 6.11 show the centre plane
section of the casting for pure aluminium with typical shrinkage cavities.
The table in Fig. 6.12 summarises the results of this set of experiments.
When R a / b = 1.0 (E7-EI2) a small cavity appears in the horizontal sec-
tion of the T-junction. The cavity is either oval in shape and positioned
at the junction (Fig. 6.llb) or flat and wide and positioned in the hor-
izontal section (Fig. 6.lla). Figs. 6.13a, b show the X-ray photographs
of these cavities. The bottom and top surfaces of the horizontal section
underneath and above the cavity are slightly concave. This could reduce
the size and position of the cavity.
The cavity at the T-junction for Ra / b = 1.67 (E13, E14) consists of two
separate parts (Fig. 6.lDe). A likely reason for that is the aspiration of
air at one of sharp corners of the junction due to a low pressure in the
metal. The same could have happened in the casting shown in Fig. 6.11b
(castings E10-E12 in Fig. 12), where the cavity has a round bubble-like
shape with smooth surface.
111
Two variations of the mould design were made in an attempt to eliminate
the defect. First, a 10 mm wide chill was inset into the mould (EI5) and
a 10 mm thick insulating pad was used in the second case (EI6, Fig.
6.14), both leading to the disappearance of the shrinkage cavity as shown
in Fig. 6.15.
(b) AI-4.5%Cu alloy {E17-E21}. When an AI-4.5% was used instead of pure
metal no shrinkage occurred in the T-junction for Ra / b = 1.0 (EI8-EI9)
but castings were subjected to some microshrinkage. The amount of
the latter depends on the hydrogen content of the melt. A small cavity
though is present at the junction in the case of R a / b = 1.67 (E20-E21).
Fig. 6.16 shows centre plane section of the alloy castings and other ex-
perimental results for the alloy are summarised in the table in Fig. 6.12.
Typical primary cavity shapes in the feeder for pure aluminium and alloy
castings are shown in Fig. 6.17.
This Section presents validation results of the full shrinkage model against two simple
test of the model first and then the experiments listed in Section 6.1. All computer
simulations are summarised in Appendix A while Fig. 6.19 gives the computational
efficiency data for some of the calculations, together with predicted solidification
times.
Constant values of material thermophysical properties are used for both for metal
and mould, though solid and liquid phase properties can differ:
112
6.2.1 2-D Test Casting Simulations
The two simple simulations presented in this Section assumed instantaneous filling
with uniform initial temperature distribution at the melting point of the metal.
A. Casting with a narrow feeding passage (SI.I, Appednix A). In this two-dimensional
pure aluminium-chill casting simulation the mould has a narrow passage through
which the left-hand part of the casting (section A) is fed at the early stage of solidi-
fication by pressure head in section B, as shown in Fig. 6.20a. This feeding pattern
prevails until the passage is blocked by fully solid metal (J.,cr = 1.0 for pure materi-
als). The direction of the feeding in section A then changes to the opposite leading
to a cavity appearing at its top (Fig. 6.20b). The critical pressure, Perl at which
internal cavities would appear in the metal, was set at -10 4 N/m 2 • The ambient
gauge pressure IS zero. As was said in Section 5.2.3, Per can have any value below
the gauge pressure.
History plots for the pressure and vertical velocity in a cell located at the centre
of section A are shown in Fig. 6.21a. First, pressure decreases continuously due
to the decrease of the pressure head in section B. Then its value drops to Per =
-10 4 NJm 2 with a simultaneous change of the vertical velocity direction when the
narrow passage solidifies.
As expected, predicted feeding flow velocities are higher at the beginning of solidi-
fication, when the cooling rate is larger, and the average is 0.8 - 1.0 mm/ s. At the
end of the process they decrease to around 0.2 - 0.4 mmJ s.
Fig. 6.22 shows the resulting porosity distribution in the solidified casting and
Figs. 6.21b - d show the iteration count for f = 0.001 s-l, the fluid volume, VI,
and the mean fluid kinetic energy, Em, evolutions. The iteration number increases
temporarily when the passages narrows at t ~ 12.0 s and feeding becomes more
difficult.
The volumetric shrinkage is 6.3%, defined by the values of the liquid and solid phase
densities: 2540 kgJm 3 and 2710 kgJm 3 , respectively. The theoretical shrinkage
volume, 6. vth, is defined by equation
113
and can be reduced by decreasing fl. Therefore, it can be concluded that the full
shrinkage model predicted the shrinkage volume practically exactly.
The mesh is uniform in both directions and the time step size is fl.t = 1.6 x 10- 2 s
throughout the calculation as dictated by the free surface stability criterion (Eq.
(2.77)). The total simulation time is 50.0 s while solidification finishes in 44.0 s.
The total CPU time 2 for this problem is tcpu = 128 s with the average iteration
count per time cycle of Nit = 15.
Fig. 6.23 shows the final cavity shape and the developed bottom-top asymmetry of
the metal and mould temperatures due to the shrinkage induced gap formation at
the top of the casting. As expected, the cavity is not located at the position of the
last-liquid-to-solidify, i.e. at the casting centre. This result is an improvement from
that of Fryer et al [121] who showed the cavity to be at the centre of the casting.
Predicted feeding flow velocities average around 0.1- 0.2 mm/ s. Due to the stricter
convergence criterion of f = 10- 6 s, fl. Vc constitutes only about 10- 5 % of the to-
tal metal volume and fl. Vp corresponds to the theoretical value with the 0.002%
accuracy.
A series of simulations of casting El described in Section 6.1 has been carried out to
test the mesh sensitivity of the shrinkage model. The filling stage is excluded from
the simulation since the casting was filled very quickly. Instead, uniform initial
temperature distribution is assumed at 90° I( superheat.
The properties of the insulating material, obtained from the supplier, are given in
1 Since l > 0 then each cell of volume A V has a net inflow or outflow of the order of l . .6. V after
the convergence criterion is satisfied. The total volumetric error is recordered at each time step as
an output variable AVe (see also Section 7.2.4).
2 All simulations were made using a Silicon Graphics Indigo 4000 computer.
114
Table 6. Initially the insulator was at 770° 1<.
Fig. 6.24 shows the computational domain for simulation S2.1 (see Appendix A).
Only half of the casting was enmeshed due to axisymmetry. The meshes were chosen
to be uniform in both x- and z-directions with ~x = ~z = ~ to ensure equal
resolution at both horizontal and vertical metal/mould interfaces. '["he following
meshes were used:
The time step size in each case was defined either by the free surface or by the
conduction stability limit, depending on whether the are liquid surface cells or
not, and in average equals to 6tS 2 .1 = 1.0 s, 6tS 2 .2 = 0.15 s, 6tS2 .3 = 0.04 sand
6ts2.4 = 0.17 s.
The total CPU time, the average CPU time per time cycle and per iteration, together
with the time step size, are given in Fig. 6.25a as functions of the total number of
cells N. The iteration number does not vary substantially between the four cases,
averaging to 3-5 iterations per time step increasing to around 10 iterations only
for the finest mesh (S2.4). Simulations were continued until the last liquid metal
solidified. Since ts is different in all four cases, tcpu in each case corresponds to a
somewhat different physical time (see Fig. 6.25b).
Predicted cavity shapes and temperature distributions at the end of solidification for
all four cases are given in Fig. 6.26 (compared with Fig. 6.3). Fig. 6.25b shows the
predicted depth of the shrinkage cavity, hili" and t, in per cent of the experimental
results as a function of N. For the finest mesh (simulation 82.4) the cavity depth is
underestimated by 14% and the solidification time is 27% above the experimental
value of 754 s.
Fig. 6.25c shows ~ Vc and the error of ~ Vp in relation to the theoretical value of the
shrinkage volume ~ vth given by Eq. (6.1). The errors are small, less than 0.3%, and
the accuracy of the shrinkage volume prediction is mainly defined by the convergence
errors ~ Vc. Both ~ Vc and ~ v" show little dependence on mesh resolution as does
the average feeding flow velocity in the bulk of the liquid metal also plotted in Fig.
115
6.25c. (~0.02 mm/s).
The progress of solidification for simulation S2.2 is given in Fig. 6.27. Figs. 6.28a, b
shows the evolution of the horizontal and vertical velocity components in a cell at
the casting geometrical centre. The absolute values of both components grow as the
solidification front approaches the cell, reaching a maximum of about 0.05 mm/ s.
The share of the horizontal feeding in the total flow in the cell also increases, although
the vertical velocity is by far the dominant one. Figs. 6.29a - d show the evolution
of the time step size, iteration count, the total fluid volume and the mean kinetic
energy.
This Section includes modelling of casting E2 and E3 described in Section 6.1. The
pouring stage is omitted since the experimental filling times are short, besides, fast
cooling rates largely define the temperature distribution in the casting.
The volumetric shrinkage is 3.85% for both the pure metal and alloy castings, defined
by the densities of the liquid and solid metal (see Tables 4 and 5).
116
6.2.3.1 Pure Iron
This simulation is referred to as simulation 53.1 (see Appendix A). The value of
the heat transfer coefficient, h = 2370.0 Wjm 2 K corresponds to the maximum heat
transfer coefficient for a bare chill [126J. Fig. 6.32b - e shows predicted shrinkage
cavities, temperature distribution (degrees Celsius) and solidification front position
at t = 50 s, t = 100 s, t = 150 s and at t = 200 s. Solidification took ta = 177 s.
Comparison with the shrinkage defects given in Fig. 6.5a shows that the predicted
shape and position of shrinkage cavities are very close to the experimental result.
Fig. 6.33 shows the plot of the maximum horizontal and vertical velocities in the
feeding flow. The flow velocity increases as the vertical liquid channel in the upper
part of the casting narrows, before it is completely solid. At its maximum it reaches
25 mm/ s which is two orders of magnitude larger than velocity in the bulk of the
liquid metal in the lower part.
• 53.2: the heat transfer coefficient reduced by a factor of ten uniformly along
the metal/chill interface, that is h = 237.0 W/m 2 K corresponding to the
minimum value for a coated chill [126J. The result is given in Fig. 6.34a;
• 53.3: the heat transfer coefficient reduced by a factor of ten only for the lower,
wider part of the chill mould (Fig. 6.34b).
• 53.4: the heat transfer coefficient reduced by a factor of ten only for the upper,
narrow part of the chill mould (Fig. 6.34c);
Solidification times in these cases are 457 s, 208 sand 448 s, respectively. Cavity
shapes predicted in simulations 53.1 and 53.2 are similar although in the latter the
metal does not freeze over in the narrow part of the mould. This suggests that
solidification in the lower part of the mould in 83.2 becomes somewhat faster than
in the upper part at later stages of the process, so that liquid metal in the vertical
channel is drained downwards before it has time to solidify.
There is a small centre line porosity region in casting 83.4 and the primary shrinkage
cavity is much shallower than that in casting 83.1 while casting 83.3 has a mush
larger secondary cavity. The variability of the heat transfer coefficient along the
interface appears to playa larger role in the final cavity shape than the actual value
of h.
117
6.2.3.2 Alloy Steel
The alloy solidification and shrinkage model was used here to simulate the steel
casting E3 shown in Fig. 6.5b. The lever rule was employed to describe the latent
heat release and two values of the critical fraction of solid, I.,er = 0.~7 (83.5) and
Is,er = 1.0 (83.6), were tested.
Fig. 6.35 shows the distribution of the fraction of solid function at intermediate
stages of solidification for Is,er = 1.0 and the final predicted shape of shrinkage
cavities is shown in Fig. 6.36. Note that for Is,er = 0.67 part of the shrinkage volume
goes into the uniform 1% microporosity in the final casting (the total volumetric
shrinkage is the same for both simulations). Comparison of the result with the
experimental casting E3 given in Fig. 6.5b shows clearly that IIJ,cr = 1.0 gives a
better approximation. This may be explained by the relatively narrow mushy zone
in the casting during solidification (Fig. 6.35) hence the alloy flow behaviour is close
to that of a pure metal. For the same reason the shape of the shrinkage cavity in
this case appears to depend little on the latent heat release model. Fig. 6.36c shows
predicted cavities for the linear mode for is,er = 1.0 (simulation 83.7).
The solidification times are (83.5) ts = 190.1 s, (83.6) t. = 192.5 s, and (83.7)
ts = 193.6 s. The average ts is 8.5% larger than that for the pure metal simulation
83.1 since the alloy solidus temperature is lower than the melting point of the pure
metal.
There is a slight increase in CPU time in the calculations for alloy modelling because
temperature in every cell is found iteratively from the value of the enthalpy (see
8ection 4.2.1). The number of iterations required is normally 2 or 3 leading to
insignificant increase in computational effort: in average 10% per time step.
Casting E4 of 8ection 6.1 has been simulated in three dimensions (simulation 84.1,
Appendix A). As before, the pouring stage was substituted here by the instantaneous
filling assumption at the uniform 10° /{ superheat. The top of the casting is open to
atmosphere and insulated3 • Figure 6.37 shows the enmeshed part of the mould. As
before, geometrical symmetry was used to simulate the solidification in only half of
the casting. The total number of cells is 9361, of which 4,645 are open to the flow,
i.e. VF < 1.
3Free surface is an adiabatic boundary in all simulations in the present work.
118
The predicted solidification time is 556.1 s which is 13.5% larger than the experi-
mental value. Fig. 6.38 shows a sequence of 2-D centre plane sections of the casting
and four 3-D views are given in Fig. 6.39. The predicted shape of the primary cavity
in the feeder is in good agreement with the experimental result (Fig. 6.7). There is
a small secondary cavity in the casting itself, indicating that the feeder is too small.
The experimental casting, however, is sound.
Feeding velocities in the liquid bulk do not exceed 0.4 mm/ s. The evolution of the
metal volume, mean kinetic energy, iteration count and time step size are given in
Fig. 6.40 and other calculation results for this simulation are presented in the table
in Fig. 6.19.
Only the bottom part of the pouring basin was included. A pressure boundary
condition, linearly decreasing with time to approximate the change of pressure head
in the basin, and a constant pouring temperature, obtained from the experiment,
was set. The filling was assumed complete when metal in the feeder was at the
level of the sprue entrance, as was done during the experiments. At that time the
boundary pressure was balanced by the pressure head in the feeder and practically
no liquid entered the casting.
The filling time is controlled by the size of the choke, dch , at the bottom of the sprue
in the same way as was done in the experimental procedure. The following are the
results for vertical to horizontal section ratio, R a / b = 1.0;
Figure 6.43 shows a 3-D view of the filling process for dch = 15.0 mm. Sloshing
and free surface break up are clearly visible in the horizontal section of the casting.
Figure 6.44 gives a 2-D view of the filling with metal and mould temperature contours
for the 80° I< superheat case (simulation S5.1). The evolution of the mean kinetic
energy and total metal volume in the casting during filling are given in Fig. 6.45.
The former serves as an indicator of the quality of filling. In the present case the
kinetic energy has two peaks: one at the time when the metal hits the bottom of
the sprue and the other, smaller, when the metal falls down the lower section of the
T-junction.
The effect of the filling time on the temperature distribution in the casting at the
end of the filling is given in Fig. 6.46 where results for d ch = 15.0 mm and d ch =
8.0 mm with 80° I< of the initial superheat are compared. At the end of the latter,
slower filling, the T-junction has already started to solidify since the minimum
temperature is 6° I< below the melting point. For dch = 15.0 mm there is still
at least 25° I< superheat left in the metal after the pouring stopped. Thermocouple
readings are also indicated for comparison for dch = 15.0 mm. The average deviation
of the predicted temperatures from the experimental results is 9.2% of the maximum
temperature difference in the casting (60° I<). The total solidification time appears
to be virtually constant for all filling rates, including the instantaneous filling and
equal to 475 s (Fig. 6.19).
Solidification progress, feeding and shrinkage formation of simulation S5.2 are shown
in Fig. 6.47. The initial metal superheat was 30° I<. A small liquid region became
isolated at the T-junction leading to a shrinkage cavity there .
. Figure 6.48 shows that the size of the cavity in the T-junction is practically inde-
pendent of the superheat.
Figure 6.49. shows a 3-D view on the predicted shrinkage cavities for the three filling
times and for an instantaneous filling, all with 80° I< superheat (simulations S5.3,
120
S5.1, S5.4 and S5.9, respectively). The shape of cavities in the junction is thin in
the vertical direction and wide in the horizontal plane, similar to the experimental
results (Figs. 6.lOb and 6.lla and castings E7-E9 in Fig. 6.12). The cavities lie
mostly in the horizontal section slightly protruding into the vertical section. That
the cavity size increases as the filling time decreases. There is no cavity when the
filling time is 22.5 s (d ch = S.O mm, Fig. 6.49a), and the cavity size-is maximum
when the filling is excluded form the simulations (Fig. 6.49d).
The general trend is that the faster the filling the larger is the cavity in the T-
junction. This may be explained by the fact that a slower filling results in a larger
temperature gradient in the direction of feeder-junction at the end of the filling,
thus promoting a directional solidification from the junction towards the feeder.
• At the end of the filling stage, at t = 7 s, metal at T-junction has the lowest
temperature and metal in the sprue is the hottest.
• Metal at the bottom of the sprue freezes in around 100 8 after the beginning
of the pouring.
• Metal in the T-junction and at the bottom of the feeder freeze practically
simultaneously arnoud 250 seconds after the beginning of the filling.
At the end of the filling stage, the measured temperatures are 726° C, 736° C and
744°C, and the predicted ones are 726°C, 735°C and 745°C, respectively. At t =
500 s, the largest error in the predicted temperatures, 80°, is at the bottom of the
sprure, where the overall temperature variation was 350°. This may be explained by
the lack of resolution in the sprue. The smallest error, 10°C, is at the bottom of the
feeder where the temperature varied by 150°Cin 500 s. The error and the variation
in temperature at the third location are 30°C and 250 0 , respectively.
The filling stage was simulated with f = 0.005 8- 1 . At the end of the filling cal-
culations were stopped and restarted from the same time with f = 0.00025 8- 1 .
A smaller f was required during solidification because velocities are significantly
smaller at this stage than during filling. The distribution of the variables at the end
of the filling are used as the initial conditions for the restart calculations. The table
121
in Fig. 6.19 contains information on the mesh sizes, run times, time step sizes and
dreezing times for this calculation.
The influence of the geometry on the shrinkage in the junction is simulated for the
fixed choke size dch = 15.0 mm with the two values of the pouring temperatures
used in Section 6.S.5A. As in S5.1 to S5.3, each simulation here included filling,
heat transfer, solidification and the use of the full shrinkage model M1.
For R a / b = 0.6 the castings appears to be sound irrespective of the pouring temper-
ature, in agreement with the experiment (Figs. 6.lOa and 6.51a, b). The pouring
temperature also does not significantly affect the predicted defects in the T-junction
for R a / b = 1.67, where a large cavity is predicted in the vertical section of the T-
junction, as shown in Fig. 6.51c - J. This is in agreement with experimental result
shown in Fig. 6.lOc. The lower part of the predicted cavity did not show the exper-
imental separate shrinkage defect in the T-junction. The difference may be due to
air aspiration during solidification in the experiment, which was not accounted for
in the simulation.
In this Section the simplified shrinkage model M2 is tested against the experimental
and M1 model results. Further modelling is carried out of the AI-4.5%Cu alloy
T -shaped casting. In each case the geometry, mesh and initial conditions are kept
the same as those in the corresponding full model calculations of Section 6.2. The
only parameter needed to be changed in the input files was ISHR, which specifies
the choice of a model (Section 5.3.1).
The four meshes described in Section 6.2.2 are employed here to investigate accu-
racy of predictions and calculation speed of the simplified model. Simulations S2.11,
S2.12, S2.13 and S2.14, correspond to simulations S2.1, S2.2, S2.3 and S2.4, respec-
tively. Figure 6.53 shows the final predicted cavities for the four cases and these are
compared with those shown in Fig. 6.26.
122
The accuracy of the predicted final values of hsh are presented in Fig. 6.25b in per
cent of the experimental result. The % errors of depth predictions are almost twice
as large as those given by the full shrinkage model. This is despite the fact that
~ Vp is practically exactly equal to the theoretical value. The volumetric error ~ Vc
is zero in the M2 model since there are no volume errors associated with fluid flow
an free surface advection.
The overall cooling rate and metal volume evolutions for simulations 82.12 with the
M2 model and 82.2 with the Ml model are compared in Fig. 6.54. In both cases
the differences between the two model predictions are negligible compared with the
total variation of the parameters.
The total CPU time and the CPU time per time cycle for simulations 82.1-82.4, using
the Ml model, and for simulations 82.11-82.14, using the M2 model, are plotted in
Fig. 6.55. Due to the fast convergence of the Ml model solution for this casting, the
CPU time per time step for the M2 model constitutes only 30 (simluation 82.11)
to 80% (simulation 82.14) of that of the Ml model. Additional gain in speed for
the M2 model is due to larger time step size, so that, for example, for the S2.14
simulation the total CPU time is more than 10 times smaller than that for the S2.4
simulation.
A. Pure Iron. Simulation S3.1 was repeated with the simplified shrinkage model
(S3.11). Figure 6.56 shows the progress of solidification with temperature contours
and solidification front position for pure iron solidification. The final shapes of the
cavities predicted by M1 (Fig. 6.32e) and M2 models are similar. The solidification
times predicted by the two models differ by 2.5 s which constitutes 1.4% of the
average time.
The time step size was tl.t = 1.36 s, which is 186 times larger than that for the M2
model simulation. The total CPU time is only 13.4 s, that is 730 times smaller than
that for the Ml model. The predicted change of volume ~ v" constitutes 99.997%
of tl. vth, achieving a better accuracy than that of the M1 model.
B. Alloy Steel. 8imulations S3.5 (is,er = 0.67) and S3.6 (isler 1. 0) were also
123
repeated using the M2 model (S3.15 and S3.16, respectively). Figure 6.57a, b shows
results of the shrinkage cavity predictions which are close to those shown in Fig.
6.36a, b of the M1 model. The total solidification times for cases S3.5 and S3.15 are
190.1 sand 189.0 s, respectively, differing by less than 0.6%.
The metal volume evolution for the pure iron and the alloy (Is,er = 1.8), calculated
by the M2 model, are given in Fig. 6.58. It takes 11 s (6.2 %) longer for the alloy
to solidify because pure metal freezes at a higher temperature.
The M2 model alloy simulations took 270 times less CPU time than those by the
M1 model (Fig. 6.19).
The simplified shrinkage model was employed to repeat simulation 84.1 of Section
6.2.4 (simulation 84.11). Figure 6.59 shows the progress of the solidification and
shrinkage formation in the centre plane. Following a tendency of the simplified
model to give shallower cavities for solidification in sand moulds than the full model,
the shrinkage cavity in Fig. 6.59f does not penetrate as far into the casting as does
the one in Fig. 6.38f. The difference between the predictions of the two model in
this case is in favour of the M2 model since its result is closer to the experiment
(Fig. 6.7).
The cooling rate and volume evolution are given in Fig. 6.60 in comparison to the
Ml model results. The simplified shrinkage model gives the exact shrinkage volume
while the full model underestimates it by 6.7% of the total shrinkage volume due to
convergence errors.
The M2 model took 7 times less CPU time to complete the calculations and the
average time step was two times larger than those in the M1 model simulation.
124
6.3.4 T-Shaped Casting: dch = 15.0 mm
Shrinkage simulations of the T-shaped casting using the M2 model were initiated
by restarting the program from the moment the filling stage ended, i.e. the filling
stage and the starting conditions are exactly the same for both shrinkage models.
However, all velocities were set to zero when the simplified model was-switched on.
Simulations were made only for the pouring temperature of T = 1013 0 J( (80 0 J(
superheat) since the full model (Section 6.2.5) and the experimental results showed
little variation of the shrinkage in the T-junction with the initial superheat. The
choke size was also fixed at dch = 15.0 mm.
If an instantaneous filling is assumed then the cavity in the junction is larger, posi-
tioned higher and further away from the feeder suggesting that the horizontal section
of the junction freezes earlier (Fig. 6.63). This result is predicted by both models
(simulations S5.1O and S5.12).
When the filling stage, which is the same for both models, is included, the gain in
speed for the simplified model is by the factor of 4.3. If only the solidification stage
is considered, then the increase in speed is by a factor of 6.1.
Given that the results of the two models differ little, the M2 model can be used
to carry out efficiently a number simulations with varied physical and/or numerical
parameters. In this particular case the following four variation fo the mould design
were made in an attempt to eliminate the cavity in the T-junction:
1. S5.13: a chill block placed along the vertical section of the T-junction (Fig.
125
6.64a).
3. 85.15: an insulating pad, 10 X 100 x 100 mm3 , placed underneath the horizontal
section (Fig. 6.64c)).
The cavity in the junction disappeared in cases 1 and 4, in agreement with the
experimental result (Figs. 6.65 and 6.15). No experimental data is available for
comparison for the other two cases. In case 2 the cavity became even larger since so-
lidification in the junction progressed slower and a larger liquid volume was trapped
there when the horizontal junction froze to the right of the insulator. In case 3 the
cavity, though very small in size, moved away from the junction but remained in
the horizontal section.
A. Ra/b = 1.0. For R a/ b = 1.0 there was no macroshrinkage in the T-junction for
is,er = 0.67 (85.16) and for is,er = 1.0 (85.17), as shown in the experment (Figs.
6.66b, c and 6.16a). Figures 6.66b, c show distributed microporosity of less than 2%
exists for the lower value of ill,er and zero microporosity is predicted for f.,er = 1.0.
The predicted total solidification time for alloy is 920 s which is substantially longer
than for the pure aluminium, 473 s. The experimental time is 870 s. Figure 6.67
shows calculated and measured cooling curves in the metal at the T-junction and
in the middle of the horizontal section for R a / b = 1.0 (simulation 85.16). Predicted
126
cooling curves are similar to the experimental ones. The maximum deviations be-
tween the computed and measured temperatures are at the end of the simulation,
t = 1000 s, and are 35°C in the T-junction and 15°C in the horizontal section.
B. R a/ b =1.67. Four values of is ,er=0.67, 0.8, 0.9 and 1.0, were tested in the
simulation of the casting with R a / b = 1.67. Figure 6.68 shows predicted porosity
distributions. A macro-cavity in the T-junction can be distinguished only for the
first two cases. It is interesting to note that the position of the cavity is the same
in both cases, though the volume is larger for iSler = 0.67 that for is,er = 0.8.
Experimental results show a cavity in the junction located lower than predicted in
these two case (Fig. 6.16b).
The degree of the distributed microporosity varies in the four cases, increasing for
smaller values of iSler, never exceeding experimentally measured porosities of 1 to
2.5%, if Is,er ~ 0.67. For iSl er = 1.0 the casting is completely sound which is never
the case in the experimental castings due to the presence of dissolved gasses in the
melt. Predicted solidification times vary with the value of is,C'r: from ts = 840 s for
iSler = 0.67 to ts = 800.0 s for iSl er = 1.0. The experimental solidification time is
770 s.
127
pouring direclty insulator
into the mould
(a) casting
sand mould
- - - -
.I
-1-':
31=50 mm
32=250 mm
33=100 mm
34=50 mm
4
(b) d2=2 mm
d3=9 mm
d4=18.5 mm
5
d5=10 mm
6
7 d6=20.5 mm
d7=29 mm
insulator
sand base
.I
84
-t ~
L" .
300
136
r"' !~
300
150
!
I
Figure 6.4. overal geometry (a) and dimensions (b) of the chill mould
cylindrical casting in mm (castings E2 and E3).
(a) (b)
(a)
sand mould
~7~
f
I'" 15
- I
14.4
7
(b)
1 ,... 22 ., 1
f
7
i
I'" 15 ~I
Figure 6.6. The overall g e ome try (a) and dime nsions in mm (b) of the
boot-shaped casting (E4).
Figure 6.7. Shrinkage cavity in aluminium-sand
casting E4.
5
a
7
----......
c
10
(a) (b)
Figure 6.8. (a) Geometry and dimensions and (b) pouring method
for the T-shaped casting.
Figure 6.9. General view of the aluminium-sand
T-shaped casting ES.
(a) (b) (c)
,,
EI2 15 1013 1.0 6.8
AI-4.5%Cu aUoy
(1))
lOmm
(b)
Figure 6.14. The chill block (a) and the insulating pad (b) experimental
setup for the T-shaped casting.
(a) (b)
100.0 nun
(a)
0.0
- _:
vrnax=O.16 mrn/s
0.0 fraction of solid 1.0
100.0 nun
(b)
0.0
pressure. velocity,
mrn/s
(a)
-10,000 - 0.43
time, s 50.0
60
iteration
count
(b)
0
8,700 time. s 50.0
fluid
volume.
(c)
mm3
8,140
9.0 time. s 50.0
mean
kinetic
energy.
0 .0
time. s 50.0
Figure 6.21. (a ) Press u re and verti cal vel o c i t y ( red line ) h istories
i n the middle of s e ct ion A in f i g ure 6.2 0 ; (b) iterat ion
cou nt , (c ) metal v olume a n d (d ) mean k inetic ene rgy evo-
l u tions f o r simulation 81. 1 .
0.0 porosity. mm3 40.0
fluid fraction
(a)
0.68
90.0 mm 933.0
temperature
oK
(b)
0.0 900.0
0.0 9O.0mm
375.0
temperature
(c)
300.0
insulator
metal
temperature,
casting axis
sand mould
0 .0 293.0
0.0 100.0 mm
... ,,'
............
,.- .. ' , . '
",
.,'
(a)
1--- tine step size -+- .aalCPU - - CPU per tine Slep
Calculation Accuracy I
sinJaions S21·S24, S211·S214
Ci\.....
[~ ~ ~+
" "'::::~_.:::_D-...--. ..------------------
_. . . ... . . .... ... .-.----------~~~~~~~~--~~.-~
..... . . .....
(b)
~ -
--
'000 ,...
Calculation Accuracy II
sinJaions S21·S24
,~------------------------------------------------------rl.M
§ 0.1
! ~ ..
(c) Ii
&i • ~ ..
..................·8........
--
·:t .. ·.............. ·.... ·........ ·...... ·........ ·· .......... ·.... ·.......................... ..
...
... '000
.aal rurDer a eels
lOGO
- •
-
I comergence erra - rna! volrne erra - - average velclcty
temperature,
oK
(a) (b)
900.0
(c) (d)
293.0
x-vc)ocity,
mmls
(a)
0.0 - - - - - - - - - - - - - - - - - "-----~
0.01
mm/s
(b)
Figure 6.28. Simulation 52.2: horizontal (a) and vertical (b) velocities
in the centre of the casting.
0.3
(a)
.6.t, S
0.1
35 time, S 1,300
(b) N
o
18,800 1.300
17,560
1.5 time, S 1,300
(d) Em,
10- 10 J/kg
O.O~ _ _ _ _ _ _ _ _ _ _ _ _ _--..I.._ _ _ _ __
time, S 1,300
Figure 6,29. Simulation S2. 2: (a) time step size, (b) iteration
count, (c) metal volume and (c) metal mean kinetic
energy.
869.6
temperature,
oK
865.4
(a) (b)
&40
. . ·. . ·. ·+·......·. ·. ·T···· .... · · · t~· . · ·~ ·: L.·..::::::::::J!::::::::::::::·.:·i!:.·.....·..·. . .
..............·i···············i.. ··········· ··1- .. ..... .
·.· .·.--. ,·.1·. ·..·.·.·.· ....·:·.·.·..·.·,-.'.·..... ·.·..... ..............
of
i ii : )( : : ! !
ii L..............L...............~...... ......
600 .......•.... . . . J.. .......... . .. j - 1---··· •••• · • • •• 1· •• ·•·•••· • • • •• 4- • • •••• •• • • • • • •
: : : : )( i. '.; .'
t 660 ...............1. .............. 1...............1. ..............1...............;....
: : : : :
..... -~.· .......·.··..!·.··.· .........t ..............
i !
(a) 620
i : : i i : ..... i................ ~:. . . ......... ..-....
480
Simulation SZ5
th'c 2. 3 rn::I 4
6501,-----~------~----~~--~------~----~----~----~~----.
(b) t
1000
1i'ne, S
Simulation S2.5
th'c 5, 6lrtd 7
650
:· ilr ; ~ r~r] : ·~ l~ ·;
600
450
U 350 ···············1·················
of
300 ............... i .... ·
ii : )(
(c)
t 250
200
1000
lire, s
t
e
rn
p
e
r
a
(a) (b) (c)
t
u
r
e
°c
V max=2 mm/s
0.0 34
0.0 75mm
(d) (e)
Figure 6.32. Pure iron solidification simulation S3.1: (a) The enmeshed
domain. Temperature distributions and cavity development at
(b) t =50 Sf (c) t =lOO Sf (d) t =150 S and (e) t = 200 s.
Simulation S3.1
maximum velocity components
20
!.
--EE
III
12
i-
"0 \, 1 ~ 1 l l ~
.•• •,. .•• .••.•• ••••• f' .•• •• •• •• ••.•• •• f- • •••••• ••• •••• t- . ... . . . . . ... ...... . . ....................... . . .
0
Qj
e
>
. 1" "'--1. I I I I
""""""",' ... .......,.................... "~>r~·'·~·:·:T:·:·~·~·~·~·~l·~·~·~·~·~·~·L~······· "··
o'+-----~--_+-----r----+---~-----+----~----~~~ -'.
o 20 ~ 00 00 100 140 160 180
time, s
fluid
fraction
0.5
solid
fraction
0.0
fluid
fraction
0.5
Figure 6.36. Predicted cavity shapes and porosity for alloy steel
simulations S3. 5- 7 : (a) fscr= O. 67, (b) fs cr =1. 0 (lever
rule) and (c) fscr=l.O (linear latent heat release ) .
There is a 1% microporosity in case (a).
75.0 mm
(a)
0.0
0.0 300.0 mm
140.0 mm
(b)
0.0
0.0 75.0 mm
140.0 mm
(c)
0 .0
0.0 300.0 mm
(f!)
...... 4,
'.
\\
,,
\
(a) (b)
, ,
,
~-:r::<"=~..
- ~\
,
,
'.
\\
__11__~.
--- -- -
'.
\\
- -- \ ................. '.\
(c) (d)
.1t. s
(a)
iteration
count
(b)
0.0
5.93 610.0
metal
volume.
(c)
105 mm3
5.55 L _ -_ _- _ -_ _- ___--==:;::::=_.
time. s 610.0
6.2
mean
(d) kinetic
energy.
10-7 J/kg
time. s 610.0
Figure 6.40. Boot-shaped casting imulation S4.1: (a) time step size,
(b) iteration count, (c) metal volume and (d) metal mea n
kinetic energy.
• :.,~ •• ..: "'... ~ '. : l , J :-. t t f.. ~f ~ ~' ii:\ 'Ii ~: 1";-:: ~ ~ a ;.. \\) ~ ~ ~ .;.' ~; " .r i ~ .:,t i'\ :"," ""-1":' ':l.::s~~;:·\j;!;.al :)'~ .• ;. ;,;it~·_.:.1~1t!,..'·~:.\J:....·..;1t.~~l':.O,:,~\,{ \ .... t{t.'io~
"·i~I·"\""\!Ii}·'t!~ i"tl?it{\o:~ft?l=~\:)-'-~~~"l.ijf~i<~t'~''''\:· I";I,J,N~ ........:,l: ~:,:: (l ~H , ,~ if t 1\\3:,';'. ~ (.;.:~. ~i:~ '.... ':~:~ )'J~';' .t" ta •..::c. ~ Ii:, tl ;..'.1 ....'.-
.... ~ ~ .. ~ ~; II :\ '. ; lJ:'I'iJ:Cl.t:w .,.~~~~c:: ; ~ .~':JIJ.' I;\. ~ a I.:-l. " .-""'Q; ~ •. -~"! ~ '!: ! ~!; !;,! ~ ~ ~ ~ ~::~:~: ~~ :; ~: ~~ .~ ~ ~ :' : ~ ;~ ~ : ~ ~
.......... ,oJQ " ' ~ " ~~.,;"; ... ~. . . ,. C' • . , ... N>' ...
.,:.Ct ..."~ . .1 "'~9,",,"I"\.:o.-,,,~
...
.... ~.
;W\-_ ....... 'r . . ,....................., .. _,"\."" .. .. ....... ,.. ,.. ... C'. .; ~ ~ ~ ~~ ~,.~
. . . . O.J ... ..:taWO."-=t.< .. n. . .... a; •
:==
-l, J .'" ~4l"f' . lClC .. " ..... , ..:~ ,_~.&o, ~ ~ ~ ~ I~ . . . . _ _ •
................-----_. __...... . . . -
...... . •••• I • • • • • • • • • • _ _ _ _ _ _ • • '!'I<, ... ... "' ..._ ... ............... r
"~.·r'l"
•••• • •••• _ _ _ _ _ _ _ _ _ _ _ _ _ _ ..
'\I . . "
• ~ • •• "
• • ,, ' I . , \ ., ...
- ." , ._
,n . .,
r~.- ','.~ I !~~~ ! ~-:~~::_.~:~~:':!~:':-~~:! :.~~~~:,~~::"~: ~ .. :r -..-..,
.;t. .. _.'" . . . ..
~·"'l
.
_
.
,
~
I"~' J..w~,. ~\I, " " \ : .... '_ .. " ,
"'" ..... " <". •••• t · t . · ..·.: ... ..,.", c.' •• _ _ _ _ _ _ . . .. . . . . . . ..... ' .... ...
.... ".' " 11111: ....
~
t,- _•••••
_.,-_ ""."....
~ _.,.. ~ MoI~'· ~.
.....- •. , .... -of • • • • , • • • • • • • • _ _ _ _ _ _ _ _ _ _ _ _ . - . . .,.
.... ............. •••• •• • • • • • • • • • _ _ _ _ _ _ _ _ _ "' ...... "" r . t ...............'" ,..)00,.''; _..: . . ' • • • • , . . . . . . . . ,., . . . . . . . . . . . ··uo "'a- I ~._
.. ~ .-c.: _ ·c .. _•••• I ~ " ~ o;4'lU' r .\,0 'II.' 0 • ..,.,- oc. 11.... • ~~.. ,"'II , - ,
,",,' .... ,........... . . . . . . . . . . . . . _ • • __ • • ___ . ......... ~.' .. I . . ., ...........
C""'I:''Clo''".;. • • • • • I '.'~ ..:e fll'l~c 'l; 0','\' 11' • • • • l"'\\(,JL II:, I~u'''''n .
. _~ ... .. \...- •••• • • • • • • • • • • • _ _ _ _ • __ . ~ ................. I ' ~ ,. . , ..~ .... .
•~ " -'"\ , ... . . . . . . l' \ - ;".' . f ~ ...... ,,: '-' ..... .... ,• . D· .'" '?P. r, ~ lil: ... OJ • ,;\ "I I ....... ,.,.. n ...·
~ .•• Q,._.,.;.,., 1 • • ~ .II',:..:. r~~'!:'!::"'C~~~~ _.,. ..~._."" .. ;. •• oI.\. ,,.;..,), ... ;1... ...... . , ..... , \ ... ' " \ 1 "Ie j)I " ...... ~ ... "'" _ tr .t" W " ., . ...... ., III" I .... -.. ,... ._
'.'1."
..1 ... ..:.~ ~ ~- , .
..... :;:,;,'0, 'j
~ .. J
;:,~;\: , I ,~~~
eC"",c:.,..::t .. ~ , .. .do.:;lrp,....ol c.1f\- \11.... .~ ...
1i;.t: .t.';:fr:;.1..-,;. S\~~
I! ~ J •.• 01 1' \
~',;.!' t'.!:I.'
l':\r'i )o.\ .....
, ... U
::t~ e-'
'i,r>o
Q.O;; """',.,..~
~ ~,' ~
.. "
Q .
. , t
,·" t \to\'
• ~"" ..... "" , "" ....... w: -". __ ...... ,', .......
l a; (I eo,
\". I>'''~ 'l." ., .
,,, . . . . . . . '. . . .
~
(c)
(d )
Figure 6.41. The enmeshed doma i n for the T-shaped casting simulat i on :
(a ), (b ). x - y plane views ; (c ) x - z plane , (d) a 3 - D view
of t he mo uld cavity and the intial position of the metal
highlighted with red colour .
average
R=a/b choke diameter experimental simulated
(mm) time (sec) time (sec)
(b) (e)
(c) (0
(a) (b)
(c)
mean
(a) kinetic
energy.
J/kg
0.0 L_--o---------===:::::::::::::=O::::==""'-_
0.0 8.0
time. s
2.0
metal
volume.
(b)
0.0.J---+---+--+---+----+----------+----+----<
0.0 8.0
time. s
Figure 6.45. Simulation S5.1: (a) mean metal kinetic energy and
(b) metal volume evolution during filling.
temperature,
vmax=O.14 m/s
oK
(a)
981.0
972.0
963.0
949
946
954.0
vmax=0.07 rn/s
945 .0
936.0
927.0
(b)
(a)
temperature.
oK
vmax=1.9 mm/s
(b)
840.0
porosity.
(c) mm 3
0.0
p
o
r
o
s
(a) (b)
t
0 .0
(c) (d)
750
.............
............. ,.. .. .........
()
...............
...<i
a...
'"
Q)
Co
E
(a) ~
710
700
6
2 3 4 6 6 7 B 9 10 11 12 13 14
Ume,s
Simulation S5.1
temperature histories
7
720
6BO
----~ ...... ...... ---
() 640 ""
(b)
~ "" ......, .....................
600 "
~... ""
Q)
Co
560
'I., .................
E "'"
~ 520
480
" '. ' ..... '"
'" ................
440 ............... -.- ....
400
a 50 100 150 200 250 300 350 400 450 600
tlme, S
p
o
r
o
s
t
y
mm3
0.0
170.0
0 .0
(f)
(e)
Figure 6.5L Simulat i ons S5. 4- 7. Porosity di s tribu tion : (a) R=O. 6,
superheat 300K~ (b) R=O.6 , s uperheat SOa K , (c ), (e) R=1.67 ,
superheat 30 0 K and (d) , (f) R=1 . 67, superheat SO OK.
905.5
t
e
m
p
e
r
a
t
u
r
(a) ( b) e
oK
898:0
(c) (d)
metal
volume.
(a) mm3
17.560 ~------------~-------
time. s 1.300.0
7.5
total
metal
(b) enthalpy.
104 J
4 .4 +-----------------_
time. s 1.300.0
F~re 6~. (a) Metal volume and (b) total enthalpy evolutions
for cylindrical casting simulations 82.2 (full model,
black) and 82.12 (simplified model, red).
Calculation Efficiency
simulations S2.1-52.4. S2.11-52.14
~.=-----------~--------------------------~~10
total CPU: M1
-....
CPU/tlme step : M1
total CPU: M2
(I)
-.:- .
<Ii
CPU/tlme step: M2 ~
:::I
CI..
U
~
01
2
O!+-~~~====~----.------r-----r-----r-----+o~
o 500 1(XX) 1500 2OCO 2500 3CXXl 3500
total number of cells
Figure 6.55. Total CPU time and CPU time per time step as
functions of the uniform mesh cell number for
the full shrinkage model HI (simulations 52.1-
52.4) and for the simplified shrinkage model
H2 (simulations 52.11-52.14).
1536
t
e
m
p
e
(a) (b) r
,.
r :\'.
~~
\\
a
t
u
t
!:
"
i
."
r
t " e
~ .,,
.,
'.~ °C
!
".
:~'
(,.
(c) (d)
f f
r
u a
i c
d
o
n
0.5
(a) (b )
F~re 657. Predi cted cavity shapes and porosity for alloy steel
simulations S3.15 and S3 . 16 using the M2 model : (a )
fs cr = O. 67 , (b) fs cr =1.0 . The r e is a 1% microporosity
in case (a).
1.48
metal
volume,
10-4 m3
1.42
time,s 210.0
Figure 6.58. Metal volume evolut ion: pure iron (black) and alloy
steel (red) simulated by the simplified model.
890.0 remperarure, OJ( 933.0
~t '
~t~~'~l,
(a) (b)
(c) (d)
(e)
(0
metal
volume,
(a)
10-5 m3
5.53
time, S 610.0
2.3
total
metal
enthalpy,
(b) 106 J
1.43
time, S 610.0
Figure 6.60. (a) Metal volume and (b) total enthalpy evolutions
for boot - shaped casting simulations S4.1 (full model,
black) and S4.11 (simplified model, red).
t
e
m
p
e
r
a
t
u
r
e
1.0
720
f f
1 r
u a
c
d
o
n
0.5
(c) (d)
(e)
Figure 6.62. Simulation S5. 11: 3-D plots of the overall fluid
(black) and liquid (red) regions. (a) t=lOO s, (b)
t=200 s, (c) t=300 s, (d) t=400 s and (e) t=500 s.
0.5 fluid fraction 1.0
(a)
(b)
lOmm lOmm
(c) (d)
Figure 6.64. (a) The chill block in simulation S5.12, (b) insulator i l
in simulation S5.13, (c) insulator i2 in simulation S5.14
and (d) insulator i3 in simulation S5.15 for the T-shaped
c asting .
- 880
e
m
p
e
<>J<.
300
1.0
u a
c
d
o
n
0.5
(c)
f f
r
u a
i c
d
o
n
(c)
0.5
700 .........................................................................................................
525 ......................................................................................... .
500 ....................................................................................................... .
475 .. ........................................................... .
450+---~----~---.----.----.----.---~----,----,--~
o 1 00 200 300 400 500 600 700 800 900 1000
time. s
(a) (b)
f f
I r
u a
c
d
(c)
o
(d)
D
0.87
Discussion
Simulations of filling and feeding flows in a variety of castings demonstrated the ad-
equacy of the use of the Navier-Stokes model of an incompressible fluid. Predicted
filling times for the T-shaped castings matched the experimental results with an
average error of only 7%. In the full model volumetric changes due to phase trans-
formation were described by introducing a volumetric source term in the continuity
equation at the liquid/solid interface. Simulations have predicted shrinkage volume
to 1% of the theoretical value defined by the values of the solid and liquid phase
densities, Eq. (6.1).
The model did not include gas evolution in the liquid metal. Experiments showed
that this assumption is valid for thoroughly degassed melts. A significant gas content
led to a substantial micro-porosity in the solidified casting which was not predicted
by the shrinkage model.
The Reynolds number in the predicted feeding flows varied in the range between
0.0 and 100, indicating the importance of viscous stresses in the balance of forces
acting on the luquid metal. Predicted feeding flow velocities are of the order of
0.01- 0.1 mm/ s for sand castings and 0.1-1.0 mm/ s for chill mould castings. This
is in agreement with experimental results and other computer simulations. Tsai [133]
showed that feeding flow velocities can be of the same order of magnitude as those
occurring during metal thermal and solutal convection and predicted velocities of
around 0.1 - 0.7 mm/s for a 1% Cr steel. Maples and Poirirer simulated solutal
convection in a AI-4.5%Cu alloy with velocities of up to 0.03 mm/ s [132].
128
Natural convection effects were not included in simulations to simplify the analysis.
Thermal and solutal convection in feeders would increase heat exchange in the ver-
tical direction resulting, in general, in a deeper primary cavity. The depth of the
latter was largely underestimated in simulations of the aluminum-sand cylindrical
casting and the T-shaped casting where a large feeder was present. generally, the
primary cavity depth in a feeder is defined by the ratio of the heat flux qh through
the sides and the top of the feeder to the heat flux qll at the feeder's bottom. The
smaller the ratio qh/911 is the smaller is the depth of the cavity. In the extreme case
of % = 0 (a fully insulated feeder) the surface of the solidified metal in the feeder is
flat and aligned with the horiztontal plane.
The error in predicted depth of the primary cavity was more significant in alloy
T-shaped casting simulations (Figs. 6.17b and 6.67). This indicates that the largest
source of the inaccuracy in alloy modelling was the solidification model. For the
AI-4.5%Cu alloy the solute is heavier than the solvent so that liquid with smaller
solute content would accumulate at the top of the feeder leading to an increase in
the local values of T/ and Ts thus accelerating the solidification process at the feeder
top. The latter would increase the effective value of qh leading to a deeper cavity.
The employed alloy solidification model did not include segregation and diffusion of
the alloy components. Instead, a simplified approach was taken by using the lever
rule. Further errors were introduced by the numerical implementation of the lever
rule (section 2.8).
Free surface was treated as an adiabatic boundary neglecting radiative and convec-
tive heat losses into atmosphere. To study the effect of these losses on the primary
cavity depth a test simulation was carried out for the cylindrical casting shown in
Fig 6.3 in which the insulator was replaced by sand. (simulation 82.45 with the
non-uniform mesh shown in Fig. 6.30a). This modification was introduced to ap-
proximate an additional heat loss due to radiation since the latter could not be
modelled directly. The final cavity shape is shown in Fig. 7.1 and its depth differs
from the experimental result by only 5% compared to the 21 % error when the insula-
tor was used (Fig. 6.30b). This result indicates that an inclusion of a non-adiabatic
boundary condition and free surface may imporove the accuracy of the shrinkage
model.
In all simulations material properties of metal solid and liquid phases and those
of mould were assumed to be temperature independent. This assumption might
be a source of additional errors in temperature predictions shown Figs. 6.31, 6.50
and 6.67. It is well known that thermal properties of sand are especially sensitive
to temperature variations [134]. However, small differences between predicted and
129
measured temperatures in Fig. 6.31, less than 10% in average, indicate that these
errors are insignificant.
Metal/mould heat transfer coefficient in sand castings was defined assuming good
contact between the two materials. This worked well while the metal was liquid as
can be seen in Fig. 6.50 where temperatures at three locations for pure metal T-
shaped casting are shown. After metal solidified in the T-junction and in the sprue,
the difference between predicted and measured temperatures increase substantially.
The measured cooling rates are smaller than the predicted suggesting that heat
transfer coefficients at these locations decreased as metal solidified, perhaps, due to
the formation of air gaps between the solidified metal and the mould. Predicted
temperature at the bottom of the feeder stays close to the measured one until the
end of the simulation. This indicates that no gap formed at this location, possibly
due to the shear weight of the feeder. In the case of the alloy, differences in the
measured and predicted temperatures are smaller and do not change abruptly after
the metal solidified (Fig. 6.67). As in the case of the pure metal, it can be seen
that predicted cooling rates at these locations are higher than the measures ones
indicating that smaller air gaps, if any, have formed in the alloy casting, perhaps,
due to the semisolid nature of the mushy zone during solidification.
Observations of the experimental T-shaped castings showed that the shape and
volume of the internal cavity in the T-junction might be changed due to deformation
of the solid metal (Fig. 6.11b). The deformation occurs because ofthe low pressure in
the trapped liquid metal as it solidifies and shrinks. This deformation led to tearing
of the solid metal and penetration of air into the metal through the permiable mould
resulting in a cavity with charateristically smooth edges (Fig. 6.13b). These effects
could not be predicted by the developed models since neither solid phase deformation
nor air flow through the mould walls were included in the modelling. However, in
most of the experimental castings these effects were small or were not present at all.
The main part of the numerical description of liquid metal flow is the VOF method
of tracking free surfaces. It is also employed in the shrinkage models to describe
shrinkage defects in solidifying castings. Comparisons with the experimental resutls
show that the VOF method is an efficient and accurate method for modelling both
filling and shrinkage cavity formation.
Upwind differencing method for the approximation of the advection terms in Navier-
130
Storkes and energy equations, central differencing for the approximation of viscous
terms and the SOR iteration technique for the solution of the coupled pressure-
velocity equations provided sufficiently accuarate solutions for the liquid metal flow.
This can be seen from the comparison of the predicted and measured filling times for
the T-shaped castings (Fig. 6.42). The simplicity of the numerical implementation of
the SOR iteration method allowed the source terms, associated with the volumteric
shrinkage, to be easily incorporated into the solution algorithm. In general, the
development of the shrinkage models on the basis of an existing general purpose
eFn code proved to be highly efficient.
The accuracy of the predictions of shrinkage volume were defined mainly by the value
of the convergence criterion f. However, in the case of the sand mould the depth of
the primary cavity was always underestimated when compared with experimental
results.
As was mentioned in Section 7.1, the primary cavity depth in a feeder is proportional
to the ratio of the heat fluxes qh and qu' Several additional factors may have led to
an underestimation of qh/ qv in the simulations:
1. The interfacial heat flux at the vertical sides of the feeder may be underesti-
mated due to the lack of mesh resolution. It was shown in Section 4.1 that the
flux is dependent on mesh size and the related value of the open volume VF in
the interfacial cells (Figs. 4.3-4.5). That might apply to the T-shaped casting
since the resolution in the casting is better than that in the feeder, as can be
seen in Fig. 6.41. However, this dependence was removed in the cylindrical
sand-mould casting simulations S2.1-S2.4 by employing uniform meshes. In
this case all interfacial cells have the same size, therefore the heat transfer
coefficient, Eq. (4.13), and truncation errors in the metal/mould heat flux
calculation, Eq.(4.14), are the same in these cells. Here, even with the finest
mesh (S2.4) the cavity depth was still 14% smaller than the experimental result
(Figs. 6.25 and 6.26).
2. Finite mesh resolution introduces error in representing the free surface of the
solidifying cell adjacent to the wall. This error affects the shape of the feeder
cavity. Figure 7.2 shows that in solidifying surface cell, the fluid level drops
lower than it would if there was a finer resolution in the horizontal direction
to resolve the curvature of the free surface in the interfacial cell (shaded area).
131
The discrepancy is greater for pure metals since the liquid flow seizes only
when fs = 1 in the cell. If the free surface in the first cell freezes at 8h I lower
than the correct position, then the metal volume 8h I Al is equally distributed
between the rest of the liquid surface cells, where Al is the cell area in the
horizontal plane. As a result the liquid level in these cells is higher than the
correct value by
8h = 8h l Al
Ao
where Ao is the total area of the remaining liquid surface cells in the horizontal
plane. Further error is introduced when the second cell solidifies and so on.
Since the total cavity volume is predicted accurately, the final cell to solidify,
which actually defines the cavity depth, accumulates all the errors. The error
becomes larger for cylindrical feeders where, as the freezing front moves to-
wards the feeder axis, Ao decreases as r decreases, therefore, 8h grows as l/r,
where r is the distance from the axis.
The mesh related factors can be decreased by refining the mesh and using a
uniformly spaced mesh.
3. An error could result from assuming conduction heat transfer between the
centre of liquid interfacial cells and the interface (Eqs. (4.9), (4.10) and (4.22)-
(4.24)). Hoadley et al [131] showed that convective heat transfer resulting from
liquid flow may be significant. The ratio of the heat transfer by convection to
conduction is described by Rayleigh Number [4]
The position and size of the secondary internal cavities were predicted with good
accuracy, especially for pure metals. These are less dependent on the accuracy of
free surface description but an accurate calculation of the interfacial heat fluxes is
also crucial.
132
An important feature of the secondary internal cavities in the T-junction is their
dependence on the filling time as shown in Fig. 6.49. The slowest filling time of
22.5 s ensures a sound casting. This result supports the general idea that the final
casting quality depends on the filling parameters. However, this dependence was
not observed in the experiments where variation in the size of the secondary cavity
for the same filling rate was comparable, or even larger, to the predicted variation
for different filling rates.
A solution to the mesh dependence problem could be the use of analytical functions
to describe the heat flux at the interfaces [66,93,95,96]. However, this remains
restricted by the small number of available analytical solutions.
The deficiencies of the heat transfer algorithm are clearly shown in the comparisons
of the temperature histories (Figs. 6.31 and 6.50). The solidification time for the
aluminium-sand cylindrical casting is about 26% above the experimental value of
754 s for the finest mesh (simulation S2.4). Additional errors are introduced by
the use of constant metal and mould thermal properties and constant interface
thermal resistances. However, these may not be crucial for accurate shrinkage defect
predictions as has been shown for the chill mould casting (simulations S3.1-S3.4 in
Fig. 6.34). These show that it is not the absolute value of the heat flux, but its
variation along the interface that largely defines the position and size of the cavities
in the solidified casting. The latter also indicates the necessity of modelling the
formation of a gap between the mould and the metal. A gap of 0.05 mm can lead to
an increase in the interfacial thermal resistance by up to two orders of magnitudes as
shown by Hou and Pehlke [24]. These authors also demonstrated that temperature
prediction may be improved if an experimentally measured time-dependent heat
transfer coefficient is used in the calculations.
A problem for alloy solidification modelling was posed by the use of the critical
fraction of solid as the feeding criterion. It was shown in simulations S5.18-S5.21 in
Fig. 6.67 that the size and even the occurrence of the internal macro-porosity in the
T-junction was sensitive to the value of f.,cr' These simulations were expected to
133
be less accurate since the AI-4.5%Cu alloy has a wide freezing range and the whole
casting was in a mushy state during most of the solidification time. This condition
magnified the deficiencies of applying the model to alloy solidification. In addition,
the use of the lever and Scheil's model in the computational cells in the form of Eqs.
(2.64) and (2.65) was not consistent with the assumptions of the two solidification
models. -
A more accurate description of feeding in the mushy zone is the use of the critical
value of the solid fraction gradient, \l is,er [112]. It takes into account that feeding
is easier along the direction of the dendrite growth than across it. However, the
value of \l is,er may vary between alloys, mould geometries and even cooling rates.
Clearly, a more sophisticated, detailed model of alloy solidification and feeding is
required.
As was shown in Section 7.1, predicted velocity magnitudes in the feeding a flow
are consistent with existing numerical and experimental data. The magnitude of
the average feeding flow velocities are defined mainly by the solidification time and
the total volumetric shrinkage (or the ratio PI! Ps). Since the latter is fixed for a
particular metal, the accuracy of the velocity predictions depends on the accuracy of
the solidification rate prediction. The mould geometry also plays a part in defining
the velocities. For the chill mould casting S3.1, for example, velocities went up to
10 mml s for a short period during solidification when feeding occurred through a
narrow passage in the metal (Fig. 6.33).
134
7.2.2 Drag Force
The drag force in the form of Eqs. (3.10) and (4.51) or (4.52) does not describe a
real physical force but arises from the fact that numerical equations are applied to
cells which contain a solid-liquid mixture and the solid phase does not_flow. Strictly
speaking this approach is valid only for pure metal solidification since for alloys
the solid phase crystals can float at low solid fraction values, and furthermore, the
frictional drag force in the mushy zone is also important. The latter is not taken
into account by the model to avoid convergence problems in the high drag regions.
The use of the drag force in calculating flow losses due to solidification is not satis-
factory also for other reasons:
• The drag coefficient, Eq. (4.51), does not ensure that the solid phase remains
frozen in space. For example, I< = 0 when dls/dt = 0 even if Is =f:. O.
The enthalpy advection algorithm is closely related to maintaining a zero solid
phase velocity (see Section 7.2.2). In the M1 model only the liquid enthalpy
content is advected. No account is taken, however, of the amount of liquid
phase left in a cell. Thus, more liquid enthalpy can be advected out of a cell
than it may provide!.
• Pressure, gravity and viscous forces in the momentum Eqs. (3.7)-(3.9) are
applied to the full cell volume while they should affect only the liquid phase
and their influence on the cell momentum should be reduced accordingly. The
effect of the drag force is thus diminished, especially when the Is value is close
to unity.
• There may be situations where the drag coefficient in nearly solid cells becomes
very high leading to excessive iterating (see Section 7.2.3).
• Finally, even for pure metals there may be significant flow losses due to shear
srtesses at liquid/solid interface in the regions of high values of the solid frac-
tion. The expression for the drag force could be modified to include the viscous
friction effects but it is not clear how to do it efficiently in a general way in
1 A test could be introduced for the amount of enthalpy to be advected through a cell face,
6Hadv, in surface cells:
6Hadv = min(u . A· H, . ~t, H, . Vj)
where u is the velocity, A the face area, H, the liquid phase enthalpy and Vj the volume of the
liquid left in the cell. It is inappropriate though for a full, internal cell due to the fluid continuity
and incompressibility conditions. These imply that 6Ha dv = u . A . H, . at must be advected
between neighbouring full cells.
135
three dimensions, even if the solidification front is assumed to be a plane in-
terface in each cell.
The situation is even more complicated for alloy solidification where the solidification
front cannot generally be represented by a surface.
A way to remove the deficiencies of the employed drag force formulation may be
found in extending the FAVOR method to take into account cell blockages by solid
phase in a manner similarly to that for the mould. This can be accomplished by the
following transformation of the fractional open volume and face areas:
-A 2 (1 - Fa fav,a)(l - Fb fav,b)
A ab,m - a.b (7.1 )
2 - Fa. fav,a. - Fb fsv,b
where Aab is the open area fraction at the common face of neighbouring cells a and b,
fsv is the volumetric solid fraction function and index m means a modified variable.
Similarly, a 'modified' fluid fraction Fm can be introduced
(7.2)
to describe the remaining liquid in the cell which can flow. The modified variables
VF,m, Aa.b,m and Fm have the same meaning as the original variables but relate to
the liquid phase only while the solid phase is treated as a part of the solid structure
of the mould.
Eqs. (7.1) and (7.2) should be applied to advection terms. The original variables
must still be used for heat conduction and transfer, and the enthalpy still refers to the
full, liquid+solid, fluid volume in a cell. In other words, a formal implementation of
this method to describe the effect of the solidification on the fluid flow would require
a combined use of the original and modified fluid, volume and area fractions.
• There is no need for the atrificial drag force in the form of Eq. (4.52) because
cell velocity in the modified method is not the mixture velocity U m but the
actual liquid phase velocity U/ (Fig. 4.8a). Therefore, possible convergence
problems in the high drag regions are removed. The momentum equations
become more consistent since all the terms refer to the liquid phase only.
• Since velocities and enthalpy advection terms are calculated with respect to
the amount of the liquid phase left in the cells, the solid phase is ensured to
be frozen in space.
136
• For pure material solidification, when the solid/liquid interface can be repre-
sented by a surface, stresses between the two phases can be calculated by the
same algorithm as for the shear stresses between the liquid and the mould
walls. Therefore, there is no need for a frictional drag force. For alloys an
additional drag force, accounting for flow losses in the mushy zone, should be
developed and introduced to momentum equations.
There will be implications for the time step size limitation since the latter depends
on the values of VF and A (Eq. (3.14)). However, the time step can only decrease
by half, as follows from Eq. (7.1):
A disadvantage of Eqs. (7.1) and (7.2) is the assumption that the solid phase it
attached to the mould walls and cannot move.
The fluctuations arise from the fact that the melting temperature Tm in a solidifying
cell is assigned to the cell centre instead of the solidification front in the cell, thus
distorting temperature gradients in the vicinity of the cell. This error may be
minimised by employing a linear interpolation (not extrapolation) of temperatures
between the solidification front and the centre of a neighbouring cell to obtain the
cell-centred va.lue of the temperature in the cell that contains the front (Fig. 7.3).
This is akin to setting pressures in the surface cells using Eq. (2.49). However, it
may require substantial numerical effort since the orientation of the solidification
front in the cell must be defined in three dimensions.
The use of such interpolation will alter the unique correspondence of the cell-centred
temperature to the cell enthalpy, thus complicating the procedure of defining the
137
cell-centred temperature from the value of the enthalpy. This difficulty could be
overcome by treating the two phases as separate fluids so that two temperatures and
enthalpies can be used in one cell. Substantial programming effort and computer
memory would then be required to introduce this modification as one cell could then
contain four materials: liquid and solid phases, mould and void.
The advection of only the liquid phase is ensured by using the enthalpy equation in
the form of Eq. (5.20). An alternative form of this equation is
aF-
PI . - H = -PI (~
- PI Hm Vn da + q'E ) - -1 -ap PI F H (7.3)
at P Va 'E P at
in which the mixture cell-centred enthalpy Hm is present in the advection terms
rather than the liquid phase enthalpy HI· Eq. (7.3) was used to simulate the pure
iron casting (S3.21) with the result shown in Fig. 7.5. It can be seen that the solid
phase is 'stripped' by the metal flow from the mould walls and carried into the bulk
of the casting which is clearly is a non-physical result.
7.2.4 Convergence
Several convergence 'failures' have occurred cylindrical casting simulations using the
Ml model, S3.1-S3.4, as shown on a diagnostics print-out in Fig. 7.6 2 • A general
increase in iterations per cycle began at t ~ 30 s (lines> 25). Most of the iteration
failures occurred when liquid in the upper half of the casting froze to leave a narrow
passage through which the rest of the liquid was fed (Fig. 6.32). Factors leading to
an increase in iteration can be summarised as follows:
1. Velocities in the liquid cells along the passage are two orders of magnitude
larger than those in the liquid bulk due to the incompressibility of the liquid (rv
10 mm/s compared to rv 0.1 mm/s, the difference magnified by the cylindrical
geometry). Therefore, the value of f may appear to be too small for these cells
to achieve reasonable convergence rate (see also Section 7.2.4).
2. A liquid cell may change from 'internal' cell to 'surface' cell accompanied by
an abrupt change of pressure in it. Consider a partially solid cell A in Fig.
7.7a. The neighbour cell B is fully liquid and cell C is completely solid. As
solidification proceeds the fluid level in cell B goes down, so does its pressure.
2 As mentioned in the footnote in Section 5.2.3, 'iteration failures' are diagnosed when the iter-
ation number in a cycle reaches a predefined maximum. These failures mean that the convergence
is slow, not that the solution does not converge.
138
By the time cell B is nearly empty at time step n pressure in it is, according
to Eq. (2.49),
(7.4)
where po is the void pressure, Gz gravity in the z-direction and ~z the cell
size in the vertical direction. Pressure in cell A is found iteratively even if
there is little fluid left in it since cell A does not have empty neighbours and is
treated as an incompressible internal cell. Its pressure then should not exceed
PH because cell B feeds cell A:
PnA <
-
pnB (7.5)
Suppose cell B is emptied during the next time cycle n + 1 (Fig. 7.7b). Then
cell A becomes a surface cell with vertically orientated free surface since cell
B is the only empty neighbour (see Section 3.3.4) and
n+l _ pn+l _ p
PA - B - 0
(7.6)
139
of the simulation S4.21 are shown in Fig. 7.8. The mean kinetic energy plot
has fewer spikes and the maximum value is 8.5 X 10- 8 J / kg which almost ten
times smaller than that in the simulation S4.1. Convergence rate also improved
with the average iteration number of around 60. A possible explanation to the
difference in the results is that when Eq. (7.3) was used then the metal in cell
A in Fig. 7.8 was more likely to be fully solid by the time cell13 was empty.
In this case the drag force in cell A was effectively infinite and no velocity
fluctuations occurred.
3. A similar effect on the fluid flow is produced when an 'internal cavity' cell,
which contains both luquid and void and has a fixed pressure, P = Per, becomes
fully empty (cell A in Fig. 7.9) during time step n + 1. According to Section
5.2.3, when cell A turns from an 'internal cavity' cell into an empty cell its
pressure does not change
nA
P = PAn+l = Per
At t = t n pressure in cell B is
PB ~ Per + PI . Gz . !:l.z
At t = t n +1 cell B becomes a surface cell and its pressure is calculated as
If FB+1 ~ 1.0, which is usually the case, then PB decreases by a value close
to 0.5p I G z !:l.z in one time step. This generates a sudden upward flow which
can put some fluid back into cell A. The latter will cause a reverse change
of pressure in cell B. These fluctuations can proceed for several time steps
requiring additional iteration effort. This can be seen in Fig. 7.6 (lines 42-52).
Numerical problems described in items 3 and 4 are the main causes of slow conver-
gence rate during shrinkage modelling using the Ml model.
To test the influence of the convergence criterion value, €, on the accuracy of the
solution and to define its optimum value, a series of simulations was carried out us-
ing the full shrinkage model for the cylindrical casting 82.2 in which the insulator is
removed so that a free surface is present throughout the simulation (82.21-82.26)3.
3 As was described in Section 2.3.2, at each time step pressure iterations are considered converged
if the modulus of the velocity divergence, divv, is smaller than ( in every mesh cell.
140
The value of f was varied between 0.0001 8- 1 and 1.08- 1 . Fig. 7.10 shows the volu-
metric error and the average iteration number as functions of f. The optimum value
of the convergence criterion in these simulations appears to be close to 0.001 8- 1
when the volumetric error is less than 0.23% of the total volume and the average
iteration number is only 7. Fig. 7.11 shows predicted distributions of the fluid
fraction function. A lack of convergence for large values of f resulted m a 'porosity'
region along the casting axis. A reasonable result, both in terms of the volumetric
error « 0.7%) and the free surface shape, wass given with the value of f up to 0.01
S-I.
For comparison, the results of the same simulations but with no shrinkage, i.e. PI =
PII, are shown in Fig. 7.12 (simulations 82.31-82.36). In contrast to the shrinkage
simulations, the volumetric error here is above 65% for f = 0.1 8- 1 , resulting in
a virtually complete emptying of the mould. It can be seen from Fig. 7.10 that
in the absence of the shrinkage induced flow a substantially smaller value of the
convergence criterion is required to obtain a solution which accurate in terms of the
conservation of the fluid volume.
For an internal cell of volume ~V the volumetric error after the convergence has
been reached is
~ Verr = div v . ~ V ~ f· ~V (7.7)
Generally, the sign of divv varies from cell to cell reducing the net volumetric error
in the domain.
For transient flow problems, e.g. mould filling, an approximate rule of thumb for
choosing the value of f is [83]
f = 0.0002· -UL (7.8)
where U is the magnitude of the velocity variation across the domain and L is the
length scale on which this variation occurs. Generally, the value of f also depends
on the time step size and geometrical features of the flow. The magnitude of the
velocity error, DUerr, due to iteration in a cell of size ~x can be estimated as
~x
DUerr ~ 6.X· f = 0.0002· L U (7.9)
or, since 6.x / L is approximately equal to the inverse of the cell number N across
the domain,
U
DUerr ~ 0.0002· N (7.10)
If for example N = 10 than the velocity error at each time step is within 0.002% of
the overall velocity variation in the domain.
141
If Eq. (7.8) is applied to the cylindrical casting simulation S2.21, where U ~
0.05 mm/s and L = 150 mm (the initial depth of the melt), then f = 7.0.10- 8 s-l.
This value is much too small for an efficient calculation as can be seen from Fig.
7.10: for f < 0.01 S-1 the gain in accuracy in terms of the volumetric error, as
compared to the results for f = 0.1 s-\ is negligible while the average iteration
number increases steeply. The reason for this is that it is not the velocity variation
with time that defines the total iteration. Feeding velocities are mainly defined by
the solidification rate and the latter varies very slowly. Pressures in the casting,
however, do vary substantially as the free surface moves downwards to reduce the
effective pressure head. For example, for casting S2.45 the pressure head changed by
h = 50 mm (h is the predicted depth of the shrinkage cavity) during solidification
and the corresponding change in pressure, !:1p, is
Ue !!
J¥
= -
PI
p
= 0.7 m/s (7.11)
Substituting Ue!! into Eq. (7.8) gives f ~ 0.0009 s-1 which is more in line with the
results of the €-test in Figs. 7.10 and 7.11.
where !:1p is the characteristic pressure change, the convergence criterion can be
estimated as
f = 0.0002 2-L -) /lp
PI
= 0.00022-L . ,/kG (7.12)
For simulation S3.1 with h = 0.25 m and L = 0.6 m, Eq. (7.12) gives f =
0.00053 S-1. This is close to the actual value used in the simulation, € = 0.001 S-1.
Eq. (7.12) also explains the increase in the iteration count when pressure fluctua-
tions occurred at the free surface in this simulation, as described in Section 7.2.3. For
example, the instantaneous pressure change given by Eq. (7.6) with !:1z = 12 mm
corresponds to Ue / / ~ 0.25 m/ s requiring the value of € = 0.0042 s-1 which is 4
times larger than the employed value.
142
7.2.6 Time Step Size Limitations
Whilst the heat transfer and the free surface time step size limitations dominate
the solidification stage, the CFL criterion, Eq. (3.14) controls b.t during the filling
stage. Figs. 7.14a, b show the variation of the time step size during the 7.0 and 22.5
seconds of filling for the T-shaped casting (S5.1 and S5.3), respectively~ The average
velocities during the fillings were 0.6 m/s and 0.2 m/s (the maximum velocity was
at the choke at 2 m / s in both cases), respectively, while the minimum cell size is
5mm.
The time step restriction during filling is a significant factor: for the T -shaped
casting simulation S5.1 the filling time was almost 60 times smaller than the solidi-
fication time (7.5 sand 475 s) but the CPU time for the filling stage simulation was
only 5.2 times smaller than that for the solidification and shrinkage (Fig. 6.19).
The shrinkage model was not used during filling since there was no significant so-
lidification at this stage. When a shrinkage model was employed the increase in the
CPU time per time step due to the additional calculations was less than 5%.
The M1 model is fundamentally more correct than the simplified M2 model as the
latter does not account for fluid flow. The full shrinkage model is also more general
since fluid flow and free surface shape are calculated dynamically without any apriori
assumptions. Besides, the mathematical and numerical formulation of the M1 model
make possible its further development to include, for example, gas evolution.
However, differences in shrinkage cavity predictions of the two models in the casting
cases that were considered in the present work are small. The similarity of the
resutls may be explained, in the first place, by the validity of the assumptions
made in the M2 model listed in Section 5.3, that is that the formation of shrinkage
cavities is mainly governed by gravity. Secondly, redistribution of heat in the melt
was mainly due to metal/mould heat transfer and conduction while advective heat
transfer played a minor role. Therefore neglecting the latter by the M2 model did
not introduce significant errors in shrinkage cavity predictions, cooling rates and
solidification times. Moreover, cavities predicted by the two models for the chill
mould casting are almost identical which may be explained by the fast cooling rate
in the chill mould so that the relative importance of the advection of heat was even
143
smaller than that in the sand castings.
The effect of the absence of the advective heat transfer can be seen in a shallower pri-
mary cavity depth predicted by the simplified model in the sand castings simulations
(e. g. Figs. 6.26 and 6.53).
The differences are greater in the numerical efficiency of the two models. The gain
in speed of the M2 model is first of all due to the use of a larger time step size
which is defined by the heat transfer process (Eq. (2.76)). For the M1 model the
time step is defined by the free surface stability limit {Eq. (2. 77)) during most of
the simulation time. The ratio of the typical time steps sizes used in M2 and M1
models varies from case to case depending on the mesh and properties of the metal
and mould. The maximum ratio was 540 for the chill casting simulations. Fig. 7.13
shows a comparison of time step sizes used by the two models for simulations 85.5
and 85.22.
Finally, since no free surface advection is involved in the M2 model and the liquid
metal free surface is always assumed to be horizontal, the predicted free surface
shape is smoother than those given by the full model as was shown, for example, in
the boot-shaped casting simulations S4.1 and S4.11 (Figs. 6.38 and 6.59).
144
mould
void
- --
metal
f f
r
u a
c
d
o
n
0.5
metal
(8) volume,
10 6 mm 3
1.42
0.0 time, s 210.0
5.3
mean
kinetic
(b)
energy,
10-7 Jlkg
0.0 .l--!.--.!.-......:......::....-.::::~-....::;..,-...::~--- _ __
0.0 time, s 210.0
Figure 7.4. (a) Metal volume and (b) mean kinetic energy
evolutions for simulation 53.1. The curves
are not smooth due to the use of the enthalpy
method.
vmax=4.5 nun/s vmax=2.5 mm/s
1.0
s f
o r
1 a
i c
d
o
n
0.0
(a) (b)
1
t
0.0001+00
-- 0-- 0 ------
7.2901-03/f. 7.2901-03 1.4181+02 1.2801+00 1.0001-03
2 pre ••ure iteretion failure at t- 7.2901-03 cyale- 1 iter-l000 delt- 7.2101-03 DOCOO- 1
3 pre ••ure iteration failure .t t- 3.2081-01 cyale- 44 iter-l000 delt- 7.2101-03 Docon- 2
4 2.0051+00 275 224 7.2901-03/fa 7.2901-03 1.4731+02 2.8021+02 1.0001-'"63
5 4.0101+00 550 36 7.2901-03/fa 7.2901-03 1.4681+02 4.7211+02 1.0001-03
6 pre ••ure iteration failure at t- 5.9711+00 cyal_ 819 iter-l000 delt- 7.2901-03 nocon- 3
7 pre ••ure iteration failure at t- 5.n81+00 cyale- 820 iter-l000 delt- 7.2901-03 Dooon- 4
1 6.0141+00 825 12 7.2I01-03/fa 7.2901-03 1.4&31+02 7.6161+02 1.0001-03
9 8.0UI+00 1100 22 7.2901-03/fa 7.21101-03 1.4621+02 1.6121+02 1.0001-03
10 1.0021+01 1375 28 7.2901l-03/fa 7.21101-03 1.46111+02 9.55111+02 1.0001-03
11 1.2031+01 1650 30 7.2901l-03/fa 7.2901-03 1.45911+02 1.0621+03 1.0001-03
1.4031+01 1125 t1 7.2901l-03/fa 7.2901-03 1.45711+02 1.2431+03 1.0001-03
12
13 1.6041+01 2200 42 7.2II01l-03/fa 7.2901-03 1.45411+02 1.4601+03 1.0001-03
14 pre ••ure iteration failure at t- 1."21+01 cyale- 22 .. iter-lOOO delt- 7.29011-03 nocon- 5
15 pre •• ure iteration failure at t- 1.6731+01 cyal_ 2295 iter-l000 delt- 7.2901-03 DOCOD- 6
pre ••ure iteration failure at t- 1.67411+01 cyale- 2296 iter-l000 delt- 7.:n01-03 Qocon- 7
U
17 pre ••ure iteration failure at t- 1.7161+01 cyale- 2354 iter-lOOO delt- 7.29011-03 Docon- 8
1.8041+01 2475 32 7.2901-03/fa 7.2901-03 1.4521+02 2.07311+03 1.0001-03
18 26., iter-l000 delt-
11 pre •• ure iteration failure at t- 1.9311+01 cyale- 7.2901-03 DOOOQ- 9
20 pre •• ure iteration failure at t- 1.9331+01 cycle· 2651 iter-l000 delt- 7.2101-03 Docon- 10
2.0051+01 2750 18 7.2901l-03/fa 7.29011-03 1.4501+02 2.20111+03 1.0001-03
21 1.45011+02 2.2701+03 1.00011-03
22 2.2051+01 3025 18 7.2901l-03/fa 7.29011-03
23 2.4061+01 3300 28 7.2901-03/fa 7.29011-03 1.44911+02 2.35211+03 1.00011-03
24 2.&061+01 3575 22 7.29011-03/ fa 7.2901-03 1.44811+02 2.4631+03 1.0001-03
25 2.1071+01 3850 19 7.2901-03/fa 7.29011-03 1.44611+02 2.58911+03 1. 00011-03
26 3.0071+01 4125 268 7.2901l-03/fa 7.2901-03 1.44511+02 2.80311+03 1.00011-03
27 3.2081+01 4404 52 7.2901-03/ fa 7.2901-03 1.44311+02 3.0098+03 1.0001-03
28 3.4081+01 4679 53 7.2901l-03/fa 7.29011-03 1.4C2II+02 3.1388+03 1.0001-03
3.6091+01 4954 371 7.2901l-03/t. 7.29011-03 1.4Cll1+02 3.37511+03 1.00011-03
2t cyale- 5061 iter-l000 delt- 7.29011-03
30 pre •• ure iteration failure at t- 3.687B+Ol Docon- 11
pre •• ure iteration failure at t- 3.6881+01 cyale- SOU iter-l000 delt- 7.29011-03 nocon- 12
31 5063 iter-1000 delt-
32 pre ••ure iteration failure at t- 3.688B+01 cycle- 7.2901-03 Docon- 13
3.8091+01 5229 34 7.290B-03/fa 7.211011-03 1.4C01I+02 3.6508+03 1.0001-03
33 7.29011-03 1.4C01I+02 3.77111+03 1. 00011-03
34 4.0101+01 5504 266 7.2901l-03/fa
pre •• ure iteration failure at t- 4.08911+01 cycle- 5613 iter-lOOO delt- 7.2901-03 Docon- 14
35 5615 iter-l000 delt- 7.29011-03
:u pre ••ure iteration failure at t- 4.0911+01 cycle- Docon- 15
4.2101+01 5779 233 7.2901-03/fa 7.2901-03 1.4391+02 4.02511+03 1.00011-03
37 6015 iter-l000 delt- 7.2901-03
38 pre ••ure iteration failure at t- 4.3821+01 cycle- nocon- 16
pre ••ure iteration failure at t- 4.3841+01 cycle- 6017 iter-1000 delt- 7.2'011-03 nocon- 17
3t 1.43511+02 4.33311+03 1.0001-03
4.4111+01 6054 521 7.290B-03/fa 7.2901-03
40 7.2901-03 1.4381+02 4.60411+03 1.0001-03
U 4.&111+01 U251 48 7.2901-03/fa
pre ••ure iteratioD failure at t- 4.65011+01 cyale- &382 iter-l000 delt- 7.2901-03 Docon- 18
42 4.652B+01 cycle- &385 iter-I 0 00 delt- 7.29011-03 Docon- 19
U pre •• ure iteratioD failure .t t-
4.8121+01 6604 34 7.2901l-031t. 7.29011-03 1.43711+02 5.02111+03 1.0001-03
44 4.9101+01 cycle- 67351 iter-l000 delt- 7.29011-03 nocon- 20
45 pre ••ure iteration failure .t t-
pre •• ure iteration failure at t- 4.91211+01 cycle- 6742 iter-lOOO delt- 7.29011-03 Docon- 21
"47
U
5.0121+01 6879
pre ••ure iteration
511 7.2901l-03/fa
f.ilure .t t-
7.21011-03
5.10011+01 cycle-
1.4361+02 5.36811+03
7000 iter-lOaD delt-
1.00011-03
7.29011-03 nocon- 22
pre •• ure iteration failure .t t- 5.10211+01 cycle- 7002 iter-1000 delt- 7.29011-03 Docon- 23
U
pre ••ure iteration failure et t- 5.10311+01 cycle- 7003 iter-1000 delt- 7.29011-03 nocon- 24
50 7.25101-03 1.43&11+02 5.81211+03 1.0001-03
51 5.2131+01 7154 54 7.HOI-03/t.
52 pre ••ure iteration failure .t t- 5.25111+01 cycle- 7207 iter-lOOO delt- 7.29011-03 nocon- 25
5.4541+01 7448 2 4.91611-01/c" 3 .1511-02 1.4351+02 4.95511+02 1.00011-03 54
53 4.515911+02 1.00011-03
5.4&01+01 7450 3 7 .25101l-03/fa 7.25101-03 1.4351+02
54
5.5121+01 7481 9 7.2901l-03/t. 7.2901-03 1.4341+02 6.40111+02 1.0001-03
55
void
(0) 1 B
liquid
void
B
(b)
(a)
0.9
mean
kinetic
(b) energy.
10-7 J/kg
0.0
0.0 610
time. s
(a) 1
liquid
void
(b)
liquid
---
.•..
Ml : Iterations
CD
Ml : volume error E
:J
"- --- (5
>
-~
CD
.0 no shrk: Iterations
E -e-
:J
c:: 0
c:: no shrk: vol. error
o if.
..:
~CD 0
t:
CD
.
.t:::
CD
20
E
-----
:J
................ 10 ~
5
....
--------------;,;.--:.;-:.:-...-. ..:.::',;;.--' ..
............ :::-; ~-,,~~ -.-......-----------
.....
-10
0.0001 0.001 0.01 0.1 1
epsllon,l/S
f f
r
(a) (b) u a
i c
d
0
n
0.5
(c) (d)
(e) (0
f f
r
(8) (b) u a
c
d
o
n
0.5
(c) (d)
(e) (0
At. s
(a)
At, S
0.0
0.0 time, S B.O
0.005
(b)
At, s
0.0
0.0 time, s 23
Conclusions and
Recommendations
8.1 Conclusions
Two numerical models were developed to describe volumetric changes during solid-
ification in metals and to predict shrinkage defects in the solidified castings. The
first model, Ml, is based on the full system of the continuity, Navier-Stokes and
enthalpy equations. Volumetric changes are described by introducing a source term
in the continuity equation which is a function of the local phase transformation rate.
The source term can include both volumetric shrinkage and expansion.
The second, simplified shrinkage (M2) model involves the solution of only the en-
thalpy equation. Simplifying assumptions that the feeding flow is governed only by
gravity and solidification rate and that phase transformation pro cedes only from
liquid to solid allowed the fluid flow equations to be excluded from consideration.
• Inclusion of the source term in the continuity equation and the equation for
the VOF function, F.
145
• Modification of the SOR pressure-velocity iteration scheme to allow internal
cavities to open and close .
In addition to these changes to the code which were required to incorporate the new
shrinkage models, a few modifications were also made to improve the accuracy of
the metal/mould heat trasnfer and solidification algorithms. These include:
• Development of a new expression for the drag coefficient for the drag force in
momentum equations which accounts for the fact that a computational cell
may contain a mixture of solid and liquid phases.
• Addition of the lever rule and Scheil equation to describe solidification III
alloys.
The critical value of the solid fraction was used as the feeding criterion in both
models.
146
As part of the development of the upwind differencing advection algorithm used
in the simulations, the Leith's method was incorporated into the standard two-
dimensional SOLA code. It is shown that the resulting scheme is unconditionally
stable despite being explicit.
1. In general, the results of the full shrinkage model are only marginally closer
to the experimental ones than those of the simplified model. The simplified
model underestimates the depth of the primary cavities even more than the
full model and predicts an earlier separation of internal liquid region from the
feeder. The latter led to larger secondary cavities in the T-shaped casting.
These differences are explained by the fact that the simplified model does not
describe advective heat transfer in the feeding flow. Nevertheless, in the case
of the pure iron and steel cylinder casting the two models gave practically
identical cavity predictions close to the experimental results.
2. Predictions of the size and position of secondary cavities are close to the ex-
perimental results, especially for pure metals.
3. Pure metal simulations for the T-shaped casting demonstrated the importance
of modelling of the mould filling stage. It is shown that the size and occur-
rence of these shrinkage defects are dependent on the filling rate. For the
T-shaped casting it was possible to eliminate the shrinkage in the T-junction
by increasing the filling time by a factor of three.
6. Cooling rates and total solidification times simulated by the two models are
very similar in most cases (less than 1% difference). The small difference
arises from a more accurate conservation of the total fluid volume by the M2
model since it does not have errors associated with finite iteration residuals
and advection of the fluid fraction function, F.
147
7. Both shrinkage models can be applied after the filling stage simulation has
been completed. The simplified shrinkage model requires less storage and
computer time, therefore, it can be used efficiently in cases where the metal
flow during solidification is not important and the cooling process is dominated
by conduction within the metal and heat transfer at metal/mould interface.
On average the simplified model runs 10 to 20 times faster than the full model.
8. The SOR iteration method, together with the VOF technique for free surface
description, are flexibile and accurate for application for shrinkage modelling.
However, when a large volume of liquid free surface is resolved only by one
or two cells, which is a likely situation during shrinkage modelling, iteration
problems arise if the free surface boundary pressure changes abruptly due to
a change in the slope. Such pressure pulses magnify long wavelength residuals
which decrease slowly in SOR iterations.
9. The time step size in the Ml shrinkage model is controlled by the free surface
stability limit during most of the simulation time, which is a consequence of
assuming hydrostatic pressure distribution in surface cells. In the simplified
model the time step size is limited by the heat transfer stability criterion. This
is usually several times larger than the free surface time step limit. The largest
deficiency of the full shrinkage model remains its computational cost due to
time step limitations and additional effort to achieve iteration convergence.
10. The shape of the cavities in solidified castings is mesh dependent though the
shrinkage volume is not. The mesh dependence arises from truncation errors
in free surface representation and temperature gradient calculations. However,
these errors decrease as the mesh is refined.
11. A value of unity for the critical fraction of solid fSl er , used as a feeding criterion,
gives good results for pure metals. For alloys a value above 0.8 appears to be
most reasonable. For fast solidification and a narrow freezing range alloy,
fSl er = 1.0 also gives best results.
12. Approximation of the temperature profile in the interfacial cells between the
cell centre and the interface by a linear function gives a second order accu-
racy with respect to the cell size ~x in the direction normal to the interface.
Assuming a constant temperature profile between these points reduces the ac-
curacy to the first order, besides, the coefficient h8T/8xl z =o of the first order
term is large for low diffusivity materials such as sand, especially at the start
of the simulation. In the latter case the first order term can be substantially
larger than the zeroth order term.
148
13. The second order heat transfer algorithm, in conjunction with the FAVOR
method for representing mould walls, substantially reduces mesh dependence
of the interfacial heat flux, compared with the first order method.
14. The solidification drag which arises from the fact that momentum equations
are applied to the solid/liquid mixture, proves to be an adequat'e method for
describing momentum changes due to phase transformation. In particular,
fully solid computational cells are marked with an infinitely large value of the
drag that is used to define the boundaries of liquid regions in the casting.
15. Since velocities in a feeding flow are usually small, between 0.1 to 10 mm/s,
but pressure in the liquid metal bulk changes significantly as the free surface
level drops, the feeding flow is pressure dominated, that is
In that case the value of the convergence criterion f required for an efficient
calculation is several orders of magnitude higher than that estimated from the
velocity scale alone. An appropriate formula to estimate f was worked out to
complement the existing one for velocity (or inertia) dominated flows. Feeding
flow stabilises numerical solution in the sense that it allows one to use a much
larger value of € than in the absence of shrinkage.
16. It has been shown that the Leith's method for approximating advection terms
in momentum equations can be used in the standard SOLA algorithm, instead
of the upwind differencing. The advantage of the Leith's method, also called
'semi-Lagrangian' method, is in its unconditional numerical stability, so that
the CFL restriction for the time step size is removed. However, preliminary
results show no gain in the efficientcy of the calculations since an increase in
the time steps size causes an increase in iteration. This indicates that large
time steps reduce the accuracy of the first-guess velocity approximation. An
application of the method in a colocated grid would improve its efficiency.
149
of flow in the mushy zone and feeding, more sophisticated than the critical
fraction of solid criterion, should be formulated. Lack of accuracy in alloy
solidification remains the largest problem in the physical model used in the
present work.
4. The full shrinkage model can be used as the basis for developing more sophis-
ticated models to include gas evolution and the formation of microporosity.
5. The truncation accuracy of the numerical equation for the metal/mould inter-
face heat transfer should be improved to minimise mesh dependence.
6. The simplified model does not employ the upwind differencing nor the donor-
acceptor advection scheme for the free surface which constitute the core of the
150
FLOW-3D. Therefore the shrinkage model could separated be from FLOW-
3D to create a smaller software package which can be used as a solidification
modelling tool on smaller computers.
7. Further work should be carried out to investigate ways of applying the Leith's
method to momentum, enthalpy and VOF equations to increase-the efficiency
of mould filling simulations in which the CFL stability criterion poses a severe
restriction on the time step size when the upwind and donor-acceptor advection
schemes are used.
151
Bibliography
[9] F.M.White, "Viscous Fluid Flow", McGraw-Hill Book Company, New York,
1974.
152
[11] S.Sundarraj, V.R.Voller, "A Discussion on Results of Recent Microsegregation
Models", in "EPD Congress 1992", eds. J.P.Hager, TMS Publication, pp 461-
479, 1991.
[16] S.G.R.Brown and J.A.Spittle, "A 2-D Implicit Finite Difference Model to Sim-
ulate the Columnar to Equiaxed Zone Transition", in "Modeling of Casting
Welding and Advanced Solidification Processes V", eds. M.Rappaz, M.R.Ozgu
and K.W.Mahin, TMS Publication, pp 395-402, 1991.
153
[20] A.Sommerfield, "Partial Differential Equations in Physics", Academic Press,
New-York, 1964.
[30] Ge-Cheng Zha, Dao-Zhi Liu and Tie-You Ma, "An Efficient Up-
wind/Relaxation Algorithm for the Euler and Navier-Stokes Equations", Int.
J. Num. Meth. Fluids, 9, pp 517-529, 1989.
154
[34] O.C.Zienkiewicz, Y.K.Cheung, "Finite Elements in the Solution of Field Prob-
lems", The Engineer, pp 507-510, 1965.
[36] J.R. Howell , ?, in "Modelling and Control of Casting and Welding Processes,
Vol. 3" , eds. S.Kou and R.Mehrabian, AIME, p. 377, 1986.
[38] G.Bonacina, G.Comini, "On the Solution of the Non-linear Heat Conduction
Equations by Numerical Methods", Int. J. Mass Transfer, 16, pp 581-589,
1983.
[43] G.P.Backer, "Finite Element Free Surface Flow Analysis: A New Tool for
Foundry Engineers", in "Modeling of Casting, Welding and Advanced Solid-
ification Processes VI", eds. T.S.Piwonka, V.Voller and L. Katgerman, TMS
Publication, pp 405-412, 1993.
155
[44] A.M.Hines, J.S.Tu, "Modeling of Fluid Flow and Heat Transfer in the Solidifi-
cation of Superalloy Investment Castings" , in "Modeling of Casting, Welding
and Advanced Solidification Processes VI", eds. T.S.Piwonka, V.Voller and L.
Katgerman, TMS Publication, pp 461-468, 1993.
[45] A.Shapiro, W.Stein and P.Rabain, "Casting Process Modeling Using Pro-
CAST and CAST2D" , in "Modeling of Casting, Welding and Advanced Solid-
ification Processes VI", eds. T.S.Piwonka, V.Voller and L. Katgerman, TMS
Publication, pp 493-500, 1993.
[47] J.Clark, "Free and Moving Boundary Problems", Clarendon Press, Oxford,
1984.
[49] R.A.Stoehr, C.Wang, "Advances in Fluid Flow, Reat Transfer and Solidifica-
tion Modeling and Applications to Actual Foundry Problems", in "Modeling
of Casting Welding and Advanced Solidification Processes V", eds. M.Rappaz,
M.R.Ozgu and K.W.Mahin, TMS Publication, pp 725-732, 1991.
[51] C.W.Hirt, B.D.Nichols, "Volume of Fluid (VOF) Method for the Dynamics of
Free Boundaries", J.Comput. Phys., 39, pp 201-225, 1981.
[52] R.A.Gentry, R.E.Martin and B.J.Daly, "An Eulerian Differencing Method for
Unsteady Compressible Flow Problems", J. Comput. Phys., 1, pp 87-118,
1966.
156
[55] K.Torrance, R.Davis, K.Eike, P.Gill, D.Gutman, A.Hsui, S.Lyons and H.Zien,
"Cavity Flows Driven by Buyoancy and Shear", J. Fluid Mechanics, 51, part
2, pp 221-231, 1972.
[60] C.Wang, "Computer Modeling of Fluid Flow and Heat Transfer in Metal Cast-
ings", Ph.D. Thesis, The University of Pittsburgh, Pittsburgh, 15261, 1990.
[65] H.-J. Liu, W.-S.Hwang, "Combined Fluid Flow and Heat Transfer Analysis
for the filling of Castings", AFS Transactions, 144, pp 447-458, 1988.
[66] R.A.Stoehr, C.Wang, "Coupled Heat Transfer and Fluid Flow in the Filling
of Castings", AFS Transactions, 176, pp 733-740,1988.
157
[67] H.J.Lin, H.L.Tsai, "Numerical Simulation of of an Integrated filling-
Solidification Casting System", in "Modeling of Casting, Welding and Ad-
vanced Solidification Processes VI", eds. T.S.Piwonka, V.Voller and L. Kat-
german, TMS Publication, pp 381-388, 1993.
[68] A.Brandt, J.E.Dendy, Jr. and H.Ruppel, "The Multigrid Method for Semi-
Implicit Hydrodynamics Codes", J. Comput. Phys., 34, pp 348-370, 1980.
[70] B.D.Nichols, C.W.Hirt, "Improved Free Surface Boundary Conditions for Nu-
merical Incompressible Flow Simulations", J. Comput. Phys., 8, pp 434-448,
1971.
[71] H.-J. Liu, W.-S.Hwang, "Three-Dimensional Fluid Flow Simulation for Mold
Filling", AFS Transactions 171, pp 855-862, 1989.
[72] C.R.Swaminathan, V.R.VoUer, "An 'Enthalpy Type' Formulation for the Nu-
merical Modeling of Mold Filling", in "Modeling of Casting, Welding and
Advanced Solidification Processes VI", eds. T.S.Piwonka, V.Voller and L. Kat-
german, TMS Publication, pp 365-372, 1993.
[77] G.K.Reziak, P.Gilotte and R.Hamar, "Modeling of Fluid Flow During Mold
Filling", in "Modeling of Casting, Welding and Advanced Solidification Pro-
cesses VI", eds. T.S.Piwonka, V.Voller and L. Katgerman, TMS Publication,
pp 435-442, 1993.
158
[78] A.Bourg, A.Latrobe, P.Large, P.Laty and C.Rigant, "SIMULOR and Fast
Solutions for Mold Filling Simulations", Bristol, September 23, 1991.
[79] C.W.Hirt, J.L.Cook and T.D.Butler, "A Lagrangian Method for Calculat-
ing the Dynamics of an Incompressible Fluid with Free Surface", J. Comput.
Phys., 5, No.1, pp 103-124, 1970.
[81] J.B.Bugg, R.D.Rowe, "Modelling the Initial Motion of Large Cylindrical and
Spherical Bubbles", Int. J. Num. Meth. Fluids, 13, pp 109-129, 1991.
[83] FLOW-3D Manual, FLOW SCIENCE Inc., Los Alamos, NM, 1991.
[87] B.van Leer, "Towards the Ultimate Conservative Difference Scheme, IV: A
New Approach to Numerical Convection", J. Comput. Phys., 23, pp 263-275,
1977.
159
[91] J.F.Evans, J.Beech and D.H.Kirkwood, "The Determination of Metal/Mold
Interfacial Heat Transfer Coefficients and the Prediction of Gross Shrink-
age Cavities in Chill Mold Castings", in "Modeling of Casting Welding
and Advanced Solidification Processes V", eds. M.Rappaz, M.R.Ozgu and
K.W.Mahin, TMS Publication, pp 531-538, 1991.
[92] L.E.Smiley, "Use of a Personal Computer to Predict Casting Heat Flow and
Solidification", AFS Transactions, 67, pp 689-696, 1988.
[94] H.Huang, J.L.Hill and J.T.Berry, "A Free Thermal Contraction Method for
Modelling the Heat-Transfer Coefficient at the Casting/Mould Interface" , Cast
Metals, 5, No.4, pp 212-216, 1993.
[95] J.A.Dantzig, S.C.Lu, "Modeling of Heat Flow in Sand Castings: Part I. The
Boundary Curvature Method", Metallurgical Transactions B, l6B, pp 195-
202, 1985.
[96] J.A.Dantzig, J.W.Wiese, "Modeling of Heat Flow in Sand Castings: Part II.
Applications of the Boundary Curvature Method", Metallurgical Transactions
B, l6B, pp 203-209, 1985.
[100] C.R.Su, H.L.Tsai, "A Direct Method to Include Latent Heat Effect for Mod-
elling Casting Solidification" , AFS Transactions, 82, pp 781-789, 1991.
160
[102] I. Ohnaka, T.Fukusako, "Calculation of Solidification of Casting by a Matrix
Method", ISIJ Transactions, 17, pp 410-418, 1977.
[105] J.H.Chen and H.L.Tsai, "Comparison of Different Modes of Latent Heat Re-
lease for Modeling Casting Solidification", AFS Transactions, pp 539-546,
1990.
[111] Llmafuku, K.Chijiiwa, "A Mathematical Model for Shrinkage Cavity Predic-
tion in Steel Castings" , AFS Transactions, 10, pp 527-540, 1983.
[113] V. de L.Davies, ?, AFS Cast Metals Res. J., 11, p. 33, 1975.
161
[116] E.Niyama, T.Uchida, M.Morikawa and S.Saito, "A Method of Shrinkage Pre-
diction and Its Application to Steel Casting Practice", AFS Int. Cast Metals
Res. J., 9, pp 52-63, 1982.
[118] V.Jesko, J.Zajac, "A New Method for Shrinkage Prediction in Casting Walls",
AFS Transactions, 48, pp 681-684, 1991.
[123] C.W.Hirt, J.M.Sicilian, "A Porosity Technique for the Definition of the Obsta-
cles in Rectangular Cell Meshes", Proc. Fourth Int. Conference Ship. Hydro.,
National Academy of Science, Washington, DC, September, pp 212-235, 1985.
[124] M.lshii and T.C.Chawla, "Local Drag Laws in Dispersed Two-Phase Flow",
Argonne Nat. Lab. report NUREG/CR-1230, 1979.
[125] V.R.Voller, "An Overview of the Modelling of Heat and Fluid Flow in Solidifi-
cation Systems" , in "Modeling of Casting Welding and Advanced Solidification
Processes V", eds. M.Rappaz, M.R.Ozgu and K.W.Mahin, TMS Publication,
pp 661-674, 1991.
[126] J.F.Evans, "Central Unsoundness in Cast Steel Rolls", Ph.D. Thesis, The
University of Sheffield, 1990.
162
[127] H.You, "Computer and Physical Modelling of Heat and Fluid Flow in Relation
to the Solidification of Alloys", Ph.D. Thesis, The University of Sheffield, in
preparation.
[128] CASM ACME Review Report No.3, The University of Sheffield, 1993.
163
Appendix A
Summary of Computer
Simulations
164
Appendix B
Semi-Lagrangian (Leith's)
Method for SOLA
The application of the Leith's method to the discretisation of the advection terms
in the momentum equation in two dimensions in the absence of free surfaces is
described in this Appendix. This advection scheme is explicit and unconditionally
stable. It posseses the same properties as the upwind differencing scheme, that is
it is transportive and conservative. It can be shown that compared to the upwind
differencing scheme the Leith's method introduces less numerical diffusion. This
because first order truncation errors in the discretisation of the advection terms in
the Leith's method do not contain terms proportional to (for x-component)
where S(t, x) comprises pressure, body force and viscous terms, may be written
using the substantive derivative rather than the partial:
du =S (B.2)
dt
Integration along the trajectory x(t) between t and to gives
Equation (B.3) is solved at each node with to = t n and t = t n +1 and x'""{t n +1 ) defined
as
Yn+! = Yi,i
Equation (BA) is solved numerically using the n-Ievel velocities to obtain the posi-
tion of the x( t n ) point. In general, Xn and Yn will not fall on a computational node.
A linear interpolation from four (eight in 3-D) nearest nodes is used to obtain the
value of u( tn, Xn , Yn). If the point (Xn' Yn) lies between these nodes (in other words
interpolation, and not extrapolation, is used) then the unconditional numerical
stability of the method with respect to advection is ensured [135].
Roache [27] mistakenly concludes that the Leith's method is only conditionally stable
with the same time step restriction as the explicit upwind differencing method. The
mistake arises from the implicit assumption by Roache that u(tn, Xn, Yn) is always
interpolated from the same node points. In that case, of course, the point (xn' Yn)
has to fall between these nodes to maintain stability and that is the origin of the
time step restriction.
(B.5)
The representation of the source term given in Eq. (B.5) is equivalent to an explicit
formulation and is first order accurate in time. It may introduce a stability time
step limit, e.g. that given by Eq. (2.76) if S includes viscous stresses.
If linear interpolation is used for evaluating u(tn, In, Yn), then a numerical diffusion
and a phase error are introduced [27]. The latter means that perturbations of
different wavelength propagate at different speeds in the numerical solution even for
166
a uniform velocity field. This is of course of no concern if a steady state solution
is sought. The diffusion and phase error effects can be reduced or even removed
completely in two ways. One is to obtain a more accurate solution Xn and Yn of
the trajectory Eq. (BA) by splitting tlt into several smaller time steps. Solving
equations at each of the intermediate step will involve also an estimation of the
velocities at each intermediate trajectory point by interpolation. The other way
is to use a higher order interpolation to estimate u(tn, X n , Yn). Normally it would
involve more nodes around (xn' Yn) but a linear interpolation will still be necessary
at the boundaries.
Apart from unconditional stability for the advection terms, the method also pos-
sesses the conservative and transportive properties [27]. It is interesting no note
that, unlike most other higher-order schemes, the Leith's method preserves the
transportive property even in its higher-order accuracy formulation.
The Leith's method has been well developed and extensively applied in the at-
mospheric circulation modelling [135,136]. The absence of solid boundaries in the
horizontal direction removes the difficulty of resolving the trajectories near walls as
well as allows one to use a higher-order accuracy version of the method throughout
the computational domain.
It appears possible to use this method successfully for an incompressible fluid flow
in the presence of walls. Here the method is incorporated into SOLA program [53]
and applied to a two-dimensional flow problem. In a staggered grid (Fig. 2.1a)
the calculation of the first guess for u-velocity in cell (i,j) at t = t n +! proceeds as
follows:
1. coordinates Xao and Yao are chosen to be at the nodal point where velocity Ui,i
is located (point ao in Fig. B.la) so that
Uao = ui,i
2. the y-direction velocity component Vao at point a is estimated by linear inter-
polation which in a uniform grid is
1
vao = 4(Vi,i-l + Vi,i + Vi+l,i-l + Vi+!,i)
3. Equation (B.4) are solved by dividing the time step into m equal parts. For
each step 1 with t:,.t m = t:,.t/m, coordinates Xal and Yal of an intermediate point
a,are found by
(B.6)
167
velocities U a, and Val at each point a, are found by linear interpolation from
the neighbouring nodes.
U':1'tl
t,J
= U
am
+!It· ST}·
1,J (B.7)
Suppose that the space x ::; 0 is blocked by an obstacle and there is a computational
cell covering a segment 0 < x < 1. One face of the cell (x = 0) is blocked by the wall
and the velocity component Uo on the opposite face (x = 1) points in the direction
away from the wall, Uo > 0 (Fig. B.lb). The velocity between the two faces is
defined by linear interpolation
u = x· Uo
In this case the solution of equation dx/dt = u is
This equation means that whatever the integration time tl - t2 is, the analytical
solution for X2 is always positive and outside the obstacle (though it can be infinitely
close to it if Uo > 0). In terms of the numerical solution it means that if the time
period between t2 and tl is divided into sufficiently small sub-time steps, then the
numerical result should be close to the analytical one, that is, also outside the wall.
The modified SOLA program is given in appendix C. The Leith's method was used
here to simulate the inviscid flow in a 'water-cooled reactor' described by Hirt et
al [53]. Figure B.2 shows the geometry and boundary conditions setup for the
168
problem. The velocity at the entrance was constant and set equal to u = -0.01 m/ s.
A 20 x 30 mesh, uniform in each direction, was used with ~x = 5 mm and ~y =
6.5 mm. With the upwind differencing scheme the time step size was chosen to
be ~t = 0.05 s to maintain numerical stability. Figure B.3 shows the result for
~t = 0.07 s where an instability occurred at t = 1. 7 s of simulation due to a violation
of the CFL criterion. The total simulation time was 8.0 s. The conver~nce criterion
was f = 0.005 S-1 in all simulations. Figure B.4 shows velocity plots at four different
times produced by the unmodified SOLA. The total CPU time was tcpu = 21.0 s
with the average iteration count equal Nit = 21. To simplify comparisons, the CPU
times and average iteration counts were normalised to those of the unmodified SOLA
simulation. Therefore, for the latter tcpu = 1 and Nit = 1.
For the Leith's method four time step sizes were used: ~tl = 0.05 s, ~t2 = 0.1 s,
~t3 = 0.2 s and ~t4 = 0.5 s. Figure B.5 displays velocity plots for Ilt = 0.05 sand
m = 1 at the same times as in Fig. B.4 showing close agreement with the upwind
method predictions. In this case tepu = 1.11 and Nit = 1.17.
Figure B.8 shows velocity plots at t = 6.0 s for the four time step sizes with m = 1.
It is clear that the solution was stable in all cases though its accuracy deteriorated
as the time step size increased. This also lead to slower convergence as shown in
Figure B.9 where the dependence of average iteration count on the time step size is
shown. Fig. B.lO gives predicted flows for ~t = 0.5 sand m = 1, m = 2, m = 10
and m = 20. As in Fig. B.6, little improvement can be seen with the increase of
m. This can be explained by the crude approximation for the pressure gradient
terms, given by Eq. (8.5). The larger the time step size the less accurate is the
assumption that the work of the pressure forces along the fluid particle trajectory
can be approximated by the work of the pressure forces at only one point of the
trajectory, made in Eq. (B.5).
Finally, Fig. 8.11 summarises CPU times of all simulations in this Appendix as
functions of the time step size and the value of m. It can be seen that none of
the simulations employing the Leith's method was more efficient than the original
method.
169
The potential of the Leith's method can be realised more efficiently if the integral
on the right-hand side of Eq. (B.3) is approximated more accurately. This will give
a better estimate for the first-guess velocities, therefore, accelerating convergence
an improving the overall accuracy of the solution.
It is unclear as to the best means of applying this method to transieftt free surface
flow. In the SOLA-VOF algorithm velocities in the void are not known, therefore
Eq. (B.4) cannot be solved in empty cells which are being filled with fluid during
the current time step. A possible means is to use the variable density method
(Section 2.4) in which velocities are estimated both in the fluid and in the air. This
will require a higher order version of the Leith's method to avoid a diffusion of the
fluid/ air interface.
The use of a staggered grid for the Leith's method reduces the efficiency of the
algorithm because in every computational cell velocity components are defined at
different locations. Therefore Eq. (B.4) has to be solved separately for each velocity
component. The starting and finishing points of the solution of Eq. (BA) are
different for each component and interpolation coefficients are also different. In
a colocated grid all three velocity components are located at the cell centre. In
that case Eq. (B.4) would have to be solved only once with the cell centre as the
starting point. The interpolation coefficients would also be computed just once for
the finishing point of the solution of Eq. (BA).
170
....... ............. . v
ao ./'
j ....,.~;.....--1f----~ uao
(a) -V·llt ..........
I .......••
al.~·
f d
(b)
____________~----~--------~~~ x
o 1
-.---
f
.. u=lO mm/s
13 mm
1
19.5 mm
..-
2mm
-~ ......
If-----10-m-m----... ~~1
I I
I~
I
I
I I
I I
, I I I I
, I I I I
II \I I I I'
I I
I \I \I I I
III \III I I'
II I \I 1 I I I'
1\ \I 11\ 1 I I I I I I
\I I I \ \\ \ I I 1\ III
111111\\\\\\ 1 1 \I II I , I .
" "11111'\\\ II 1111\\\1\11
""'IIII\\"~ II 11\\\\\\\\\:/
,, ,'I 1\\\\\\\\\'-/
"\ \ \\\"\"'--/
_/~ y,'/
I I , , ,
, , , , I \ \ \ \ \ \ , .....
" .... , ........ - - " ' , / 1 I"" " " " " __ /111
........... --~"",'
, ... " ............ - - - ; / 1 /
.. --- ... ',"
, ............
.............. - - ....... , I ,
(a) (b)
1 1 //-
1 1
, '~
11 , 1
1
1
, 1
1 I
1\
1 I \
I "I"
I " I I "
I " I I , 1/.
1 /I I / I I , .
II /1111".
II 111111.,
II I 1/ I I I • \
1\\\\1111\ ' \/
1\\\\\\\\\-1/
I I I \ \ \ \ \ \ \ , ..... -...-
I I \ I\\\\\"--/~Hq
\ \ \ I \ \ \ " ' -...... / , , / ,
" \\\""'---//~~I
" , ................ ----.-./111
(c)
I~
I I I I I I I I I I
I I I I
I I I I I I I
I I I I I
I I I I I I
I I
I I I I I I I I I I I I
1111 1111 IIII
I I I
1\ I I I I II I 1\ I I
I I \ , , , , I 1\ I ,
I I \ I \ I
I I I \ \ \ 111111111/111
I I I I \ I ""'" I I I \I
I I I \ I \ II , I I , , I I "
I I I I \ I \ 111111111111,
\ I I II \I \ 11111111111"
I I I II \I \ 11111111111, ..
I I 11111111\ 111111/1/1, •. .
I I 111111\\\1 1111111111,.
\I 1111111//1, .,
,
I I I I I \ \ \ \ \ \ \
I I , , 1\\\\\1\\\ I 111/11.-'\
I I , , ,\ \ \ \ \ \ \ \ I (111.'\\
1\ \ , I\\\\\\"~ I 11 I • ' \ \
I I I I '\ \ I
II \ \\\\\ .I I II
\\\\\\\\\\,~'///,
I I \ \ \ \\ \ \ ,,_ .... / 11,
\\\\\\\\',-..../Il Z
\\ \\ \\ \\\, ""',------..":' //'/. ~Zr.I I
" . . . . . .__:::.'/'/ I
,............................... :,....."////1
.. _------- --///1'
(b)
(a)
I~
r .. .. ~ •
(c) (d)
I
'ff 11111111111
11111111111
11111111111
11111111111
I I I I 11111111111
I I I I 11111111111
I I I I 111',,,,",
I I I I 11111"'111
I I I I 11111111111 ,
I I I I 11/111111'".
\II II 1111111111".
\I
I I I I 1111111111 •.. -
\I
I II I' 1111111111,.
I \I
I I \I I 1111111111,."
\1\1\1\\1 111//"111, ' I
1111\1\\"
111111\\\\\1
111\1\\\\\\\\
III'II[/~/I':::
I I If : . , \ I
I \I ,,
"'I\\'\"'~
I : I , ,.
(b)
(a)
111111111111'11 ./~
111111111111111
111/1111'111111
11111111111""
1111111111/'"1
111111111111110
111111111111 11 ,
11111111111''' ,
1I1I1I'lfll'"
/1111111" .. .
1111111" . ..
'Zlllll" .. ..
l l l l l , ... ..
Z
dIll, ... "
Il
d~~:::~::
1/"''11\
111""'1
IJ, • , I , \ I
II, '" I I
~~~~~~::~l);
\'~\'\\~~~:~~;;}~ 11/1 I
,\ \ \"-////11/ 11 7'
\ \ \"----.. . .///'l'l~ I •
1\\\" .....--...........////1 1 ' .
1\" ..........' . _ ....../ / / / 1 ' · ·
, , , ......... -------,.,,.. / , , . . . ....... _------ ... , .
(c) (d)
I
ZI~~~~:::~: :
I ; •.• \ ,-
I
I.IZ~~;::::::
111',\11, I~~:'''II
111"'111 II • .• , 1I I
III .• I I , I I I, . , 1"
II
~~\\"'IIII/
I "1',I
\" \ I , .• I I , I
\ I.
.
\\\,"'//1'
n~\, ~"1/11
'/ I , I
~~~~~2~;:~t~~
~~~
~
\ \ \ \ " - - - -.....
1
I '
\\\\~~~:.:;~~/~ ~ ,
\ ~ ~~~~~.::~~~~ ~ ~ :
----,//hl
1\ " , '.......
I \ " '......~·
I' . 1 \ ' " ,.....::::::::-////11'
I" .... --~////II~
.
-"....////1;·
, " ..... , - - - -.- --,.,,/; , .. ·
"-""'/;1'
.
..... - - - - - - - - - , , . #
(b)
(a)
11111111111""
111111""'1'" ./~ 111/111'1"", '/~
I ~~Y.':~
111111""'"''
1111111/"',", lit(, I \I I \I /I / I I " "
1/ I I / " , / / "" ,
I I 1// " "
1111"""/1111
'" 1 " /!J~ 111111111/'""
11111/111"""
11111111111//1' '"111111/,,,,.
,,'1111111"01, 1111111'"' . . . .
11111111111 . . . . I / l I l I f l l l f I , ..
111/1111/11 ••. . 1111/1111.
1 I I I I I I u ... . 1111111, '"
1111111; .. . .
~/lII/;, .. . ,
HII//I,:
ll/",.
:: :
d ll ",,.. \ 1 Zl
"ZII!/", .• ,
l.d lll
" .. "
/, I. ! I I , •• \ I
I
/~~;::::::
I.IZ~~;::::::
111"\11\
I , , . . \ \ 1I
I I , , • \I 1 I 1/ I II I I
• •
II • .• I 1" III • I I, I
\ \. ' I "1
II 1 'I I , I
\" . " / I I I
\ I \ . • I I I I •
~\"""'I/I
\
\\"""11/1
~~\\"""II//, n\\\"~~"IIf/,
\ " ... _~////~
\\\\\\\~~~:':':':~~~l/II 7 ~
.11
(c) (d)
5
c::
1.
E
:J
c:
c:: c:
o o
~ ~
~ 1. I 26 ~
.. ...
I
"0 \ I "0
.'
CD / CD
, !!
!!
"'E 1.1
o
23 "'oE
c:: c:
1+---.--.---r---r--r---r--r--~--._-_+~
o 2 4 6 8 10 12 14 16 18 ~
number of steps. m
'~
I I I , I I I I I I I II , ,
11111111/111"
11111111/11'" Il / I I I I I i / I I ' "
IIIIllIfllI"" 1 I ' " I I I i / I Il "
I I , , II I I I I I I " , 1111111111 " " ,
1111111111"" 11/11111111", .
111111111"'00 11/111111"10,
111//11111/100 I f l l l I l l l l " ..
/1111111" .. 1111111110 . .
/111111" . .. 1111111".
l l l l / I I , .... "11111, . '.
Z
1,"/111, .. "
II/III •. ' "
I.UIIII .. : ::-
1.1~)i~~~::::::
'
/
11/"\\\1
lIZ~~~: = ~:::
I I • • , \ \ I
I I I •• \ \ \ I I, ' I \ II
II, . • , \ I I fI , ~ , II
\
1\, ' I I I I
\ \ \ • , I I /I '
\' \ \ \ .
,"
' I
' I
I \
I I
I /
I '
n\\\"'IIII/ \,n\\\\:~' 11/1/
~~~~~2~;~~~t ~
.\1\ \\, '''1 11 1,
I I \\\\\\,,:.:. . . ///t ZZ
\ \ \ \ \"---///1.1, Y. z
\.,----_~~~ ~
\ \ \ \ I 1\ \ \ \ ' ,_ _ _......../ / / / ~ I
1\ " " '_____///1 / I'
I \ \ ' " , .....- . . . : ; - -/ ' / / / / I I
I \ " ,......- - - - / / / / I" 1"' ..........- :"';////11
", ........ - - . . . - / / /
- ------ -
. .. . - - ... ;
/'
~
~
. ....... _------- ~, -.
(a) (b)
"" ..........
-=-.. ..,,, /
I I I I
_------*"'.,.;"
- - , . ' / / I"
. .. ----------, .. ..~--------------~;
(c) (d)
....
II>
---
.c
E
:J
c:
c:
o
~
~
"'0
II>
I/)
:;;
E
o
c:
I II 11111 I , ,
I II II II I I • , .
I II I I I I I , • :
•
'11111,,'"
1 I 1111, , ' • ,
I I I I I " • I
111/ • . I I
·
I
:
1111 I , • I I " I , •• • , 1
1111, , • II I , . 'II
1111, . , : I' \ II I . • , I
II I I •. \ \ \ \ ' . , • /I
II 1\ 'II \\1\,"''''/
\\I\\~~I/I/
\\\\\\\,,_II//. 1,
I \\\\\\,, __ ////.1.1.
1\1\\\"' ___ ~///Z1.
I I \\\\,,
\
I \\\\\",,--'///ZI.J
1\\\ ' ,"•' 1 1 "1
\\\\" . I//'/.Y,%
I '/ Z
I,
,,_=. . . . . /
\ II\\\"'---,~~~~~, \\\\"," '//////
\ \ \ \ ",,----,,// I / I I \ \ \ "",..:::::.-..... / / / 11II
,
,\ _________
\\""------~~//II
...... ,,; '" '" " ,'\ " ........ _ - - : - - / / / 1 / I
... - --_._------
--*!'""'"
~ ~
,
,
I
.
(b)
(a)
111111111111111
I I I 1 I I I I I I / I I " ,,,___ I I III I I I I I , I I " 1/'--
I I I / I / I I , , 1IIIIIIIIIIIIIt lit(
I I 1 I / / I , " II11111/ " " , " /,~
I I I I I I I , , " , ' 11/11111/11""
1111111111/1'" 111I1111t11l1l
1111111111111 •
",,"/1/11""
, " " ' 1 / 1 1 11 " I /I I I I I I I I I I
'"1/111/111" 11,'///11//111 : . .
11111111/1" .. 1111111,
1/111111111". 1 11/11111," .
/1/111111" " 1 /111111,
I I / l f l l I , ... , , I I I I ,
II/III, , . ' '" I I , : : ' I
IIII I , , •: I III I • 'I
II I , , .• , I
I111 I , , ( I, , • I t I
1 111 , , II I I • . . • ,
I'III, . ..
I, , .
'I
III, . " "
III,. " .,
II I I • . I
.
• ,
\ I \1\, . '" I
\ \ 1
I , , . "
\ \ \ \ , . , , '/ \
I \ ~ ~ ~ ~ ~ . : : ~!I'/,t 'I
I ... . "/
\ \ ' · · · · " 11
\\\\\"""'!~'/. I \\ --"/ I.
\ 1\\\",---/
\ \ \ \, ,, __- / ~ 1/ Z \\\\,~~~:--/~~~~!
'
\\\\\""--.. . -//~~Jj
\ 1 \ \ " ' , , - -.......... ////11
I I \ \ \ " , _ ' -...../ / I I I
" \ " , , - _ ............... //111
, \ , ..... .......... _ ---..-" ;I' / I I I ,
\ , \ " " - - - - - - _..... / / / / I I
, " , .... _ _ _ _ _ _ _ _ ..... ,. / , I ,
• " . __ -----,,,,,, I
-- - .,. ,.
--------, .
;
~
(d)
(c)
E
.,
:::l
1.6
1.5 .. -.-
Leith's: m=1
Leith's: m=2
:::l 1.4 -[3-
Q...
U Leith's: m=5
.,
"C 1.3 -H-
!l Leith's: m=10
r; 1.2
§
0
.....
~
1.1 Leight's: m=20
0.9
-
a81----.--~r---~--_.--~----r_--~--_r--_,----r_--~
o 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55
tlme step size, s
Figure Bll Total CPU time as function of the time step size
for the original SOLA upwind differencing method
(filled square) and the Leith's (or semi-Lagrangian)
method for a range of m values.
Appendix C
171
program sola
c
c Modified SOLA program with linear Leith's advection
c term approximation instead of the original
c upwind-differencing scheme. The Leith's
c method is implemented in subroutines
c EQTN and GINTER
c
dimension u(36,36),v(36,36),
* un(36,36),vn(36,36),p(36,36),xput(24)
real nu, pn(36,36)
real*8 cpul,cpu2
real*4 tarray(2)
integer cycle,wl,wr,wt,wb
c
integer nf(36,36)
c
common /dat/ nf, neqn, imax, jmax
common /mesh/ delt, delx, rdx, xmin, xmax, dely,
* rdy, ymin, ymax, rrxy
common /vel/ un,vn
call etime(tarray)
cpul-tarray(1)+tarray(2)
open(unit-10,file-"semiin",status-"old",form-"formatted")
open(unit-9,file-"semiout",status-"old",form-"formatted")
c
c read and print initial input data
c
xmin-O.O
xmax-l.O
ymin-O.O
ymax-l.O
read(lO,26) num
read(lO,25) (xput(i),i-l,num)
ibar-xput(l)
jbar - xput(2)
delx - xput(3)
dely - xput(4)
delt - xput(5)
nu - xput(6)
cyl .. xput(7)
epsi - xput(8)
dzro - xput(9)
gx I: xput(lO)
gy = xput(ll)
ui .. xput(12)
vi .. xput(13)
velmx • xput(14)
twfin - xput(15)
cwprt - xput(16)
cwplt - xput(17)
omg - xput(18)
alpha .. xput(19)
wl .. xput(20)
wr -= xput(21)
wt = xput(22)
wb - xput(23)
neqn=xput(24)
write(9,SO)(xput(i),i-l,num)
25 format(4(6x,e12.5))
26
format(6x,i2)
27
format(lh ,l8x,lOa8,lx,alO,2(lx,a8))
35
format(lhl)
44
format(6x,7hcycle- ,is,8x,4htd- ,lpel2.S,8x,
1 4ht2= ,el2.5,9x,6hiter- ,is)
45 format(lOa8)
71 format(lx,i4,e10.3,i4,/,lOO(3elO.3))
72 format(lx,i4,5x,3(6x,lpelO.3))
46 format(80x,3ht- ,1pelO.3,4x,7hcycle- ,i4)
47 format(6x,lhi,7x,lhj,l2x,lhu,17x,lhv,18x,lhp)
48 format(4x,i3,Sx,i3,3(6x,lpe12.S))
49 format(6x,6hiter- ,is,lOx,6htime- ,lpe12.S,lOx,7hcycle- ,i4)
50 format(lx,el2.5)
c 50 format(lh ,Sx,6hibar- ,lpel2.s/6x,6hjbar- ,el2.s/
c 16x,6hdelx- ,e12.5/6x,6hdely. ,el2.S/6x,6hdelt- ,
c 2e12.5/8x,4hnu- ,el2.S/7x,5hcyl- ,e12.S/6x,6hepsi- ,
c 3e12.5/6x,6hdzro- ,el2.5/8x,4hgx- ,el2.S/8x,4hgy- ,
c 4el2.S/8x,4hui- ,e12.5/8x,4hvi- ,
c Sel2.S/5x,7hvelmx- ,e12.S/5x,7htwfin- ,el2.s/sx,
c 67hcwprt- ,e12.S/Sx,7hcwplt- ,e12.5/7x,Shomg- ,
c 7el2.5/Sx,7halpha- ,e12.S/8x,4hwl- ,e12.S/8x,
c 84hwr- ,el2.S/8x,4hwt- ,el2.S/8x,4hwb- ,el2.S)
cc
c compute constant terms and initialize necessary variables
c
imax=ibar+2
jmax-jbar+2
iml-imax-l
jml=jmax-l
rdx=l.O/delx
rdy-l.O/dely
rrxy-rdx*rdy
delx2-delx**2
dely2-dely**2
rdx2-rdx**2
rdy2-rdy**2
dxdy-delx*dely
dtdx-delt*rdx
dtdy-delt*rdy
jm2-jmax-2
im2=imax-2
c
c * * NF is used as a flag for blocked and open cells
c * * NF=l - open cell, NF-O - blocked cells
c
do 10 i=2,im1
do 10 j=2,jml
10 nf(i,j)-l
c
taO.
iter-O
cycle-O
twprt-O.
twplt .. O.
beta- omg/(2.*delt*(rdx**2+rdy**2))
c
c special input data
c
c
c set initial velocity field into u and v arrays
c
do 560 i=2,im1
do 560 j-2, jm1
u(i,j)- ui
v(i,j)- vi
560 continue
assign 5000 to kret
go to 2000
c
c start cycle
c
1000 continue
iter==O
flg-l.
assign 3000 to kret
c
c compute temporary u and v
c******************************************************************
c * * temporary u and v are computed in EQTN * *
c*******************************************************************
do 1111 i ... 2,im1
xxu-xmin+delx*(i-1)
xxv-xxu-0.5*delx
do 1111 j ... 2,jm1
yyv-ymin+dely*(j-1)
yyu-yyv-0.5*dely
if(nf(i,j).eq.0.or.nf(i+1,j).eq.0) goto 1112
call eqtn(xxu,yyu,uint,vtest)
visx- nu*((un(i+l,j)-2.*un(i,j)+un(i-l,j))*rdx2+
1 (un(i,j+l)-2.*un(i,j)+un(i,j-l))*rdy2
2 +cyl*((un(i+1,j)-un(i-1,j))/(2.*delx2*float(i-1))
3 -un(i,j)/(delx2*float(i-1)**2)))
u(i,j)- uint+delt*((p(i,j)-p(i+1,j))*rdx + gx+visx)
1112 continue
if(nf(i,j).eq.0.or.nf(i,j+1).eq.0) goto 1111
call eqtn(xxv,yyv,utest,vint)
visy= nu*((vn(i+1,j)-2.*vn(i,j)+vn(i-1,j))*rdx2+
1 (vn(i,j+1)-2.*vn(i,j)+vn(i,j-l))*rdy2
2 +cyl*(vn(i+1,j)-vn(i-1,j))/(2.*delx2*(float(i)-1.5)))
v(i,j)= vint+delt*((p(i,j)-p(i,j+1))*rdy + gy+visy)
1111 continue
c********************************************************************
c********************************************************************
c
c set boundary conditions
2000 continue
do 2200 j-1,jmax
go to(2020,2040,2060,2080),wl
2020 u(l,j)-O.O
v ( 1 , j ) -v ( 2 , j )
c
nf(l,j)-O
c
go to 2100
2040 u(l,j)-O.O
v ( 1 , j ) a-V ( 2, j )
c
nf(l,j)=O
c
go to 2100
2060 if(iter.gt.O .and. flg.gt.0.5)go to 2100
u ( 1 , j ) "'U ( 2 , j )
v ( 1 , j ) -v ( 2 , j )
c
nf(1,j)-1
c
go to 2100
2080 u(1,j)cu(im2,j)
v ( 1 , j ) -v ( i m2 , j )
c
nf(1,j)=1
c
2100 go to (2120,2140,2160,2180),wr
2120 u(im1,j)-0.0
v(irnax,j)-v(irn1,j)
c
nf(imax,j)-O
c
go to 2200
2140 u(irn1,j)-0.0
v(irnax,j)--v(im1,j)
c
nf(imax,j)-O
c
go to 2200
2160 if(iter.gt.O .and. f1g.gt.O.5) go to 2200
u(im1,j)- u(im2,j)*(float(im2-1)/float(im2)*cyl+(1.0-cyl))
v(irnax,j)-v(im1,j)
c
nf(imax,j)-1
c
go to 2200
2180 u(irn1,j)-u(2,j)
v ( i rn1 , j ) -v ( 2 , j )
p(irn1,j)-p(2,j)
v(imax,j)-v(3,j)
c
nf(imax,j)-1
c
2200 continue
do 2500 i-1,irnax
go to (2320,2340,2360,2380),wt
2320 v(i,jrn1)-0.0
u(i,jrnax)-u(i,jrn1)
c
nf(i,jmax)-O
c
go to 2400
2340 v(i,jm1)-0.0
u(i,jmax)--u(i,jm1)
c
nf(i,jmax)-O
c
go to 2400
2360 if(iter.gt.O .and. flg.gt.0.5) go to 2400
v(i,jm1)-v(i,jm2)
u(i,jmax)-u(i,jm1)
c
nf(i,jmax)-1
c
go to 2400
2380 v(i,jrn1)=v(i,2)
u(i,jrn1)-u(i,2)
p(i,jrn1)=p(i,2)
u(i,jrnax)=u(i,3)
c
n f ( i , j rna x ) -1
c
2400 go to (2420,2440,2460,2480),wb
2420 v(i,l)-O.O
u ( i , 1 ) -u ( i , 2 )
c
nf(i,l)=O
c
go to 2500
2440 v(i,l)=O.O
u(i,1)--u(i,2)
c
nf(i,l).O
c
go to 2500
2460 if(iter.gt.O .and. flg.gt.0.5) go to 2500
v ( i , 1 ) -v ( i , 2 )
u ( i , 1 ) -u ( i , 2 )
c
nf(i,l)-l
c
go to 2500
2480 v(i,1)-v(i,jrn2)
u(i,1)-u(i,jrn2)
c
nf(i,l)=l
c
2500 continue
c special boundary conditions
c
do 2801 j-12,jrnl
c
nf(17,j)-0.0
c
u(16,j)-0.0
u(17,j).0.0
v(17,j).0.0
2801 continue
v(17,11).0.0
u(2l,29)--1.0
u(2l,30)--1.0
u(21,31)--1.0
do 2800 i-18,irn1
v(i,jrnl)-O.O
2800 continue
do 2802 j=2,28
u(21,j)-O.O
c
nf(22,j)-0.0
c
2802 continue
c
go to kret(3000,5000)
3000 continue
c
c has convergence been reached
c
if(flg.eq.O.)go to 4000
iter-iter+l
if(iter.lt.5000) go to 3050
c if(iter.lt.500) go to 3050
c if(cycle.lt.lO) go to 4000
if(cycle.1t.100) go to 4000
t- 1.e+10
go to 5000
3050 f1g-0.0
c
c compute updated cell pressure and velocities
c
do 3500 j-2,jml
do 3500 i ... 2,iml
c
if(nf(i,j).eq.O) goto 3500
c
d-rdx*(u(i,j)-u(i-1,j))+rdy*(v(i,j)-v(i,j-1))+cyl*(u(i,j)
1+u(i-1,j))/(2.*delx*(float(i)-1.5))
if(abs(d/dzro).ge.epsi)flg-1.0
delp. -beta*d
p(i,j)-p(i,j)+delp
u(i,j)-u(i,j)+dtdx*delp
u(i-1,j)-u(i-1,j)-dtdx*delp
v(i,j)-v(i,j)+dtdy*delp
v(i,j-1)-v(i,j-1)-dtdy*delp
3500 continue
go to 2000
4000 continue
c
c print and plot
c
5000 continue
if(t.gt.O.)go to 5030
print 50,(xput(i),i-1,num)
5030 continue
if(cycle.le.O) go to 5100
if(t+1.e-6 .It. twplt) go to 5600
twplt=twplt+cwplt*delt
5100 continue
print 49, iter,t,cycle
print 46,t,cycle
c
c list velocity and pressure fields
c
5600 continue
if(cycle.le.O) go to 5800
if(t+l.e-6.1t.twprt) go to 6000
twprt-twprt+cwprt*delt
5800 continue
do 5900 i- 2,im1
do 5900 j=2, jml
write(9,48) i,j,u(i,j),v(i,j),p(i,j)
5900 continue
c
c set the advance time velocities u and v into the un and vn arrays
c
6000 continue
do 6100 i-1,imax
do 6100 j-1,jmax
un ( i , j ) -u ( i , j )
vn ( i , j ) -v ( i , j )
pn ( i , j ) -p ( i , j )
6100 continue
c
c advance time t-t+delt
c
t-t+delt
if(t.gt.twfin) go to 6500
cycle-cycle+1
go to 1000
6500 continue
call etime(tarray)
cpu2-tarray(1)+tarray(2)
cpu2-cpu2-cpu1
print*, , total CPU time: ',cpu2
call exit(l)
end
c *****************************************************************
c
subroutine eqtn(x,y,uint,vint)
c
c solve trajectory equations dx/dt--u and dy/dt/--v
c and get interpolated velocities at the
c newly computed (x,y) location
c
real u(36,36), v(36,36)
c
integer nf(36,36)
c
common /vel/ u,v
common /mesh/ dt,dx,rdx,xmin,xmax,dy,rdy,ymin,ymax,dd
common /dat/ nf, neqn, imax, jmax
c
c interpolate velocities to (x,y)
c
call ginter(x,y,uii,vii)
c
c NEQN - the number of sub-time steps to solve
c the trajectory equations
c
ddt-dt/neqn
n-neqn
xO-x
yO-y
uO-uii
yO-vii
c
do 10 k-1,n
xl-xO-uO*ddt
y1-yO-vO*ddt
i1-(xO-xmin)*rdx+2.0
jl-(yO-ymin)*rdy+2.0
if(i1.le.O) i1-1
if(i1.gt.imax) i1-imax
if(j1.le.O) j1-1
if(j1.gt.jmax) j1-jmax
ddt1-ddt
m-1
key-O
if(nf(i1,j1).eq.1) go to 40
c
c subdivide the time step further if
c the trajectory hits a blocked cell
c
30 continue
x2-xO
y2-yO
u2-uO
v2-vO
ddt1-ddtl*0.5
m-2*m
do 20 l-l,m
x1-x2-u2*ddt1
y1-y2-v2*ddtl
il-(x2-xmin)*rdx+2.0
j1-(y2-ymin)*rdy+2.0
if(i1.le.0) i1-1
if(il.gt.imax) il-imax
if(jl.le.O) jl-l
if(j1.gt.jmax) j1-jmax
if(nf(i1,j1).eq.1) goto 50
key-key+1
if(key.gt.20) then
uint-u2
vint-v2
goto 100
endif
go to 30
50 continue
if(l.eq.m) goto 40
x2-x1
y2-yl
call ginter(x2,y2,u2,v2)
20 continue
40 continue
xO-xl
yO-y1
c
c interpolate velocities to (xO,yO)
c
call ginter(xO,yO,uO,vO)
10 continue
uint-uO
vint-vO
100 continue
return
end
c**************************************************
c
subroutine ginter(x,y,uint,vint)
c
c interpolated velocities to the point (x,y)
c
real u(36,36), v(36,36)
integer nf(36,36)
common jvelj u,v
common jmeshj dt,dx,rdx,xmin,xmax,dy,rdy,ymin,ymax,dd
common /dat/ nf, ne, imax, jmax
c
c define the indices of the cells which will be
c used for interpolation
c
il-(x-xmin)*rdx
if(il.ge.l) goto 11
il-1
11 continue
if(il.1e.imax) goto 12
il-imax
12 continue
x1-x-(i1-1.0)*dx
if(i2.ge.l) goto 13
i2-1
13 continue
if(i2.1e.imax) goto 14
i2-imax
14 continue
x2-(i2-1.0)*dx-x
j1-(y-ymin)*rdy
if(jl.ge.l) goto 15
jl-l
15 continue
if(jl.le.jmax) goto 16
jl-jmax
16 continue
j2-jl+1
y1-y-(j1-1.0)*dy
if(j2.ge.l) go to 17
j2-1
17 continue
if(j2.le.jmax) goto 18
j2-jmax
18 continue
y2-(j2-1.0)*dy-y
c
iO-2*x1*rdx
jO-2*yl*rdy
i3-il+2*iO
j3-jl+2*jO
c
c linearly interpolate velocities
c
c-abs(-y1+0.5*dy)
d-abs(-yl+dy*(2*jO-0.5))
uint-(d*(u(il,j2)*x2+u(i2,j2)*x1)+
* c*(u(il,j3)*x2+u(i2,j3)*xl))*dd
c-abs(-xl+dx*O.5)
d-abs(-xl+dx*(2*iO-0.S))
vint-(d*(v(i2,jl)*y2+v(i2,j2)*yl)+
* c*(v(i3,jl)*y2+v(i3,j2)*yl))*dd
return
end