Numerical Simulation of the Filling and Curing Stages in Reaction Injection Moulding


Master's Thesis, 2007
120 Pages

Free online reading

Universidade de Aveiro
2007
Departamento de Engenharia Mecânica
Rui Alexandre
Aires da Trindade
Igreja
Numerical Simulation of the Filling and Curing
Stages in Reaction Injection Moulding, using CFX
Simulação Numérica das Fases de Enchimento e
Cura em Reaction Injection Moulding, usando o CFX

Universidade de Aveiro
2007
Departamento de Engenharia Mecânica
Rui Alexandre
Aires da Trindade
Igreja
Numerical Simulation of the Filling and Curing
Stages in Reaction Injection Moulding, using CFX
Simulação Numérica das Fases de Enchimento e
Cura em Reaction Injection Moulding, usando o CFX
Dissertação apresentada à Universidade de Aveiro para cumprimento dos
requisitos necessários à obtenção do grau de Mestre em Engenharia
Mecânica
O trabalho conducente a esta dissertação foi realizado no âmbito do
Projecto de Investigação POCTI/EME/39247/2001, financiado pela
Fundação para a Ciência e Tecnologia e co-financiado pelo FEDER.

palavras-chave
Reaction Injection Moulding (RIM), Enchimento de Moldes, Simulação
Numérica, Escoamentos Bifásicos, (ANSYS) CFX.
resumo
Os métodos habitualmente utilizados para a simulação de injecção em
moldes envolvem um número considerável de simplificações, originando
reduções significativas do esforço computacional mas, nalguns casos
também limitações. Neste trabalho são efectuadas simulações de Reaction
Injection Moulding (RIM) com o mínimo de simplificações, através da
utilização do software de CFD multi-objectivos CFX, concebido para a
simulação numérica de escoamentos e transferência de calor e massa.
Verifica-se que o modelo homogéneo para escoamentos multifásicos do
CFX, geralmente considerado o apropriado para a modelação de
escoamentos de superfície livre em que as fases estão completamente
estratificadas, é incapaz de modelar correctamente o processo de
enchimento. Este problema é ultrapassado através da implementação do
modelo não homogéneo juntamente com a condição de fronteira de
escorregamento livre para o ar.
A reacção de cura é implementada no código como uma equação de
transporte para uma variável escalar adicional, com um termo fonte. São
testados vários esquemas transitórios e advectivos, com vista ao
reconhecimentos de quais os que produzem os resultados mais precisos.
Finalmente, as equações de conservação de massa, quantidade de
movimento, cura e energia são implementadas conjuntamente para simular
os processos simultâneos de enchimento e cura presentes no processo
RIM. Os resultados numéricos obtidos reproduzem com boa fidelidade
outros resultados numéricos e experimentais disponíveis, sendo
necessários no entanto tempos de computação consideravelmente longos
para efectuar as simulações.

keywords
Reaction Injection Moulding (RIM), Mould Filling, Numerical Simulation,
Two-Phase Flow, (ANSYS) CFX.
abstract
Commonly used methods for injection moulding simulation involve a
considerable number of simplifications, leading to a significant reduction of
the computational effort but, in some cases also to limitations. In this work,
Reaction Injection Moulding (RIM) simulations are performed with a
minimum of simplifications, by using the general purpose CFD software
package CFX, designed for numerical simulation of fluid flow and heat and
mass transfer.
The CFX's homogeneous multiphase flow model, which is generally
considered to be the appropriate choice for modelling free surface flows
where the phases are completely stratified and the interface is well defined,
is shown to be unable to model the filling process correctly. This problem is
overcome through the implementation of the inhomogeneous model in
combination with the free-slip boundary condition for the air phase.
The cure reaction is implemented in the code as a transport equation for an
additional scalar variable, with a source term. Various transient and
advection schemes are tested to determine which ones produce the most
accurate results.
Finally, the mass conservation, momentum, cure and energy equations are
implemented all together to simulate the simultaneous filling and curing
processes present in the RIM process. The obtained numerical results
show a good global accuracy when compared with other available
numerical and experimental results, though considerably long computation
times are required to perform the simulations.

i
Table
of
contents
Resumo
Abstract
Table of contents
i
Nomenclature
iii
1 Introduction
1
1.1 Overview of the RIM process
1
1.2 Literature review
4
1.2.1 Injection moulding
4
1.2.2 Reaction Injection Moulding
5
1.3 The 2½D Hele-Shaw flow approach
6
1.4 Motivation and objectives
8
2 Simulation of the filling process
11
2.1 The homogeneous model
11
2.1.1 Physical model
11
2.1.2 Numerical model
12
2.1.3 Case study
16
2.2 The inhomogeneous model
28
2.2.1 Physical model
28
2.2.2 Numerical model
31
2.2.3 Case study
33
2.3 Conclusions
41

Table of contents
ii
3 Simulation of the curing process
43
3.1 The cure reaction
43
3.2 Physical model
44
3.3 Numerical models
44
3.4 Case study
46
3.5 Conclusions
51
4 Simulation of the RIM filling and curing stages
53
4.1 The energy equation
53
4.1.1 Physical model
53
4.1.2 Numerical model
54
4.2 Case study 1
55
4.2.1 Case description
55
4.2.2 Numerical issues
57
4.2.3 Results
61
4.3 Case study 2
73
4.3.1 Case description and numerical issues
73
4.3.2 Results
75
4.4 Case studies 3 and 4
80
4.4.1 Cases description
80
4.4.2 Numerical issues
81
4.4.3 Results
83
4.5 Conclusions
91
5 Conclusions
93
References
99

iii
Nomenclature
Symbols:
A Area
[m
2
]
A
Constant in the viscosity equation
[­]
A
DE
Interfacial area per unit volume
[m
-1
]
A
P
Pre-exponential factor in the viscosity equation
[kg/(m
s)]
A
1
Pre-exponential factor in the cure equation
[s
-1
]
A
2
Pre-exponential factor in the cure equation
[s
-1
]
B
Constant in the viscosity equation
[­]
cs Control
volume
cs
Surface of the control volume
C
Degree of cure
[­]
C
D
Drag
coefficient
[­]
C
g
Solidification (gel) point
[­]
Cp
Specific heat at constant pressure
[J/(kg
K)]
Cr Courant
number
[­]
Cv
Specific heat at constant volume
[J/(kg
K)]
d
DE
Interface length scale
[m]
Dh Hydraulic
diameter
[m]
e
Specific internal energy
[J/kg]
E
P
Viscosity activation energy
[J/mol]
E
1
Reaction activation energy
[J/mol]
E
2
Reaction activation energy
[J/mol]
f Face
F Force
[N]
Fg
JJG
Gravity force vector
[N]
Fw
Wall shear force
[N]

Nomenclature
iv
g
G
Gravity acceleration vector
[m/s
2
]
h
Mould half thickness
[m]
h Specific
enthalpy
[J/kg]
htot
Specific total enthalpy
[J/kg]
H Mould
thickness
[m]
k Thermal
conductivity
[W/(m
K)]
k
1
Parameter in the cure equation
[s
-1
]
k
2
Parameter in the cure equation
[s
-1
]
L Mould
length
[m]
Lf
Position of the flow front
[m]
m
Constant in the cure reaction
[­]
m Mass
[kg]
m
Mass flow
[kg/s]
M
DE
JJJJJG
Volumetric interface momentum transfer
[N/m
3
]
n
Constant in the cure reaction
[­]
n
G
Outward normal surface vector
[­]
N Shape
function
P, p
Pressure
[Pa]
P Perimeter
[m]
Q
Volumetric flow rate
[m
3
/s]
R
Q
Volumetric heat generation rate
[W/m
3
]
Q
t
Total volumetric heat of reaction
[J/m
3
]
r Volume
fraction
[­]
R
Universal gas constant
[J/(mol
K)]
R
JG
Vector from the upwind node to ip
[m]
Re Reynolds
number
[­]
S Cross
section
S Fluidity
[m
s]
S Surface

Nomenclature
v
S
C
Rate of cure
[s
-1
]
S
Cp
Linear source coefficient in the cure equation
[s
-1
]
S
Cv
Source value in the cure equation
[s
-1
]
S
E
Volumetric heat source
[W/m
3
]
S
Ep
Linear source coefficient in the energy equation
[W/(m
3
K)]
S
Ev
Source value in the energy equation
[W/m
3
]
M
S
JJJG
Volumetric momentum source
[N/m
3
]
t Time
[s]
t
f
Filling
time
[s]
t
R
Time of residence
[s]
T Period
[s]
T Temperature
[K]
Ti Initial
temperature
[K]
u
Velocity x component
[m/s]
U
JG
Velocity vector
[m/s]
v
Velocity y component
[m/s]
V Volume
[m
3
]
w
Velocity z component
[m/s]
W Mould
width
[m]
W
v
Volumetric viscous work rate
[W/m
3
]
x x
direction
[m]
y y
direction
[m]
z z
direction
[m]
Greek symbols:
E
Blend factor for the advection term discretization
[­]
J
Shear rate
[s
-1
]
G
Identity matrix

Nomenclature
vi
't
Time step
[s]
'T
ad
Adiabatic temperature rise
[K]
'x
Length of the mesh elements
[m]
T
Quotient between the old and the new time step
[­]
P
Viscosity
[kg/(m
s)]
U
Density [kg/m
3
]
W
w
Wall shear stress
[N/m
2
]
I
General variable
Subscripts:
air Air
in Inlet
ip Integration
point
max Maximum
value
n Node
out Outlet
ref Reference
resin Resin
up Upwind
node
w Wall
x
x component of the vector
y
y component of the vector
z
z component of the vector
D Phase
D
E Phase
E

Nomenclature
vii
Superscripts:
CFX
Results from CFX
hom Homogeneous
model
inhom
Inhomogeneous model
0
Old time level
00
Time before the old time level
Operators:
x
Inner product
Tensor product
Nabla vector
T
X
Transpose of X
Others:
X
Average value of X

Nomenclature
viii

1
1 Introduction
1.1 Overview of the RIM process
The application of synthetic polymeric materials, which are nowadays commonly used in
all of the major market sectors, experienced a dramatic expansion in the second half of the
20
th
century. Due to the continuous progress made in polymeric engineering, these
materials can be synthesized to meet a wide range of mechanical, chemical, optical and
electrical properties. Their low density, resulting in a relatively high specific strength and
stiffness, their aptitude to be manufactured into complex shaped parts and, more important,
their ability to be integrated in automatic mass production processes, which together with
the relatively low cost of raw material, lead to the low cost of the final product, make these
materials very attractive, specially from an economical point of view.
Amongst the various polymer processing techniques, injection moulding is
certainly one of the most important, representing approximately one third of all
manufactured polymeric parts [1, 2]. Extrusion represents approximately another third
(32% and 36% of weight, respectively [1]).
Synthetic polymers can be classified in two major categories: thermoplastic
polymers and thermosetting polymers. For thermoplastics, by far the largest volume (about
80% of the raw material used in injection moulding [3]), the long molecules are not
chemically joined, and thus they can be melted by heating (typically 200-350 ºC [3]),
solidified by cooling and remelted repeatedly. Thermosets are initially made of short chain
molecules, and upon heating or mixing with appropriate reagents, they undergo an
irreversible chemical reaction which causes the short chain molecules to bond. This
process leads to the formation of a rigid three-dimensional structure, which once formed
will not remelt by heating [2, 4].
In Thermoplastics Injection Moulding (TIM), the hot polymer melt is pushed at
high pressure into a cold cavity where it undergoes solidification by cooling down below
the glass transition temperature. Although this process is extensively used for non-specific
and undemanding applications, it presents some weaknesses which limit its use for more
technical parts. Thermoplastics usually have lower mechanical properties than
thermosetting polymers. Their relatively high viscosity (usually exceeding 100 kg/(m
s)
[2]) requires high injection pressures (typically 50-200 MPa [3]), limiting the dimensions
of parts to typically 1 m
2
. High viscosity also limits the use of reinforcements, which are
necessary to meet more demanding requirements, and the production of parts with complex

Introduction
2
shapes. In order to overcome these limitations, reactive moulding techniques have been
developed for thermosetting polymers [2].
Reactive moulding is quite different from conventional TIM, since it uses
polymerization in the mould, instead of cooling, to form the stiff solid parts.
Reaction Injection Moulding (RIM) is a process for rapid production of complex
parts through the mixing and chemical reaction of two or more components. The liquid
components, usually isocyanate and polyol, held in separated temperature controlled tanks,
feed to metering units. When injection begins, the valves open and the components flow at
moderate to high pressures, typically between 10 and 20 MPa, into a mixing chamber. The
streams are intensively mixed by high velocity impingement and due to the shape of the
mixing chamber, and the mixture begins to polymerize, or cure, as it flows into the mould
cavity. As the mixture is initially at low viscosity, low pressures, less than 1 MPa, are
needed to fill the mould, typically in less than five seconds. Inside the mould, cure occurs,
forming the polymer, solidifying and building up enough stiffness and strength such that
the mould can be opened and the part removed, often in less than one minute. Postcuring
may be necessary [5, 6]. Fig. 1 [7] shows the diagram of a typical RIM process.
Fig. 1: Diagram of a typical RIM process [7].

Introduction
3
During RIM's early years, its main use was for high density rigid polyurethane
foam parts for the automotive industry, such as bumpers and fascia. Over the years it has
developed uses in many different areas. The end use applications for RIM are nowadays so
varied that polyurethane can be found in forms such as protective coatings, flexible foams,
rigid foams and elastomers. Although polyurethanes comprise the majority in RIM
processing, other types of chemical systems can be processed, for instance polyureas,
nylons, dicyclopentadienes, polyesters, epoxies, acrylics and hybrid urethane systems [8].
Due to the low viscosity RIM liquids, ranging from 0.1 to 1 kg/(m
s), large parts
can be manufactured with relatively small metering machines, complex shapes, with
multiple inserts, can be fabricated and low pressures can be used to fill the moulds, and
these may be constructed from light-weight materials often at lower costs than for TIM.
Mould clamping forces are much lower, requiring lower cost presses. This opens up short
production runs, and even prototype applications [5]. Low viscosity also opens options in
reinforcement as Structural Reaction Injection Moulding (SRIM), and Reinforced Reaction
Injection Moulding (RRIM) [8].
RIM temperatures are typically lower than those for TIM, with less energy
demands [5]. However, handling of reactive, and often hazardous liquids, requires special
equipments and procedures. As some components freeze at room temperature, a
temperature controlled environment is required for their shipping and storage, increasing
the costs. Due to the low viscosity, it is difficult to seal moulds, increasing leakage and
flash. As low viscosity liquids penetrate mould surfaces, there is the need for release
agents, which has been a major problem for high RIM production.
If flow into the mould is too rapid, air may be entrained and large bubbles appear in
the final part, which is perhaps the greatest cause of scrap. Also, due to low pressure
during filling, it is difficult to remove air from pockets formed behind inserts and from
corners. Thus, vent location becomes extremely important. Some of these problems can be
prevented if moulds are filled slowly, but this can lead to short shots.
As asserted by [9], "Successful moulders must use RIM long enough to learn all the
peculiarities of chemical manipulation, mould manufacture, and processing parameters.
These typically differ for each project, resulting in a long learning curve for the few
companies that choose to offer RIM-produced products". This is why a better
understanding of the injection and curing processes is important. Numerical simulations of
these processes may be helpful to choose the chemical components, to design the mould,
its position and vents location, and to design the process itself: shot time and inlet and
mould temperatures.

Introduction
4
1.2 Literature review
1.2.1 Injection moulding
As mentioned in [3] and [10], the first attempts to study the filling stage in injection
moulding were reported by Spencer and Gilmore [11]. They visually studied the filling of
the mould and derived an empirical equation to determine the filling time. Other important
early flow visualization and tracer studies include the contributions from [12] and [13].
Numerical simulations applied to injection moulding have basically started in the
early 1970s. According to [2] and [14], these first developments were applied to the filling
stage in simple geometries [15­23]. Only tubular, circular and rectangular shapes were
considered, allowing the flow to be accurately assumed as unidirectional. The temperature
field was two-dimensional, one coordinate in the flow direction and the other in the
thickness direction, leading to the so-called 1½D approach. The injected polymer was
assumed to be a Newtonian fluid, and finite difference techniques were used to numerically
solve the set of balance equations.
In [15], one-dimensional flow analysis was coupled with a heat balance equation
for a rectangular cavity. The work presented in [16, 17] dealt with the filling of a disc
mould, using the assumption of a radial creeping flow, and [21] modelled one-dimensional
tubular flow of polymer melts.
In order to expand the previous approaches to more realistic geometries, conformal
mapping [24­26] or decomposition of complex shape cavities in a number of simple
elements [27, 28] were used to extend the 1½D approach to more complex flow situations.
However, these methods lack sufficient generality to be satisfactory and the solution
accuracy strongly depends on how the geometry is partitioned, requiring astute judgment
from the user.
The real breakthrough came with the development of a general 2½D approach,
originally proposed in [29], combining finite elements along the midsurface of the cavity
with finite differences along the thickness direction. The pressure field was solved in two
dimensions by finite element method and the temperature and velocity fields were solved
in three dimensions by means of a mixed finite element / finite difference method.
However, based on the Hele-Shaw approximation, the 2½D approach was unable to
represent the complex flow kinematics of the flow front region, the so-called fountain
flow, see Fig. 2 [30], first reported in [31]. The description of this phenomenon was
addressed by many authors by means of experimental [32­36], analytical [37­39] and
numerical [33, 40í45] methods, leading to approximate models able to capture its basic
flow kinematics without resolving the complex 3D flow details.

Introduction
5
(a) (b)
Fig. 2: Kinematics of fountain flow.
(a) Reference frame of mould; (b) Reference frame of the moving flow front. [30]
Different techniques have been used to handle the time-dependency of the flow
domain during filling. One solution consisted in the use of the control volume method
[46-48], while alternative solutions included the use of boundary fitted coordinates [49, 50]
or the use of a front tracking and remeshing techniques [51­53].
To date several commercial and research three-dimensional simulation programs
for injection moulding have been developed. In particular, [54] developed a three-
dimensional finite-element code for predicting the velocity and pressure fields governed by
the generalized Navier-Stokes equations, [55] analysed the three-dimensional mould filling
of an incompressible fluid and the shape of the fountain flow front, [56, 57] incorporated
the polymer compressibility, by treating its density as a function of temperature and
pressure, in a three-dimensional mould filling process.
1.2.2 Reaction Injection Moulding
As both thermoplastics and reaction injection moulding are basically governed by the same
flow kinematics, most of the developments made in the simulation of the filling stage of
thermoplastics are applicable to reaction injection moulding.
Early studies were dedicated to the static analysis of heat transfer and cure, on the
assumption that the curing stage could be decoupled from the filling stage [58, 59]. More
realistic models were obtained by extending the 1½D approach to reaction injection
moulding [39, 60­66]. The work of Castro and Macosko [39], presenting experimental and
numerical results was of particular importance, it was considered by many authors as a
benchmark solution (it is used also in this work, in Section 4.2 and Section 4.3, as a
benchmark solution).
Extensions to the 2½D approach were only reported more recently [2, 67­71], and
3D simulations have already been reported in some works [72­76].

Introduction
6
1.3 The 2½D Hele-Shaw flow approach
Despite many commercial codes' ability to perform 3D injection moulding simulations, the
2½-dimensional Hele-Shaw approach remains the standard numerical framework for
simulation of injection moulding. Therefore, it becomes important to succinctly describe
its assumptions and formulation, and to discuss its advantages and limitations relatively to
a full three-dimensional approach.
Injection moulding is generally characterized by the part thickness being much
smaller than its overall dimensions and by the high viscosity of the polymer (the latter is
not entirely valid for RIM resins), resulting in low ratios of the inertia force to the viscous
forces (characterized by low Reynolds numbers). In the Hele-Shaw flow formulation, the
inertia effect and the velocity component and thermal convection in the gap-wise direction
are neglected. Moreover, due to small thickness of the part, the velocity gradient in the
gap-wise direction is considered to be much larger than in the other directions. Thus,
denoting the planar coordinates by x and y, the gap-wise coordinate by z, and the velocity
components by u, v and w, respectively, the continuity and momentum equations are
reduced to [77, 78]:
u
v
0
t
x
y
w U
w U
wU
w
w
w
(1)
p
u
0
x
z
z
§
·
w
w
w
P
¨
¸
w
w
w
©
¹
(2)
p
v
0
y
z
z
§
·
w
w
w
P
¨
¸
w
w
w
©
¹
(3)
The boundary conditions for u and v are [77]:
u
v
0
at z
h
on mold wall
(
)
r
(4)
u
v
0
at z
0
on middle plane
z
z
w
w
w
w
(
) (5)
As the pressure is independent of the z direction, by integrating twice equations
(2) and (3), and taking into account the boundary conditions, one arrives to:
h
z
p
z
u
dz
x
w
w
P
³
(6)
h
z
p
z
v
dz
y
w
w
P
³
(7)

Introduction
7
By integrating equation (1) between z = 0 and z = h, the continuity and momentum
equations are merged into a single Poisson-like equation in terms of pressure and fluidity
[14, 77]:
h
0
p
p
dz
S
S
0
t
x
x
y
y
§
·
§
·
w
w
w
w
w
§
·
U
¨
¸
¨
¸
¨
¸
w
w
w
w
w
©
¹
©
¹
©
¹
³
(8)
where S is the fluidity, expressed as [77]:
2
h
0
z
S
dz
U
P
³
(9)
The two above equations could be even further simplified if assuming constant
polymer density.
The gap-wise direction averaged velocities may be obtained as:
2
h
0
z
dz
p
u
x
h
w
P
w
³
(10)
2
h
0
z
dz
p
v
y
h
w
P
w
³
(11)
Besides the velocity and heat convection in the gap-wise direction being omitted,
the heat conduction in the flow directions is assumed to be negligible when compared with
that in the gap-wise direction [14, 78]. The energy equation is then simplified to [77]:
2
R
T
T
T
T
Cp
u
v
k
Q
t
x
y
z
z
§
·
§
·
w
w
w
w
w
U
P J
¨
¸
¨
¸
w
w
w
w
w
©
¹
©
¹
(12)
The term (
P J
2
) represents the heat generated by viscous dissipation, and
R
Q the rate of
heat generation due to chemical reaction.
J is the shear rate, which, according to the
assumptions mentioned above, is expressed as [79]:
2
2
u
v
z
z
§
·
§
·
w
w
J
¨
¸
¨
¸
w
w
©
¹
©
¹
(13)
or:
2
2
z
p
p
x
y
§
·
§
·
w
w
J
¨
¸
¨
¸
P
w
w
©
¹
©
¹
(14)

Introduction
8
However, for typical RIM situations, the viscous dissipation term in the energy equation
can be neglected when compared to the heat generation term due to the cure reaction [39,
80].
Equations (8) and (12) for pressure and temperature, are the basic equations of the
Hele-Shaw approximation for mould filling in thin cavities [78]. Equation (8) is commonly
solved by the finite-element method and equation (12) by the finite-difference method
[14].
Due to the Hele-Shaw formulation assumptions and simplifications mentioned
above, the usage of computational storage and CPU time can be considerably reduced
when compared with the case of a full three-dimensional (unsteady) formulation. But,
although its successful applications to the injection moulding process, it has its limitations
mostly owing to those assumptions [14].
The inertia and three-dimensional effects may become significant enough to
influence the flow, especially in complex thick parts, in situations of branching flow,
where the part thickness changes abruptly, or in regions around special features such as
bosses, corners and ribs. For the RIM process, because of resins small viscosities, the fluid
inertia and the gravity force cannot be omitted [14].
At the filling front, the fluid moves away from the centre, spilling out like a
fountain to the walls [5] (thus the designation "fountain flow"), and therefore the flow is as
important in the transverse as in the planar directions. Simple fountain flow
approximations have been implemented with the Hele-Shaw formulation. However, in
view of its importance in RIM, since it determines the path and the location of each
reactive fluid element during filling, a three-dimensional formulation is expected to
provide more detailed and accurate information [14].
1.4 Motivation and objectives
The core objective of this work is to assess the potential of the commercial general-purpose
CFD software, CFX-5 (recently renamed ANSYS CFX), a fully implicit, unstructured,
node centred, finite-volume based code, to simulate the complete three-dimensional
behaviour of the RIM simultaneous filling and curing processes.
Because of the inherent RIM process and parts characteristics, such as resins' low
viscosities, common complex geometries possibly with the presence of inserts, together
with the importance to accurately determine the flow behaviour at the flow front region, a
full 3D approach seems to be much more valuable, relatively to the 2½D approach, than in
simulations of the TIM process.

Introduction
9
Despite the existence of commercial injection moulding software packages able to
perform 3D mould filling analysis, the 2½D approach, as previously mentioned, is still the
standard analysis in injection moulding simulations and, extremely little has been made
regarding the simulation of this important kind of industrial process making use of general-
purpose CFD software packages.
Although three-dimensional simulations for injection moulding are still notorious
for their excessive computation time requirements, as computer performances increase
they are anticipated to increase in the near future [14]. It is reported in [81] that in 2003 the
same CFD problem could be run over 100 times faster than in 1996 (7 years before), a
factor of 30 due to increase in processors speed (double every 18 months, as stated by the
modern and most popular formulation of the Moore's law [82]), and the remaining due to
enhancements in numerical methods and faster solver algorithms. Additionally, the
possibility to run calculations on a network of PCs, something unavailable not too long
ago, and improvements in parallel computing algorithms, will tend towards faster time
solutions to many CFD problems [81] at relatively low hardware costs.
Commercial multi-purpose CFD software packages are usually cheaper than
commercial specialized injection moulding packages, and are able to perform numerical
simulations of many other varieties of processes, making them a much more versatile and
valuable tool than the specialized packages usual in injection moulding.

Introduction
10

11
2 Simulation of the filling process
2.1 The homogeneous model
2.1.1 Physical model
The process of an enclosed volume being filled with a fluid A and displacing the original
fluid B, usually air, is common in engineering practice. The filling of a mould is an
example of it. This process on its own, even not considering the energy and the chemical
reaction, is complex, as it is a transient two-phase free surface flow.
For this type of flows, where the phases are completely stratified and the interface
is well defined, it makes sense to assume that both phases share a single velocity field [83].
Thus, choosing the homogeneous model of CFX, the flow is described by [83]:
the equation of conservation of total mass:
U
0
t
(
)
wU x U
w
JJG
(15)
the momentum equations:
x
direction:
x
2
M
u
u
u v
u w
t
x
y
z
p
u
u
v
u
w
2
S
x
x
x
y
y
x
z
z
x
(
)
(
)
(
)
(
)
w U
w U
w U
w U
w
w
w
w
ª
º
ª
º
§
·
§
·
§
·
w
w
w
w
w
w
w
w
w
P
P
P
«
»
«
»
¨
¸
¨
¸
¨
¸
w
w
w
w
w
w
w
w
w
©
¹
©
¹
©
¹
¬
¼
¬
¼
(16)
y
direction:
y
2
M
v
u v
v
v w
t
x
y
z
p
v
u
v
v
w
2
S
y
x
x
y
y
y
z
z
y
(
)
(
)
(
)
(
)
w U
w U
w U
w U
w
w
w
w
ª
º
ª
º
§
·
§
·
§
·
w
w
w
w
w
w
w
w
w
P
P
P
«
»
«
»
¨
¸
¨
¸
¨
¸
w
w
w
w
w
w
w
w
w
©
¹
©
¹
©
¹
¬
¼
¬
¼
(17)
z
direction:
z
2
M
(
w)
(
u w)
(
v w)
(
w )
t
x
y
z
p
w
u
w
v
w
2
S
z
x
x
z
y
y
z
z
z
w U
w U
w U
w U
w
w
w
w
ª
º
ª
º
§
·
§
·
§
·
w
w
w
w
w
w
w
w
w
«
»
«
»
¨
¸
¨
¸
¨
¸
w
w
w
w
w
w
w
w
w
©
¹
©
¹
©
¹
¬
¼
¬
¼
(18)

Simulation of the filling process
12
or, in general vector form:
T
M
U
U
U
p
U
U
S
t
(
)
(
)
(
)
w U
ª
º
x U
x G P
«
»
¬
¼
w
JG
JG JG
JG
JG
JJJG
(19)
the conservation of mass of phase
D:
r
r
U
0
t
(
)
(
)
D
D
D
D
w
U
x
U
w
JG
(20)
and the constraint that the two phases completely fill up the available volume:
r
r
1
D
E
(21)
Therefore, the flow is characterized by 6 equations (5 of them partial differential
equations) and 6 unknowns (p, u, v, w, r
D
and r
E
).
In these equations, and are the volume fraction weighted mixture density and
viscosity evaluated as [83]:
r
r
D D
E E
U U U (22)
r
r
D D
E E
P P P (23)
S
M
represents the source of momentum due to internal forces, the gravity in this case, and
is given by [83]:
M
ref
S
g
U U
JJJG
G
(24)
where
U
ref
is the reference density.
2.1.2 Numerical model
CFX is based on the Finite Volume Method (FVM), and each node in the mesh is at the
centre of a finite control volume, Fig. 3.
The partial differential equations are integrated over all the control volumes, using
Gauss' Divergence Theorem to convert volume integrals on to surface integrals:
V
S
dV
n
dS
I
I
³
³
G
(25)
where n
G
represents the outward normal surface vector.

Simulation of the filling process
13
Fig. 3: Representation of the control volume associated with each mesh node.
As the control volume surface is composed by a series of surface segments, or
faces, the surface integrals can be transformed into a sum of integrals over the faces [84]:
f
S
f
n
dS
n
dS
§
·
¨
¸
I
I
¨
¸
©
¹
¦
³
³
G
G
(26)
These face integrals are then discretely represented at integration points, located at the
centre of each face:
ip
ip
ip
f
ip
f
n
dS
n
A
§
·
¨
¸
I
|
I
¨
¸
©
¹
¦
¦
³
G
JJJG
(27)
A
ip
being the area of the face associated with the integration point.
According to this, the discrete forms of the governing equations are [83]:
the equation of conservation of total mass:
ip
ip
V
dV
U n A
0
t
§
·
w ¨
¸
U
U x
¨
¸
w ©
¹
¦
³
JJG G
(28)
the momentum equations:
x
direction:
x
ip
ip
x
ip
ip
ip
V
x
y
z
M
ip
ip
u dV
m
u
n
A p
t
u
u
v
u
w
2
n
n
n
A
S
V
x
y
x
z
x
§
·
w ¨
¸
U
¨
¸
w ©
¹
½
ª
º
§
·
§
·
w
w
w
w
w
°
°
P
®
¾
«
»
¨
¸
¨
¸
w
w
w
w
w
©
¹
©
¹
°
°
¬
¼
¯
¿
¦
¦
³
¦
(29)

Simulation of the filling process
14
y
direction:
y
ip
ip
y
ip
ip
ip
V
x
y
z
M
ip
ip
v dV
m
v
n
A p
t
v
u
v
v
w
n
2
n
n
A
S
V
x
y
y
z
y
§
·
w ¨
¸
U
¨
¸
w ©
¹
½
ª
º
§
·
§
·
w
w
w
w
w
°
°
P
®
¾
«
»
¨
¸
¨
¸
w
w
w
w
w
©
¹
©
¹
°
°
¬
¼
¯
¿
¦
¦
³
¦
(30)
z
direction:
z
ip
ip
z
ip
ip
ip
V
x
y
z
M
ip
ip
w dV
m
w
n
A p
t
w
u
w
v
w
n
n
2
n
A
S
V
x
z
y
z
z
§
·
w ¨
¸
U
¨
¸
w ©
¹
½
ª
º
§
·
§
·
w
w
w
w
w
°
°
P
®
¾
«
»
¨
¸
¨
¸
w
w
w
w
w
©
¹
©
¹
°
°
¬
¼
¯
¿
¦
¦
³
¦
(31)
and the mass conservation equation of phase
D [85]:
ip
ip
V
r dV
r
U n A
0
t
D D
D
D
§
·
w ¨
¸
U
U x
¨
¸
w ©
¹
¦
³
JG G
(32)
ip
m being the discrete mass flow rate through a face of the finite volume, obtained from
the preceding iteration as [86]:
ip
ip
m
U n A
U x
JJG G
(33)
n
x
, n
y
and n
z
are the Cartesian components of the outward normal surface vector, and V the
volume of the control volume.
Choosing the second order backward Euler scheme, the transient term on the
equation of conservation of total mass (28) is approximated as [87]:
2
0
00
V
V
1
dV
2
1
t
t
1
§
·
½
w
°
°
ª
º
¨
¸
U
|
T T U T U U
®
¾
«
»
¬
¼
¨
¸
w
'
T T
°
°
¯
¿
©
¹
³
(34)
and, on the momentum equations (29), (30) and (31) as:
i
V
2
0
0
00
00
i
i
i
U dV
t
V
1
2
U
1
U
U
t
1
§
·
w ¨
¸
U
|
¨
¸
w ©
¹
½
°
°
ª
º
T T U
T U
U
®
¾
«
»
¬
¼
'
T T
°
°
¯
¿
³
(35)

Simulation of the filling process
15
where U
i
represents the i component of the velocity vector,
0
denotes the old time
level,
00
the time before the old time level, and
T is the quotient between the old and the
new time step:
0
0
00
0
t
t
t
t
t
t
'
T
'
(36)
At this point, the values of the variables and of the derivatives at the integration
points must be obtained from the values of the variables stored at the mesh nodes. The
value of the pressure p at ip is evaluated using finite element shape functions (linear in
terms of parametric coordinates) [83]:
ip
n
n
ip
n
p
N
p
ª
º
«
»
¬
¼
¦
(37)
where N
n
is the shape function for node n, and p
n
the value of p at the node n. For the
diffusion terms, the derivatives at the integration points are also evaluated making use of
shape functions. For example, the derivative of
I in the x direction at ip is obtained as [83]:
n
n
ip
ip
n
N
x
x
ª
º
§
·
§
·
w
wI
«
»
I
¨
¸
¨
¸
w
w
«
»
©
¹
©
¹
¬
¼
¦
(38)
On the advection terms, a variable
I at the integration points is obtained as:
ip
up
R
I I E I x
JG
(39)
where
I
up
is the value of
I at the upwind node,
I
is the upwind gradient of
I, R
JG
is the
position vector from the upwind node to ip, and is a blend factor. Different values of
produce different advection schemes. If = 0 the scheme is a first order Upwind
Differencing Scheme (UDS), if = 1, the scheme is a Second Order Upwind (SOU) biased
scheme. Here, the advection scheme chosen for the total mass and momentum equations is
based on that of [88], the high resolution scheme, where is variable and locally computed
to be as close as possible to 1 without violating boundedness principles.
However, for free surface applications, this scheme still causes too much numerical
diffusion to the volume fraction equation. Instead, in order to keep the interface sharp, and
as the order of accuracy is not the most important when the solutions are inherently
discontinuous, when the "free surface flow" option is selected, CFX uses a compressive
differencing scheme for the advection term of the volume fraction equation. This is made
allowing 1, but reducing it as much as necessary to still maintain boundedness [85]. A

Simulation of the filling process
16
compressive transient scheme is also used for the transient term in the volume fraction
equation.
2.1.3 Case study
A simple test case, the filling of a space between two parallel plates with a resin with
constant density and viscosity, as shown in Fig. 4, was simulated. The width of the plates is
considered to be infinite, so a two-dimensional approach is performed.
Fig. 4: Geometry of the simulated test case.
As CFX does not have a 2D solver, a three-dimensional simulation, with an
one-element-thick mesh and symmetry boundary conditions imposed on the "front" and
"back" faces, has to be performed. This applies whenever a two-dimensional simulation is
mentioned in this work.
The properties of the filling fluid were chosen to be similar to those of the common
RIM resins: density 1000 kg/m
3
and viscosity 0.05 kg/(m
s). The space between the plates
is initially filled with air, whose density and viscosity, taken as constant, are 1.185 kg/m
3
and 1.831
u10
-5
kg/m
s, respectively.
The imposed boundary conditions are: a parabolic velocity profile with an average
velocity of 0.1037 m/s at the inlet, pressure equal to zero at the outlet, and the condition of
no slip on the walls. The values for the inlet average velocity and for the thickness H were
chosen to be the same used for the mould modelled in Section 4.2.
The Reynolds number may be obtained as:
h
U D
Re
U
P
(40)
where D
h
is the hydraulic diameter defined as [89]:

Simulation of the filling process
17
h
4 A
D
P
(41)
which, for a rectangular cross section with infinite width W, takes the value:
W
h
h
4 W H
D
D
2 H
2 W
H
o
f
,
(
)
(42)
Thus, the Reynolds number is approximately 13 for the resin and 43 for the air. The
flow is therefore unquestionably laminar, and the parabolic velocity profile imposed at the
inlet boundary makes sense as a fully developed laminar flow. The velocity profile of the
air at the outlet boundary is also expected to be parabolic.
Considering the control volume surrounded by the dashed line represented in
Fig. 5:
Fig. 5: Representation of the control volume and forces.
and applying the Newton's Second Law of motion in the x direction and an integral
balance, one obtains [89]:
x
cv
cs
d
d
F
m u
u
dV
u
U n dA
dt
dt
§
·
¨
¸
U
U
x
¨
¸
©
¹
¦
³³³
³³
JG G
(43)
where cv designates the control volume and cs its surface.
Being the velocity on the walls zero, and assuming that the velocity profile at the
outlet is parabolic due to the flow being laminar, as at the inlet, the second term on the
right side of equation (43) reduces to:
2
2
cs
outlet
inlet
2
resin
air
u
U n dA
u dA
u dA
6
U
H W
5
U
x
U
U
U
U
³³
³³
³³
JG G
(44)
where W indicates the width of the control volume.

Simulation of the filling process
18
Assuming a flat flow front, the density is constant on each cross section. Thus, the
first term on the right side of equation (43) may be rewritten as:
L
cv
0
S
d
d
u
dV
u dA dx
dt
dt
ª
º
§
·
§
·
«
»
¨
¸
¨
¸
U
U
¨
¸
¨
¸
«
»
©
¹
©
¹
¬
¼
³³³
³ ³³
(45)
where S indicates a cross section. But, because of conservation of mass, and as the fluids
are considered incompressible, the volume flow rate is constant:
in
in
S
u dA
U
A
U H W
constant
(
)
³³
(46)
Then, equation (45) may be simplified as:
^
`
resin
air
cv
d
d
u
dV
H W
U
Lf
L
Lf
dt
dt
§
·
¨
¸
ª
º
U
U
U
¬
¼
¨
¸
©
¹
³³³
(47)
As the plates are parallel and therefore the cross section area is constant, the
position of the flow front, Lf, may be obtained by integrating the inlet velocity over the
time:
t
0
Lf
U dt
³
(48)
leading finally to:
`
resin
air
air
cv
2
resin
air
d
dU
u
dV
H W
Lf
L
dt
dt
U
§
·
°
¨
¸
ª
º
U
U
U
U
®
¬
¼
¨
¸
°¯
©
¹
U
U
³³³
(49)
If the inlet velocity is constant, the derivative on the right side of the above equation is
zero, and the sum of forces is therefore:
N
2
x
in
in
out
out
resin
air
0
1
F
P
A
P
A
Fg
Fw
U
H W
5
U
U
¦
(50)
The gravity force, Fg, is achieved by integrating (
Ug) over the control volume,
which considering again that the density is constant on each cross section, and that the
cross section area is constant and equal to H
W, leads to:

Simulation of the filling process
19
L
resin
air
air
cv
0
Fg
g dV
g H W
dx
g H W
Lf
L
ª
º
U
U
U
U
U
¬
¼
³³³
³
(51)
The wall shear stress is obtained, for Newtonian fluids, as:
w
wall
du
dy
@
ª
º
W P
«
»
¬
¼
(52)
which, for a laminar flow with parabolic velocity profile may be reworked as:
w
6
U
H
P
W
(53)
and the walls shear force, Fw, is obtained by integrating
W
w
over the walls:
L
L
w
2
walls
0
0
resin
air
air
2
6
U
12 U H W
Fw
dA
2
W dx
dx
H
H
12 U H W
Lf
L
H
P
W
u
P
ª
º
P
P
P
¬
¼
³³
³
³
(54)
Note that, due to the fountain flow, the velocity profile at the flow front region is
not parabolic, and therefore there is an error on the above expression for Fw. However, as
the flow front region is small compared with the total length L, this error may be neglected.
Finally, inserting the expressions of equations (51) and (54) into equation (50), the
inlet pressure, P
in
, is obtained as:
in
resin
air
air
2
resin
air
air
resin
air
2
P
g
Lf
L
12 U
1
Lf
L
U
5
H
ª
º
U
U
U
¬
¼
ª
º
P
P
P
U
U
¬
¼
(55)
On this last equation, the density and viscosity of the air could be fairly neglected, as they
are considerably smaller than those of the resin:
U
air
/
U
resin
# 0.12% and P
air
/
P
resin
# 0.04%.
2
in
resin
resin
resin
2
12 U
1
P
g Lf
Lf
U
5
H
#
U
P
U
(56)
The errors caused by this simplification would be, for Lf = 0.1L, 1 % on the first term of
the right side of the equation, 0.3 % on the second term and 0.12 % on the last term. For
Lf = 0.5L, the errors would be 0.12 %, 0.04 % and 0.12 %, respectively.
On the CFX model, the buoyancy reference density,
U
ref
, was set to the air density.
This is interpreted as a momentum source on equation (19) equal to:

Simulation of the filling process
20
M
ref
air
S
g
g
(
)
(
)
U U
U U
JJJG G
G
(57)
which means that the air hydrostatic pressure effect is omitted. Consequently, to obtain the
inlet pressure according to the CFX homogeneous model, P
in
hom
, the term L
U
air
shall be
withdrawn from the first term of the right side of equation (55):
hom
in
resin
air
2
resin
air
air
resin
air
2
P
g Lf
12 U
1
Lf
L
U
5
H
U
U
ª
º
P
P
P
U
U
¬
¼
(58)
This causes an error of 1.2 % on the first term of the right side of this last equation, for
Lf = 0.1L, and 0.24 % for Lf = 0.5L. The last term on equations (55) and (58) corresponds
only to 2.1 Pa.
To calculate the value of the pressure at each position x, the velocity profile is
considered to be parabolic at each cross section, which is valid for the whole domain
except at the vicinity of the flow front. The same previous analysis is performed, except
that this time the control volume is considered to comprise only the region between a cross
section at the position x and the outlet. Thus the pressure depends on x as:
resin
air
resin
air
2
2
resin
air
air
air
2
g
Lf
x
L
Lf
12 U
Lf
x
L
Lf
for x
Lf
H
1
P x
U
5
12 U
g
L
x
L
x
for x
Lf
H
(
)
( )
(
)
ª
º
U
U
¬
¼
°
°
ª
º
P
P
°
¬
¼
°
°
U
U
®
°
°
°
° U
P
!
°
¯
(59)
and the pressure according to the CFX homogeneous model as:
resin
air
resin
air
2
2
hom
resin
air
air
2
g
Lf
x
12 U
Lf
x
L
Lf
for x
Lf
H
1
P
x
U
5
12 U
L
x
for x
Lf
H
(
)
( )
(
)
U
U
°
°
ª
º
P
P
¬
¼
°
°
°
U
U
®
°
°
°
°
P
!
°
¯
(60)

Simulation of the filling process
21
On these two last equations there is a discontinuity at x = Lf. This is due to the
assumptions of a flat flow front and of a parabolic velocity profile at every cross section,
which is not true at the flow front region. However, this discontinuity represents only
2.1 Pa.
The derivative of the pressure along x is then:
resin
resin
2
air
air
2
12 U
g
for x
Lf
H
dP x
dx
12 U
g
for x
Lf
H
(
)
( )
(
)
§
·
P
U
° ¨
¸
¨
¸
° ©
¹
°
®
° §
·
P
° U
!
¨
¸
¨
¸
° ©
¹
¯
(61)
and:
resin
resin
air
2
hom
air
2
12 U
g
for x
Lf
H
dP
x
dx
12 U
for x
Lf
H
(
)
( )
(
)
ª
º
P
U
U
° «
»
«
»
° ¬
¼
°
®
°
P
°
!
°
¯
(62)
Due to symmetry, only half of the geometry was modelled. The simulations were
performed with four different meshes: mesh1 with 24 elements on the longitudinal
direction by 5 on the transversal direction, mesh2 with 48 by 5 elements, mesh3 with
24 by 20 elements, and mesh4 with 48 by 20 elements.
The position of the resin-air interface after 0.35 s is represented in Fig. 6a-6d. The
three contour lines denote the resin volume fractions of 0.9, 0.5 and 0.1. The vertical line
indicates the theoretical position of the flow front, obtained with equation (48), which for
t = 0.35 s results in Lf = 36.3 mm.
It is noticeable that, for the four meshes, the interface is captured within two mesh
elements, proving the efficiency of the compressive discretization schemes. Increasing the
number of elements on the longitudinal direction improves its definition in that direction,
but does not change its position. By increasing the number of elements on the transverse
direction, the interface position gets closer to its theoretical position, though it is always
ahead.
However, the most important conclusion from these simulations is that the interface
does not touch the walls, there existing a layer of air between the resin and the walls. This
is clearly not in accordance with reality. The volume of resin which is ahead of the
theoretical position of the flow front (the vertical line) shall correspond to the volume of
resin lacking near the walls.

Simulation of the filling process
22
(a)
(b)
(c)
(d)
Fig. 6: Position of the resin-air interface for t = 0.35 s. Contour lines denote the resin
volume fractions 0.9, 0.5 and 0.1. (a) mesh1: 24
u5 elements; (b) mesh2: 48u5 elements;
(c) mesh3: 24
u20 elements; (d) mesh4: 48u20 elements.

Simulation of the filling process
23
In Fig. 7a-7d, the value of the resin volume fraction along the transverse direction,
at midway between the inlet and the theoretical position of the flow front, is represented.
0.0
0.4
0.8
1.2
1.6
0.0
0.5
1.0
Vol. Fraction
y [
m
m
]
(a)
0.0
0.4
0.8
1.2
1.6
0.0
0.5
1.0
Vol. Fraction
y [
m
m
]
(b)
0.0
0.4
0.8
1.2
1.6
0.0
0.5
1.0
Vol. Fraction
y
[
mm]
(c)
0.0
0.4
0.8
1.2
1.6
0.0
0.5
1.0
Vol. Fraction
y
[
mm]
(d)
Fig. 7: Resin volume fraction along the transverse direction, at x = ½ Lf, for t = 0.35 s.
(a) mesh1, average volume fraction = 0.925; (b) mesh2, average volume fraction = 0.928;
(c) mesh3, average volume fraction = 0.960; (d) mesh4, average volume fraction = 0.962.
From this last figure, one concludes that increasing the number of mesh elements on the
longitudinal direction does not change this behaviour, while increasing the number of
elements on the transverse direction causes the average volume fraction to increase but the
volume fraction on the wall nodes to decrease.
The viscosity is obtained as the volume fraction weighted phases' viscosities by
equation (23). Because the resin volume fraction on the walls is close to zero, the
computed viscosity on the walls nodes will be much smaller than its actual value, resulting
in a higher velocity gradient on the transverse direction next to the walls and consequently
a wrong velocity profile. When solving the cure and energy equations, this will cause an
error on the advection terms of those equations, which will possibly lead to completely
wrong final results.
Also due to the lower computed values of the viscosity on the walls nodes, the wall
shear stress, which is obtained by equation (52), will be lower than its real value, leading to
a lower value of the viscous effect contribution to pressure. This is of major importance, as

Simulation of the filling process
24
injection pressure is one of the most important parameters to be predicted by numerical
simulation.
Comparing, for t = 0.35 s, the inlet pressure obtained from CFX with mesh1 (the
coarsest mesh), 536 Pa, and the theoretical pressure obtained with equations (58) or (60),
574 Pa, the pressure obtained with CFX is just 6.5% lower. However, the error caused by
the lower shear stress on the walls is masked by an increase of the gravity effect in
pressure due to the fact that the flow front is ahead of its theoretical position.
The pressure due to the gravity effect may be evaluated from the hydrostatic law
as [89]:
dP grav
g
dx
(
) U (63)
or:
hom
air
dP
grav
g
dx
(
) U U (64)
which, considering the pressure at the outlet as zero, leads to:
L
x
P grav
g dx
(
)
U
³
(65)
or:
L
hom
air
x
P
grav
g dx
(
)
U U
³
(66)
The pressure due to the viscous effect may be regarded as the difference between
the whole pressure and the pressure due to the gravity effect: P-P(grav), or P
hom
-P
hom
(grav).
This is not completely true, as there is also a small contribution from the rate of change of
linear momentum, 2.1 Pa, and thereby it is not designated here by P(visc).
Thus, the comparison will be done between P
hom
-P
hom
(grav), which is obtained as:
resin
air
2
2
hom
resin
air
hom
air
2
12 U
Lf
x
L
Lf
H
1
P
x
U
5
P
grav x
12 U
L
x
H
ª
º
P
P
°
¬
¼
°
°
§
·
U
U
°
¨
¸ ®
¨
¸ °
©
¹
°
°
P
°
¯
( )
(
)( )
for x
Lf
for x
Lf
!
(
)
(
)
(67)

Simulation of the filling process
25
and the pressure obtained from CFX, P
CFX
, minus its gravity contribution:
L
CFX
CFX
CFX
air
x
CFX
P
P
grav
P
g dx
(
)
§
·
¨
¸
U U
¨
¸
©
¹
³
(68)
In this equation, the integration is done along the centre line where the flow front has its
most advanced position.
The comparison between the results obtained from the simulations and the
theoretical ones are represented in Fig. 8a-8d. The four graphics, obtained with the four
different meshes, are quite similar.
The derivative of (P
CFX
-P
CFX
(grav)) is compared with theoretical values in
Fig. 9a-9d. It is possible to observe the error induced in (P
CFX
-P
CFX
(grav)). The derivative
should be constant for each phase, from equation (67) it should take the value -6076 Pa/m
for the resin and -2.2 Pa/m for the air. Instead, a nearly ramp function is observed.
In Table 1, the results obtained at the inlet are presented and compared with the
theoretical ones. It is interesting to notice that when increasing, the number of mesh
elements on the transverse direction by a factor of 4, from mesh1 to mesh3 and from
mesh2 to mesh4, the error of the obtained pressure, P
CFX
, which is always negative,
increases, in absolute value, from 6.7 % to 9.0 % and from 6.3 % to 8.5 %, respectively. At
a first glance this might look somehow unexpected, but the reason is that the error of the
gravity contribution to pressure, P
CFX
(grav), which is always positive, decreased from
10.9 % to 5.4 % and from 10.3 % to 5.0 %, respectively. This was already observed in
Fig. 6a-6d, the interface positions obtained with mesh3 and mesh4, though always ahead,
are closer to the theoretical position than the obtained with mesh1 and mesh2.
Table 1: Comparison between the theoretical pressures and the pressures obtained from the
simulations, at x = 0 (inlet), at t = 0.35 s. The errors relative to the theoretical values are
presented inside brackets.
P
hom
P
hom
(grav) P
hom
-P
hom
(grav)
Theory
573.7 355.3 218.4
P
CFX
P
CFX
(grav) P
CFX
-P
CFX
(grav)
535.2 394.0 141.1
mesh1:
24×5 elem.
(-6.7 %)
(+10.9 %)
(-35.4 %)
537.3 391.7 145.7
mesh2:
48×5 elem.
(-6.3 %)
(+10.3 %)
(-33.3 %)
522.1 374.3 147.8
mesh3:
24×20 elem.
(-9.0 %)
(+5.4 %)
(-32.4 %)
524.9 372.9 152.0
mesh4:
48×20 elem.
(-8.5 %)
(+5.0 %)
(-30.4 %)

Simulation of the filling process
26
(a)
0
100
200
300
400
500
600
0
10
20
30
40
50
x [mm]
P [Pa]
Phom
PCFX
PCFX(grav)
Phom(grav)
(Phom-Phom(grav))
(PCFX-PCFX(grav))
(b)
0
100
200
300
400
500
600
0
10
20
30
40
50
x [mm]
P [Pa]
Phom
PCFX
PCFX(grav)
Phom(grav)
(Phom-Phom(grav))
(PCFX-PCFX(grav))
(c)
0
100
200
300
400
500
600
0
10
20
30
40
50
x [mm]
P [Pa]
Phom
PCFX
PCFX(grav)
Phom(grav)
(Phom-Phom(grav))
(PCFX-PCFX(grav))
(d)
0
100
200
300
400
500
600
0
10
20
30
40
50
x [mm]
P [Pa]
Phom
PCFX
PCFX(grav)
Phom(grav)
(Phom-Phom(grav))
(PCFX-PCFX(grav))
Fig. 8: Comparison between the pressures obtained from simulations and theoretical
values, for t = 0.35 s. (a) mesh1; (b) mesh2; (c) mesh3; (d) mesh 4.

Simulation of the filling process
27
(a)
-8000
-6000
-4000
-2000
0
2000
0
10
20
30
40
50
x [mm]
dP/dx [Pa/m]
d(PC FX-PC FX(grav))/dx
d(Phom-Phom(grav))/dx
(b)
-8000
-6000
-4000
-2000
0
2000
0
10
20
30
40
50
x [mm]
dP/dx [Pa/m]
d(PC FX-PC FX(grav))/dx
d(Phom-Phom(grav))/dx
(c)
-8000
-6000
-4000
-2000
0
2000
0
10
20
30
40
50
x [mm]
dP/dx [Pa/m]
d(PCFX-CFX(grav))/dx
d(Phom-Phom(grav))/dx
(d)
-8000
-6000
-4000
-2000
0
2000
0
10
20
30
40
50
x [mm]
dP/dx [Pa/m]
d(PC FX-PC FX(grav))/dx
d(Phom-Phom(grav))/dx
Fig. 9: Comparison between the derivative of (P
CFX
-P
CFX
(grav)) and theoretical
values, for t = 0.35 s. (a) mesh1; (b) mesh2; (c) mesh3; (d) mesh 4.

Simulation of the filling process
28
The error of (P
CFX
-P
CFX
(grav)), in absolute value, shows only a tiny decrease when
increasing the number of mesh elements. Mesh4 has 8 times the number of mesh elements
of mesh1, and the error only decreased from 35.4 % to 30.4 %. Thus, it seems that
increasing the number of mesh elements, at least within reasonable limits, is not the cure
for the problem.
The computing time, for each mesh, to complete 0.47 s of simulation is shown in
Fig. 10. CFX was running in double precision, in a workstation PC with a Pentium
®
4
2.5 MHz processor and 512 MB RAM memory, using Windows
®
2000 operating system.
mesh1
mesh2
mesh3
mesh4
0
2
4
6
8
0
200
400
600
800
1000
nr. of mesh elements
CPU time [h]
Fig. 10: Computing time for the simulations to complete 0.47 s.
This unwanted behaviour has although been mentioned by some authors, who claim
that that in filling simulations the no-slip condition on mould walls should be imposed only
on the filled portion of the mould, as a no-slip boundary condition in the air will prevent
the flow front from touching the walls [14, 56].
2.2 The inhomogeneous model
2.2.1 Physical model
Despite the physical model described in the previous sub-section being theoretically
correct, it is unable to produce satisfactory results with the also physically correct no-slip
boundary condition on the walls.
CFX does not allow the implementation of conditional boundary conditions: it is
not possible to define a wall boundary condition as no-slip if the volume fraction of the
liquid phase is above a certain value, 0.5 for instance, and as free-slip if it is below.
Therefore, the only way to prescribe a no-slip boundary condition on the filled portion of

Simulation of the filling process
29
the mould and a free-slip boundary condition in the empty portion is by using the CFX's
inhomogeneous multiphase flow model. In the inhomogeneous model, each phase has its
own flow field and they interact via interphase transfer terms [83]. According to the
assumption made by this model, the flow is described by:
the equations of conservation of mass of each phase:
r
r
U
0
t
(
)
(
)
D
D
D
D
D
w
U
x
U
w
JJJG
(69)
r
r
U
0
t
(
)
(
)
E E
E
E
E
w U
x
U
w
JJJG
(70)
the momentum equations for phase
D:
x
direction:
x
x
2
M
r
u
r
u
r
u
v
r
u
w
t
x
y
z
u
u
v
p
r
2 r
r
x
x
x
y
y
x
u
w
r
S
M
z
z
x
(
)
(
)
(
)
(
)
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
DE
w
U
w
U
w
U
w
U
w
w
w
w
ª
º
w
w
w
§
·
§
·
w
w
w
P
P
«
»
¨
¸
¨
¸
w
w
w
w
w
w
©
¹
©
¹
¬
¼
ª
º
w
w
§
·
w
P
«
»
¨
¸
w
w
w
©
¹
¬
¼
(71)
y
direction:
y
y
2
M
r
v
r
u
v
r
v
r
v
w
t
x
y
z
v
u
v
p
r
r
2 r
y
x
x
y
y
y
v
w
r
S
M
z
z
y
(
)
(
)
(
)
(
)
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
DE
w
U
w
U
w
U
w
U
w
w
w
w
ª
º
w
w
w
§
·
§
·
w
w
w
P
P
«
»
¨
¸
¨
¸
w
w
w
w
w
w
©
¹
©
¹
¬
¼
ª
º
w
w
§
·
w
P
«
»
¨
¸
w
w
w
©
¹
¬
¼
(72)
z
direction:
z
z
2
M
r
w
r
u
w
r
v
w
r
w
t
x
y
z
w
u
w
v
p
r
r
r
z
x
x
z
y
y
z
w
2 r
S
M
x
z
(
)
(
)
(
)
(
)
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
DE
w
U
w
U
w
U
w
U
w
w
w
w
ª
º
ª
º
w
w
w
w
§
·
§
·
w
w
w
P
P
«
»
«
»
¨
¸
¨
¸
w
w
w
w
w
w
w
©
¹
©
¹
¬
¼
¬
¼
w
§
·
w
P
¨
¸
w
w
©
¹
(73)

Simulation of the filling process
30
or, in general vector form:
T
M
r
U
r
U
U
t
r
p
r
U
U
S
M
(
)
(
)
(
)
D
D
D
D
D
D
D
D
D
D
D
D
D
DE
w
U
x
U
w
ª
º
x
P
«
»
¬
¼
JJJG
JJJG JJJG
JJJG
JJJG
JJJJG JJJJJG
(74)
the momentum equations for phase
E, for which just the vector form is presented:
T
M
r
U
r
U
U
t
r
p
r
U
U
S
M
(
)
(
)
(
)
E
E
E
E E
E
E
E
E
E
E
E
E
ED
w
U
x
U
w
ª
º
x
P
«
»
¬
¼
JJJG
JJJG JJJG
JJJG
JJJG
JJJJG JJJJJG
(75)
and the constraint that the volume fractions sum to unity:
r
r
1
D
E
(76)
By assuming a different velocity field for each phase, the flow is characterized by
9 equations (8 of them partial differential equations) and 9 unknowns (p, u
D
, v
D
, w
D
, u
E
, v
E
,
w
E
, r
D
and r
E
).
As in the homogeneous model, S
M
represents the source of momentum due to the
gravity force, and is given, for phase
D and E, respectively, by:
M
ref
S
r
g
D
D
D
U U
JJJJG
G
(77)
M
ref
S
r
g
E
E
E
U U
JJJJG
G
(78)
M
DE
and M
ED
are interphase momentum transfer terms, and they represent, respectively, the
force per unit volume exerted by phase
E on phase D, and by phase D on phase E. These
interfacial forces are equal in absolute value and opposite in direction [83]:
M
M
DE
ED
JJJJJG
JJJJJG
(79)
As both phases, resin and air, are continuous, there is not any implemented model in CFX
for M
DE
. It is simply defined as [83, 90]:
D
M
C
A
U
U
U
U
DE
DE
DE
E
D
E
D
U
JJJJJG
JJJG JJJG JJJG JJJG
(80)
where C
D
is a non-dimensional drag coefficient,
U
DE
is the mixture density given by:
r
r
DE
D D
E E
U
U U (81)
and A
DE
is the interfacial area per unit volume, given by:

Simulation of the filling process
31
r
r
A
d
D E
DE
DE
(82)
where d
DE
is a user-specified interface length scale.
CFX requires the definition of C
D
and d
DE
. However, substituting the expressions of
equations (81) and (82) into equation (80), leads to:
D
C
M
r
r
r
r
U
U
U
U
d
DE
D E
D D
E E
E
D
E
D
DE
U U
JJJJJG
JJJG JJJG JJJG JJJG
(83)
from which one concludes that what matters is the quotient C
D
/d
DE
and not their individual
values. This was checked by setting in CFX different values of C
D
and d
DE
but keeping
their quotient constant, and exactly the same results were obtained.
It is interesting to notice that as the volume fraction of one of the phases
approximates zero, the term M
DE
vanishes. Therefore, as one would expect, M
DE
is only
meaningful at the interface.
2.2.2 Numerical model
The discretization of the governing equations is done in the same way as for the
homogeneous model. The discrete forms of the equations are therefore [91]:
the equations of conservation of mass of each phase:
ip
ip
V
r
dV
r
U
n A
0
t
D
D
D
D
D
§
·
w ¨
¸
U
U
x
¨
¸
w ©
¹
¦
³
JJJJG G
(84)
ip
ip
V
r
dV
r
U
n A
0
t
E
E
E
E E
§
·
w ¨
¸
U
U
x
¨
¸
w ©
¹
¦
³
JJJG G
(85)
the momentum equations for phase
D:
in x direction:
ip
ip
ip
x
x
x
ip
ip
ip
V
x
y
z
ip
ip
M
r
u dV
m
u
r
n
A p
t
u
u
v
u
w
r
2
n
n
n
A
x
y
x
z
x
S
V
M
V
D
D
D
D
D
D
D
D
D
D
D
D
D
D
DE
§
·
w ¨
¸
U
¨
¸
w ©
¹
½
ª
º
w
w
w
w
w
§
·
§
·
°
°
P
®
¾
«
»
¨
¸
¨
¸
w
w
w
w
w
©
¹
©
¹
°
°
¬
¼
¯
¿
¦
¦
³
¦
(86)

Simulation of the filling process
32
in y direction:
ip
ip
ip
y
y
y
ip
ip
ip
V
x
y
z
ip
ip
M
r
v dV
m
v
r
n
A p
t
v
u
v
v
w
r
n
2
n
n
A
x
y
x
z
y
S
V
M
V
D
D
D
D
D
D
D
D
D
D
D
D
D
D
DE
§
·
w ¨
¸
U
¨
¸
w ©
¹
½
ª
º
w
w
w
w
w
§
·
§
·
°
°
P
®
¾
«
»
¨
¸
¨
¸
w
w
w
w
w
©
¹
©
¹
°
°
¬
¼
¯
¿
¦
¦
³
¦
(87)
in z direction:
ip
ip
ip
z
z
z
ip
ip
ip
V
x
y
z
ip
ip
M
r
w dV
m
w
r
n
A p
t
w
u
w
v
w
r
n
n
2
n
A
x
z
y
z
z
S
V
M
V
D
D
D
D
D
D
D
D
D
D
D
D
D
D
DE
§
·
w ¨
¸
U
¨
¸
w ©
¹
½
ª
º
w
w
w
w
w
§
·
§
·
°
°
P
®
¾
«
»
¨
¸
¨
¸
w
w
w
w
w
©
¹
©
¹
°
°
¬
¼
¯
¿
¦
¦
³
¦
(88)
where
ip
m
D
is the discrete mass flow of phase
D through a face of the finite volume,
obtained from the previous iteration as:
ip
ip
m
r
U
n A
D
D
D
D
U
x
JJJJG G
(89)
The three momentum equations for phase
E are analogous to the three momentum
equations for phase
D presented above.
As the second order backward Euler scheme is used, the transient term on the
momentum equations (86) to (88) is approximated as:
`
i
i
i
i
V
2
0
0
0
00
00
00
V
1
r
U
dV
t
t
1
2
r
U
1
r
U
r
U
D
D
D
D
D
D
D
D
D
D
D
D
§
·
w
°
¨
¸
U
|
®
¨
¸
w
'
T T
°¯
©
¹
ª
º
T T U
T
U
U
«
»
¬
¼
³
(90)
where
T is the quotient between the old and the new time step, as already defined in
equation (36). The values of the variables and of the derivatives at the integration points
are obtained according to what was previously stated for the homogeneous model.

Simulation of the filling process
33
2.2.3 Case study
The case under study is the same one used for the homogeneous model: the filling of a
space between two parallel plates with a resin, as represented in Fig. 4. However, now the
no-slip condition on the walls is only applied to the equations referring to the resin phase,
while the free-slip condition is applied to the equations referring to the air phase.
Consequently, as the wall shear stress is zero for the air phase and, due to this, the
velocity profile of the air at the outlet is expected to be flat, the theoretical pressure defined
with these boundary conditions, P
inhom
, is slightly different than the one defined with the
boundary conditions used in the homogeneous model, P
hom
.
For these boundary conditions:
2
resin
air
cs
6
u
U n dA
U
H W
5
§
·
U
x
U
U
¨
¸
©
¹
³³
JG G
(91)
and:
Lf
resin
w
2
walls
0
12 U H W
Lf
6
U
Fw
dA
2
W dx
H
H
P
P
W
u
³³
³
(92)
Therefore, P
inhom
(x) is given by:
resin
air
2
resin
inhom
resin
2
g
Lf
x
12 U
1
Lf
x
U
P
x
5
H
0
U
U
°
P
°°
U
®
°
°
°¯
( )
for x
Lf
for x
Lf
!
(
)
(
)
(93)
As in the analysis done for the homogeneous model, the gravity effect in pressure
may be obtained from:
L
inhom
air
x
P
grav x
g dx
U U
³
(
)( )
(94)
and the pressure due to the viscous effect as:
2
resin
resin
2
inhom
inhom
12 U
1
Lf
x
U
for x
Lf
5
H
P
x
P
grav x
0
for x
Lf
P
U
°
°
§
·
°
¨
¸ ®
¨
¸ °
©
¹
°
°
!
¯
(
)
( )
(
)( )
(
)
(95)

Simulation of the filling process
34
As mentioned before, CFX requires the definition of the parameters C
D
and d
DE
to
evaluate M
DE
. But, because of the fact already stated, that what matters is their quotient
(C
D
/d
DE
) and not their individual values, d
DE
was kept constant with the default CFX value
(1
u10
-3
m), and C
D
was varied between 0.05 and 500.
The simulations were performed with mesh2 (48 by 5 elements) used for the
homogeneous model. The time step is 2
u10
-4
s, the residuals convergence target which
shall be achieved by all equations, except the volume fraction equations, is 10
-5
, up to a
maximum of 50 iterations within each time step. CFX automatically sets the convergence
criterion for the volume fraction equations, as these are usually more difficult to converge,
to 10 times the specified convergence target [83, 92, 93].
The positions of the resin-air interface, for t = 0.35 s, obtained with C
D
set to 0.05,
0.5, 5 and 50, are represented in Fig. 11a-11d.
The results obtained with C
D
set to 0.05, 0.5 and 5 are very similar. Unlike in the
homogeneous model simulations, now the interface touches the walls and, clearly upstream
of the interface there is resin and downstream of it there is air. Considering the interface
position the location where the volume fraction is 0.5, it is in very good agreement with its
theoretical position.
Although, when C
D
is increased to 50, Fig. 11d, the behaviour becomes similar to
that of the homogeneous model. The reason for this is that when M
DE
increases, the phases'
flow fields tend to equalize each other. Actually, the homogeneous model may be regarded
as a special case of the inhomogeneous model for a very large value of M
DE
.
For these simulations, the comparison between the pressures obtained from
simulation and the theoretical ones is shown in Fig. 12a-12d. As expected, because of the
results shown in Fig. 11a-11d, the graphics of Fig. 12a-12c are similar and the results
obtained from simulations nearly match the theoretical ones. As also expected, for
C
D
= 50, Fig. 12d, the pressure values obtained from simulation become similar to those
obtained with the homogeneous model.
The comparison between the derivative along x of (P
inhom
(x)-P
inhom
(grav)(x)) and
(P
CFX
-P
CFX
(grav)) is shown in Fig. 13a-13d. On the first three graphics, the results from
simulations match the theoretical values for the majority of the domain, except for the flow
front region. However for the two smallest values of C
D
, the error on this region is greater
than for C
D
= 5. For C
D
= 50, the graphic is again similar to the homogeneous model
results.
In Table 2, the results and the errors relative to theoretical values, at the inlet, for
each value of C
D
are presented. Comparing these results with the ones presented in
Table 1, obtained with the homogeneous model, one concludes that, for C
D
equal to 0.05,
0.5 and 5, the errors were considerably reduced. For the same mesh, mesh2, the whole

Simulation of the filling process
35
(a)
(b)
(c)
(d)
Fig. 11: Position of the resin-air interface for t = 0.35 s, obtained using mesh2. Contour lines
denote the resin volume fractions 0.9, 0.5 and 0.1. The vertical line represents the theoretical
position of the flow front. (a) C
D
= 0.05; (b) C
D
= 0.5; (c) C
D
= 5; (d) C
D
= 50.

Simulation of the filling process
36
(a)
0
100
200
300
400
500
600
0
10
20
30
40
50
x [mm]
P [Pa]
Pinhom
PCFX
PCFX(grav)
Pinhom(grav)
(Pinhom-Pinhom(grav))
(PCFX-PCFX(grav))
(b)
0
100
200
300
400
500
600
0
10
20
30
40
50
x [mm]
P [Pa]
Pinhom
PCFX
PCFX(grav)
Pinhom(grav)
(Pinhom-Pinhom(grav))
(PCFX-PCFX(grav))
(c)
0
100
200
300
400
500
600
0
10
20
30
40
50
x [mm]
P [Pa]
Pinhom
PCFX
PCFX(grav)
Pinhom(grav)
(Pinhom-Pinhom(grav))
(PCFX-PCFX(grav))
(d)
0
100
200
300
400
500
600
0
10
20
30
40
50
x [mm]
P [Pa]
Pinhom
PCFX
PCFX(grav)
Pinhom(grav)
(Pinhom-Pinhom(grav))
(PCFX-PCFX(grav))
Fig. 12: Comparison between the pressures obtained from simulations, with mesh2, and
theoretical values, for t = 0.35 s. (a) C
D
= 0.05; (b) C
D
= 0.5; (c) C
D
= 5; (d) C
D
= 50.

Simulation of the filling process
37
(a)
-10000
-8000
-6000
-4000
-2000
0
2000
4000
6000
8000
0
10
20
30
40
50
x [mm]
dP/dx [Pa/m]
d(PC FX-PC FX(grav))/dx
d(Pinhom-Pinhom(grav))/dx
(b)
-10000
-8000
-6000
-4000
-2000
0
2000
4000
6000
8000
0
10
20
30
40
50
x [mm]
dP/dx [Pa/m]
d(PC FX-PC FX(grav))/dx
d(Pinhom-Pinhom(grav))/dx
(c)
-10000
-8000
-6000
-4000
-2000
0
2000
4000
6000
8000
0
10
20
30
40
50
x [mm]
dP/dx [Pa/m]
d(PC FX-PC FX(grav))/dx
d(Pinhom-Pinhom(grav))/dx
(d)
-10000
-8000
-6000
-4000
-2000
0
2000
4000
6000
8000
0
10
20
30
40
50
x [mm]
dP/dx [Pa/m]
d(PC FX-PC FX(grav))/dx
d(Pinhom-Pinhom(grav))/dx
Fig. 13: Comparison between the derivative of (P
CFX
-P
CFX
(grav)), obtained with mesh2, and
theoretical values, for t = 0.35 s. (a) C
D
= 0.05; (b) C
D
= 0.5; (c) C
D
= 5; (d) C
D
= 50.

Simulation of the filling process
38
pressure error was cut down from 6.3 % to around 0.5 %, the error of the gravity
contribution from 10.3 % to around 1.5 % and the error of the viscosity contribution from
33.3 % to between 3.5 % and 4.1 %. For C
D
equal to 50 and 500 the errors become higher,
for C
D
= 500 they are close to the ones obtained with the homogenous model.
Table 2: Comparison between the theoretical pressures and the pressures obtained from
the simulations, at x = 0 (inlet), at t = 0.35 s. The errors relative to the theoretical values
are shown inside brackets.
P
inhom
P
inhom
(grav) P
inhom
-P
inhom
(grav)
Theory 573.7 355.3
218.4
P
CFX
P
CFX
(grav) P
CFX
-P
CFX
(grav)
570.2 360.7
209.5
C
D
= 0.05
(-0.6 %)
(+1.5 %)
(-4.1 %)
571.1 360.3
210.8
C
D
= 0.5
(-0.4 %)
(+1.4 %)
(-3.5 %)
570.6 361.0
209.5
C
D
= 5
(-0.5 %)
(+1.6 %)
(-4.1 %)
555.4 368.7
186.8
C
D
= 50
(-3.2 %)
(+3.8 %)
(-14.5 %)
537.4 376.1
161.4
C
D
= 500
(-6.3 %)
(+5.9 %)
(-26.1 %)
The evolutions with time of the inlet pressure obtained from simulations
(P
CFX
@Inlet) and of its theoretical value (P
inhom
@Inlet) obtained from equation (93), are
shown in Fig. 14a-14d. The error of the theoretical value according to the CFX
inhomogeneous model, where the air hydrostatic pressure and the air viscous effects are
neglected, relative to the full theoretical value (P
in
) obtained by equation (55), defined as:
inhom
in
in
P
P
Inlet
error1
P
@
(96)
and the error of the inlet pressure obtained from simulation, relative to the inhomogeneous
model theoretical value, defined as:
inhom
CFX
inhom
P
Inlet
P
Inlet
error2
P
Inlet
@
@
@
(97)
are also presented.

Simulation of the filling process
39
(a)
0.01%
0.10%
1.00%
10.00%
100.00%
0.0
0.1
0.2
0.3
0.4
t [s]
error
0
200
400
600
800
P@Inlet [Pa]
error2
error1
PCFX@Inlet
Pinhom@Inlet
(b)
0.01%
0.10%
1.00%
10.00%
100.00%
0.0
0.1
0.2
0.3
0.4
t [s]
error
0
200
400
600
800
P@Inlet [Pa]
error2
error1
PCFX@Inlet
Pinhom@Inlet
(c)
0.01%
0.10%
1.00%
10.00%
100.00%
0.0
0.1
0.2
0.3
0.4
t [s]
error
0
200
400
600
800
P@Inlet [Pa]
error3
error2
error1
PCFX@Inlet
Pinhom@Inlet
(d)
0.01%
0.10%
1.00%
10.00%
100.00%
0.0
0.1
0.2
0.3
0.4
t [s]
error
0
200
400
600
800
P@Inlet [Pa]
error1
error2
PCFX@Inlet
Pinhom@Inlet
Fig. 14: Evolution with time of the inlet pressure obtained from simulation and of its theoretical value.
(a) C
D
= 0.05; (b) C
D
= 0.5; (c) C
D
= 5; (d) C
D
= 50.

Simulation of the filling process
40
For
C
D
= 0.05 and C
D
= 0.5, it is possible to observe that, although the pressure
obtained from CFX follows the theoretical value, it oscillates around the theoretical value.
The oscillations of the curve error2 are well visible in Fig. 14a-14b.
These time oscillations are very probably related to the spatial oscillations at the
flow front region shown in Fig. 13a-13d. Indeed, lower values of C
D
show higher
oscillations both on Fig. 13a-13d and on Fig. 14a-14d.
For
C
D
= 5, the pressure obtained from CFX and the theoretical one are nearly
superposed. For t = 0.032 s error2 becomes smaller than 5 %, and for t = 0.15 s smaller
than 1 %.
In Fig. 14c the error of the inlet pressure obtained from CFX relative to the full
theoretical model value:
CFX
in
in
P
P
Inlet
error3
error1 error2
error1 error2
P
u
@
(98)
is also shown. However, the curve relative to error2 is nearly the same as the one relative
to error3, meaning that the error introduced by the inhomogeneous model may be regarded
as irrelevant.
For C
D
= 5, the oscillation of the curve error2 is not as evident as for C
D
= 0.05 or
C
D
= 0.5, but it is still present. And, it is interesting to note that the period of this
oscillation, T(error2), corresponds to the period at which the flow front crosses the mesh
elements.
x
T error2
U
(
)
'
(99)
where
'x represents the x direction length of the mesh elements.
In Fig. 15, the computation time to complete 0.47 s of simulation, for the four
different values of C
D
, is shown:
0
2
4
6
8
10
12
0.01
0.1
1
10
100
C
D
CPU tim e [h]
CPU tim e
Fig. 15: Computation time for the simulations to complete 0.47 s.

Simulation of the filling process
41
For C
D
= 0.05, the necessary computing time is nearly 2.6 times longer than for
C
D
= 5. This occurs because, due to the observed oscillations in pressure for the smallest
values of C
D
, the number of iterations to achieve the convergence target within each time
step is larger than for C
D
= 5.
From all these presented results, it becomes clear that the value of C
D
= 5, among
the five analysed values, is the most suitable for this filling process.
2.3 Conclusions
The homogeneous multiphase flow model, which is the appropriate one to simulate free
surface flows, would supposedly be the appropriate one to simulate filling processes,
where the two phases are completely stratified and the interface is well defined. However,
it was shown that it is unable to give acceptable results for the case under analysis, the
filling of a thin space, which is the typical situation in any kind of injection moulding
process.
The problem was overcome with the implementation of the inhomogeneous model,
together with the no-slip boundary condition for the resin phase and the free-slip condition
for the air phase. As it was shown, the fact that the air shear stress on the walls is not taken
into account causes an insignificant error, as the air is considerably less viscous than the
RIM resins.
Five simulations with five different values for the drag coefficient, C
D
, were
performed, of which C
D
= 5 revealed to be the most suitable one. The error of the obtained
inlet pressure with this value, for t = 0.35 s, is merely 0.5 % relatively to theory.
Nevertheless, an important question remains unanswered: would this be also the
most suitable drag coefficient value for different meshes (different elements size and
different elements aspect ratio), different time steps, different thicknesses, different inlet
velocities or different material properties? However a very large amount of work would be
necessary to show how this "ideal" value varies with all these variables.

Simulation of the filling process
42

43
3 Simulation of the curing process
3.1 The cure reaction
The key issue of the RIM process is the chemical reaction between the two or more liquid
components, causing the mixture to polymerize. This chemical reaction, the so called cure
of the resin, is described by a Kamal-type equation [94]:
n
m
C
1
2
dC
S
k
k
C
1 C
dt
(100)
where k
1
and k
2
have an Arrhenius dependence on temperature, as follows:
E1
R T
1
1
k
A
e
(101)
E2
R T
2
2
k
A
e
(102)
and where C is the degree of cure, T the absolute temperature, R the universal gas constant
and A
1
, A
2
, E
1
, E
2
, m and n are model parameters.
The cure reaction leads to an increase of the mixture viscosity, which is expressed
by [5, 39]:
E
R T
A B C
g
g
C
A
e
C
C
P
§
·
¨
¸
©
¹
P
§
·
P
¨
¸
¨
¸
©
¹
(103)
where C
g
is the gel, or solidification, point of the mixture, and A
P
, E
P
, A and B are model
parameters. Obviously, equation (103) is meaningless for C
[C
g
, 1].
The reaction is exothermic, and the heat rate released per unit volume is
proportional to the rate of cure, and given by:
R
t
C
Q
Q
S
(104)
where Q
t
is total volumetric amount of heat produced by the complete reaction.
The degree and the rate of cure play, therefore, an important role in the whole RIM
process, as exemplified in Fig. 16. Consequently, numerical errors involving the degree of
cure will give rise to errors involving all the other variables, even if these are numerically
obtained with accuracy. CFX is a code designed to deal mainly with fluid flow and heat

Simulation of the curing process
44
transfer, and does not have any implemented model of the cure reaction; the degree of cure
has to be implemented as an additional variable.
Fig. 16: Influence of the degree and rate of cure in the RIM process.
3.2 Physical model
As the molecular diffusion for the RIM materials is very low, ~ 10
-11
m
2
/s [39], the cure is
dominated by the reaction and convection terms:
C
r
C
r
C U
r
S
t
(
)
(
)
D
D
D
D
D
D
w
U
x
U
U
w
JG
(105)
where S
C
is the rate of cure given by equation (100), and
D refers to the liquid (resin)
phase.
3.3 Numerical models
Equation (105) is integrated over each control volume, and discretely represented as:
ip
ip
C
ip
V
r
C dV
m
C
r
S
V
t
D
D
D
D
D
§
·
w ¨
¸
U
U
¨
¸
w ©
¹
¦
³
(106)
where m
Dip
is the discrete mass flow of phase
D through a face of the finite volume, given
by equation (89).

Simulation of the curing process
45
The transient term may be approximated by the first order Euler transient scheme
as:
0
0
0
V
r
C
r
C
r
C dV
V
t
t
D
D
D
D
D
D
§
·
U
U
w ¨
¸
U
|
¨
¸
w
'
©
¹
³
(107)
or by the second order Euler transient scheme as:
`
V
2
0
0
0
00
00
00
V
1
r
C dV
t
t
1
2
r
C
1
r
C
r
C
D
D
D
D
D
D
D
D
§
·
w
°
¨
¸
U
|
®
¨
¸
w
'
T T
°¯
©
¹
ª
º
T T U T
U
U
«
»
¬
¼
³
(108)
where
T is the quotient between the old and the new time step, defined by equation (36).
The degree of cure, C, at the integration points is obtained by:
ip
up
C
C
C R
E x
JG
(109)
where C
up
is the value of the degree of cure at the upwind node,
C
is its gradient, R
JG
is
the vector from the upwind node to ip, and is a blend factor. The high resolution scheme
computes locally to be as close as possible to 1 without violating boundedness principles
[83]. The compressive differencing scheme allows to be greater than 1, but reducing it as
much as necessary to still maintain boundedness [85].
In CFX a source term S
I
is linearised as:
v
p
n 1
S
S
S
I
I
I
I I
(110)
where S
Iv
is the source value and S
Ip
the linear source coefficient, both obtained from the
previous iteration, and
I
n-1
the value of
I also from the previous iteration. As convergence
is approached, (
I-I
n-1
) approaches zero and the value of S
Ip
becomes irrelevant and does
not affect accuracy of the converged solution. However, it may improve convergence and
stability. It should always be less or equal to zero and be set to the value [83]:
v
p
S
S
I
I
w
w I
(111)
If this linearization coefficient is positive, even if the neighbour coefficients are
positive, the centre point coefficient can become negative. Computationally, it is essential
to keep S
Ip
negative so that instabilities and physically unrealistic solutions do not arise.
Physically, a positive S
Ip
implies that, as
I increases, the source term increases, which, if a
removal mechanism is not present, will lead to an increase of
I, and so on [95].

Simulation of the curing process
46
3.4 Case study
The case under study is the one-dimensional (and isothermal) filling of a space initially
filled with air, with "fictitious" water entering with an inlet velocity of 0.1 m/s. The
densities of the water and air are 997 kg/m
3
and 1.185 kg/m
3
, respectively. The rate of cure
is chosen to be given simply by:
C
S
k 1 C
(
)
(112)
The geometry, mesh and boundary and initial conditions are shown in Fig. 17.
Fig. 17: Geometry, mesh and boundary and initial conditions.
To simulate the one-dimensional filling process, the top and bottom boundaries are set as
free-slip walls.
Equation (105) may also be written as:
C
d r
C
r
S
dt
D
D
D
D
U
U
(113)
and taking into account the mass conservation equation, equation (69), leads to:
C
dC
dt
S
(114)
Integrating equation (114) between t = 0 and t = t
R
, t
R
representing the time of residence,
with the condition C = 0 for t = 0, leads to:
R
k t
C
1 e
(115)
The time of residence is given, for the liquid phase, by:
R
x
t
u
(116)
and for the air phase, by:

Simulation of the curing process
47
R
t
t
(117)
The value of k was chosen to be 0.921, so that when the flow front reaches the
position x = 0.5 m, the value of C is 0.99.
For t = 4 s, the value of C as function of x, obtained with equations (115), (116) and
(117), is shown in Fig. 18. The vertical dashed line indicates the flow front position.
t = 4 s
0.0
0.2
0.4
0.6
0.8
1.0
0
0.1
0.2
0.3
0.4
0.5
x [m]
C
Fig. 18: Analytical value of C for t = 4 s.
The simulations were performed with three different Courant numbers, which is
expressed as:
u
t
Cr
x
'
'
(118)
and with the combinations of the first order Euler transient scheme (TrSc1) or the second
order Euler transient scheme (TrSc2) and the high resolution advection scheme (HR) or the
compressive advection scheme (Comp), for the discretization of the cure equation. The
definition of the compressive advection scheme for the cure equations has to be done by
editing the CFX Command Language (CCL) file.
The residuals convergence target which shall be achieved by all equations is 10
-5
,
up to a maximum of 50 iterations within each time step.
Initially there were problems in obtaining acceptable values of the degree of
cure, C. After consulting the CFX technical services, it was found that the cause was an
inconsistency in the discretization of the volume fraction multiplier in the source and the
transient terms in the cure equation, in CFX-5.6. The problem was overcome with a custom
solver provided by the CFX technical services. This problem has been fixed in the most
recent versions of the code.
The degree of cure, C, obtained from simulations is represented, for x between
0.3 m and 0.5 m, in Fig. 19a-19c.

Simulation of the curing process
48
(a)
t = 4 s ; Cr = 0.1
0.925
0.950
0.975
1.000
0.30
0.35
0.40
0.45
0.50
x [m]
C
Theoretical
TrSc1 Comp
TrSc1 HR
TrSc2 Comp
TrSc2 HR
(b)
t = 4 s ; Cr = 0.5
0.925
0.950
0.975
1.000
0.30
0.35
0.40
0.45
0.50
x [m]
C
Theoretical
TrSc1 Comp
TrSc1 HR
TrSc2 Comp
TrSc2 HR
(c)
t = 4 s ; Cr = 1.0
0.925
0.950
0.975
1.000
0.30
0.35
0.40
0.45
0.50
x [m]
C
Theoretical
TrSc1 Comp
TrSc1 HR
TrSc2 Comp
TrSc2 HR
Fig. 19: Values of C, for x between 0.3 and 0.5, obtained from simulations.
(a) Cr = 0.1; (b) Cr = 0.5; (c) Cr = 1.0.

Simulation of the curing process
49
The error defined as:
max
C
Ctheor
Error
Ctheor
(119)
where C designates the values from numerical simulations, Ctheor the theoretical values
obtained with equations (115), (116) and (117), and Ctheor
max
the maximum theoretical
value of C at the instant t, is represented in Fig. 21a-21c.
From these figures, Fig. 21a-21c, one observes that: the compressive advection
scheme creates artificial gradients on the values of C; for Cr = 0.1 the results obtained with
the first and the second order transient schemes become analogous and; for all the three
Courant numbers the best results are obtained with the combination of the second order
Euler transient scheme with the high resolution advection scheme.
With these discretization schemes, simulations were performed in another mesh,
consisting of tetrahedral mesh elements instead of hexahedral elements, as shown in
Fig. 20b.
(a)
(b)
Fig. 20: Meshes used for simulations.
(a) meshA, formed by hexahedral elements; (b) meshB, formed by tetrahedral elements.
MeshA is the one used for the previous simulations, and meshB has exactly the same nodes
as meshA, but it is formed by tetrahedral elements.
Table 3: Characteristics of each mesh.
Nr. of mesh
nodes
Nr. of mesh
elements
Type of
mesh elements
meshA 612
250 Hexahedral
meshB 612
1250 Tetrahedral

Simulation of the curing process
50
(a)
Error = (C - Ctheor)/Ctheor
max
t = 4 s ; Cr = 0.1
-4%
-2%
0%
2%
0.0
0.1
0.2
0.3
0.4
0.5
x [m]
Error
TrSc1 Comp
TrSc1 HR
TrSc2 Comp
TrSc2 HR
(b)
Error = (C - Ctheor)/Ctheor
max
t = 4 s ; Cr = 0.5
-4%
-2%
0%
2%
0.0
0.1
0.2
0.3
0.4
0.5
x [m]
Error
TrSc1 Comp
TrSc1 HR
TrSc2 Comp
TrSc2 HR
(c)
Error = (C - Ctheor)/Ctheor
max
t = 4 s ; Cr = 1.0
-4%
-2%
0%
2%
0.0
0.1
0.2
0.3
0.4
0.5
x [m]
Error
TrSc1 Comp
TrSc1 HR
TrSc2 Comp
TrSc2 HR
Fig. 21: Error of the values of C obtained from numerical simulations.
(a) Cr = 0.1; (b) Cr = 0.5; (c) Cr = 1.0.

Simulation of the curing process
51
Errors obtained with both meshes and Cr = 0.5, are represented in Fig. 22. Despite
the fact that the two meshes have exactly the same nodes, because the mesh elements are
different, the integration points are not the same and therefore the results differ. The error
from meshA, where the hexahedral elements are aligned with the velocity vectors, is lower
than the error from meshB.
Error = (C - Ctheor)/Ctheor
max
t = 4 s ; Cr = 0.5
-4%
-2%
0%
2%
0.0
0.1
0.2
0.3
0.4
0.5
x [m]
Error
meshA (hexa)
meshB (tetra)
Fig. 22: Errors obtained from both meshes.
3.5 Conclusions
CFX does not have any implemented model to solve the cure equation, which may,
although, be modelled as a transient convection-source equation for an additional variable
representing the degree of cure, C. However, as it is not a standard equation, various
discretization schemes were tested, and the combination of the second order Euler transient
scheme with the high resolution advection scheme revealed to be the most accurate to
approximate the values of C at the integration points.
It was also shown that the results obtained with a mesh made up of hexahedral
elements aligned with the velocity vector, are considerably more accurate than those
obtained with a mesh containing exactly the same nodes but formed by tetrahedral
elements.

Simulation of the curing process
52

53
4 Simulation of the RIM filling and curing stages
4.1 The energy equation
As mentioned at the end of Section 1, the objective of this work is to simulate the RIM
filling and curing stages using the CFX software package. The filling stage is characterized
by the simultaneous motion, cure and energy equations, after which, when the injection of
resin has stopped and the material velocities inside the mould fall to zero, the pure curing
stage, characterized by only the cure and energy equations, takes place.
Up to this point, only the motion and cure equations were focused, but in this
section also the energy equation is taken into consideration.
4.1.1 Physical model
It has been shown in this work, that only the inhomogeneous model produces accurate
results for the filling of thin spaces. Therefore it will be used to model the filling stage.
The conservation of mass and momentum equations were already described by
equations (69) to (76), in Section 2. The cure is described by equation (105), in Section 3.
The energy equation for phase
D (resin) may be obtained by applying the First Law of
Thermodynamics to an infinitesimal control volume, leading to [92]:
v
E
r
htot
p
r
r
htot
U
t
t
k
T
W
S
D
D
D
D
D
D
D
D
D
D
D
D
w
U
w
x
U
w
w
x
JJJG
(120)
where k
D
represents the resin thermal conductivity, W
v
D
the viscous work, S
E
D
heat
sources, and htot
D
the total enthalpy defined as:
p
htot
h
0 5 U
U
u
0 5 U
U
.
.
D
D
D
D
D
D
D
D
x
x
U
JJJG JJJG
JJJG JJJG
(121)
and u
D
is the internal energy.
As for multiphase situations, which is the case, CFX-5.6 does not model the
viscous and pressure work neither the kinetic energy effects, the full energy equation
becomes simply the thermal energy equation [83]:

Simulation of the RIM filling and curing stages
54
E
r
h
r
h
U
k
T
S
t
D
D
D
D
D
D
D
D
D
D
w
U
x
U
x
w
JJJG
(122)
For incompressible fluids, the enthalpy is expressed solely in terms of temperature [92]:
ref
ref
h
h
c
T
T
D
D
D
(123)
where c
D
is the specific heat, and ref indicates a reference state for enthalpy evaluation (the
default reference state in CFX is: h
ref
= 0 J/kg for T
ref
= 0 K [83, 92]). No distinction is
made between Cp
D
and Cv
D
, as for incompressible fluids:
Cp
Cv
c
D
D
D
(124)
The heat source, S
E
D
, is given by:
E
R
t
C
S
r
Q
r
Q
S
D
D
D
(125)
An energy equation for phase
E (air), similar to equation (122) but without a heat
source term, could also be modelled, what would allow the implementation of a symmetric
extra term in each energy equation, the interphase heat transfer, representing the heat
transfer between the phases at their interface. However, mainly because, except for the first
instants, the interface area is significantly small compared to the contact area between the
resin and the walls, but also because the air thermal conductivity is typically smaller than
that of RIM resins, the heat transfer between the resin and the air will be insignificant
compared to the heat transfer between the resin and the walls, and consequently the
interphase heat transfer term can be omitted on the energy equation of the resin without a
foreseen loss of accuracy. Due to this, and as the air temperature is not an important issue,
there is no need to model an energy equation for the air.
In the most recent versions of CFX, CFX-5.7 and ANSYS CFX 10.0, it is already
possible to use the full energy equation for multiphase flows [92, 93]. Although, as
previously mentioned, for typical RIM situations the viscous work can be neglected [39,
80], and the kinetic energy effects are only important for high speed flows. Therefore
equation (122) is perfectly valid to model the energy balance in RIM.
4.1.2 Numerical model
As all the other equations, the energy equation for phase is integrated over each control
volume. It is then discretely represented as [83, 96]:

Simulation of the RIM filling and curing stages
55
ip
ip
ip
V
x
y
z
R
ip
ip
r
h
dV
m
h
t
T
T
T
k
n
n
n
A
r
Q
V
x
y
z
D
D
D
D
D
D
D
§
·
w ¨
¸
U
¨
¸
w ©
¹
ª
º
§
·
w
w
w
«
»
¨
¸
w
w
w
©
¹
¬
¼
¦
³
¦
(126)
The transient term, as in all the equations except the volume fraction equations, is
approximated by the second order Euler transient scheme as:
`
V
2
0
0
0
00
00
00
V
1
r
h
dV
t
t
1
2
r
h
1
r
h
r
h
D
D
D
D
D
D
D
D
D
D
D
D
§
·
w
°
¨
¸
U
|
®
¨
¸
w
'
T T
°¯
©
¹
ª
º
T T U
T
U
U
«
»
¬
¼
³
(127)
where
T, the quotient between the old and the new time steps, is obtained from equation
(36).
On the advection term, the discrete mass flow of phase
D, m
Dip
, is obtained from
equation (89), and h
Dip
, from equation (39), being the blend factor
E computed according to
the high resolution advection scheme. On the diffusion term, the temperature derivatives at
the integration points are evaluated making use of shape functions as indicated in
equation (38).
4.2 Case study 1
4.2.1 Case description
The work presented in [39], as already mentioned in Section 1.2, has been considered by
many authors as a benchmark solution. It is used also in this work, in this section (Case
study 1) and in Section 4.3 (Case study 2), as a benchmark solution to validate the results
produced by CFX when simulating the whole RIM model, consisting of the mass,
momentum, volume fractions, cure and energy equations.
In this Case study 1, the experimental system 9/21/1 conducted by [39] is used.
This system was also numerically modelled in [39, 72].
The mould is a simple rectangularly shaped mould, placed vertically, with a full
gate at the bottom, as shown in Fig. 23. Its dimensions, together with the filling conditions,
are indicated in Table 4.

Simulation of the RIM filling and curing stages
56
Fig. 23: Mould geometry.
Table 4: Mould dimensions and filling conditions [39].
L Length 0.505
[m]
W Width 0.101
[m]
H Thickness 3.2
u10
-3
[m]
Q
in
Inlet flow rate
3.35
u10
-5
[m
3
/s]
T
in
Inlet
temperature
55.3
[ºC]
T
w
Walls
temperature
82.0
[ºC]
The chemical system used was a RIM type polyurethane, whose density and
thermal properties, assumed being constant, are specified in Table 5.
Table 5: Density and thermal properties of the chemical system [39].
U
Density 1000
[kg/m
3
]
Cp Specific
heat 1880
[J/(kg
K)]
k Thermal
conductivity 0.17 [W/(m
K)]
A second order cure kinetics is assumed to describe the chemical reaction. The
kinetic parameters are indicated in Table 6.
Table 6: Kinetic parameters of the chemical system [39].
k
o
Pre-exponential factor in cure equation
36
u10
3
[m
3
/(mol
s)]
E
Reaction activation energy
57.8
u10
3
[J/mol]
C
OHo
Initial concentration of reactive species
2.6
u10
3
[mol/m
3
]
H
R
Heat of reaction
83
u10
3
[J/mol]

Simulation of the RIM filling and curing stages
57
Therefore, for this resin the rate of cure, S
C
, is given by:
E
R T
2
C
o
OHo
dC
S
k
C
e
1 C
dt
(128)
and the released heat rate per unit volume of resin, Q
R
, by:
E
R T
2
2
R
t
C
R
OHo
C
R
o
OHo
Q
Q
S
H
C
S
H
k
C
e
1 C
(129)
The rheological parameters, necessary for the viscosity equation (103) are given in
Table 7.
Table 7: Rheological parameters of the chemical system [39].
A
P
Pre-exponential factor in viscosity equation
4.1
u10
-8
[kg/(m
s)]
E
P
Viscosity activation energy
38.3
u10
3
[J/mol]
C
g
Solidification (gel) point
0.85
[ ]
A
Constant in viscosity equation
4.0
[ ]
B
Constant in viscosity equation
-2.0
[ ]
The filling time, t
f
, is the mould volume divided by the volumetric inlet flow rate:
f
in
H W L
t
4 87 s
Q
.
(130)
4.2.2 Numerical issues
As the mould width is much larger than its thickness (W/H
| 32), the side walls effect may
be neglected, and a two-dimensional simulation may be performed without important loss
of accuracy.
Due to symmetry, only half of the geometry was modelled. As convergence
problems arise when the flow front gets close to the outlet boundary, the calculation
domain length was extended to 0.52 m instead of the 0.505 m mould length. The problem
was run with two different meshes, both with 208 elements on the longitudinal direction,
each 2.5 mm long, by 5 elements on the transverse direction in Mesh 1, Fig. 24a, and
10 elements in Mesh 2, Fig. 24b.
The inhomogeneous model is used, together with the walls no-slip boundary
condition for the liquid (resin) phase equations and the free-slip condition for the air phase
equations.

Simulation of the RIM filling and curing stages
58
(a)
(b)
Fig. 24: Detail of the meshes used to model the process. (a) Mesh 1: five mesh elements on the
transverse direction; (b) Mesh 2: ten mesh elements on the transverse direction.
At the inlet a parabolic velocity profile is imposed, with an average velocity of:
in
in
Q
U
0 1037 m s
H W
.
/
|
(131)
resulting in a Reynolds number of 13, according to equations (40) and (41) and based on
the resin viscosity at the inlet temperature and zero conversion, 5.05
u10
-2
kg/(m
s). At the
outlet, the boundary condition is zero pressure.
The reference density is that of the air, 1.185 kg/m
3
. The resin volume fraction is 1
at the inlet and zero at the outlet. It is assumed that the mixture has not reacted before it
enters the mould, and therefore the condition of zero degree of cure is imposed at the inlet.
For the energy equation, the boundary conditions are the temperatures indicated in
Table 4. In CFX, the energy flow by diffusion at an inlet boundary is assumed to be
negligible when compared to advection, and equated to zero [83].
The imposed initial conditions all over the domain are fluids at rest, pressure equal
to zero, temperature equal to T
in
and, due to the mould being initially filled with air, the
initial resin and air volume fractions are 0 and 1, respectively.
At the end of the filling stage the inlet velocity is reduced to zero in a time interval
of 0.01 s, corresponding only to 0.2% of the filling time, via a third-order polynomial as
represented in Fig. 25. This is more realistic than a sharp step or a ramp function, and
guarantees the temporal continuity of the term
wu/wt in the momentum equations. A lower
order polynomial would cause a discontinuity in this term, leading to a discontinuity of
wp/wx and therefore a discontinuity of the inlet pressure which is physically unrealistic.

Simulation of the RIM filling and curing stages
59
0
t [s]
U [m/s]
'
t = 0.01 s
t
f
t
f
-
'
t / 2
t
f
+
'
t / 2
U
in
Fig. 25: Inlet velocity reduced to zero, via a third-order polynomial, at the end of the filling stage.
According to equations (110) and (111), the source term in the cure equation was
linearised as:
C
Cv
Cp
n 1
S
S
S
C
C
(132)
where the source value, S
Cv
, was set to:
E
R T
2
Cv
o
OHo
S
k
C
e
1 C
(133)
and the linear source coefficient, S
Cp
, which is always negative, to:
E
R T
Cv
Cp
o
OHo
S
S
2 k
C
e
1 C
C
w
w
(134)
The source term in the energy equation could also be linearised relatively to
temperature in a similar way:
E
Ev
Ep
n 1
S
S
S
T
T
(135)
where the source value, S
Ev
, would be:
E
R T
2
2
Ev
R
o
OHo
S
r
H
k
C
e
1 C
D
(136)
and the linear source coefficient, S
Ep
:
E
R T
2
2
Ev
Ep
R
o
OHo
2
S
E
S
r
H
k
C
e
1 C
T
R T
D
w
w
(137)

Simulation of the RIM filling and curing stages
60
However, as the source coefficient, S
Ep
, would always be positive, the source term was not
linearised and only the source value, S
Ev
, was used.
The time step is 4
u10
-4
s, and the maximum number of iterations within each time
step is 50. Three simulations were run with Mesh 1, with the residuals convergence target
set to 10
-3
, 10
-4
and 10
-5
. With Mesh 2 only one simulation was performed, with the
residuals convergence target set to 10
-5
. Note that, as previously mentioned, the
convergence criterion for the volume fraction equations is automatically set by CFX to 10
times the specified convergence target.
All the governing equations kept being solved up to t = 4.95 s, and after this instant
the fluids are assumed to be completely at rest and only the energy and the cure equations
are solved. Note that when the velocity at the inlet boundary is zero, the energy flows by
diffusion and by advection are equated to zero, corresponding consequently to an adiabatic
boundary.
The computation times, in hours, to perform the four simulations are shown in
Fig. 26. The steep lines correspond to the filling stage when all the equations are solved,
while the low slope lines correspond to the curing stage when only the cure and the energy
equations are solved and the convergence target is achieved in only one or two iterations
within each time step.
0
24
48
72
96
120
144
0
2
4
6
8
10
12
Simulation time [s]
CP
U t
im
e
[
h
]
Mesh 2, Residuals Target = 1e-5
Mesh 1, Residuals Target = 1e-5
Mesh 1, Residuals Target = 1e-4
Mesh 1, Residuals Target = 1e-3
Fig. 26: Computation time versus simulation time.

Simulation of the RIM filling and curing stages
61
4.2.3 Results
All the presented results, in this and in the next case studies, were obtained with
conservative variable values, which are, actually, the results produced by the CFX-Solver.
That is, a variable value on a boundary node is not precisely equal to the specified
boundary condition (the boundary conditions are imposed on the integration points laying
on the boundary), and instead it represents the average variable value in the control volume
surrounding that node [83], Fig. 27. This seems to be more correct than using hybrid
variable values, obtained by taking the results produced by the CFX-Solver and
over-writing them on the boundary nodes with the specified boundary conditions [83].
Fig. 27: Two-dimensional representation of a boundary control volume.
All the results, except when the contrary is mentioned, were obtained with Mesh 1
and with the residuals convergence target set to 10
-5
.
In Fig. 28a-28e the degree of cure, and in Fig. 29a-29e the temperature, are
represented at five different instants: at t = 3 s (during the filling stage), at t = 4.87 s (the
end of the filling stage), at t = 6.8 s (when the maximum temperature is observed), at
t = 9 s and at t = 15 s.
Note that because of the length of the mould being much bigger than its thickness,
the longitudinal and transverse directions are represented at different scales in order to be
possible to properly visualize the results.
At the end of the filling stage, the maximum value of the degree of cure is 0.79,
which is close to the degree of cure at which the resin solidifies (gel or solidification
point), C
g
= 0.85, and the highest temperature is 151.6 ºC. At t = 6.8 s, the maximum
temperature during the whole process, 177.9 ºC, is achieved, and the resin which ended up

Simulation of the RIM filling and curing stages
62
(a)
X [m]
Y [mm]
0.4
0.3
0.2
0.1
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(b)
X [m]
Y [mm]
0.
1
0.2
0.3
0.4
0.5
0.6
0.7
0.7
9
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(c)
X [m]
Y [mm]
0
.2
0
.3
0
.4
0.5
0.
6
0
.7
0
.8
0.
85
0.
9
0.
96
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(gel
poi
nt)
(d)
X [m]
Y [mm]
0.7
0.8
0.85
0.9
0.
98
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(gel point)
(e)
X [m]
Y [mm]
0.8
0.85
0.9
0.98
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(gel point)
Fig. 28: Degree of cure. (a) t = 3 s; (b) t = 4.87 s, when the filling stage ends;
(c) t = 6.8 s; (d) t = 9 s; (e) t = 15 s. (Figure not to scale).

Simulation of the RIM filling and curing stages
63
(a)
X [m]
Y [mm]
60
70
80
90
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(b)
X [m]
Y [mm]
6
0
70
80
90
100
110
120
1
3
0
1
4
0
1
5
0
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(c)
X [m]
Y [mm]
8
0
9
0
90
100
110
120
130
140
150
160
170
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(d)
X [m]
Y [mm]
160
150
140
130
120
110
100
90
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(e)
X [m]
Y [mm]
14
0
130
120
110
100
90
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
Fig. 29: Temperature [ºC]. (a) t = 3 s; (b) t = 4.87 s, when the filling stage ends; (c) t = 6.8 s, when the
maximum temperature (177.9 ºC) is observed; (d) t = 9 s; (e) t = 15 s. (Figure not to scale).

Simulation of the RIM filling and curing stages
64
at the forward centre part of the mould (about 42 %) has already solidified. At t = 9 s and
at t = 15 s, the temperature contour lines are nearly parallel to the mould walls, indicating
that at these instants the most important energetic process is heat conduction from the
mould centre to the walls. About 81 % of the resin has already solidified at t = 9 s, and
93 % at t = 15 s, remaining just a thin layer next to the mould walls yet to solidify.
In Fig. 29b, the use of conservative temperature values is well visible. The resin
temperature next to the wall near to the flow front, ~140 ºC, is quite above the wall
temperature, 82 ºC. This is due to the fountain flow: at the flow front the fluid on the centre
part of the mould at high temperature moves towards the walls. This behaviour is
mentioned in [2] as "thermal shock", and it also happens at the beginning of the filling
stage when the resin temperature is lower than the mould walls temperature. It is
characteristic of various injection moulding processes, and results in high fluid temperature
gradients next to the mould walls.
The minimum value of the degree of cure at each instant is presented in Fig. 30. It
corresponds, for t below 8.1 s to the degree of cure at the centre position of the inlet
(x, y = 0, 0 [mm]), and for t above 8.1 s to the degree of cure at the location of the inlet in
contact with the mould walls (x, y = 0, 1.6 [mm]). This can also be observed in
Fig. 28b-28d (t = 6.8 s, t = 9 s and t = 15 s), and it occurs because the temperature at the
mould walls is always 82 ºC and therefore the cure reaction is isothermal, while at the inlet
centre the cure reaction is nearly adiabatic and the temperature is 55.3 ºC during the filling
stage but increasing during the curing stage up to 167 ºC at t = 9.7 s.
0.0
0.2
0.4
0.6
0.8
1.0
0
5
10
15
20
25
30
t [s]
D
e
gr
e
e
of
c
u
re
Degree of cure at (x,y)= (0,0) [mm]
Degree of cure at (x,y)=(0,1.6) [mm]
Minimum value of the degree of cure
(gel point )
Fig. 30: Minimum value of the degree of cure along the time.
According to the simulation, the totality of the resin has solidified 22.8 s after the
injection of resin has started.

Simulation of the RIM filling and curing stages
65
In Fig. 31, the degree of cure at those points is compared with the degree of cure
obtained by an adiabatic cure reaction with an initial temperature equal to T
in
, 55.3 ºC, and
by an isothermal cure reaction at a temperature equal to T
w
, 82 ºC.
0.0
0.2
0.4
0.6
0.8
1.0
0
5
10
15
20
25
30
t [s]
D
e
gr
e
e
of
c
ur
e
Degree of cure at (x,y)= (0,0) [mm]
Adiabatic cure (Ti = 55.3 ºC)
Degree of cure at (x,y)=(0,1.6) [mm]
Isothermal cure at T = 82 ºC
(gel point )
Fig. 31: Comparison of the degree of cure at (x, y) = (0, 0) [mm] and at (x, y) = (0, 1.6) [mm] with the
degree of cure obtained by an adiabatic cure reaction and by an isothermal cure reaction.
The curve for the isothermal cure reaction is obtained by integrating equation
(128) for a constant temperature, with the initial condition C = 0 for t = t
f
, which leads to:
E
R T
o
OHo
f
1
C
1
1 k
C
e
t
t
(138)
where T is equal to 355.15 K (82 ºC), and t
f
is the filling time, 4.87 s.
The curve corresponding to the adiabatic cure reaction is obtained by integrating
the energy equation (122) without the convection and diffusion terms, with the initial
condition T = Ti for C = 0, where Ti is the initial temperature, which gives:
r
OHo
H
C
T
C
Ti
Cp
U
(139)
The above expression for T is introduced in equation (128), leading to:
^
`
H C
r
OHo
Cp
E
R Ti
C
2
o
OHo
C
K
C
1 C
e
t
U
ª
º
«
»
¬
¼
w
w
(140)
and this last equation is numerically integrated, also with the initial condition C = 0 for
t = t
f
. Ti is equal to T
in
, 328.45 K (55.3 ºC).

Simulation of the RIM filling and curing stages
66
The term H
r
C
OHo
/
UCp in equation (139) is known as the adiabatic temperature
rise,
'T
ad
[5, 39], as it corresponds to the temperature rise due to a complete adiabatic cure
reaction. For this particular resin,
'T
ad
= 114.8 ºC.
It is very interesting to notice that for an initial temperature of 55.3 ºC, a complete
adiabatic cure reaction would lead to a resin temperature of 170.1 ºC (55.3 ºC + 114.8 ºC),
which is below the maximum temperature observed during the whole process, 177.9 ºC.
This means that, although the walls temperature is significantly below these values, the
part of the resin which reached that temperature had to have a positive energy transfer
balance with its surroundings along its path. This may be explained partially by the walls
temperature being above the resin inlet temperature, and partially by the resin distribution
due to the fountain flow effect, as may be seen in Fig. 29b, where there is a gradient of
temperature from "older" resin, located between the centre and the wall, to "younger"
resin, located at the mould centre.
This clearly demonstrates that the adiabatic temperature rise cannot be taken as an
upper limit for estimating the maximum temperature observed during the RIM process, and
it attests the importance of the accurate simulation of the fountain flow behaviour and the
complexity of the simultaneous flow motion, cure reaction and heat transfer processes.
From equations (138) and (140), one concludes that, for the adiabatic and for the
isothermal cure reactions, the time to reach a given degree of cure is determined, besides
the resin properties, by the filling time, t
f
, and by the initial temperature for the adiabatic
cure, or by the constant temperature for the isothermal cure.
For the resin under analysis, the time (t-t
f
) to reach the solidification point
(C = C
g
= 0.85) and the time to reach the degree of cure C = 0.9, for an isothermal cure
reaction are represented in Fig. 32a as function of the temperature, and for an adiabatic
cure reaction in Fig. 32b as function of the initial temperature.
The curves in Fig. 32a are obtained from:
E
R T
f
o
OHo
C
t
t
1 C
K
C
e
(141)
derived from equation (138), while the curves in Fig. 32b are obtained by solving
numerically equation (140).
The part can be removed from the mould when the totality of the resin has reached,
typically, the degree of cure C = 0.9 [97, 98], such that enough molecular weight or
network structure has been built. Therefore the cycle time will be determined by the
longest between the isothermal and the adiabatic cure reactions times to reach that degree
of cure, which, for common RIM operation temperatures, is very likely to be the one
corresponding to the isothermal reaction. That is, unless the mould walls temperature, T
w
,
is much higher than the inlet temperature, T
in
, the last portion of resin to solidify and to

Simulation of the RIM filling and curing stages
67
reach the degree of cure such that the part may be removed, will be that in contact with the
inlet and the mould walls.
t - t
f
0
10
20
30
40
50
60
70
80
90
100
110
120
T [ºC]
t [s]
time to reach C = 0.9
time to reach C = Cg
(a)
t
- t
f
0
2.5
5
7.5
10
12.5
40
50
60
70
Ti [ºC]
t [s]
time to reach C = 0.9
time to reach C = Cg
(b)
Fig. 32: Time to reach the solidification (gel) point (C = Cg = 0.85) and to reach C = 0.9.
(a) Isothermal cure reaction at temperature T; (b) Adiabatic cure reaction with initial temperature Ti.
For the conditions under analysis, the part may be removed from the mould
approximately 35 s (4.87 s + 30.47 s) after the injection of resin has started. This time may
be considerably reduced, as can be observed in Fig. 32a, by increasing T
w
. However,
increasing the mould walls temperature will cause higher temperatures at the part centre,
which must not exceed the polymer degradation point. A temperature of 200 ºC is probably
an upper limit for polyurethanes [5]. The maximum temperature during the process may
only be predicted by numerical simulation.
The viscosity at the end of the filling stage is shown in Fig. 33. As the viscosity is
solely a function of the degree of cure and of the temperature, as expressed by equation
(103), its distribution could be foreseen by analysis of Fig. 28b and Fig. 29b: the highest
values of viscosity occur where the degree of cure is higher and the temperature lower.
X [m]
Y [mm]
0.05
0.05
0.1
0.3
0.5
0.7
1
1.5
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
Fig. 33: Viscosity [kg/(m
s)] at t = 4.87 s, when the filling stage ends. (Figure not to scale).

Simulation of the RIM filling and curing stages
68
This viscosity distribution causes the pressure distribution along the longitudinal
direction and its gradient represented in Fig. 34. The highest absolute values of the
pressure gradient correspond obviously to the position where the highest values of
viscosity are located, dropping to zero after the flow front.
0.0E+00
5.0E+03
1.0E+04
1.5E+04
2.0E+04
2.5E+04
0
0.1
0.2
0.3
0.4
0.5
x [m]
P [Pa]
-2.0E+05
-1.6E+05
-1.2E+05
-8.0E+04
-4.0E+04
0.0E+00
dP/dx [Pa/m]
P
dP/dX
Fig. 34: Pressure distribution along the x direction, just before the end of the filling stage.
The evolution of the inlet pressure with time obtained with CFX, together with the
numerical and experimental results obtained by Castro and Macosko [39], and the pressure
for the case of constant viscosity, are represented in Fig. 35.
0.0E+00
5.0E+03
1.0E+04
1.5E+04
2.0E+04
2.5E+04
0
1
2
3
4
5
t [s]
P [Pa]
CFX
[39] Numerical
[39] Experimental
Press(visc0)
Fig. 35: Evolution of the inlet pressure with time.

Simulation of the RIM filling and curing stages
69
The pressure for the case of constant viscosity is given by:
2
in
in
in
2
12 U
P visc0
inlet
t
g U
t
H
(
) @
P
U
(142)
where
P
in
indicates the initial mixture viscosity:
2
in
C
0 T
Tin
5 05 10
kg m s
(
;
)
.
/(
)
P P
|
u
(143)
In Fig. 35, the results obtained with CFX are in good agreement with the
experimental and numerical results obtained in [39].
The evolution of the inlet pressure obtained with the three residuals target values is
shown in Fig. 36. The results obtained with 10
-5
and 10
-4
are practically superposed, while
the results obtained with 10
-3
are somewhat different. The error of the inlet pressure
obtained with the residuals target set to 10
-4
relatively to that obtained with 10
-5
is always
below 1 %, and most of the time below 0.1 %, while the error of the inlet pressure obtained
with 10
-3
rises above 40 % near the end of the filling stage
0.0E+00
5.0E+03
1.0E+04
1.5E+04
2.0E+04
2.5E+04
3.0E+04
0
1
2
3
4
5
t [s]
P [Pa]
Residuals Target = 1e-3
Residuals Target = 1e-4
Residuals Target = 1e-5
Fig. 36: Inlet pressure obtained with the different residuals target values.
The reason why for the first 3 seconds the inlet pressure follows a straight line, is
that as the cure reaction takes place, the degree of cure increases, but due to the released
heat also the temperature increases, leading to equilibrium in the viscosity equation. This
behaviour can be seen in Fig. 37, where the degree of cure and the viscosity are
represented for the case of adiabatic cure reaction with initial temperature equal to T
in
,
55.3 ºC.

Simulation of the RIM filling and curing stages
70
0.0
0.2
0.4
0.6
0.8
1.0
0
1
2
3
4
5
t [s]
P
[kg/(m
s)]
0.0
0.2
0.4
0.6
0.8
1.0
C
Viscosity
Degree of cure
Fig. 37: Viscosity and degree of cure along the time for an adiabatic cure reaction.
The viscosity along the time, for the case of adiabatic cure reaction, with various
initial temperatures, is represented in Fig. 38. The curves were obtained from equation
(103), together with equation (139) and the numerical integration of equation (140).
The initial viscosity is as lower as higher are the values of Ti. However, the time
when a rapid increase of viscosity is observed is as shorter as higher are the values of Ti.
0.0
0.2
0.4
0.6
0.8
1.0
0
1
2
3
4
5
6
7
t [s]
P [kg/(ms)]
Ti = 65ºC
Ti = 60ºC
Ti = 55.3ºC
Ti = 50ºC
Ti = 45ºC
Fig. 38: Viscosity along the time for adiabatic cure reactions with various initial temperatures.
From Fig. 38, and for the mould under analysis, which is filled in 4.87 s, one may
expect for an inlet temperature of 50 ºC, a nearly constant viscosity (
P
in
=
6.36×10
-2
kg/(m
s)) during the filling stage, and hence a linear evolution of the inlet
pressure along the time resulting in 8.8 KPa at the end of the filling stage, which is almost
one third of the predicted inlet pressure for T
in
= 55.3 ºC.

Simulation of the RIM filling and curing stages
71
From the numerical integration of equation (140) or from Fig. 32b one concludes
that, for Ti = 50 ºC, the time t-tf for an adiabatic cure reaction to reach C = 0.9 is 6.45 s.
From equation (141) or from Fig. 32a one concludes that the necessary temperature for an
isothermal cure reaction to reach that degree of cure in approximately the same time is
112 ºC (for T = 112 ºC, t-t
f
= 6.63 s). Therefore it is expected that, if the inlet and the
mould walls temperatures were modified to 50 ºC and to 112 ºC, respectively, the time for
the part may be removed would be reduced from 35 s to 11.1 s (4.47 s + 6.63 s)
Fig. 39 shows the temperature on a line between the mould centre (y = 0) and the
mould wall (y = H/2), at x = 0.53L, during the filling and curing stages. The grid represents
the temperature evolution along the time at points spaced by 0.1
u(H/2), and temperature
profiles on the y-direction at every 0.5 s.
Fig. 39: Temperatures during filling and curing, at x = 0.53L
| 0.268 m. The white line
represents the instant when the filling stage ends (t = 4.87 s).
Due to the fountain flow, the temperature at the flow front is nearly uniform, and
therefore when it reaches x = 0.53L, at t = 2.58 s, the temperature is almost constant along
the mould thickness. Between this instant and the end of the filling stage (represented by a
white line in the figure), the lowest temperature is verified at the point located in the centre
plane, where the velocity and hence the energy convection term have their highest values.
The maximum temperature is verified at y = 0.7
u(H/2).

Simulation of the RIM filling and curing stages
72
The temperature profile at the end of the filling stage shows clearly the difficulty of
capturing the temperature curvature with only five mesh elements on the thickness
direction.
After the filling stage has ended and the convection terms turned into zero, the
temperatures at the centre points rapidly increase, reaching the maximum approximately at
t = 7 s. After that, the heat lost by conduction to the mould walls is higher than the heat
released by the chemical reaction, and the temperatures start to decrease.
Fig. 40 shows the comparison between the results obtained with CFX, numerical
and experimental results obtained by Castro and Macosko [39], and numerical results
obtained by Lo [72]. The temperatures obtained with CFX, with Mesh 1, are, for most of
the time, slightly higher than all the other results, quite probably due to the five mesh
elements in the mould thickness direction being unable to capture very accurately the
temperature gradients on that direction, leading to underestimation of the heat conduction
from the mould interior to the mould walls. This conclusion may be strengthened by the
observed lower temperatures obtained with Mesh 2 (10 mesh elements on the thickness
direction), which are closer to the numerical results obtained in [39]. However, as Mesh 2
has twice the mesh elements and 11/6 times the mesh nodes than Mesh 1, the computation
time is nearly the double, as could be observed in Fig. 26.
50
70
90
110
130
150
170
2
4
6
8
10
12
t [s]
T [ºC]
CFX (Mesh 1)
CFX (Mesh 2)
[39] Numerical
[39] Experimental
[72] Numerical
Center
0.6
0.7
0.9
Tw
Tin
Fig. 40: Temperatures at x = 0.53L. Numerical results obtained with CFX, with Mesh 1 and
Mesh 2, at the centre plane, at y/(H/2) = 0.6, 0.7 and 0.9. Numerical results obtained in [39] at
the centre plane, at y/(H/2) = 0.6, 0.7 and 0.9. Experimental results obtained in [39] at
y/(H/2) = 0.6. Numerical results obtained in [72] at y/(H/2) = 0.6.

Simulation of the RIM filling and curing stages
73
Nevertheless, from Fig. 40 one may conclude that there is a general good
agreement between the results obtained with CFX, even those obtained with Mesh 1, and
the other numerical results. When compared with the experimental results at y/(H/2) = 0.6,
the results obtained with CFX seem to be as good as the other numerical results.
The comparison between the temperatures obtained with the three residuals target
values is shown in Fig. 41. The temperatures obtained with 10
-5
and 10
-4
are nearly
superposed, while the ones obtained with 10
-3
are somewhat different.
70
80
90
100
110
120
2.5
3
3.5
4
4.5
5
t [s]
T [ºC]
Residuals Target = 1e-3
Residuals Target = 1e-4
Residuals Target = 1e-5
Center
0.6
0.9
Tw
Fig. 41: Temperatures at x = 0.53L, at the centre plane, at y/(H/2) = 0.6 and at
y/(H/2) = 0.9, obtained with the three different residuals target values.
4.3 Case study 2
4.3.1 Case description and numerical issues
This case is the experimental system 9/21/2 conducted and numerically modelled by [39].
The mould, as in Case 1, is a simple rectangularly shaped mould with a full gate at the
bottom. Its dimensions according to Fig. 23, and filling conditions are indicated in Table 8.
Table 8: Mould dimensions and filling conditions [39].
L Length 0.487
[m]
W Width 0.101
[m]
H Thickness 3.2
u10
-3
[m]
Q
in
Inlet flow rate
2.75
u10
-5
[m
3
/s]
T
in
Inlet
temperature
54.0
[ºC]
T
w
Walls
temperature
70.0
[ºC]

Simulation of the RIM filling and curing stages
74
The resin is the same as in Case 1, whose properties were reported in Tables 5-7.
The filling time, t
f
, is equal to:
f
in
H W L
t
5 72 s
Q
.
(144)
Because of the same reasons reported in Case 1, a two-dimensional simulation is
performed and only the top half of the mould is modelled. The calculation mesh is Mesh 1
used in Case 1, with 5 mesh elements on the transverse direction, represented in Fig. 24a.
A parabolic velocity profile, with an average value of:
in
in
Q
U
0 0851 m s
H W
.
/
|
(145)
is imposed at the inlet boundary. This results, according to equations (40) and (41), in a
Reynolds number of 10, based on the resin viscosity at the inlet temperature and zero
conversion, 5.34
u10
-2
kg/(m
s).
The major differences between this Case 2 and Case 1 are a lower inlet flow rate,
leading to a longer filling time, and the 12 ºC lower walls temperature.
At t = t
f
= 5.72 s the inlet velocity is set to zero via a third-order polynomial, in a
time interval of 0.01 s, as in the previous case, but all the equations keep being solved up
to t = 5.85 s. After that the fluids are assumed to be at rest and only the energy and the cure
equations are solved.
All the other numerical issues, including time step, are the same as employed in
Case 1. Furthermore, as in Case 1, three simulations were performed with the residuals
convergence target set to 10
-3
, 10
-4
and 10
-5
.
The computation time to perform the simulation, for each of the residuals target
value, is shown in Fig. 42, where the difference between when all the equations are solved
and when only the cure and the energy are solved is clear, as in Fig. 26.
0
24
48
72
96
0
2
4
6
8
10
12
Simulation time [s]
CP
U t
im
e
[
h
]
Residuals Target = 1e-5
Residuals Target = 1e-4
Residuals Target = 1e-3
Fig. 42: Computation time versus simulation time.

Simulation of the RIM filling and curing stages
75
4.3.2 Results
All the following presented results were obtained with the residuals convergence target set
to 10
5
, except when the contrary is mentioned.
The degree of cure, in Fig. 43a-43e, and the temperature in Fig. 44a-44e, are
represented at five different instants: at t = 4 s (during the filling stage), at t = 5.72 s (at the
end of the filling stage), at t = 7.6 s (when the maximum temperature is observed), at
t = 10 s and at t = 15 s. The longitudinal and transverse directions are represented at
different scales in order to allow a proper visualization of the results.
According to the results obtained from the simulation, at the end of the filling stage
a small portion (~1 %) of the resin has already solidified (the degree of cure is above 0.85),
Fig. 43b. This simulation was only possible to perform by setting the maximum resin
viscosity to 10
3
kg/(m
s). At this instant, the highest value of temperature is 159.3 ºC. At
t = 7.6 s, the maximum temperature, 179.0 ºC, is observed, and 29 % of the resin has
already solidified. At t = 10 s and t = 15 s, the temperature contour lines are nearly parallel
to the mould walls, due to the heat conduction from the mould centre to the walls, and
67 % and 84 %, respectively, of the resin has solidified.
It is curious to notice that, although the inlet and the mould walls temperatures are
1.3 ºC and 12 ºC, respectively, lower than in Case 1, the maximum observed temperature is
1.1 ºC higher. This is probably due to the higher degree of cure and temperature achieved
at the end of the filling stage. As in Case 1, but now the difference being bigger, the
maximum temperature, 179 ºC, is above the temperature caused by a complete adiabatic
cure reaction, 168.8 ºC (54.0 ºC + 114.8 ºC).
The inlet pressure along the time obtained with CFX is represented in Fig. 45, and
shows a good agreement with the experimental and numerical results obtained in [39] also
shown in Fig. 45.
In Fig. 45, the dashed line indicated by Lengthening Reactor represents numerical
results obtained in [39] neglecting the fountain flow effect, while the other dashed line
represents results obtained in [39] with a fountain flow model. As CFX solves the full
three-dimensional equations there is no need for any fountain flow model, but it is
important to recognize that CFX is capable to model its effect even when the mesh
elements longitudinal dimensions are considerably larger than their transverse dimensions
as shown in Fig. 24a.
Fig. 46 shows the inlet pressure obtained with the three residuals target values. As
in previous comparisons, the results obtained with 10
-5
and 10
-4
are practically superposed,
while the results obtained with 10
-3
are different. But in this case the simulation performed
with 10
-3
, after the premature nearly asymptotic increase of pressure, diverged.

Simulation of the RIM filling and curing stages
76
(a)
X [m]
Y [mm]
0.1
0.2
0.3
0.4
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(b)
X [m]
Y [mm]
0.1
0.2
0.3
0.4
0.5
0.
6
0.
7
0
.8
0
.8
5
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(c)
X [m]
Y [mm]
0
.2
0
.3
0
.4
0
.5
0
.6
0.
7
0.
8
0.8
5 (g
el p
oin
t)
0.9
0
.9
6
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(d)
X [m]
Y [mm]
0.5
0.6
0.7
0.8
0.85 (gel p
oint)
0.9
0
.9
8
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(e)
X [m]
Y [mm]
0.7
0.8
0.85 (gel point)
0.9
0.98
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
Fig. 43: Degree of cure. (a) t = 4 s; (b) t = 5.72 s, when the filling stage ends;
(c) t = 7.6 s; (d) t = 10 s; (e) t = 15 s. (Figure not to scale).

Simulation of the RIM filling and curing stages
77
(a)
X [m]
Y [mm]
60
70
80
90
1
0
0
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(b)
X [m]
Y [mm]
60
70
80
90
10
0
1
1
0
1
2
0
1
3
0
1
4
0
1
5
0
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(c)
X [m]
Y [mm]
8
0
17
0
16
0
150
140
130
120
110
100
90
80
0