(Thesis in pdf available here: http://goo.gl/EZzlT6 and here: http://hdl.handle.net/10773/2419)
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 Ansys CFX, designed for numerical simulation of fluid flow and heat and mass transfer.
The Ansys 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.
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 [1523]. 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 [2426] 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 [3236], analytical [3739] 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 [5153].
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, 6066]. 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, 6771], and
3D simulations have already been reported in some works [7276].
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
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
(d)
X [m]
Y [mm]
1
6
0
15
0
150
140
130
120
110
100
90
80
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
80
1
3
0
0
0.1
0.2
0.3
0.4
0.5
0
0.4
0.8
1.2
1.6
Fig. 44: Temperature [ºC]. (a) t = 4 s; (b) t = 5.72 s, when the filling stage ends; (c) t = 7.6 s, when the
maximum temperature (179.0 ºC) is observed; (d) t = 10 s; (e) t = 15 s. (Figure not to scale).
Simulation of the RIM filling and curing stages
78
0.E+00
2.E+04
4.E+04
6.E+04
8.E+04
1.E+05
0
1
2
3
4
5
6
t [s]
P [Pa]
CFX
[39] Numerical
[39] Numerical - Lengthing Reactor
[39] Experimental
Fig. 45: Inlet pressure along the time.
0.E+00
2.E+04
4.E+04
6.E+04
8.E+04
1.E+05
0
1
2
3
4
5
6
t [s]
P [Pa]
Residuals Target = 1e-3
Residuals Target = 1e-4
Residuals Target = 1e-5
Fig. 46: Inlet pressure obtained with the different residuals target values.
From these comparisons between the results obtained with the three different
values of the residuals target, one observes that the results obtained with 10
-4
are nearly the
same as those obtained with 10
-5
but the computation time is more than 40 % lower.
Therefore it seems completely reasonable to accept the value 10
-4
as good compromise
between accuracy and computation time.
Fig. 47 and Fig. 48 show, respectively, the predicted degree of cure and
temperature at the end of the filling stage versus the vertical position y* = y/(H/2) for some
longitudinal positions x* = x/L, obtained with CFX and obtained in [39].
Simulation of the RIM filling and curing stages
79
C
y*
x*=
0.
05
x*
=0
.5
8
x
*=
1
x*
=0
.9
2
0
0.2
0.4
0.6
0.8
1.0
0
0.2
0.4
0.6
0.8
1.0
Fig. 47: Degree of cure, C, at the end of the filling stage. x* = x/L, y* = y/(H/2).
Solid lines: results obtained with CFX. Dashed lines: numerical results obtained in [39].
T[ºC]
y*
60
80
100
120
140
160
0
0.2
0.4
0.6
0.8
1.0
x*
=
0
x
*
=
0
.0
5
x
*
=
0
.5
8
x
* =
0.
92
x*
=
1
Fig. 48: Temperature [ºC] at the end of the filling stage. x* = x/L, y* = y/(H/2).
Solid lines: results obtained with CFX. Dashed lines: numerical results obtained in [39].
The results obtained with CFX are once again in good agreement with the results
obtained in [39].
Simulation of the RIM filling and curing stages
80
4.4 Case studies 3 and 4
4.4.1 Cases description
In both previous cases, because of the moulds large ratio width/thickness, two-dimensional
simulations could be performed without significant loss of accuracy.
Now, in the next two cases, geometries where all the three dimensions effects are
important and where the flow front splits and remerges during the filling stage, were
selected with the purpose of testing the ability of the software to model such situations.
The geometry used for Case 3, represented in Fig. 49, is the same used in [99] to model the
isothermal mould filling with a high viscous polymer melt. The geometry used for Case 4,
represented in Fig. 50, is just a slight modification, the rectangular hole was substituted by
a circular one.
The resin used in these two cases is again the same RIM type polyurethane used in
[39] and in both previous cases, whose properties were described in Tables 5-7. The inlet
and walls temperatures are the same as in Case 1, 55.3 ºC and 82 º C, respectively.
The filling time, t
f
, was chosen to be the same as in [99], 2 s, corresponding to the
inlet flow rates, Q
in
, and inlet velocities, U
in
, indicated in Table 9.
Fig. 49: Mould geometry, and dimensions in millimetres, used for Case 3.
Simulation of the RIM filling and curing stages
81
Fig. 50: Mould geometry, and dimensions in millimetres, used for Case 4.
Table 9: Filling conditions for Case 3 and Case 4.
Case 3
Case 4
V
Mould volume
2.960
u10
-5
m
3
2.869
u10
-5
m
3
t
f
Filling time
2 s
2 s
Q
in
= V / t
f
Inlet flow rate
1.480
u10
-5
m
3
/s 1.434
u10
-5
m
3
/s
U
in
= Q
in
/ A
in
Inlet velocity
0.370 m/s
0.359 m/s
T
in
Inlet temperature
55.3 ºC
55.3 ºC
T
w
Mould temperature
82.0 ºC
82.0 ºC
4.4.2 Numerical issues
Due to symmetry of the geometries and boundary conditions, only one fourth of the
geometries were modelled.
As previously mentioned, convergence difficulties arise when the flow front
approximates the outlet boundary. For this reason the calculation domain length was
extended from 0.18 m (mould length) to 0.19 m. The calculation meshes used to model
these cases are shown in Fig. 51 and Fig. 52. As in the previous cases, the mould half
thickness is divided into 5 elements, resulting in 6320 elements and 8190 nodes in Case 3,
and 6400 elements and 8304 nodes in Case 4. The extra 0.01 m length corresponds to four
rows of mesh elements, which has been shown to be enough to overcome the convergence
problems when the resin flow front comes close to the outlet boundary.
Simulation of the RIM filling and curing stages
82
Fig. 51: Mesh for Case 3 (6320 elements and 8190 nodes).
Fig. 52: Mesh for Case 4 (6400 elements and 8304 nodes).
The highest values of the Reynolds number are verified at the mould entrance zone,
where the cross section area is smaller and consequently the velocity is higher, which for
both cases, and based on the resin viscosity at the inlet temperature and zero conversion,
have approximately the value of 26.
As in both previous cases, the inlet velocity is set to zero at the end of the filling
stage via a third-order polynomial, as it was represented in Fig. 25. All the equations keep
being solved until t = 2.05 s, and after that the fluids are assumed to be at rest and only the
energy and cure equations are solved.
In these two cases the residuals convergence target was set to 10
-4
, since in the
previous cases it was shown to produce identical results to 10
-5
and reducing considerably
Simulation of the RIM filling and curing stages
83
the computation time, which becomes more and more important as the number of mesh
elements and nodes increase.
All the other numerical issues are the same employed in the previous cases. The
only exception is the time step in Case 4 after 2.05 s (when only the cure and energy
equations are solved), which was increased from 4
u10
-4
s to 4
u10
-3
s. This value is small
enough to obtain convergence on this phase of the process simulation and allows reducing
slightly the computation time, shown in Fig. 53.
0
24
48
72
96
120
0
1
2
3
4
5
6
7
8
Simulation time [s]
C
P
U
ti
m
e
[h
]
Case 3
Case 4
Fig. 53: Computation time versus simulation time.
4.4.3 Results
Fig. 54a-54i and Fig. 55a-55i show, for Case 3 and Case 4, respectively, the obtained resin
volume fraction distribution at the midplane (left side of the images) and at the mould wall
(right side of the images) during the filling stage.
The flow front is well captured at the midplane at every instant. However, at the
mould wall, when the flow front velocity is higher its definition becomes less accurate,
Fig. 54a and Fig. 55a, and Fig. 54d-54e and Fig. 55d-55e next to the side wall. This
behaviour is due to the value for the drag coefficient obtained in Section 2.2, C
D
= 5, used
in the inhomogenous momentum equations, being inadequate for these conditions. This
"ideal" value was obtained based on a mould thickness of 3.2 mm and an average velocity
of 0.1 m/s, while in these two last cases the mould thickness is 2 mm and in the entrance
zone the average velocity is around 0.36 m/s. This shows that the "ideal" value for C
D
is
not constant, and it varies at least with the velocity. A lower value should be used, where
and when the flow front velocity is higher, in order to decrease the coupling between the
phases' velocities and allow the air to leave the region to be occupied by the resin.
Simulation of the RIM filling and curing stages
84
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
Fig. 54: Volume fraction of resin at the midplane (left side) and at the mould wall (right side), during
the filling stage, for Case 3. The three contour lines indicate resin volume fractions of 0.1, 0.5 and 0.9.
(a) t = 0.1 s; (b) t = 0.2 s; (c) t = 0.35 s; (d) t = 0.45 s; (e) t = 0.55 s; (f) t = 0.8 s;
(g) t = 1.1s; (h) t = 1.4 s; (i) t = 2 s (end of the filling stage).
Simulation of the RIM filling and curing stages
85
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
Fig. 55: Volume fraction of resin at the midplane (left side) and at the mould wall (right side), during
the filling stage, for Case 4. The three contour lines indicate resin volume fractions of 0.1, 0.5 and 0.9.
(a) t = 0.1 s; (b) t = 0.2 s; (c) t = 0.35 s; (d) t = 0.45 s; (e) t = 0.55 s; (f) t = 0.8 s;
(g) t = 1.05s; (h) t = 1.4 s; (i) t = 2 s (end of the filling stage).
Simulation of the RIM filling and curing stages
86
However, this observed wrong behaviour is not expected to have a significant
influence on the solution. As may be seen in the figures, the regions where the resin
volume fraction is below 0.9 are small compared to the mould size.
Fig. 56a-56c and Fig. 57a-57c for Case 3, and Fig. 58a-58c and Fig. 59a-59c for
Case 4, show the degree of cure and the temperature at the midplane, at the bottom wall, at
some longitudinal and transverse planes, at the side walls, and at the flow front (defined as
the surface where the resin volume fraction is 0.5), obtained at three instants: at t = 1.05 s
(when the flow front remerges), at t = 2 s (end of the filling stage) and at t = 4.2 s (when
the maximum temperature is observed). To allow a better visualization of results, the
thickness direction was scaled up 8 times relatively to the other directions.
In both cases, when the flow fronts remerge, at t = 1.05 s, the degree of cure of the
streams meeting each other is around 0.05, Fig. 56a and Fig. 58a. If this value is above a
critical value, typically 0.40 [5], the material will not interpenetrate and a knit line will be
visible and weak. At this instant the highest values of the degree of cure are 0.22, Fig. 56a
and Fig. 58a, and the highest temperatures are 84.9 ºC, Fig. 57a and Fig. 59a, just slightly
above the walls temperature.
At the end of the filling stage, t = 2 s, the resin with higher degree of cure, Fig. 56b
and Fig. 58b, and higher temperature, Fig. 57b and Fig. 59b, is located next to the walls,
contrasting to the end of the of the filling stage in Case 1 (with the same inlet and wall
temperatures, but a longer filling time) where, at the end of the filling stage, the resin with
higher degree of cure and temperature was located between the mould centre and the walls,
Fig. 28b and Fig. 29b. At this instant, the highest degrees of cure are 0.37, for both cases,
and the highest temperatures 92.5 ºC and 92.7 ºC, for Case 3 and Case 4, respectively.
At t = 4.2 s the situation has inverted, the resin at the mould centre has higher
degree of cure, Fig. 56c and Fig. 58c, and its temperature is higher than at the walls,
Fig. 57c and Fig. 59c. At this instant the maximum temperatures during the whole process,
155.4 ºC, were observed in both cases. This temperature is 22.5 ºC below the maximum
temperature observed in Case 1, despite the fact that the resin and the inlet and the walls
temperatures are the same. This may be explained by the mould thickness being smaller,
increasing the heat conduction from the resin to the walls, but also by the shorter filling
time, and consequently lower temperatures at its end. At this instant the highest values of
the degree of cure are 0.90 and 0.89, for Case 3 and Case 4, respectively, but only 3.0 %
and 2.8 % of the resin has solidified, Fig. 56c and Fig. 58c.
At t = 5 s the percentage of resin that has solidified increases to 36 %, and at
t = 10 s to 78 %. At t = 20.7 s the last portion of resin reaches the gel point, C
g
= 0.85, and
at t = 32 s the totality of the resin has reached the degree of cure of 0.9 and the parts may
be safely removed from the moulds.
Simulation of the RIM filling and curing stages
87
(a)
(b)
(c)
Fig. 56: Degree of cure for Case 3. (a) t = 1.05 s, when the flow front remerges;
(b) t = 2 s, at the end of the filling stage; (c) t = 4.2 s (the black curves indicate C = C
g
).
(The thickness direction is scaled up 8 times relatively to the other directions).
Simulation of the RIM filling and curing stages
88
(a)
(b)
(c)
Fig. 57: Temperature [ºC] for Case 3. (a) t = 1.05 s, when the flow front remerges;
(b) t = 2 s, at the end of the filling stage; (c) t = 4.2 s, when the maximum temperature is observed.
(The thickness direction is scaled up 8 times relatively to the other directions).
Simulation of the RIM filling and curing stages
89
(a)
(b)
(c)
Fig. 58: Degree of cure for Case 4. (a) t = 1.05 s, when the flow front remerges;
(b) t = 2 s, at the end of the filling stage; (c) t = 4.2 s (the black curves indicate C = C
g
).
(The thickness direction is scaled up 8 times relatively to the other directions).
Simulation of the RIM filling and curing stages
90
(a)
(b)
(c)
Fig. 59: Temperature [ºC] for Case 4. (a) t = 1.05 s, when the flow front remerges;
(b) t = 2 s, at the end of the filling stage; (c) t = 4.2 s, when the maximum temperature is observed.
(The thickness direction is scaled up 8 times relatively to the other directions).
Simulation of the RIM filling and curing stages
91
From equation (141), according to the theoretical analysis previously presented in
Section 4.2, these times would be 21.2 s and 32. 5 s. The differences occur because the last
portions of resin to reach C = C
g
= 0.85 and C = 0.9 are in contact with the inlet and the
walls, and as CFX calculates the source terms based on the nodes values, which, for
boundary nodes, as previously mentioned, are not precisely equal to the imposed boundary
conditions but represent the average of the boundary control volume, the source term on
the cure equation is calculated based in a temperature slightly above the imposed wall
temperature. But the differences, being small (2.4 % and 1.5 %), shall not constitute a
concern, because in practice it will be impossible to keep the mould walls temperature
constant, and this fact will probably cause a higher source of error.
The obtained inlet pressures, for both cases, are shown in Fig. 60. The initial steep
parts of the curves, before approximately 0.1 s, correspond to when the flow front is at the
entrance zone moving at a higher velocity. Due to the slightly lower inlet velocity, the inlet
pressure in Case 3 before 0.65 s is a little lower than in Case 4, but after this instant,
because of the smaller cross section around the circular hole, it becomes slightly higher.
0.E+00
1.E+03
2.E+03
3.E+03
4.E+03
5.E+03
0
0.5
1
1.5
2
t [s]
P [Pa]
Case 3
Case 4
Fig. 60: Inlet pressure along the time.
4.5 Conclusions
Ten partial differential equations are necessary to model the RIM simultaneous filling and
curing processes. A small time step value and a large number of iterations within each time
step are required to achieve convergence, leading to long computing times, even for
Simulation of the RIM filling and curing stages
92
two-dimensional analyses and simple geometries with relatively small number of mesh
elements.
Nevertheless, two cases were simulated and the results obtained with CFX were
compared with available numerical and experimental results, showing a general good
agreement with them, what demonstrates that it is possible to simulate the complex RIM
simultaneous filling and curing processes using CFX, with a relatively good accuracy, even
for situations unlikely to happen in reality, where the gel point is nearly reached in one
case, and actually reached in the other, during the filling stage.
Simulations were performed with three different values for the residuals
convergence target. The ones performed with 10
-4
provided analogous results to the ones
with 10
-5
, though requiring a computation time around 40 % shorter, allowing one to
consider 10
-4
as a fair reasonable compromise between accuracy and computation
performance.
Two other cases, where all the three-dimensional effects are important and where
the flow fronts splits and remerges, were simulated. The flow front was always well
captured at the midplane, however, at the wall its definition becomes less accurate when its
velocity is higher. This occurs because the drag coefficient value used (previously obtained
based on other conditions) is too high for these conditions.
In order to try to reduce the computation time, these cases were run in parallel in
two PCs, but the computation times were longer than when run in serial. It is mentioned in
the CFX documentation [83, 92, 93] that, typically, a minimum of 30
u10
3
nodes, for
tetrahedral meshes, or 75
u10
3
, for hexahedral meshes, should be used per partition to see
improvements in the computational performance.
93
5 Conclusions
Synthetic polymeric materials are nowadays commonly used in all of the major market
sectors. They are very attractive from an economical point of view, but also by their
aptitude to be manufactured into complex parts and their relatively high weight specific
strength and stiffness,
Synthetic polymers can be classified in two major categories: thermoplastics and
thermosets. For thermoplastics, by far the largest volume, as the long molecules are not
chemically joined, they can be melted, solidified and remelted again. Thermosets upon
heating or mixing with appropriate reagents undergo an irreversible chemical reaction,
causing the short chain molecules to bound, and leading to the formation of a rigid
structure.
Amongst the various polymer processing techniques, injection moulding is one of
the most important, representing about one third of all manufactured polymer parts. In
TIM, the hot polymer is pushed at high pressure into a cold cavity where it undergoes
solidification by cooling. Reactive moulding is quite different, it uses polymerization,
instead of cooling, to form the solid part.
RIM is a process for rapid production of complex parts through the mixing and
chemical reaction of two or more components. When injection begins, the liquid
components, held in separated tanks, flow into a mixing chamber, the streams are
intensively mixed, and the mixture begins to polymerize, or cure, as it flows into the mould
cavity.
The end use applications for RIM can, at present, be found in a large variety of
forms, and although polyurethanes comprise the majority in RIM processing, there are also
other types of chemical systems that can be processed.
Due to RIM liquids' low viscosity, large parts can be produced with relatively
small machines, complex shapes with multiple inserts can be fabricated, 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, and mould clamping forces are much lower requiring lower
cost presses, opening up short runs applications and even prototype applications. Low
viscosities also open options in the use of reinforcements. RIM's temperatures are typically
lower than those for TIM, with less energy demands.
However, handling of reactive, and often hazardous liquids, requires special
equipments and procedures, and as some components freeze at room temperature, a
temperature controlled environment is required, increasing the costs. Because of the low
viscosity, moulds are difficult to seal, increasing flash, and low viscosity liquids penetrate
Conclusions
94
mould surfaces, requiring the use of release agents. Due to low pressure, it is difficult to
remove air from behind inserts and from corners, rendering vent locations extremely
important. If flow into the mould is too rapid, air may be entrained and large bubbles
appear in the final part, but moulds filled too slowly may lead to short shots.
The first attempts to study the filling stage in injection moulding took place in the
beginning of the 1950s with flow visualizations and tracer studies.
Numerical simulations basically started in the early 1970s. The first developments
were applied to the filling stage of simple tubular, circular and rectangular geometries,
allowing the flow to be assumed as unidirectional. The temperature was considered as two-
dimensional, one coordinate in the flow direction and the other in the thickness direction,
leading to the so-called 1½D approach, and finite difference techniques were used to solve
the set of equations.
The real breakthrough in injection moulding simulations came, in the beginning of
the 1980s, with the development of a general 2½D approach, combining finite elements
along the midsurface of the cavity with finite differences along the thickness direction.
However, based on the Hele-Shaw approximation, this approach was not able to represent
the complex flow kinematics of the flow front region, the fountain flow. The description of
this phenomenon was addressed by means of experimental, analytical and numeric
methods, leading to approximate models able to capture its kinematics without resolving
its complex 3D details.
To date, several commercial and research three-dimensional programs for injection
moulding simulations have been developed.
Early studies in RIM simulation were dedicated to the static analysis of heat
transfer and cure, assuming that the curing stage could be decoupled from the filling stage.
More realistic models were obtained by extending the 1½D approach to RIM. Extensions
to the 2½D approach were only reported more recently, in the 1990s, and some 3D
simulations have already been reported in some works.
Although many commercial codes are able to perform 3D injection moulding
simulations, the 2½D Hele-Shaw approach is still the standard numerical framework for
simulation of injection moulding. Because typically in injection moulding the part
thickness is much smaller than its overall dimensions and the polymer viscosity is high
(resulting in low Reynolds numbers), in the 2½D formulation the inertia effect, the velocity
component and thermal convection in the thickness direction, and, the velocity gradients
and the heat conduction in the flow directions, are neglected. Based on these
approximations, the continuity and momentum equations are simplified and merged into a
simple two-dimensional equation in terms of pressure and fluidity, and a simplified energy
equation is obtained. Due to those assumptions and simplifications, the computation time
can be considerably reduced relatively to full three-dimensional formulations. But although
its successful applications, the 2½D formulation has its limitations: in some cases the
Conclusions
95
inertia and three-dimensional effects may become significant enough to influence the flow,
for the RIM process, because of resins low viscosity, the inertia and gravity effects cannot
be omitted, and due to the importance of the fountain flow in RIM, a three-dimensional
formulation is expected to provide more accurate information than simple fountain flow
approximations.
A very simple case, the filling of a space between two parallel plates with a liquid
with constant density and viscosity, was studied. Applying the Newton's Second Law of
motion, the theoretical values of the pressure, and the theoretical values of the pressure
according to the CFX Homogenous Model, which omits the air hydrostatic pressure effect,
were derived. The case was simulated, with four different meshes, using the CFX
Homogeneous Model. This model, which assumes that both liquid and gaseous phases
share the same velocity field, should be the appropriate for this type of process, as the two
phases are completely stratified and the interface is well defined. In the simulations, the
resin-air interface was well captured within two mesh elements, proving the efficiency of
the advection and transient compressive discretization schemes. However, the interface is
always ahead of its theoretical position, and, much more important, it does not touch the
walls (there is always a layer of air between the resin and the walls) what obviously is not
in accordance with reality. Moreover, reducing the mesh elements' size, either on the
longitudinal direction or on the transverse direction, does not contribute to improvement of
this incorrect behaviour. Because of the resin volume fraction next to the walls being close
to zero, the computed viscosity is much smaller than its actual value, leading to a wrong
velocity profile and to an underestimation of the viscous effect contribution to pressure.
This is of major importance, as the injection pressure is one of the key parameters to be
predicted from numerical simulations, but also because, when solving the cure and the
energy equations, errors on the advection terms may lead to completely wrong final results.
The gravity and the viscous effects contribution to pressure, as well as the whole pressure,
obtained from the simulations were compared with theoretical values. The error of the inlet
pressure due to the viscous effect, obtained with the four different meshes, was comprised
between 30 % and 35 %, showing only a tiny decrease with the increase of the number of
mesh elements, which does not seem to be a solution for the problem. This behaviour has
already been mentioned by some authors, the physically correct no-slip boundary condition
applied to the air phase prevents the flow from touching the walls, and they claim that the
no-slip condition should be imposed only on the filled portion of the mould.
As CFX does not allow the implementation of conditional boundary conditions, the
only way to prescribe the no-slip boundary condition only on the filled portion of the
mould, and the free-slip condition on the empty portion, is by employing the
inhomogeneous model, which assumes that each phase has its own velocity field and they
interact via interphase transfer terms. The no-slip condition was applied to the resin
equations and the free-slip condition to the air equations. The interphase momentum
Conclusions
96
transfer term, which links the two sets of momentum equations, is function of the phases'
relative velocities, their density and volume fractions (it tends to zero as the volume
fraction of one of the phases tends to zero, being meaningful only at the interface), and two
parameters that must be defined by the user: a non-dimensional drag coefficient, C
D
, and
an interface length scale, d
DE
. However, as it was shown, the individual values of these
parameters do not matter, but their quotient does. Therefore, the parameter d
DE
was kept
constant (with the CFX default value) and simulations were performed with different
values of C
D
, from 0.05 to 500.
Contrasting with the simulations employing the homogenous model, in the
simulations performed with C
D
set to 0.05, 0.5 and 5, the resin-air interface touches the
walls and is in very good agreement with its theoretical position. But when C
D
increases to
50, the behaviour becomes similar to that of the homogenous model, what is explained by
the increase of the interphase momentum transfer term causing the phases' flow fields to
equalize. As it was done for the simulations employing the homogeneous model, the
results obtained from the simulations employing the inhomogeneous model were compared
with theoretical values. For C
D
equal to 0.05, 0.5 and 5, the errors were considerably
reduced relatively to the ones obtained with the homogenous model (the error of the
viscous contribution to pressure decreased from 33 % to 4 %). The evolution of the
obtained inlet pressures with time was compared with theory and, despite following very
closely the theoretical pressure, the pressures obtained with the three lowest values of C
D
exhibit oscillations with the same period at which the interface crosses the mesh elements.
The amplitudes of these oscillations are bigger as smaller is the value of C
D
, and for
C
D
= 5, although still present, the oscillations nearly vanish. The shortest computation time
was obtained with C
D
equal to 5, what together with the presented results, allows one to
conclude that, among the five different analyzed values, C
D
= 5 is the most suitable for the
studied case. However, this leaves an unanswered question: would this value also be the
most suitable one for other physical and/or numerical conditions?
The cure of the resin can be modelled in CFX as a transient convection-source
equation for an additional variable, the degree of cure. But, as this is not a standard
equation in CFX, and as errors involving the degree of cure will give rise to errors
involving the other variables, various schemes for the discretization of the transient and
advection terms of the cure equation were tested. A one-dimensional filling case with a
first order rate of cure (allowing one to obtain the analytical solution) was studied.
Simulations were performed with combinations of the first or the second order Euler
transient schemes and the high resolution or the compressive advection schemes. The
results were compared with theory, showing that the most accurate ones were obtained
with the combination of the second-order Euler transient scheme with the high resolution
advection scheme.
Conclusions
97
As CFX-5.6, for multiphase flows, does not model the viscous and pressure works
neither the kinetic energy effects, the energy equation solved for the resin is simply the
thermal energy equation. Although in the most recent versions of CFX it is possible to use
the full energy equation for multiphase flows, because the viscous work is neglectable for
typical RIM situations and because the kinetic energy effects are only important for high
speed flows, the thermal energy equation is perfectly valid to model the energy balance in
RIM.
Two case studies, experimentally conducted by [39], were used to validate the
whole model in CFX, consisting of the mass, momentum, volume fraction, energy and cure
equations. Both moulds are simple rectangularly shaped cavities and filled with a RIM type
polyurethane. In both cases, as the mould width is considerably larger than the thickness,
two-dimensional simulations were performed, and due to symmetry only half of the
geometry was modelled.
Both cases were run on a mesh with 5 elements on the transverse direction, and one
of them also on a mesh with 10 elements on the transverse direction. The inhomogeneous
model was employed, together with the walls no-slip boundary condition for the resin
equations and the free-slip condition for the air equations. All the governing equations kept
being solved slightly longer (~ 2 %) than the filling time, and after that the fluids were
assumed to be completely at rest and only the energy and the cure equations were solved.
Because of the very small time step required to obtain convergence, the computation times
were considerable. In both cases high degrees of cure were achieved at the end of the
filling stage, resulting in a significant raise of viscosity and consequently the increase of
the inlet pressure. The highest temperature observed during the whole process was, in both
cases, above the temperature resulting from a complete adiabatic cure reaction, meaning
that, although the walls temperature are significantly lower than those highest temperature
values, the part of the resin that reached that temperature had to have o positive energy
balance, along its path, with its surroundings. This clearly shows that the adiabatic
temperature rise cannot be used directly as an upper limit for the maximum resin
temperature during the RIM process, and attests the complexity of the simultaneous flow
motion, cure reaction and heat transfer processes. Some results obtained from the
simulations, such as the evolution of the inlet pressure and of temperatures with time, and
the degree of cure and the temperature distributions at the end of the filling stage, were
compared with available numerical and experimental results from [39] and numerical
results from [72], showing a general good agreement with those.
Two other cases, where all the three dimensions effects are important and where the
flow front splits and remerges, were simulated. At the midplane the flow front was always
well captured, but at the wall its definition becomes less accurate when its velocity is
higher. The times, obtained from the simulations, for the resin to have completely
Conclusions
98
solidified and to have reached the degree of cure such that the parts may be removed from
the mould, showed only a small difference relatively to the theoretical analysis presented.
It was demonstrated that CFX is capable to simulate the complexities of the RIM
simultaneous three dimensional filling and curing processes with good accuracy. However,
even for simple geometries and meshes with very few number of elements (in the context
of what is nowadays common in CFD), due to the very small time steps necessary to obtain
convergence, the computation times are considerably long. More complex geometries,
because of details, will require finer meshes and very probably also smaller time steps,
which will drive the computation times to some impractical values. Even taking cases 3
and 4 as a reference, which took around 5 days of CPU time to complete 8 s of simulation,
if the parts were slightly different and did not have the two planes of symmetry and the
whole geometry had to be modelled, the computation time, making a rough extrapolation,
would rise to 20 days. But, in order to increase accuracy, if one wanted to double, from 5
to 10, the number of elements on the thickness direction (in the 2½D approach the
thickness direction is typically divided in 8-20 layers [14]), keeping the same aspect ratio
of the mesh elements, the computation time would jump to 160 days (20 days
u 2
3
).
In all the simulated cases, the moulds and the 2D geometries were filled upwards
and they had a full outlet at the top, placed at least 4 mesh elements ahead of the final
position of the flow front. In RIM, because of the resins low viscosities, the moulds are
commonly filled upwards to facilitate the exit of the air and to avoid the formation of air
pockets (in TIM, because of the polymers high viscosity, the gravity effect, and therefore
the mould orientation, are irrelevant), but surely real moulds do not have full outlets.
Instead, they contain vents, which shall be ideally located so that the totality of the air can
leave the mould. But, because of their tiny dimensions, vents cannot be meshed. In fact
they could, but that would require the mould to be meshed with elements with the same
order of size as the vents, what would mean a completely unrealistic computation effort.
However, and leaving a suggestion for future work, the parts of the mould walls containing
vents may possibly be modelled as regions resistant to the advance of the resin but free to
the passage of the air, that is, as subdomains with resistive source terms on the resin
momentum equations.
99
References
[1] D. V. Rosato, D. V. Rosato and M. G. Rosato, Injection Molding Handbook
3
rd
Ed., Kluwer Academic Publishers, Boston, 2000.
[2]
O.
Mal,
Modelling and Numerical Simulation of Reaction Injection Moulding
Processes, PhD Thesis, Univertsité catholique de Louvain, Louvain-la-Neuve,
March 1999.
[3] J. J. van der Werf and A. H. M. Boshouwers, INJECT-3, A Simulation Code for the
Filling Stage of the Injection Moulding Process of Thermoplastics, PhD Thesis,
Technische Universiteit Eindhoven, Eindhoven, May 1988.
[4] J. Vlachopoulos and D. Strutt, Polymer processing, Materials Science Technology,
Vol. 19, No. 9, 2003, p. 1161-1169.
[5] C. W. Macosko, RIM - Fundamentals of Reaction Injection Molding, Hanser
Publishers, Munich, 1989.
[6] Bayer Polymers LLC. (2003 http://www.rimmolding.com)
[7] C. D. Snyder, Materials for Reaction Injection Molding (RIM) Processing, 2001
Annual Conference of the Composites Fabricators Association (COMPOSITES
2001), Tampa, Florida, USA, October 2001.
[8] Kawin International, Inc. (2003 http://www.kawincomposites.com)
[9] F. T. Wickis Jr., Reaction injection molding, Medical Device Diagnostic Industry
Magazine, April 1996.
(2003 http://www.devicelink.com/mddi/archive/96/04/014.html)
[10]
I. Barros, Modelação do Comportamento Térmico de Moldes de Injecção,
PhD Thesis, Universidade do Minho, Guimarães, June 2004.
[11] R. S. Spencer and G. D. Gilmore, Some flow phenomenon in the injection moulding
of polystyrene, Journal of Colloid Science, Vol. 6, No. 2, 1951, p. 118-132.
References
100
[12] R. L. Ballman, T. Shusman and H. L. Toot, Injection molding: flow of a molten
polymer into a cold cavity, Industrial and Engineering Chemistry, Vol. 51, No. 7,
1959, p. 847-850.
[13] F. J. Rielly and W. F. Price, Plastic flow in injection molds, Society of Plastic
Engineers Journal, 1961, p. 1097-1101.
[14] S.-W. Kim and L.-S. Turng, Developments of three-dimensional computer-aided
engineering simulations for injection moulding, Modelling and Simulations in
Materials Science and Engineering, Vol. 12, No. 3, 2004, p. S151-S173.
[15] D. H. Harry and R. G. Parrot, Numerical simulation of injection mold filling,
Polymer Engineering and Science, Vol. 10, No. 4, 1970, p. 209-214.
[16] M. R. Kamal and R. Kening, The injection molding of thermoplastics. Part I:
Theoretical model, Polymer Engineering and Science, Vol. 12, No. 4, 1972,
p. 294-301.
[17] M. R. Kamal and R. Kening, The injection molding of thermoplastics. Part II:
Experimental test of the model, Polymer Engineering and Science, Vol. 12, No. 4,
1972, p. 302-308.
[18] J. L Berger and C. G. Gogos, A numerical simulation of the cavity filling process
with PVC in injection molding, Polymer Engineering and Science, Vol. 13, No. 2,
1973, p. 102-112.
[19] P.-C. Wu, C. F. Huang, C. G. Gogos, Simulation of the mold filling process,
Polymer Engineering and Science, Vol. 14, No. 3, 1974, p. 223-230.
[20] C. Gutfinger, E. Broyer and Z. Tadmor, Melt solidification in polymer processing,
Polymer Engineering and Science, Vol. 15, No. 7, 1975, p. 515-524.
[21] G. Williams and H. A. Lord, Mold-filling studies for the injection molding of
thermoplastic materials. Part I: The flow of plastic materials in hot and cold walled
circular chanels, Polymer Engineering and Science, Vol. 15, No. 8, 1975,
p. 553-568.
References
101
[22] H. A. Lord and G. Williams, Mold-filling studies for the injection molding of
thermoplastic materials. Part II: The transient flow of plastic materials in the
cavities of injection-molding dies, Polymer Engineering and Science, Vol. 15,
No. 8, 1975, p. 569-582.
[23] J. R. A. Pearson, Mechanics of Polymer Processing, Elsevier, Amsterdam, 1985.
[24] S. Richardson, Hele-Shaw flow with a free boundary produced by the injection of
fluid into a narrow channel, Journal of Fluid Mechanics, Vol. 56, No. 4, 1972,
p. 609-618.
[25] Y. Kuo and M. R. Kamal, The fluid mechanics and heat transfer of injection mold
filling of thermoplastic materials, AIChE Journal, Vol. 22, No. 4, 1976, p. 661-669.
[26] M. E. Ryan and T.-S. Chung, Conformal mapping analysis of injection mold filling,
Polymer Engineering and Science, Vol. 20, No. 9, 1980, p. 642-651.
[27] S. M. Richardson, H. J. Pearson and J. R. A. Pearson, Simulation of injection
moulding, Journal of Fluid Mechanics, Vol. 56, No. 4, 1980,
p. 609-618.
[28] S. M. Richardson, Simplified geometry models, in Fundamentals of Computer
Modeling for Polymer Processing, Hanser Publishers, Munich, Chapter 4, 1989.
[29] C. A. Hieber and S. F. Shen, A finite element / finite difference simulation of the
injection molding filling process, Journal of Non-Newtonian Fluid Mechanics,
Vol. 7, No. 1, 1980, p. 1-32.
[30] A. C. B. Bogaerds, Stability Analysis of Viscoelastic Flows, PhD Thesis, Technische
Universiteit Eindhoven, Eindhoven, December 2002.
[31] W. Rose, Fluid-fluid interfaces in steady motion, Nature, Vol. 191, No. 4785, 1961,
p. 242-243.
[32] L. R. Schmidt, A special mold and tracer technique for studying shear and extension
flows in a mold cavity during injection molding, Polymer Engineering and Science,
Vol. 14, No. 11, 1974, p. 797-800.
References
102
[33] C. G. Gogos, C.-F. Huang and L. R. Schmidt, The process of cavity filling including
the fountain flow in injection molding, Polymer Engineering and Science, Vol. 26,
No. 20, 1986, p. 1457-1466.
[34] D. J. Coyle, J. W. Blake and C. W. Macosko, The kinematics of fountain flow in
mold-filling, AIChE Journal, Vol. 33, No. 7, 1987, p. 1168-1177.
[35] E. Vos, H. E. H. Meijer and G. W. Peters, Multilayer injection molding,
International Polymer Processing, Vol. 6, No. 1, 1991, p. 42-50.
[36]
W. Zoetelief, Multi-component Injection Moulding, PhD Thesis, Technische
Universiteit Eindhoven, Eindhoven, April 1995.
[37] S. Bhattacharji and P. Savic, Real and apparent non-Newtonian behaviour in viscous
pipe flow of suspensions driven by a fluid piston, Proceedings of the 1965 Heat
Transfer and Fluid Mechanics Institute, Stanford University Press, 1965,
p. 248-262.
[38] Z. Tadmor, Molecular orientation in injection molding, Journal of Applied Polymer
Science, Vol. 18, No. 6, 1974, p. 1753-1772.
[39] J. M. Castro and C. W. Macosko, Studies of mold filling and curing in the reaction
injection molding process, AIChE Journal, Vol. 28, No. 2, 1982, p. 250-260.
[40] M. R. Kamal, E. Chu, P. G. Lafleur and M. E. Ryan, Computer simulation of
injection mold filling for viscoelastic melts with fountain flow, Polymer
Engineering and Science, Vol. 26, No. 3, 1986, p. 190-196.
[41] H. Mavridis, A. N. Hrymak and J. Vlachopoulos, Finite element simulation of
fountain flow in injection molding, Polymer Engineering and Science, Vol. 26,
No. 7, 1986, p. 449-454.
[42] R. A. Behrens, M. J. Crochet, C. D. Denson and A. B. Metzner, Transient free
surface flows: Motion of a fluid advancing in a tube, AIChE Journal, Vol. 33,
No. 7, 1987, p. 1178-1186.
[43] M. R. Kamal, S. K. Goyal and E. Chu, Simulation of injection mold filling of
viscoelastic polymer with fountain flow, AIChE Journal, Vol. 34, No. 1, 1988,
p. 94-106.
References
103
[44] H. Mavridis, A. N. Hrymak and J. Vlachopoulos, The effect of fountain flow on
molecular orientation in injection molding, Journal of Rheology, Vol. 32, No. 6,
1988, p. 639-663.
[45] H. Mavridis, A. N. Hrymak and J. Vlachopoulos, Transient free surface flows in
injection mold filling, AIChE Journal, Vol. 34, No. 3, 1988, p. 403-410.
[46] K. K. Wang and V. W. Wang, Computer aided mold design and manufacturing, in
Injection and Compression Molding Fundamentals, Marcel Dekker, Inc., New
York, Chapter 9, 1987, p. 607-669.
[47] H. P. Wang and H. S. Lee, Numerical techniques for free and moving boundary
problems, in Fundamentals of Computer Modeling for Polymer Processing, Hanser
Publishers, Munich, Chapter 8, 1989, p. 340-401.
[48]
P.
Kennedy,
Flow Analysis of Injection Molds, Hanser Publishers, Munich, 1995.
[49] S. I. Güçeri, Finite difference solution of field problem, in Fundamentals of
Computer Modeling for Polymer Processing, Hanser Publishers, Munich,
Chapter 5, 1989, p. 141-236.
[50] S. Subbiah, D. L. Trafford and S. I. Güçeri, Non-isothermal flow of polymers into
two-dimensional, thin cavity molds: A numerical grid generation approach,
International Journal of Heat and Mass Transfer, Vol. 32, No. 3, 1989, p. 415-434.
[51] A. Couniot and M. J. Crochet, Finite elements for the numerical simulation of
injection molding, Proceedings of The 2
nd
International Conference on Numerical
Methods in Industrial Forming Processes (Numiform'86), Balkema Publishers,
Rotterdam, 1986, p. 165-170.
[52]
A.
Couniot,
Développement d'un Algorithme de Remaillage Automatique pour la
Simulation du Moulage par Injection, PhD Thesis, Université catholique de
Louvain, Louvain-la-Neuve, 1991.
[53] F. Dupret, A.Couniot, O. Mal, L. Vanderschuren and O. Verhoyen, Modeling and
simulation of injection molding, in Advances in the Flow and Rheology of
Non-Newtonian Fluids, Elsevier, Amsterdam, Chapter 7, 1999, p. 939-1010.
References
104
[54] J. F. Hétu, D. M. Gao, A. Garcia-Rejon and G. Salloum, 3D finite element method
for the simulation of the filling stage in injection molding, Polymer Engineering
and Science, Vol. 38, No. 2, 1998, p.223-236.
[55] E. Pichelin and T. Coupez, Finite element solution of the 3D mold filling problem
for viscous incompressible fluid, Computer Methods in Applied Mechanics and
Engineering, Vol. 163, No. 1, 1998, p. 359-371.
[56] G. A. A. V. Haagh and F. N. van de Vosse, Simulation of three-dimensional
polymer mould filling process using a pseudo-concentration method, International
Journal for Numerical Methods in Fluids, Vol. 28, No. 9, 1998, p. 1355-1369.
[57] G. A. A. V. Haagh, G. W. M. Peters, F. N. van de Vosse and H. E. H. Meijer, A 3-D
finite element model for gas-assisted injection molding: Simulations and
experiments, Polymer Engineering and Science, Vol. 41, No. 3, 2001, p. 449-465.
[58] E. Broyer and C. W. Macosko, Heat transfer and curing in polymer reaction
molding, AIChE Journal, Vol. 22, No. 2, 1976, p. 268-276.
[59] L. J. Lee and C. W. Macosko, Heat transfer in polymer reaction injection molding,
International Journal of Heat and Mass Transfer, Vol. 23, No. 11, 1980,
p. 1479-1492.
[60] Z. Tadmor, C. G. Gogos, Principles of Polymer Processing, John Wiley
Sons, Ltd., New York, 1979.
[61] J. D. Domine and C. G. Gogos, Simulation of reactive injection molding, Polymer
Engineering and Science, Vol. 20, No. 13, 1987, p. 847-858.
[62] C. N. Lekakou and S. M. Richardson, Simulation of reacting flow during filling in
reaction injection moulding (RIM), Polymer Engineering and Science, Vol. 26,
No. 18, 1986, p. 1264-1275.
[63] I. Manas-Zloczower, J. W. Blake and C. W. Macosko, Space-time distribution in
filling a mold, Polymer Engineering and Science, Vol. 27, No. 16, 1987,
p. 1229-1235.
References
105
[64] M. R. Kamal and M. E Ryan, Thermoset Injection Molding, in Injection and
Compression Molding Fundamentals, Marcel Dekker, Inc., New York, Chapter 9,
1987, p. 329-376.
[65] C. D. Lack and C. A. Silebi, Numerical simulation of reactive injection molding in a
radial flow geometry, Polymer Engineering and Science, Vol. 28, No. 7, 1988,
p. 434-443.
[66] D. S. Kim and S. C. Kim, Rubber modified resin. III: Reaction injection molding
process, Polymer Engineering and Science, Vol. 35, No. 7, 1995, p. 564-576.
[67] M. A. Garcia, C. W. Macosko, S. Subbiah and S. I. Güçeri, Modelling of reactive
filling in complex cavities, International Polymer Processing, Vol. 6, No. 1, 1991,
p. 73-82.
[68] R. E. Hayes, H. H. Dannelongue and P. A. Tanguy, Numerical simulation of mold
filling in reaction injection molding, Polymer Engineering and Science, Vol. 31,
No. 11, 1991, p. 842-848.
[69] R.-Y. Chang and C.-H. Chen, Simulation of filling and curing processes in epoxy
reactive molding, Journal of the Chinese Institute of Chemical Engineers, Vol. 26,
No. 3, 1995, p. 183-194.
[70] G. A. A. V. Haagh, G. W. M. Peters and H. E. H Meijer, Reaction injection
molding: Analyzing the filling stage of a complex product with a highly viscous
thermoset, Polymer Engineering and Science, Vol. 26, No. 20, 1996, p. 2579-2588.
[71] D. Seo, J. R. Youn and C. L Tucker, Numerical simulation of mold filling in foam
reaction injection molding, International Journal for Numerical Methods in Fluids,
Vol. 42, No
.
10, 2003, p. 1105-1134.
[72] Y.-W. Lo, D. D. Reible, J. R. Collier and C.-H. Chen, Three-dimensional modeling
of reaction injection molding. I, Polymer Engineering and Science, Vol. 34, No. 18,
1994, p. 1393-1400.
[73] Y.-W. Lo, D. D. Reible, J. R. Collier and C.-H. Chen, Three-dimensional modeling
of reaction injection molding. II: Application, Polymer Engineering and Science,
Vol. 34, No. 18, 1994, p. 1401-1405.
References
106
[74] R.-Y. Chang, C.-C. Hung and W.-H. Yang, Three-dimensional CAE Analysis of
underfill flow of flip-chips, ANTEC 2002 Proceedings, May 2002, p. 588-592.
[75] R. Sekula, P. Saj, T. Nowak, K. Kaczmarek, K. Forsman, A. Rautiainen and
J. Grindling, 3-D modeling reactive molding processes: From tool development to
industrial application, Advances in Polymer Technology, Vol. 22, No. 1, 2003,
p. 42-55.
[76] D. Seo and J. R. Youn, Numerical analysis on reaction injection molding of
polyurethane foam by using a finite volume method, Polymer, Vol. 46, No. 17,
2005, p. 6482-6493.
[77] Y. K. Shen, S. H. Chen and H. C. Lee, Analysis of the mold filling process on flip
chip package, Journal of Reinforced Plastics and Composites, Vol. 23, No. 4, 2004,
p. 407-428.
[78] R.-Y. Chang, W.-H. Yang, S.-J. Hwang and F. Su, Three-dimensional modeling of
mold filling in microelectronics encapsulation process, IEEE Transactions on
Components and Packaging Technologies, Vol. 27, No. 1, 2004, p. 200-209.
[79] C. v. L. Kietzmann, J. P. Van Der Walt and Y. S. Morsi, A free-front tracking
algorithm for a control-volume-based Hele-Shaw method, International Journal for
Numerical Methods in Engineering, Vol. 41, No. 2, 1998, p. 253-269.
[80] J. M. Castro, M. Cabrera Ríos and C. A. Mount-Campbell, Modelling and
simulation in reactive polymer processing, Modelling and Simulations in Materials
Science and Engineering, Vol. 12, No. 3, 2004, p. S121-S149.
[81] A. H. Haidari and B. Matthews, Future trends for computational fluid dynamics in
the process industry, Third International Conference on CFD in the Minerals and
Process Industries, Melbourne, Australia, December 2003, p. 331-338.
[82] Wikipedia contributors, Moore's Law, Wikipedia, The Free Encyclopedia,
October 7
th
2006.
(http://en.wikipedia.org/w/index.php?title=Moore%27s_Lawoldid=80010669)
[83]
ANSYS,
Inc.,
CFX-5 Version 5.6 Documentation, 2003.
References
107
[84]
O.
Ubbink,
Numerical Prediction of Two Fluid Systems with Sharp Interfaces,
PhD Thesis, Dep. of Mech. Engineering, Imperial College of Science, Technology
and Medicine, London, January 1997.
[85] P. J. Zwart, M. Scheuerer and M. Bogner, Free surface flow modelling of an
impinging jet,
ASTAR International Workshop on Advanced Numerical Methods
for Multidimensional Simulation of Two-Phase Flow, Garching, Germany,
September 2003, p. 1-12.
[86] ANSYS Canada Ltd, CFX-5 Supplementary Theory, 2004.
[87] AEA Technology Engineering Software, Ltd., CFX-TASCFlow Computational Fluid
Dynamics Software Version 2.12 Theory Documentation, 2002.
[88] T. J. Barth and D. C. Jesperson, The design and application of upwind schemes on
unstructured meshes, AIAA Paper 89-0366, 1989.
[89]
F.
M.
White,
Fluid Mechanics 2
nd
Ed., McGraw-Hill International Editions,
Singapore, 1988.
[90] T. Franck, Advances in computational fluid dynamics (CFD) of 3-dimensional
gas-liquid multiphase flows, NAFEMS Seminar: Simulation of Complex Flows
(CFD) Applications and Trends, Wiesbaden, Germany, April 2005, p. 1-18.
[91] A. Burns, A. Splawski, S. Lo and C. Guetari, Application of coupled solver
technology to CFD modelling of multiphase flows with CFX, First International
Conference on Computational Methods in Multiphase Flow, Orlando, U.S.A,
March 2001.
[92]
ANSYS,
Inc.,
CFX-5 Version 5.7 Documentation, 2004.
[93]
ANSYS,
Inc.,
ANSYS CFX Release 10.0 Documentation, 2005.
[94] M. R. Kamal, S. Sourour and M. Ryan, Integrated thermo-rheological analysis of the
cure of thermosets, SPE Technical Paper, 19, 1973, p. 187-191.
[95]
S. V. Patankar, Numerical Heat Transfer and Fluid Flow, Taylor Francis,
Philadelphia, 1980.
References
108
[96] F. R. Menter, P. F. Galpin, T. Esch, M. Kuntz and C. Berner, CFD simulations of
aerodynamic flows with a pressure-based method, 24
th
International Congress of
the Aeronautical Sciences (ICAS 2004), Yokohama, Japan, August-September
2004.
[97] S. R. Estevez and J. M. Castro, Applications of a reaction injection molding process
in the analysis of premature gelling, demold time and maximum temperature rise,
Polymer Engineering and Science, Vol. 24, No. 6, 1984, p. 428-434.
[98] J. M. Castro and C. W. Macosko, Reaction injection molding, in Encyclopedia of
Materials Science and Engineering, Vol. 6, Pergamon Press, Oxford, 1986,
p. 4085-4090.
[99] R.-Y. Chang and W.-H. Yang, Numerical simulation of mold filling in injection
molding using a three-dimension finite volume approach, International Journal for
Numerical Methods in Fluids, Vol. 37, No. 2, 2001, p. 125-148.
Frequently asked questions
What is the main topic of this document?
The document focuses on the numerical simulation of the filling and curing stages in Reaction Injection Moulding (RIM), using CFX (now ANSYS CFX) software. It includes analyses of the filling process using both homogeneous and inhomogeneous models, simulation of the curing process, and a combined simulation of filling and curing.
What is Reaction Injection Moulding (RIM)?
RIM is a process for rapid production of complex parts through the mixing and chemical reaction of two or more components, typically isocyanate and polyol. The liquid components are injected into a mixing chamber and then flow into a mold cavity where they polymerize and cure.
What are the key differences between Thermoplastics Injection Moulding (TIM) and Reaction Injection Moulding (RIM)?
TIM involves injecting hot polymer melt into a cold cavity where it solidifies by cooling. RIM involves injecting liquid components into a mold where they undergo chemical reaction (curing) to form a solid part. TIM uses thermoplastics, which can be re-melted, while RIM uses thermosets, which undergo an irreversible chemical reaction.
What is the homogeneous model used for in the document?
The homogeneous model assumes that both phases (resin and air) in the filling process share a single velocity field. It's typically used for modeling free surface flows where phases are completely stratified.
What are the limitations of using the homogeneous model for RIM simulation, according to the document?
The document finds that the homogeneous model, with the usual no-slip boundary condition, is unable to model the filling process correctly, leading to an incorrect interface that doesn't touch the walls.
What is the inhomogeneous model and how does it address the limitations of the homogeneous model?
The inhomogeneous model uses separate flow fields for each phase (resin and air) and accounts for interphase momentum transfer. It overcomes the limitations of the homogeneous model by implementing a free-slip boundary condition for the air phase at the walls, allowing for a more accurate representation of the filling process.
How is the cure reaction implemented in the numerical simulations?
The cure reaction is implemented as a transport equation for an additional scalar variable representing the degree of cure, with a source term representing the rate of cure. Various transient and advection schemes are tested for accuracy.
What are the key findings related to the energy equation in the context of RIM simulation?
The viscous dissipation term can often be neglected compared to the heat generated by the cure reaction. A high temperature gradient can be observed on the mould walls when the flow front arrives.
What is the significance of the Courant number (Cr) in the simulations?
The Courant number is related to numerical stability. The advection scheme chosen for the volume fraction equation automatically allows high Courant numbers locally with boundedness principles to keep the interface sharp.
What are the key software and numerical methods used in this work?
The simulations are performed using the commercial CFD software package ANSYS CFX. Numerical methods include the finite volume method (FVM), different transient and advection schemes, and various two-phase flow models.
What are key themes explored in this document based on the table of contents?
Key themes explore the filling process in the RIM, including using the homogenous and inhomogenous model. There is a key focus on the curing process and simulation of the RIM filling and curing stages.
Excerpt out of 120 pages
- scroll top
- Quote paper
- Rui Igreja (Author), 2007, Numerical Simulation of the Filling and Curing Stages in Reaction Injection Moulding, Munich, GRIN Verlag, https://www.grin.com/document/273166
Look inside the ebook
