This article was downloaded by: [CMERI Central Mechanical Engineering Res. Inst.

]
On: 14 February 2013, At: 07:41
Publisher: Taylor & Francis
Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered
office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK

Numerical Heat Transfer, Part A:

Applications: An International Journal of
Computation and Methodology
Publication details, including instructions for authors and
subscription information:
[Link]

A Fast Numerical Simulation for

Modeling Simultaneous Metal Flow and
Solidification in Thin Cavities Using the
Lubrication Approximation
a b
Alexandre Reikher & Krishna M. Pillai
a
Albany–Chicago, LLC, Fox Point, Wisconsin, USA
b
Department of Mechanical Engineering, University of
Wisconsin–Milwaukee, Milwaukee, Wisconsin, USA
Version of record first published: 19 Nov 2012.

To cite this article: Alexandre Reikher & Krishna M. Pillai (2013): A Fast Numerical Simulation
for Modeling Simultaneous Metal Flow and Solidification in Thin Cavities Using the Lubrication
Approximation, Numerical Heat Transfer, Part A: Applications: An International Journal of
Computation and Methodology, 63:2, 75-100

To link to this article: [Link]

PLEASE SCROLL DOWN FOR ARTICLE

Full terms and conditions of use: [Link]

This article may be used for research, teaching, and private study purposes. Any
substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing,
systematic supply, or distribution in any form to anyone is expressly forbidden.

The publisher does not give any warranty express or implied or make any representation
that the contents will be complete or accurate or up to date. The accuracy of any
instructions, formulae, and drug doses should be independently verified with primary
sources. The publisher shall not be liable for any loss, actions, claims, proceedings,
demand, or costs or damages whatsoever or howsoever caused arising directly or
indirectly in connection with or arising out of the use of this material.
 Numerical Heat Transfer, Part A, 63: 75–100, 2013
Copyright # Taylor & Francis Group, LLC
ISSN: 1040-7782 print=1521-0634 online
DOI: 10.1080/10407782.2012.724331

A FAST NUMERICAL SIMULATION FOR MODELING

SIMULTANEOUS METAL FLOW AND SOLIDIFICATION
IN THIN CAVITIES USING THE LUBRICATION
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

APPROXIMATION

Alexandre Reikher1 and Krishna M. Pillai2

1
Albany–Chicago, LLC, Fox Point, Wisconsin, USA
2
Department of Mechanical Engineering, University of Wisconsin–
Milwaukee, Milwaukee, Wisconsin, USA

A numerical algorithm for modelling steady flow of liquid metal accompanied by solidifi-
cation in a thin cavity is presented. The problem is closely related to a die cast process
and in particular to the metal flow phenomenon observed in thin ventilation channels. Using
the fact that the rate of metal flow in the channel is much higher than the rate of solidifi-
cation, a numerical algorithm is developed by treating the metal flow as steady in a given
time-step while treating the heat transfer in the thickness direction as transient. The flow
in the thin cavity is treated as two dimensional after integrating the momentum and conti-
nuity equations over the thickness of the channel, while the heat transfer is modelled as a
one-dimensional phenomenon in the thickness direction. The presence of a moving solid-
liquid interface introduces non-linearity in the resulting set of equations, and which are solved
iteratively. The location and shape of the solid-liquid interface are found as a part of the sol-
ution. The staggered grid arrangement is used to discretize the flow governing equations and
the resulting set of partial differential equations is solved using the SIMPLE algorithm. The
thickness direction heat-transfer problem accompanied by phase change is solved using a
control volume formulation. The results are compared with the predictions of the commercial
software FLOW3D1 which solves the full three-dimensional set of flow and heat transfer
equations accompanied with solidification. The Reynolds’s lubrication equations accom-
panied by the through-the-thickness heat loss and solidification model can be successfully
implemented to analyze flow and solidification of liquid metals in thin channel during the
die cast process. The results were obtained with significant savings in CPU time.

INTRODUCTION
Global competition for manufacturing superiority has entered a new stage. As
economists predicted for quite some time, there is no a single country or region
which can claim absolute world dominance in manufacturing capabilities. Wide-
spread use of numerical analysis software and free exchange of information allow
engineers around the world to design, analyze, and manufacture new products in
record times. The cast industry is not an exception. Flow, thermal, and distortion
analyses are an integral part of developing die cast process parameters as well as

Received 3 March 2012; accepted 6 August 2012.

Address correspondence to Krishna M. Pillai, Department of Mechanical Emgineering, University
of Wisconsion-Milwaukee, P.O. Box 784, Milwaukee, WI 53211, USA. E-mail: krishna@[Link]

75
 76 A. REIKHER AND K. M. PILLAI

NOMENCLATURE
A surface area of the control volume, m2 ^
w dimensionless velocity in z direction,


C specific heat, J=kg K ¼ W
Lw
H0
H0 height of the channel, m


W reference velocity in z direction,
h fluid height, m  
^
 ¼V
H0
L
h dimensionless fluid height, ¼ Hh0 ^
x dimensionless channel length in x
k coefficient of heat transfer, J=ms K  
direction, ¼ Lx
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

L length of the channel, m ^

dimensionless channel width in y
y  
Lf latent heat of fusion, J=kg
_
direction, ¼ Ly
n normal vector ^
z dimensionless channel height in z
p pressure, Pa

^
 2  direction, ¼ Hz0
p H0 L
dimensionless pressure, L mV
p

q
~ conductive heat flux vector Greek symbols


Re Reynolds number, ¼ qVm L a artificial compressibility
D solidification parameter
S solidification rate, ¼ ks Tq mL
T w
f Ho
m=s g enthalpy of formation, J=kg
s

s location of the solid–liquid interface, m q density of metal, kg=m3

^
s dimensionless location solid–liquid m viscosity, kg= m s


interface, ¼ Hso s reference time, s
d ¼ e2Re Remodified Reynolds number
T temperature of metal, C proposed in reference [1].
t time, s  
e cavity aspect ratio, ¼ HL0
u velocity in x direction, m=s
^
u dimensionless velocity in x direction, Subscripts


¼ Vu
c cold metal
v velocity in y direction, m=s CS controlled surface
^
dimensionless velocity in y direction, CV controlled volume
v
 m molten metal
¼ V
v melt melting
V reference velocity in x and y l liquid
directions, m=s o inlet
w velocity in z direction, m=s s solid

die-cast die design. However, due to an increase in complexity in part design, it takes
longer to go through the complete numerical analyses cycle; in many cases, it takes
several iterations to achieve the desired results.
With the development of faster computers as well as more efficient and accu-
rate numerical approximations, engineers can examine more design options and
achieve better results in a much shorter time. However, in spite of the latest advances
in numerical simulations, detailed examinations of the flow and solidification inside
thin channels remain challenging.
Liquid flow and solidification in channels is a complex phenomenon which has
gained much attention of researchers in the past few decades. Complexity of the
fluid-flow physics and solidification, as well as changes in the flow regime along
the length of the channel, create quite a few challenges in the development of numeri-
cal algorithms to predict the location and shape of the liquid-solid interface, as well
as velocity and temperature distributions in the channel. Detailed descriptions of the
fluid flow, heat transfer, and solidification in the straight channels was conducted by
 MODELLING METAL FLOW AND SOLIDIFICATION 77

Epstein and Chung [2]. The numerical analysis of fluid flow and solidification in
channels requires the solution of the 3-D Navier–Stokes equations. The thin cavities
with high length-to-thickness aspect ratios require a large number of computational
cells in order to achieve accuracy and convergence.
Many numerical algorithms were developed to analyze flow and solidification
between two parallel plates. In order to simplify the 3-D problem, it is reduced into a
2-D one, where the original governing equations are converted from the Cartesian coor-
dinate system into the curvilinear coordinates. The numerical model developed by B.
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

Weigand et al. [3] successfully solved the two-dimensional Navier-Stokes equation

coupled with the heat transfer equation. The numerical analysis of heat transfer during
solidification in a channel, in most cases, requires the conversion of the heat conduction
equation from the Cartesian coordinates into the curvilinear coordinates as well [4].
Though several numerical methods have been proposed to model solidification
of materials in the recent past, not all of them are useful for modeling the flow and
solidification in thin channels. For example, a generalized finite difference method
was shown to be an efficient technique to model the solidification of metals in refer-
ence [5]. However, the model was developed for stationary liquid metals and requires
a numerical procedure for identification and generation of nodes throughout the
computational domain; such a method will be difficult to implement in thin channels
as it will require to generate a large number of computational nodes. Similarly, a
cellular automaton model used for microscopic modeling of heat transfer and
coupled with the finite volume method for macroscopic modeling of solidification
process was introduced by Yao et al. [6]. The model allows for accurate prediction
of the solidification parameters in both macro and micro scales. However, the model
is implemented for stationary fluids only, and requires finite volume descritization in
the direction of solidification. A numerical algorithm for modeling two-phase flow
was proposed in reference [7], where forcing terms are added to the Navier-Stokes
equations to account for the properties variation between the two fluids. The method
doesn’t account for the inertia driven flows and does not predict solidification, and
requires computational descritization of the domain in the direction of the interface
between two fluids during advection. In the last few years, several finite-difference
and finite-volume based methods have been developed and implemented to model
solidification of metals during casting [8–10]. These models concentrate on develop-
ing accurate and efficient algorithms to predict temperature distribution during the
solidification process. Such models require a large number of computational cells to
be generated in the direction of solidification-front motion, and do not account for
the effects of the inertia driven flow on temperature distribution in the liquid phase.
Another approach is to convert the original 3-D governing equations into the
depth-averaged equations [11,12]; this approach is widely used in the shallow-depth
fluid-flow models. However, these algorithms still require the solution of the turbu-
lence models. Several turbulence models have been developed to be used with the
depth-averaged governing equations [12–15].
Although the above-mentioned approaches simplify the governing equations,
the inclusion of turbulence and other details in numerical models requires substantial
amount of CPU time to achieve converged solution. Moreover, conversion of the
governing equations into curvilinear coordinates in references [12–15] creates added
complexities in the development of the numerical algorithm.
 78 A. REIKHER AND K. M. PILLAI

In this article, we propose a simpler, less-computationally expensive approach

where the three-dimensional problem of flow and solidification in a thin channel is
reduced to a two-dimensional one based on the Hele-Shaw approximation [16]. This
approach is based on the Reynolds lubrication theory. The fundamental assumption
of the lubrication equation is that in thin cavities with the flow at relatively small
Reynolds numbers, the inertia forces are much smaller than the viscous forces and
can be neglected. In such a situation, the three-dimensional Navier-Stokes equation
can be reduced to a Reynold lubrication equation and used to analyze flows in thin
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

cavities [17]. In spite of its limitations, the Reynold lubrication formulation remains
the foundation of the numerical analysis in thin cavities.
Owing to the high-speed nature of the die cast process [18], inertia effects in the
metal flow cannot be neglected. Some attempts were made to include the influence of
inertia in the lubrication equation. For example, validity of integration of the gov-
erning equation over a cavity thickness after assuming a parabolic distribution of
the velocity was experimentally confirmed [19]. Similarly, the inertia effects in
thin-channel flows were included in the lubrication equation and validity of the
altered lubrication equation for a wide range of Reynolds numbers was established
[20,21].
In this article, a numerical solution of flow in a thin cavity using the lubrication
approximation along with a control-volume based solidification model will be pre-
sented. The staggered grid arrangement is used to discretize the governing equations.
Then, an iterative SIMPLE algorithm is used to solve the discretized equations for
momentum in the centerline 2-D plane within the channel, while another iterative
scheme is used to model the out-of-plane solidification.
The present article is organized as follows. In section 1, the problem formu-
lation is presented. Section 2 gives a detailed description of the steps required to
reduce the three-dimensional governing equations for flow within the channel into
the two-dimensional governing equations for centerline flow after integrating along
the thickness. The solution procedure is outlined in section 3. Results of the simula-
tion and a numerical validation using the commercial software FLOW3D1 are given
in section 4. Finally, the conclusions are presented in section 5.

An Introduction to the Die Casting Process

Die casting (also known as the die cast process) is a precision manufacturing
process in which molten metal is injected at high pressure and velocity into a perma-
nent metal mold.

Cold Chamber Process

The cold chamber die-cast process is used for alloys with a high melting point,
such as aluminium and brass.
A typical process consists of several steps (Figure 1).

1. Molten metal is ladled into the shot sleeve (Figure 1a).

2. Hydraulic cylinder applies pressure on the plunger (Figure 1b).
 MODELLING METAL FLOW AND SOLIDIFICATION 79
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

Figure 1. Process steps in the cold chamber die-cast process. a) Molten metal is ladled into the shot sleeve;
b) hydraulic cylinder applies pressure on plunger; c) plunger pushes metal from the sleeve through the gat-
ing system into the cavity; d) high pressure is maintained during solidification; e) after solidification is
complete, the die opens; and f) the part is ejected from the cavity (color figure available online).

4. The plunger pushes the metal from the sleeve through the gating system into the
cavity (Figure 1c).
5. High pressure is maintained during the solidification process (Figure 1d).
6. After solidification is complete, the die opens (Figure 1e).
7. The part is ejected from the cavity (Figure 1f).

Essential details of the die-cast die design are shown in Figure 2. Metal fills the cavity
and as metal replaces air, a ventilation channel allows the air to escape.

1. PROBLEM FORMULATION
Metal flow and solidification in a thin channel is a subject of this study. Molten
metal is fed from the left of the channel (see Figure 3a) in positive x direction. Flow
is induced due a pressure difference between the left side (inlet) and right side (outlet)
of the channel. In the present study, it is assumed that a steady flow of metal has
 80 A. REIKHER AND K. M. PILLAI
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

Figure 2. A schematic of a typical die-cast die. The injected metal enters from the left, it fills up the die-cast
die cavity, and then the extra metal spills through the ventilation channel.

Figure 3. Straight channel with a rectangular cross-section. The liquid metal enters from the left-most
section in the y-z plane, flows along the x direction, and then exits from the other end.

been established before the onset of solidification at the walls. (Such an assumption
is justified since the filling of such channels happen within a second.) After a suf-
ficient amount of heat has been extracted from the metal, a solid-liquid interface
formed next to the channel walls grows and meets at the center of the channel. Fluid
flow is assumed to feed the solidification front while the heat is being extracted.
Both the metal and channel are at a superheated temperature initially. The
channel walls are suddenly cooled to a temperature below the solidification tempera-
ture. Solidification fronts will be forming near the walls of the channel, propagating
inside the melted metal. During solidification, the metal is moving under a
pressure-driven flow with a prescribed inlet velocity.
 MODELLING METAL FLOW AND SOLIDIFICATION 81

Table 1. Properties of pure aluminium

Metal density (q) 2700, 2300 [kg m
3]

Heat conduction coefficient (ks, kl) 210 [J m
1 s
1 K
1]
Specific heat (C) 938 [J kg
1 K
1]
Liquid metal viscosity (m) 10
3 [kg m
1 s
1]
Melting point 500 [ C]
Latent heat (Lf) 3.97  105 [J kg
1]
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

Now, we present a comparison of the typical speed with which the metal soli-
difies versus the speed with which the metal passes through the channel—such a com-
parison will help us to ignore solidification during the filling of the channel. A typical
solidification rate inside the horizontal channel can be found [22] as

Tm
Tw 600
10 m
S ¼ ks ) S ¼ 94 ¼ 5:18e
2 ð1Þ
q L f Ho 2710  3:95e5  0:001 s

Using the properties of pure aluminium, shown in Table 1. Meanwhile the typi-
cal rate for metal flow in the channel1 without the solid–liquid interface being present
is 1 ms
1. The characteristic values clearly show that the rate of solidification is much
smaller then the rate of metal flow in the channel. In fact, as the solid-liquid interface
converges at the center of the channel, the rate of metal flow increases and the ratio
of the solidification rate to the flow rate further reduces and goes almost to zero.
Based on these conclusions, we are justified in developing our numerical algorithm
for transient solidification in the channel accompanied by liquid-to-solid heat
conduction, while treating the metal flow to be quasi-steady.
The proposed numerical algorithm is developed based on the assumptions that
at time greater than zero, the liquid metal is entering the channel with its temperature
above the melting point. Due to their low thermal resistance, the channel walls are
assumed to remain at a constant temperature below the melting point to induce sol-
idification. Since the variation in the solid-layer thickness with a position along the
channel length is small, quasi stable-state can be assumed for the heat conduction in
the solid. The liquid-metal temperature is taken to be a constant, while the metal
velocity at the channel entrance is considered to be fully developed and steady. All
physical properties for both the liquid and solid phases are considered constants.
The problem is applicable to the metal die-cast process involving flow in thin
cavities. A brief description of the process is given below. Thin-wall castings, flow
in ventilation channels, etc., are some examples where the proposed algorithm can
be utilized. The proposed numerical solutions can be used to reduce the number
of design iterations employing the full 3-D simulation algorithm. One of the ways
to reduce the computational time is to reduce a 3-D problem to a 2-D one. In the
present case, the die-cast mold cavity is thin and, hence, the flow in the vertical direc-
tion is neglected. Assumption of negligible inertial forces allows one to reduce the
Navier-Stokes equations to the Reynolds lubrication approximation; but, such an

1
The typical velocity corresponds to the end of the die cast process, after the main cavity is filled
and liquid metal is in the ventilation channel.
 82 A. REIKHER AND K. M. PILLAI

approximation is valid only for small Reynolds numbers. A major characteristic of

the die casting process is metal flow under high pressure and velocity, and since the
Reynolds numbers can achieve quite high values, the inertia effects cannot be
neglected. Hence, a modified lubrication approximation, after including the inertial
effect, is employed here to predict the thickness-averaged in-plane flow in the die-cast
mold.
The presented algorithm considers the thickness-averaged 2-D steady-state
flow in a die-cast mold in the in-plane directions and 1-D along-the-thickness tran-
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

sient heat-transfer. A set of nonlinear partial differential equations is developed to

solve for flow, which are then discretized using the finite difference scheme after
employing the SIMPLE algorithm. A direct control-volume based formulation is
proposed to model heat transfer and solidification along the thickness direction.
The two sets of equations are solved in an iterative manner using Matlab, and the
obtained results are validated by comparing with the ones achieved using the com-
mercial software FLOW3D2, where the control volume method is used to discretize
the governing equations while the enthalpy method is used to estimate temperature
distribution in the cavity. Simulation results show that the developed numerical
approximation can be successfully used for predicting velocity and temperature dis-
tributions in a thin cavity with acceptable accuracy and substantial savings in
computational time.

2. GOVERNING EQUATIONS
We are considering a three-dimensional flow in a straight channel with a
rectangular cross-section, as shown in Figure 3. The flow is considered to be incom-
pressible, viscous, and Newtonian. Due to the fact that the rate of metal flow in the
channel is much higher than the rate of solidification, a steady-state flow with tran-
sient heat conduction from liquid into solid is assumed. Since the inertial effects
characterized by high Reynolds numbers are dominant in the flow, so the gravi-
tational forces are neglected. The governing equations are expressed in the Cartesian
coordinate system with x coordinate in the direction of flow (along the cavity
length), y in the direction normal to the flow (along the cavity width), and z in the
direction transverse to the x-y plane (along the cavity height); u, v, and w are the cor-
responding velocities. The governing equations used are the continuity, momentum,
and energy equations in the liquid and solid phases with momentum and energy
boundary conditions specified at the channel walls, inlet, and outlet as well as at
the solid-liquid interface. Location and shape of the solid-liquid interface is found
as a part of the solution of the presented algorithm. The steady-state conservation
equations governing the transport of mass, momentum and energy are expressed
as follows.

2
FLOW3D1 [23] is a general purpose commercial CFD software which solves three-dimensional
fluid-flow and solidification problems using the finite different approximation. FLOW3D utilizes the
volume-of-fluid technique and the FAVOR method to track free surfaces as well as solid-liquid
interfaces. The two equation k-e model is used to resolve the turbulent properties of the flow. The
averaged Navier–Stokes equations coupled with the energy equation allow the software to achieve an
accurate solution for turbulent metal flow undergoing solidification.
 MODELLING METAL FLOW AND SOLIDIFICATION 83

2.1. Liquid Region

Continuity equation

qu qv qw
þ þ ¼0 ð2Þ
qx qy qz

Momentum balance equations

Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

  !
qu qu qu q2 u q2 u q 2 u qp
q u þv þw ¼m 2
þ 2þ 2
ð3aÞ
qx qy qz qx qy qz qx

  !
qv qv qv q2 v q 2 v q2 v qp
q u þv þw ¼m 2
þ 2þ 2
ð3bÞ
qx qy qz qx qy qz qy

  !
qw qw qw q2 w q2 w q2 w qp
q u þv þw ¼m 2
þ 2þ 2
ð3cÞ
qx qy qz qx qy qz qz

Energy balance equation

In differential form, the energy balance equation is expressed as follows.

qT
qC þ qC~
v  rT
r  ~
q¼0 ð4Þ
qt
Solution of the governing equations (2)–(4) presents several problems. To begin with,
the convective terms on the left-hand side of Eq. (3) are nonlinear. All equations are
coupled because velocity components are present in each equation. On comparing
the rate of solidification (see Eq. 1) to the rate of flow, one can define a solidification
parameter as follows.

S
D¼ ð5Þ
u
Based on the flow and solidification characteristic of the presented problem and the
characteristic x-direction velocity of u ¼ 1 m=s, the solidification parameter in our
case will be as follows.

0:0518
D¼ ¼ 0:0518
1:0
In the limit where D ! 0, the advection terms becomes dominant. The temperature
distribution defines the location of the solid-liquid interface, which in its term defines
the liquid domain where the momentum equations are solved.
In order to further simplify the governing equations, we conducted an
order-of-magnitude analysis to determine the importance of each term on the flow
characteristics. Under this, the dimensionless variables were defined as follows.
 84 A. REIKHER AND K. M. PILLAI

^ x ^ y ^ z ^ h T
Tm ^ t
x¼ ; y¼ ; z ¼ ; h¼ ; h¼ ; t ¼
L L Ho Ho Tw
Tm s
 2
^ u ^ v ^ w Lw ^ H0 L ð6Þ
u¼
; v ¼
; w¼
¼
; p¼
p
V V W V Ho L mV
^ s
s ¼
Ho
Owing to a small aspect ratio of the cavity height to its length and width, length and
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

width of the cavity are considered on the same order of magnitude and will be
denoted by L along both x and y. For notational convenience, the tildes were
dropped from nondimensional variables.
On nondimensionalizing the continuity equation, Eq. (2), we get the following.

V qu V qv W qw
þ þ ¼0 ð7Þ
L qx L qy Ho qz

will result in the following.
Multiplying all terms of Eq. (9) by L=V

qu qv W L qw
þ þ
qz ¼ 0 ð8Þ
qx qy Ho V

We have to define the characteristic velocity in z direction. In order to ensure that all
the terms of Eq. (8) are on the same order of magnitude, the characteristic velocity in
z direction is defined as follows.

W L Ho

¼1)W ¼V ð9Þ
Ho V L
Equation (9) indicates that the characteristic velocity in z direction is much smaller

<< V
than those in x and y directions, i.e., W
because Ho<<L. After absorbing this
conclusion, the resultant nondimensional continuity equation, Eq.(8), reduces to the
following.

qu qv qw
þ þ ¼0 ð10Þ
qx qy qz

Substitution of dimensionless variables in the momentum balance equations,

Eq. (3), leads to the following.

  !
H 2V qu qu qu q2 u H 2 q2 u q2 u qp
q o u þv þw ¼ 2 þ 2o þ
ð11aÞ
Lm qx qy qz qz L qx2 qy2 qx

  !
H 2V qv qv qv q2 v H 2 q2 v q2 v qp
q o u þv þw ¼ 2 þ 2o þ
ð11bÞ
Lm qx qy qz qz L qx2 qy2 qy
 MODELLING METAL FLOW AND SOLIDIFICATION 85

  !
H 4V qw qw qw Ho2 q2 w Ho4 q2 w q2 w qp
q o3 u þv þw ¼ 2 2þ 4 2
þ 2
ð11cÞ
L m qx qy qz L qz L qx qy qz

H2
Due to the small cavity aspect ratio, i.e., HLo < < 1, all terms on the order L2o or
higher can be neglected. Then, the in-plane momentum balance equations result in
the following.

 
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

Ho2 V qu qu qu q2 u qp
q u þv þw ¼
ð12aÞ
Lm qx qy qz qz2 qx
 
Ho2 V qv qv qv q2 v qp
q u þv þw ¼ 2
ð12bÞ
Lm qx qy qz qz qy

While the momentum equation in the direction transverse to the flow reduces to the
following.

qp
0¼ ð13Þ
qz
Equation (13) indicates that the fluid pressure is uniform in the z direction regardless
of the inertia effects in the flow and, hence, the pressure is p ¼f (x, y, t) regardless of
the high-Re character of the flow.
Previous work on high-speed flow in thin channels [20] has assumed a para-
bolic distribution of flow velocities. We also will assume a parabolic distribution
of velocity along the x and y directions for further analysis.

u ¼ Uðx; yÞðz2
zhÞ ð14aÞ

v ¼ V ðx; yÞðz2
zhÞ ð14bÞ

On being integrated over the thickness of the channel, the continuity equation, Eq.
(10), becomes the following.

qU qV
þ ¼0 ð15Þ
qx qy

Note, that based on the no-penetration boundary condition on the top and bottom
and the small cavity size in z the direction, w velocity variation is negligible and is set
to zero, i.e., qw
qz ¼ 0.
On nondimensionalizing Eq. (12), the in-plane two-dimensional momentum
equations are expressed as follows.
8

< d U qU þ V qU ¼ q2 U2
qP
qx qy qz qx

 ð16aÞ
: d U qV þ V qV ¼ q2 V2
qP
qx qy qz qy
 86 A. REIKHER AND K. M. PILLAI
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

Figure 4. A typical control volume, defined around the nodes of the mid-level x-y plane, is used to model
the z-direction heat loss and subsequent solidification in the thin cavity.

After integrating Eq. (16a) across the cavity thickness from 0 to h, the momentum
equations become the following.3
   
d 4 qU qU d 2 qh qh 3 qp
h U þV þ U þ UV h ¼ 2U
ð16bÞ
30 qx qy 6 qx qy qx
   
d 4 qV qV d qh 2 qh qp
h U þV þ UV þV h3 ¼ 2 V
ð16cÞ
30 qx qy 6 qx qy qy

For a detailed derivation of the equation set Eq. (16), see Appendix A.

2.2. Energy Equation

A general balance equation for energy is developed for the arbitrary control
volume shown in Figure 4. On integrating the conservative form of Eq. (4) over
the control volume after neglecting variations in fluid properties and velocity, we
get the following.
ZZZ  
q
ðqCT Þ þ rðqCT ~
vÞ
r  ~
q dV ¼ 0 ð17Þ
CV qt

On rewriting the last two terms of Eq. (17) as surface integrals, the energy balance
equation over the fixed control volume changes to the following.
ZZZ ZZ ZZ
q _
qCTdV ¼
v  n dA þ
qCT ~ q  ^ndA
~ ð18Þ
qt CV CS CS

The control volume can be treated as an open system that exchanges heat with its
surroundings, and where mass can flow in and out; hence, Eq. (18) represents the

3
Derivation of Eq. (16) is for two-dimensional variation of the cavity thickness h. The second term
of the equations has a denominator of 6 instead of 12 in reference [12].
 MODELLING METAL FLOW AND SOLIDIFICATION 87

energy balance that can be described as: rate of heat accumulation in control
volume ¼ net rate of heat transport into control volume (by fluid flow)–net rate of heat
transferred out of control volume to surrounding through conduction. Note, that due
to high Peclet numbers involved in this problem, the energy transfer between the fluid
metal and the channel wall, or between the fluid and solidified metal, is driven by con-
vection; the heat transfer through the liquid metal is taken to be purely convective as
well. Therefore, the conduction terms are ignored. (See reference [26] for implemen-
tation of this idea in the corresponding discretized equations.)
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

2.3. Boundary and Initial Conditions

At the initial time t ¼ 0, the same uniform temperature, T ¼ 600 C, is applied to
liquid metal lying within the computational domain. At the entrance (x ¼ 0), the fluid
temperature is set as T ¼ 600 C. At the solid-liquid interface, Ts ¼ Tl ¼ Tmelt. The
wall temperatures at z ¼ 0 and z ¼ h are set to 10 C. The flow is driven by a uniform
velocity imposed at x ¼ 0 location. The outflow boundary condition specified at the
qu
end of the cavity, x ¼ L, is qx ¼ 0. The no-slip boundary conditions are applied at the
walls. Additionally, a no-penetration boundary condition in the form of the velocity
qv
gradient in the direction normal to the wall being equal to zero (i.e., qy ¼ 0), is
applied. An additional condition is needed to express the velocity of the moving
solid-liquid interface as a function of heat transfer in both the solid and liquid
phases. This is called the Stefan condition and can be expressed as follows.

ds qTs qTl
qLf ¼ ks x¼sðtÞ
kl x¼sðtÞ ð19Þ
dt qx qx
In order to establish the validity of Eqs. (15), (16), and (18) that form the gov-
erning equations for the presented problem, they were solved numerically and the
results were compared with the solution of the incompressible Navier–Stokes equa-
tions fully-coupled with the three-dimensional energy equation during solidification
that was solved using the commercial software FLOW3D.

3. SOLUTION PROCEDURE
The system of dimensionless equations (Eqs. (15), (16), and (18)), gives a
complete mathematical formulation of the presented problem of liquid-metal flow
and solidification in a thin channel. The solution involves determination of velocity
and temperature distribution in the liquid phase, as well as the temperature distri-
bution in the solid phase of the thin channel. The governing equations in a liquid
phase are coupled through the interface (Stefan) condition, Eq. (19). The solution
of the Stefan condition gives the location of the solid-liquid interface as a function
of time and position along the length of the channel.
The problem is solved in a straight channel of rectangular cross-section, shown in
Figure 3. A uniform velocity is applied at the x ¼ 0 location to drive the flow. Constant
temperatures are specified at z ¼ 0 and z ¼ h walls, while the walls at y ¼ 0 and y ¼ ymax
are considered adiabatic. Owing to the weak coupling between the momentum and
 88 A. REIKHER AND K. M. PILLAI

energy equations, the temperature distribution within the computational domain can
be solved first. This establishes the location and shape of the solid-liquid interface,
and thus defines the boundaries of the liquid domain. Momentum equations are then
solved using the semi-implicit method for pressure-linked equations (SIMPLE) [24]
procedure, where the momentum and continuity equations are solved in a coupled
manner. The momentum equation, Eq. (16), uses the guessed pressure field and solves
for the preliminary velocities U and V. Then, the modified continuity equation, Eq.
(15), is used to calculate the corrected value of the pressure field.
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

 
qP qU qV
þ a2 þ ¼0 ð20Þ
qt qx qy

Figure 5. Flow chart for the solution algorithm.

 MODELLING METAL FLOW AND SOLIDIFICATION 89

In the usual incompressible form, the mass conservation or the continuity

equation does not have any pressure term. An artificial compressibility term a has
to be added to the modified continuity equation, Eq. (20), to allow for a solution
of P. The usual values are 0.1 > a > 1.2—it is taken to be 0.8 [25] in the present
algorithm.
The solution procedure is described in Figure 5, and can be broken down as
follows.
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

1. Guess pressure values in the first time step.

2. Solve the momentum equation to estimate the preliminary values of the u and v
velocity components.
3. Use the modified continuity equation, Eq. (20), to correct the pressure values.
4. Correct velocities using the new pressure value and continue iterating until the
continuity Eq. (15) is satisfied.

At this point of analysis, height h of the channel in the momentum equations

still remains unknown. In the absence of solidification, h is equal to the channel
height, and the momentum and continuity equations alone will allow us to estimate
the velocity and pressure distributions in the channel. In order to close the system of
governing equations, the energy balance equation has to be solved to find tempera-
ture distribution in the channel. Based on the result of the energy equation and using
the Stefan condition, Eq. (19), location of the solid-liquid interface can be estab-
lished for every particular time step. On knowing the location of the solid-liquid
interface, h value can be updated and used in the momentum equation for the next
time step. It is assumed that there is a perfect contact between the solidified metal
and walls of the channel. Walls of the channel are assumed to have high thermal
mass and conductivity, and therefore, their temperatures remain constant during
the calculation procedure.
Note, that the momentum equations were solved in dimensionless form while
the energy equation was solved in dimensional form. See reference [26] for more
details.

4. RESULTS
Governing equations were solved as indicated in section 3 (solution procedure).
The material properties used in the results presented in this section are of pure
aluminium (see Table 1).
The proposed algorithm is verified for flow and solidification in a straight
channel of a rectangular cross-section (Figure 3). At the initial time step itself, the
flow is considered fully developed. Flow is driven by a uniform axial velocity
imposed at the entrance of the cavity at x ¼ 0. At the time t ¼ 0, metal temperature
is considered to be 600 C and a uniform temperature of 10 C is applied to the top
and bottom of the cavity (z direction). At the inflow boundary, the metal tempera-
ture is set to a constant 600 C. Analyses were run for 1 s. Velocity, temperature dis-
tribution, and location of solid–liquid interface were plotted at three locations.
Velocity u ¼ 1 m=s was applied at x ¼ 0 location. The proposed algorithm was
 90 A. REIKHER AND K. M. PILLAI
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

Figure 6. Grid independence study conducted at z ¼ 0.5 plane. (The cavity width in the y direction was
nondimensionalized as y=L, Eq. (6), after using the length of the cavity as L ¼ 0.1 m. The velocity was ren-

, Eq. (6), after employing the characteristic velocity value of V
dered dimensionless as u=V
¼1 m=s).

verified against results obtain using the commercial software FLOW3D, which simu-
lated a fully-coupled three-dimensional flow analysis with solidification.
Channel (Figure 3a) dimensions are 10  1  0.1 (mm) in the x, y, and z direc-
tions, respectively. Grid independence was insured by comparing 2-D results4 of the
analysis with grid densities 100  10, 200  20, 300  30, 400  40, 500  50, 600  60,
shown in Figure 6. Since the difference between 500  50 and 600  60 results are less
than 0.1%, the analyses were conducted with originally with 500  50 grid. In order
to reduce angularity in the interface-location plots, the mesh density along the thick-
ness z-direction was later taken to be 150 grid points.
The governing equations, Eqs. (15) and (16), were solved using the algorithm
described in the last section. Convergence of the solution was judged by the maximum
change in each variable values during each iteration. The solution was considered
converged, when changes in a dimensionless variables value was less than 10
8.
To verify analyses obtained using the presented algorithm, three-dimensional
flow and solidification solutions from the commercial CFD code FLOW3D were
obtained using the same boundary and initial conditions. Presented results include
fluid velocity, temperature distribution, as well as location of the solid-liquid inter-
face. Three control points along the x direction at dimensionless locations x ¼ 0.2,
x ¼ 0.5, and x ¼ 0.9 were chosen for the plots of z-averaged velocities based on the
solidification patterned observed in the cavity.

4
The mesh densities are for solving the z-averaged velocity fields along x and y directions.
 MODELLING METAL FLOW AND SOLIDIFICATION 91
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

Figure 7. Velocity destribution. a) Velocity distribution at x ¼ 0.2; b) velocity distribution at x ¼ 0.5; and
c) velocity distribution at x ¼ 0.9. (The cavity width in the y direction was nondimensionalized as y=L, Eq.
(6), after using the length of the cavity as L ¼ 0.1 m. The velocity was rendered dimensionless as u=V
, Eq.
(6), after employing the characteristic velocity value of V
¼1 m=s).

Results presented in Figures 7a–7c show velocity variation along the cavity
length where the velocities predicted by our program are compared with the veloci-
ties predicted by FLOW3D. We observe that a fairly close flow-prediction is made by
our simulation based on the lubrication approximation. We also observe that the x-
direction velocity increases with x.
Temperature distribution shown in Figure 8a–8c are plotted at the same loca-
tions as used for Figure 7. Temperature distribution, as it falls below liquidus tem-
perature or the melting point of aluminium (Table 1), suggests the presence of
solid-liquid interface some distance away from the cavity wall. Moving solid-liquid
interface reduces cavity height, and as a result, causes an increase in the melt velocity
(Figure 7) due to conservation of mass. Evolution of the solid-liquid interface along
the channel length is shown using Figure 9a and 9b. We note that some discrepancy
exists between the lubrication approximation solution and the FLOW3D solution in
the beginning. However, we achieve a better convergence of results as the time
increases. The difference in the results may be attributed to the turbulent nature
 92 A. REIKHER AND K. M. PILLAI
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

Figure 8. Temperature distribution along the cavity thickness at time ¼ 1 s; a) Temperature distribution at
x ¼ 0.2; b) temperature distribution at x ¼ 0.5; and c) temperature distribution at x ¼ 0.9. The coordinate in
the z direction was nondimensionalized as z= Ho, Eq. (6), while using Ho ¼ 0.001 m as the cavity thickness.

of the flow employed in the FLOW3D simulation: as the channel height decreases,
turbulence is less prevalent in the flow, and the results predicted by the presented
algorithm are closer to the FLOW3D solution.
At the specified x locations, differences in the position of the solid–liquid inter-
face are significant enough to cause visible velocity differences in Figure 7. As the
above given discussion indicates, velocity changes caused by the reduction in cavity
height h corresponds nicely with the changes in velocity estimated by FLOW3D. All
results are within 10% of the solution obtained by running three-dimensional analy-
ses utilizing the commercial software FLOW3D.
A significant computational advantage is achieved through a dramatic
reduction in CPU time. Owing to the simplification of the governing equations using
the lubrication approximation, the CPU time for the proposed algorithm was
observed to be 20 s. In contrast, the CPU time for the corresponding three-
dimensional analysis with FLOW3D software was 12 min. This 36-fold reduction
in CPU time clearly demonstrates that the proposed algorithm based on reduced
physics is quite fast without a significant sacrifice in the accuracy.
 MODELLING METAL FLOW AND SOLIDIFICATION 93
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

Figure 9. Evolution of the solid–liquid interface with time for u ¼ 1. a) The interface location at t ¼ 0.5s:
and b) the interface location at t ¼ 1 s. (The cavity length in the x direction was nondimensionalized as
x=L, Eq. (6), after using the length of the cavity as L ¼ 0.1 m. The coordinate in the z direction was non-
dimensionalized as z= Ho, Eq. (6), while using Ho ¼ 0.001 m as the cavity thickness).

5. CONCLUSION
Results of the presented analyses indicate that Reynolds lubrication approach
can be successfully implemented to investigate the flow and solidification of the mol-
ten metal in thin cavities during the die cast process. The proposed 2.5-D algorithm
allows one to estimate the thickness-averaged liquid-metal velocity in the plane of
the cavity using the finite difference method. Then, a finite-volume based algorithm
allows one to estimate temperature distribution along the thickness direction as well
as the location of the solid—liquid interface. The numerical simulation based on the
algorithm is verified by comparing its predictions with the solution of the
three-dimensional Navier–Stokes equation fully coupled with a three-dimensional
energy equation, as predicted by the commercial software FLOW3D1. Results indi-
cate that the proposed simulation is fairly accurate in predicting the averaged velo-
city fields, temperatures along the thickness, and gap thicknesses inside the cavity.
Considering small error and significant savings in computational time, the proposed
algorithm can be used to reduce time on the initial stages of process development of
the die-cast process. It will allow to expedite flow analysis of the die casting process,
by using the presented algorithm in the cases where a high aspect ratio of the thin
cavity require a large number of the computational cell to achieve the converged sol-
ution. It can be especially useful in analyzing fluid flow and solidification in venti-
lation channels of the die-cast die.

REFERENCES
1. B. J. Hamrock, Fundamental of Fluid Film Lubrication, McGraw-Hill, New York, 1994.
2. M. Epstein, F. B. Cheung, Complex Freezing-Melting Interfaces in Fluid Flow, Fluid
Mech., vol. 15, pp. 293–319, 1983.
3. B. Weigard, H. Beer, Ice-Formation Phenomena for Water Flow Inside a Cooled Parallel
Channel: An Experimental and Theoretical Investigation of Wavy Ice Layers, Int. J. Heat
Mass Transfer, vol. 36, no. 3, pp. 685–693, 1993.
 94 A. REIKHER AND K. M. PILLAI

4. A. Rouboa and E. Monteiro, Heat Transfer in Multi-Block Grid during Solidification:

Performance of Finite Differences and Finite Volume Method, J. Material Processing
Tech., vol. 204, pp. 541–458, 2008.
5. B. Mochnacki, and E. Majchrzak, Numerical Modelling of Casting Solidification using
Generalized Finite Difference Method, Material Sci. Forum, vols. 638–642, pp.
2676–2681, 2010.
6. X. Yao, M. S. Dargusch, A. K. Dahle, C. J. Davidson, D. H. StJohn, Investigation Into
the Effects of Nucleation Parameters on Grain Formation during Solidification using a
Cellular Automation–Finite Control Volume Metod, J. Materials Res., vol. 23, no. 9,
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

pp. 2312–2325, 2008.

7. H. Montazeri, M. Bussmann, and J. Mostaghimi, Accurate Implementation of Forcing
Term for Two-Phase Flows into SIMPLE Algorithm, Int. J. Multiphase Flow, vol. 45,
pp. 40–52, 2012.
8. Z. Domanski, M. Ciesielski, and B. Mochnacki, Application of Control Volume Method
using the Voronoi Tessellation in Numerical Modelling of Solidification Process, AIP
Conf. Proc., vol. 1220, pp. 17–26, 2010.
9. V. Grozdanic, Finite-Difference Methods for Simulating the Solidification of Castings,
Matriali in Tehnologije, vol. 43, no. 5, pp. 233–237, 2009.
10. A. Rouboa and E. Monteiro, Heat Transfer in Multi-Block Grid during Solidification:
Performance of Finite Difference and Finite Volume Methods, J. Material Processing
Tech., vol. 204, no. 1–3, pp. 451–458, 2008.
11. C. Yang, and C. Li, Application of SIMPLE Algorithm on Nonstaggered Grid to
Numerical Simulation of the Yellow River Flow, and Sediment Transport in Shapotou
Reach, ISWREP 2011 — Proc. of 2011 Inte. Symp. on Water Resource and Environmental
Protection, vol. 1, no. 5893085, pp. 626–629, 2011.
12. M. Larmaei, T. -F. Mahdi, Analysis of SIMPLE Algorithm for Depth-Averaged Simula-
tions, Environmental Hydraulics Proc. of the 6th Int. Symp. on Environmental Hydraulics,
vol. 2, pp. 967–972, 2010.
13. A. K. Rastogi, Predictions of Heat and Mass Transfer in Open Channels, J. the Hydrau-
lics Division, ASCE, vol. 104, no. 3, pp. 397–420, 1978.
14. Y. S. Chen and S. M. Kim, Computation of Turbulent Flow using an Extended k
e
Turbulence Closure Model, NASA, CR—179204, pp. 21–25, 1987.
15. V. Yakhot, S. A. Orszag, S. Thangam, and T. B. Speziale, Development of Turbulence
Models for Shear Flow by a Double Expansion Technique, Phys. Fluids A, vol. 4, no.
9, pp. 1510–1520, 1992.
16. H. S. Hele-Shaw, The Flow of Water, Nature,vol. 58, pp. 34–38, 1898.
17. A. N. Szeri, Tribology: Friction, Lubrication, and Wear, McGraw-Hill, New York, 1980.
18. K. Venkatesan and R. Shivpuri, Numerical Simulation and Comparison with Water
Modeling Studies of the Inertia Dominated Cavity Filling in Die Casting, in S. F. Shen
and P. R. Dawson–Balkema, (eds.), Simulation of Materials Processing: Theory, Methods
and Applications, Springer, Rotterdam, pp. 1203–1211, 1995
19. V. Aleksenko, V. V. Nakoryakov, and B. G. Pokusaev, Wave Formation on a Vertical
Falling Liquid Film, AIChE J, vol. 31, pp. 1446–1458, 1985.
20. M. R. Siddique, and R. Khayat, Influence of Inertia and Topography in Thin-Cavity,
Physics of fluids, vol. 14, no. 5, pp. 1703–1719, 2002.
21. Ch. Karcher and P.H. Steen, High Reynolds Number Flow in a Narrow Gap Driven by
Solidification, Physics of fluids, vol. 13, no. 4, pp. 834–840, 2001.
22. J. K. Carpenter and P. H. Steen, Heat Transfer and Solidification in Planar-Flow Melt-
Spinning: High Wheel Speeds, Int. J. Heat Mass Transfer, vol. 40, no. 9, pp. 1993–
2007, 1997.
23. Flow Science, FLOW3D1, www.fl[Link], 2012.
 MODELLING METAL FLOW AND SOLIDIFICATION 95

24. S. V. Patankar, Numerical Heat Transfer and Fluid Flow, McGraw-Hill, New York, 1980.
25. D. Kwak and C. C. Kiris, Computation of Viscous Incompressible Flows, Springer,
2011.
26. A. Reikher, Numerical Analysis of Die-Casting Process in Thin Cavities using Lubrication
Approximation, Ph.D. thesis, University of Wisconsin Milwaukee, WI, USA, 2012.

APPENDIX A
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

The momentum equations in dimensionless form are as follows.

8

< d U qU þ V qU ¼ q2 U2
qP

qx qy
 qz qx
ðA1Þ
: d U qV þ V qV ¼ q2 V2
qP
qx qy qz qy

On integrating Eq. (A1) over the thickness of the cavity h, we get the following.

8 hR R i R R
< d h qU 2 dz þ h qUV dz þ h qPdz ¼ h q2 U2 dz
0 qx 0 qy 0 qx 0 qz
hR R i R R ðA2Þ
: d h qUV dz þ h qV 2 dz þ h qPdz ¼ h q2 V2 dz
0 qx 0 qy 0 qy 0 qz

Assume a parabolic velocity distribution as follows.

u ¼ Uðx; y; tÞðz2
zhÞ ðA3Þ

v ¼ V ðx; y; tÞðz2
zhÞ

Inserting the expressions of the velocities, Eq. (A3), into Eq. (A2) gives us the
following.

8    
> R h q U 2 ðz2
zhÞ2 R h q UV ðz2
zhÞ2 R h qp R h q2 ðU ðz2
zhÞÞ
>
<d 0 dz þ 0 dz þ 0 qx dz ¼ 0 dz
qx qy qz2
    ðA4Þ
>
> R q UV ðz2
zhÞ 2
R h q V 2 ðz2
zhÞ 2
Rh R h q2 ðV ðz2
zhÞÞ
: d 0h qx dz þ 0 qy dz þ 0 qp
qydz ¼ 0 qz2
dz

Rearranging the terms yields the following.

2
 2 
 2  3
Z h q U 2 z2
zh Z h q UV z2
zh
d4 dz þ dz5
0 qx 0 qy
Z h Z h  
qp q qUz2 qUzh
þ dz ¼
dz ðA5Þ
0 qx 0 qz qz qz
|fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl}
X
 96 A. REIKHER AND K. M. PILLAI

2
 2 
 2  3
Z h q UV z2
zh Z h q V 2 z2
zh
d4 dz þ dz5
0 qx 0 qy
Z h Z h   ðA6Þ
qp q qVz2 qVzh
þ dz ¼
dz
0 qy 0 qz qz qz
|fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl}
XI

Using the chain rule,

Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

qðf ðxÞgðxÞÞ qf ðxÞ qgðxÞ

¼g þf ðA7Þ
qx qx qx

the term X in Eq. (A5) will now be

qUz2
X¼
Uh ðA8Þ
qz

since

qUh
¼0 ðA9Þ
qz

qVh
¼0 ðA10Þ
qz
After inserting velocity values in every term of the Eqs. (A5) and (A6), the
terms X and XI in Eqs. (A5) and (A6) will, respectively, be as follows.

q2 Uz2
X¼ ðA11Þ
qz2

q2 Vz2
XI ¼ ðA12Þ
qz2
2 3
6q 2 4  q  7 qp q2 Uz2
d6 3 2 2
4qx U z
2z h þ z h þ qy UV z
2z h þ z h 5 þ qx ¼ qz
4 3 2 2 7
ðA13Þ
|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
XX XXX

2 3
6q  q 2 4  7 qp q2 Vz2
d6 4 3 2 2
4qx UV z
2z h þ z h þ qy V z
2z h þ z h 5 þ qy ¼ qz
3 2 2 7
ðA14Þ
|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
XXI XXXI

On opening brackets, using the chain rule Eq. (A7) and collecting the terms in
Eq. (A13), leads to the following simplifications.
 MODELLING METAL FLOW AND SOLIDIFICATION 97

q 2 4  qU 2 qðU 2 hÞ qðUhÞ2
XX ¼ U z
2z3 h þ z2 h2 ¼ z4
2z3 þ z2
qx qx qx qx
qU 2 4 qU 2
qh qU 2
qh 2
¼ z
2z3 h
2z3 U 2 þ z2 h2 þ z2 U 2
qx qx qx qx qx
 2 2 2
   ðA15Þ
qU 4 3 qU 2 2 qU 3 2 qh 2 2 qh
¼ z
2z h þz h þ
2z U þ 2z U h
qx qx qx qx qx
2 
qU qh  2 
¼ z4
2z3 h þ z2 h2 þ 2U 2 z h þ z3
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

qx qx

  
qV z4
2z3 h þ z2 h2 4 qV 3 qðVhÞ 2 q Vh2
XXX ¼ U ¼z U
2z U þz U
qy qy qy qy
2
qV qV qh qV qh
¼ z4 U
2z3 hU
2z3 UV þ z2 h2 U þ z2 UV
qy qy qy qy qy
    ðA16Þ
4 qV 3 qV 2 2 qV 2 qh 3 qh
¼ z U
2z hU þz h U þ 2z hUV
2z UV
qy qy qy qy qy
qV  4  qh  
¼U z
2z3 h þ z2 h2 þ 2UV z2 h
z3
qy qy

 2

q 4 3 2 2
 4 qUV 3 qðUVhÞ 2 q UVh
XXI ¼ UV z
2z h þ z h ¼ z
2z þz
qx qx qx qx
qUV 4 3 qUV 3 qh 2 2 qUV 2 qh2
¼ z
2z h
2z UV þz h þ z UV
qx qx qx   qx qx  ðA17Þ
qUV 4 3 qUV 2 2 qUV 3 qh 2 qh
¼ z
2z h þz h þ
2z UV þ 2z UVh
qx qx qx qx qx
qUV  4  qh  
¼ z
2z3 h þ z2 h2 þ 2UV z2 h þ z3
qx qx

Once again, on opening brackets, using the chain rule given in Eq. (A7) and
collecting terms in Eq. (A14) results in the following.

  
qV z4
2z3 h þ z2 h2 4 qV 3 qðVhÞ 2 q Vh2
XXXI ¼ V ¼z V
2z V þz V
qy qy qy qy
2
qV qV qh qV qh
¼ z4 V
2z3 hV
2z3 V 2 þ z2 h2 V þ z2 V 2
qy qy qy qy qy
    ðA18Þ
qV qV qV qh qh
¼ z4 V
2z3 hV þ z 2 h2 V þ 2z2 hV 2
2z3 V 2
qy qy qy qy qy

qV 4  
qh 2 
¼V z
2z3 h þ z2 h2 þ 2V 2 z h
z3
qy qy

Inserting the results of Eqs. (A15)–(A18) into Eqs. (A13) and (A14) will result
in the following.
 98 A. REIKHER AND K. M. PILLAI

8 " 2     4 #
4 3
>
>
qU
qx z
2z h þ z2 h2 þ 2U 2 qx
qh 2
z h þ z3 þ U qV z
2z3 h þ z2 h2 2
Uz2
>   qy
¼ q qz qp

qx
< d þ2UV qh z2 h þ z3
>
" qy
 4     4 #
>
>
qUV
z
2z 3
h þ z 2 2
h þ 2UV qh 2
z h þ z 3
þ V qV
z
2z 3
h þ z 2 2
h
>d qx qx qy 2
Vz2

qp
: þ2V 2 qh z2 h þ z3 
> ¼ q qz qy
qy

ðA19Þ
8 2 3
>
>


Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

>
> 6 2  4   2 7
>
> d4 qU þ U qV z
2z3 h þ z2 h2 þ 2U 2 qx qh
þ 2UV qh z h þ z3 5
>
> qx qy
|fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} qy
|fflfflfflfflfflffl{zfflfflfflfflfflffl}
>
>
>
> I II
>
> 2
>
> q Uz 2
>
> ¼ qp

qx
>
> qz
>
< |fflfflffl{zfflfflffl}
2 III 3 ðA20Þ
>
>


>
> 6 qUV qV 2  4   2 7
>
> d4 qx þ qy z
2z3 h þ z2 h2 þ 2UV qx qh
þ 2V 2 qh z h þ z3 5
>
> |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} qy
|fflfflfflfflfflffl{zfflfflfflfflfflffl}
>
>
>
> IX IIX
>
>
>
> q 2
Vz 2
>
> qp
> ¼ qz
qy
>
>
: |fflffl{zfflffl}
IIIX

Let us integrate the following terms in Eq. (A20).

Integration of the term I results in the following.
 2 Z h  2 Z h  2 Z h
qU qV 4 qU qV 3 2 qU qV
þU z dz
2h þU z dz þ h þU z2 dz
qx qy 0 qx qy 0 qx qy 0
 2   2   2 
qU qV z5 h qU qV z4 h 2 qU qV z3 h
¼ þU
2h þU þh þU
qx qy 5 0 qx qy 4 0 qx qy 3 0
 2  5  2  4  2  3
qU qV h qU qV h qU qV h
¼ þU
2h þU þ h2 þU
qx qy 5 qx qy 4 qx qy 3
 2   2   2 
qU qV h5 qU qV h5 qU qV h5
¼ þU
þU þ þU
qx qy 5 qx qy 2 qx qy 3
 2 
1 qU qV 5
¼ þU h
30 qx qy
ðA21Þ

Integration of the term IX yields the following.

 Z h  Z h  Z h
qUV qV 2 4 qUV qV 2 3 2 qUV qV 2
þ z dz
2h þ z dz þ h þ z2 dz
qx qy 0 qx qy 0 qx qy 0
     
qUV qV 2 z5 h qUV qV 2 z4 h qUV qV 2 z3 h
¼ þ
2h þ þ h2 þ
qx qy 5 0 qx qy 4 0 qx qy 3 0
 MODELLING METAL FLOW AND SOLIDIFICATION 99

     
qUV qV 2 h5 qUV qV 2 h4 qUV qV 2 h3
¼ þ
2h þ þ h2 þ
qx qy 5 qx qy 4 qx qy 3
 2
 5  2
 5  2
 5
qUV qV h qUV qV h qUV qV h
¼ þ
þ þ þ ðA22Þ
qx qy 5 qx qy 2 qx qy 3
 2

1 qUV qV
¼ þ h5
30 qx qy
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

Integration of the term II leads to the following.

 Z h  Z h Z h 
qh qh qh qh
2 U 2 þ UV ðz2 h
z3 Þdz ¼ 2 U 2 þ UV z2 hdz
z3 dz
qx qy 0 qx qy 0 0
  3 4 h
   4 4

qh qh z h z qh qh h h
¼ 2 U 2 þ UV h
¼ 2 U 2 þ UV

qx qy 3 0 4 0 qx qy 3 4
 
1 qh qh 4
¼ U 2 þ UV h
6 qx qy
ðA23Þ

Integration of the term IIX results in the following.

 Z h  Z h Z h 
qh 2 qh 2 3 qh 2 qh 2 3
2 UV þV ðz h
z Þdz ¼ 2 UV þV z hdz
z dz
qx qy 0 qx qy 0 0
  3     
qh qh z h z4 h qh qh h4 h4
¼ 2 UV þ V2 h
¼ 2 UV þ V2

qx qy 3 0 4 0 qx qy 3 4
 
1 qh qh
¼ UV þ V2 h4
6 qx qy
ðA24Þ

Integration of term III yields the following.

  Z h   
q qUz2 q qUz2 qUz2 h
¼ dz ¼ ¼ 2Uz ¼ 2Uh ðA25Þ
qz qz 0 qz qz qz 0

Integration of term IIIX results in the following.

  Z h   
q qVz2 q qVz2 qVz2 h
¼ dz ¼ ¼ 2Vz ¼ 2Vh ðA26Þ
qz qz 0 qz qz qz 0

The result of the integration of Eq. (A1) is given as follows.

8

<
d qU h4 þ d qU 2 h6 þ 1 U 2 qh h5 þ 1 UV qh h5 þ qUV h6 ¼ 2Uh2
h2 qp
6 qt qx 30 6 qx 6 qy qy 30 qx

 ðA27Þ
:
d qV h4 þ d qUV h6 þ 1 UV qh h5 þ 1 V 2 qh h5 þ qV 2 h6 ¼ 2Vh2
h2 qp
6 qt qx 30 6 qx 6 qy qy 30 qy
 100 A. REIKHER AND K. M. PILLAI

Dividing terms of Eq. (A27) by h2 results in the following.

8
2 4 
< d qU h þ 1 U 2 qh h3 þ 1 UV qh h3 þ qUV h4 ¼ 2U
qp

qx 30 6 qx 6 qy qy 30
 qx
ðA28Þ
: d qUV h4 þ 1 UV qh h3 þ 1 V 2 qh h3 þ qV 2 h4 ¼ 2 V
qp
qx 30 6 qx 6 qy qy 30 qy

The final lubrication approximation equation with inertia forces included will
be as follows.
Downloaded by [CMERI Central Mechanical Engineering Res. Inst.] at 07:41 14 February 2013

8


< d 4
U qU þ qU
þ d 2 qh
þ qh 3 qp
30 h V U UV qy h ¼ 2U
qx

qx qy
 6
qx  ðA29Þ
: d h4 U qV þ V qV þ d UV qh þ V 2 qh h3 ¼ 2 V
qp
30 qx qy 6 qx qy qy