7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Paper No
MUL TIPLE WA YCOUPLING OFAIRFIBERDROPLET PHASES IN THE POL YURETHANE
FIBERREINFORCED COMPOSITE SPRA Y MOLDING MANUFA CTURING PROCESS
P. Diffo*, P. Wulf* and M. Breuer
Department of Mechanical and Production Engineering, University of Applied Sciences Hamburg,
Berliner Tor 21, D20099 Hamburg, Germany
diffoi~rzbt.hawhambure. de, wulfi~rzbt. hawhamburg de
tDepartment of Fluid Mechanics, Institute of Mechanics, HelmutSchmidtUniversity Hamburg, Holstenhofweg 85, D22043
Hamburg, Germany
breueri~hsuhh.de
Keywords: Multiphase Flow, Fiber Dynamics, FiberDroplet Interaction, Spray Modeling
Abstract
To investigate and compute the fiber orientation and density distribution in polyurethanebased (PUR) composites manufac
tured by spray molding, a new approach is presented in this paper. This manufacturing technology is applied by spraying the
initially liquid polyurethane matrix together with laterally injected reinforcing fibers in a tool form or on a substrate. It is
presumed that the final position and orientation of a fiber on a substrate result from its dynamics and coupled interactions
with air, PURdroplets and other fibers within the spray cone. Therefore, a model of the spray molding process is built, that
computes the spray cone within the environment of a finite volume based CFD code (ANSYS Fluent) using the
EulerLagrange approach with an airdroplet2waycoupling, and a dropletdropletinteraction. Furthermore, the fibers are
modeled with a code (FIDYST), based on the dynamical KirchhoffLovetheory, that considers a fiber as a onedimensional,
inextensible elastic rod. Its dynamics is governed by the gravitational and the surrounding fluid forces. The effects of air flow
turbulence are computed by adding a stochastic linear fraction to the nonlinear deterministic aerodynamic forces which are
derived from the fluidfiberrelative velocity (airfibercoupling).
For the dropletfibercoupling an algorithm is used, that averages the timedependent dynamic data (mass and velocity) of
discrete PURdroplets present in a finitevolumecell over the cell volume to compute the effective composite fluid forces
acting on the fibers in that cell. The effective forces are considered in FIDYST in a straightforward manner by adjusting some
specific fluid flow field quantities already used for the interaction of the fibers with the continuous air phase. In this paper
two approaches of the mentioned adjustment are presented introducing the expression "homogenized fluid" (Hfluid) to
characterize the numerical / mathematical result of the airparticlecellspecific mixture.
The coupling of the two codes happens sequentially by first simulating the spray process without fibers and writing periodi
cally some specific node based fluid data, and then starting a fiber simulation using the fluid data. The final position of each
fiber's discretized node is written in a separate file, which is evaluated with a MZlATLABbased algorithm. This MlATLAB rou
tine projects the 3D fiber's node position on the 2D substrate before the mean fiber orientation and density distribution of
each quadrilateral cell of the preliminary partitioned substrate can be quantified and visualized.
1. Introduction
In polyurethanefiberreinforced composite spray molding
polyurethane (PUR) is sprayed in a form or on a substrate
while chopped short or long fibers (natural, glass) are later
ally injected in the spray cone for wetting before the
PURfibermixture hits the substrate and cures to form a
composite (see Figure la).
This manufacturing process is at present mainly used in the
sanitaryware industry for example to reinforce the rear sides
of bathtubs and shower trays. Further applications can be
found in the agricultural machinery and utility vehicle sector
and in the automotive industry. Especially the latter sector
requests low priced and high flexible technologies which
allow manufacturing automobile parts of nearly arbitrary
shape, resistance and stiffness. The PUR spraying technol
ogy may be able to fulfil some of these requests. First
parts like roof modules, fenders or cowls were produced as
composites by this technology.
For a wider utilisation and acceptance of the production
technology it is necessary to obtain information about the
average orientation and density distribution of the fibers in
the composite to predict at least some mechanical properties
of the parts. In order to compute these distributions and to
quantify the influence of some process parameters on them,
the different interaction processes of the three phases (air,
glass fibers and polyurethane droplets) within the spray jet
are to be modeled and simulated [1]. The longterm goal is
to identify the process parameters which influence the av
erage density and orientation distributions of the fibers in
the composites. This information would make the process
controllable and the distributions predictable.
The aim of this paper is to present the mathematical ap
proach and methodology for modeling and simulation of
this process. The fluid phases, both PUR and air, are mod
eled with the CFD code ANSYS Fluent. The presence of a
solid phase (fibers) in the multiphase problem requires the
use of a second software, FIDYST, developed by the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Paper No
FraunhoferInstitute for Industrial M2lathematics ITWM~l
(Kaiserslautern, Germany).
Together with the modeling of the three phases, their
5waycoupling (see Figure lb) over two software packages
will be covered in this paper with a special focus on the
fiberdroplet coupling.
E
O
t~
E
p
Y
u
~
E
me
ro
rp
a,
nE
~
I~iAWll~naA
Op~
5;~
/5
Time axis
Figure la: Longfiber injection spray process with the air
cone, the sprayed droplets and the chopped fibers.
L= gas phase
T= droplet phase
$ F = fiber phase
L, F = phase interaction
Figure lb: Interaction's model between the three phases of
the spraylongfibberinj section process.
2. Simulation Methodology
The modeling of the longfibersinjection spray process
presents a high degree of complexity. The ejected liquid
(phase 1) is first atomized by injected air (phase 2), then the
airdropletmixture encounters the added fibers (phase 3).
Due to the presence of a solid phase and the amount of solid
elements in this multiphase problem it is advisable to carry
out the modeling in segregated form. Thus two simulations
are performed sequentially. First a transient, RANS
EulerLagrange based, twoway coupled CFD simulation is
realized, followed by a fiber dynamics simulation. The latter
needs some specific data of the gas phase that are used to
compute drag and resistance forces on each fiber. In detail,
the density p, the velocity vector u, the turbulent kinetic
energy k and the turbulent dissipation rate E have to be
transferred into the second simulation,
In Figure 2 the detailed workflow of a complete simulation
shows GAMZBIT, a CADenvironment and mesher, at the
beginning of the complete simulation. This allows a
3Dreconstruction and meshing of the flow domain to be
computed.
Figure 2: Workflow of a complete process simulation.
Within the ANSYS Fluent environment the two fluids phases
are modeled using the EulerLagrangeapproach for the gas
phase and the particle phase, respectively. Both phases are
fully coupled (twowaycoupling) and require the definition
of the spray parameters [6]. Here also the particleparticle
interaction is taken into account by making use of the parcel
model. The spray injector is modelled as a point source and
can be moved through the computational domain with help
of User Defined Functions (UDFs), thus avoiding the use of
remeshing methods. The final step in the CFD computation
is to export the appropriate data for the fiber dynamics
simulation.
The at first glance apparent approach to couple the gas and
liquid phases separately to the fibers would require for the
liquid phase some extensional algorithms for the drop
letfiber collision's probability and the drop
letfiberoneway momentum transfer (see Figure 3). The
computational cost of such an approach with collision de
tection and eventdriven momentum transfer would be
prohibitive.
Thus, the authors present for the dropletfibercoupling a
new approach where the timedependent dynamic data
(mass, density and velocity) of the gas phase and liquid
droplets present in a finitevolumecell are averaged over
the cell volume. A fictive homogenized fluid is mathemati
cally created which cell specific data (p, u, k, e) at specified
time steps are homogenized properties of PURdroplets and
air (see Figure 4). This method avoids the direct calculation
of the probability of local collision events, which obviously
have to some extent a random nature, and the modeling of
the explicit momentum transfer due to a collision of a drop
let with a fiber (see Figure 4). The flow field quantities (p,
u, k, e), that are already used for the interaction of the fibers
with the continuous air phase, are adjusted to allow in addi
tion a dropletfiber interaction. The mathematical method is
called "homogenization" and the homogenized fluid is
denoted "Hfluid" by the authors. The mathematical ap
proach is presented further below in this paper with two
variations.
Paper No
Figure 3: Apparent approach for the dropletfiber coupling
FiFe ren po iin
Fiber orientation
Figure 4: Present approach for the dropletfiber coupling
combined with the airfiber coupling.
Continuing on the time axis of the workflow in Figure 2, a
fiber dynamics simulation is conducted in the next step by
making use of the stored data of the Hfluid. The properties
of the Hfluid are imported into a new domain meshed by
3Delements as well as 2DElements for visualization and
fiberwallcontact (substrate) computation. Fibers getting in
contact with a wall of the new domain are stored in a data
file for a posterior distribution analysis. This step is not
depicted in the Figure 2.
3. Governing equations and solution method
3.1.Gas Phase
The presence of air as a continuous fluid in a spray process
leads to an accelerated liquid core disintegration. This
airassisted atomization produces in general smaller droplets
as the plain jet (airless) atomization. .
The formulation of the gas phase govemmig equations em
ployed herein is in Eulerian form based on the conservation
equations for mass and momentum using the mncompressible
continuity equation and the NavierStokes equations:
pV. u = S,, (1)
p Su+ pV(uu) = Vp + V.? + pg + F (2)
The source term S, refers to the mass added to the con
tinuous phase by the atomizer or ejector, respectively. u is
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
the velocity vector, p is the gas density, p is the static pres
sure, I is the stress tensor, g is the gravitational accelera
tion and the source term F refers to external forces per unit
volume that arise from the interaction with the disperse
phase. The specific interaction force source term is updated
every iteration of the disperse phase equation and couples
the gas flow with the liquid droplets.
The modeling of statistical turbulence is based on the Rey
noldsaveraged approach using the realizable kemodel
with k as the turbulence kinetic energy and e as the turbu
lence dissipation rate (Shi et al. [15]). Therefore, two addi
tional transport equations for k and e are added to the
precedent system of equations. This model differs from the
standard model of Launder and Spalding [16] in two points:
the turbulent viscosity was reformulated
the dissipation rate transport equation was changed
Notice that the energy equation is not solved in this project.
The transient computation of the gas phase through the
domain occurs with a pressurebased solver where:
 the time is 1st order implicitly discretized
 the pressure is interpolated using the momentum equation
coefficients. This is the Standairl approach in 4NSTS Fluent
 the pressurevelocity coupling is achieved by using the
SIMPLE algorithm
 the spatial discretization of the convective terms of the
different transport equations (momentum and turbulence)
uses the 1st upwind scheme
 the gradients are reconstructed at the cell faces according
to the GreenGaussCell based approach
 the air ejected by the movable atomizer is represented by
two sources using UDFs. The air volume pouring into the
domain is defined by a mass source, whereas a momentum
source accounts for accelerating the air up to the ejector
velocity
a third UDF determines the cell in which the atomizer
point presently resides. This information is needed for the
origination of the liquid droplets.
3.2. Liquid Phase
Spray processes play an important role in many technical
and industrial applications. Crop spraying, spray painting or
spray cooling are typical examples.
The numerical simulation of sprays with help of CFD
methods has been preferentially used in industry for model
ing fuel injection processes in thermal engines [2] (motors,
turbines and rocket machines), metal coating processes [3]
and driving processes in the chemical industry.
The process of disintegration of the continuous fluid is
highly complex, not completely understood and highly
dependent on (controlled by) the nozzle geometry and pres
ence of compressedair [12]. However, the occurring
breakup regimes in the air assisted atomization process can
be divided into primary and secondary atomization and each
can be allocated to a region in the spray cone [13]. The first
type is based on the capillary or chaotic breakup of a jet or
liquid sheet ejected with high velocity from the nozzle due
to high injection pressure, whereas the second type is based
on the aerodynamic breakup caused by shear stresses at the
liquidgas interface.
The modeling of the primary atomization would mean the
do = 1.88d (1f + h)
dL iS the diameter of the liquid ligament prior to the disinte
Paper No
computation of the liquid bulk's transition from a stream to
droplets through growth of disturbance waves in the liquid.
This alone typically requires the numerical method "Volume
of Fluid (VoF)" with a highly resolved mesh of the compu
tational domain for the resolution of possible, small droplets
shedding from the stream. The attempt would require high
computer costs in order to obtain a few thousands of drop
lets (particles).
Therefore, in a lot of the CFD codes the particles's atomiza
tion is not part of the spray simulation. The reason for that is
that in the atomization shed particles are accelerated (enter a
higher relative velocity regime), deformed and divided
(shattered) in smaller, "satellite" droplets by the surrounding
gas phase and by the occurring dropletdropletcollision.
This period is the secondary atomization of the liquid.
The applied secondary droplet breakup model is the Taylor
Analogy Breakup (TAB) model based on the analogy of
Taylor [20] between a distorting, oscillating droplet and a
spring mass system. The equation governing the droplet
distortion y, meaning the deviation of the droplet equator
from its equilibrium position, is based on the differential
equation:
d2y F, p k~ C4 1~ dy
=  y (3)
dt2 CbPP r2 pp3 pp2 dt
p and p, are the continuous phase and discrete phase densi
ties, u, is the relative velocity between the gas phase and the
droplet, r is the undisturbed droplet radius, o is the droplet
surface tension, and pu is the droplet viscosity. The dimen
sionless constants Cb, d F and Ck, have been chosen 1/2,
5, 1/3 and 8, respectively, to match experiments and theory
[19].
Breakup occurs when y>1. Thus size and velocity of the
"child" droplets must satisfy the equation of energy conser
vation according to [21]. The TAB model is well suited for
low Weber number injections (We < 100) [6]. The dimen
sionless droplet Weber number relates the dynamic pressure
to the surface tension and is a characteristic measure of the
droplet breakup behaviour. It is defined as:
(4)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
gration and Oh is the Ohnesorge number, a dimensionless
number which divides the square root of the Weber number
by the Reynolds number.
The shedded droplet diameters are statistically distributed
around the mean diameter do by a stochastic function enti
tled "RosinRammler distribution" [6]. The diameter do is
assumed to be the most probable droplet size of the
twoparameter droplet diameter RosinRammler distribution
function. That implies the existence of an exponential rela
tionship between the droplet diameter, d, and Y,, the mass
fraction of droplets with diameter greater than d:
Yd ("/do n (6)
do is the mean diameter and n is the spread parameter. Ac
cording to the work of Schmidt et al. [17] the applied spread
parameter is 3.5.
For computing the trajectories of the dispersed phase (drop
lets) the forces acting on droplets are integrated according to
the Newton's 2nd law given in a Lagrangian reference frame:
dup p
at FD (H Hp) + g PPP+ F1 F2
(7)
u, and p, are the velocity and the density of a spray droplet
respectively while u and p denote the velocity and the den
sity of a the gas phase. FD(H Hp) TepreSents the drag
force per unit mass and:
18p, CdRe
F ppd: 24 8
p;l is the dynamic fluid viscosity, p, is the density of the
particle, and d, is the particle diameter. Re is the relative
Reynolds number, which is defined as:
pdp u, u
e, (9)
Cd is the drag coefficient computed dynamically by consid
ering the Weber number dependent droplet shape distortion
during the droplet motion through the gas. The drag coeffi
cient varies linearly from the drag coefficient of a sphere
(Cd,sphere) 10 the coefficient of a disc by evaluating the dif
ferential equation of the droplet distortion y in equation (3)
[18].
The equation of drag variation reads:
C, = Cd,sphere(1 + 2.632y) (10)
pdpur2
W
e = 
d, is the particle diameter,
To solve the problem of computing simultaneously millions
of particles' trajectories and possible intercollisions, the
concept of "parcels" is used [31]. Parcels are statistical
representations of a number of individual droplets. The
representative parcel has the same properties as the particles
which it represents and carries on top the number of repre
sented particles as supplementary information. Nevertheless,
the reduction of the number of particle trajectories to be
computed has to be balanced by modeling the effects of the
"missing" particles in the coupling algorithms with the gas
phase and with other parcels (parcelparcelcollisions).
At the begin of the simulation a number of parcels is prede
fined. A parcel mean diameter, do, is computed from the
predefined atomizer parameters (nozzle diameter, liquid
mass flow, spray cone angle):
Re > 1000
Re <; 1000
Cd~sphere =~ 6.21 2 3
and y is the droplet distortion governed by equation (3).
The second term of equation (7) corresponds to the gravity
term. Fl describes the "virtual mass" force which has to be
obeyed as the fluid surrounding the particle is also acceler
ated:
1p d
F, = (u u,) (12)
2 pp dt
( 5)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Paper No
Finally, F2 describes an additional force arising from the
pressure gradient in the fluid generating a force on the parti
cle due to different pressures acting on the particles surface:
F2 UpYH (13)
The tivowaycoupling of air and particles is realized by the
determination of momentum exchange between the gas
phase and the parcels, according to the previously presented
equations. When a parcel passes through each control vol
ume (cell) of the computational mesh the momentum trans
fer is initially computed by evaluating the momentum
change of the parcel in the present time step and adding a
sink term to the gas phase momentum equation (2). The
subsequent solution of equation (2) yields to a new gas
velocity which is then reintroduced in the parcel momentum
equation (5). This operation is repeated until convergence is
reached by both velocities,
The computation of the interaction between particles can be
divided in three steps: The search for collision partners, the
probability of the collision event and the outcome of the
collision. The applied particleparticlecollision is based on
the O'Rourke algorithm [22] which assumes that two par
cels collide only if they reside in the same cell of computa
tional mesh during the present time step. Furthermore, the
probability of the smaller parcel being within the collision
volume of the bigger parcel is evaluated. In this concept the
exact location of the two colliding droplets is not required,
The collision's probability is computed statistically based on
the local distribution of the droplet presence. The probabil
ity distribution of the number of collisions between the
bigger parcel and all the particles of the smaller parcel is a
Poisson distribution function,
Once it is decided that two parcels collide, the O'Rourke
algorithm considers only two collision's regimes, the coa
lescence and the bouncing of parcels, depending on the
parcel's Weber number We (equation (4) with d, as the mean
diameter of the two parcels) and on the impact angle.
As the outcome of the collision the post collision velocities
of either one coalesced droplet or two bounced off droplets
is computed based on the conservation of momentum and
kinetic energy, considering losses to viscous dissipation and
angular momentum generation,
The turbulence dispersion of particles is modeled using the
stochastic tracking model (random walk) which uses sto
chastic methods to include the effects of instantaneous
turbulent velocity fluctuations of the gas phase on the parti
cles. This interaction of a particle with a turbulent eddy
characterized by a Gaussian distributed random velocity
fluctuation V' and the characteristic life time re, prevail
during the eddy life time. The values of V' are predicted by
making use of the quantities from the ke model under the
presumption of isotropy. In addition, the values of V' obey a
Gaussian probability distribution:
V' = ( (14)
where 5 describes a random number normally distributed
and k the turbulence kinetic energy.
The computation of the eddy life time occurs through a
random calculation:
ze = T log(r) ,
(15)
where r is a uniform random number between 0 and I and
TL is the integration time scale which describes the time of
a particle spent in a turbulent eddy along the particle path:
T = 0.15 (16)
and E is the turbulent dissipation rate.
The transient computation of the gasliquid spray mixture in
the domain is realized by the .4NSYS Fhient builtin "air
blast atomizer model". This is in brief an injection point
from which the gas and all the particles or parcels, respec
tively, originate. The injection point can be moved in the
domain during the computation to simulate a robot driven
path of the original PUR spray nozzle. The atomizer model
requires:
 a number of parcels per inj section time step
 a (time dependent) position vector for the injector
 a relative velocity vector of the born parcels
 a spray cone angle
 a liquid mass flow rate
 a nozzle diameter
 a sheet and a ligament constant
 an atomizer dispersion angle.
The equation of motion of the parcels is solved by using an
automated tracking scheme selection which switches in an
automated fashion between a chosen lower order scheme,
here the implicit scheme, and a higher order scheme, here
the trapezoidal scheme. The automatic change depends on
the numerical distance of a parcel from its hydrodynamic
equilibrium. When the distance is small the lower order
scheme is applied and vice versa.
3.3. Fiber Phase and Fiber Dynamics
The modeling of fibers within a fluid field has been sub
jected to intensive research in the past century. In the year
1922 Jeffreys already published a model for the prediction
of the motion of an ellipsoid in a fluid flow [23]. With the
appearance of short fiber reinforced polymers and their
possible FEMcomputing (design) some new semiempirical
models were developed in the 80's for FEM based filling
and curing studies. These models considered all the fibers as
one continuum and not as separate elements within the flow
[2427]. However, these models were limited to short fiber
computations. First numerical simulations of single long
fibers in a single flow were conducted at the beginning of
the 90's with a particle simulation model which built the
fibers as a chain of particles according to the beadmodel
[28]. Lately, the Fraunhofer institute of industrial M2lathe
inatics ITWIMl has developed a simulation tool (FIDYST) to
compute the dynamics of endless fibers for the
meltspinning process of nonwoven materials [3638]. Dur
ing this process hundreds of individual endless fibers are
continuously extruded, stretched and entangled by highly
turbulent air.
To model long but finite fibers in a spray flow an adaption
of the FIDYST software algorithm has been made.The tra
jectory and deformation of a fiber is based on the dynamical
Paper No
KirchhoffLovetheory for the timedependent behavior of a
Cosserat rod (Cosserat continuum). This theory allows large,
geometrically nonlinear deformations ([10], [14]). Under
presumption of the domination of elastic (extensile) effects
the BemnoulliEuler material law is applied while consider
ing a fiber as a onedimensional inextensible elastic rod
parameterized by its arc length s and time t. A single fiber
can be specified mathematically by the momentum conser
vation equation and the length conservation equation:
pAzra, =Ta`) E,I~ 94 a+'Y sray (
PFFat2 aS\ Of S4
~2= 1 (18)
r(s,t) is the three dimensional position of the fiber, T is a
modified normal force, p, is the fiber's density, .4, is the
fiber cross section, E, is the fiber's elasticity modulus, I is
the momentum of inertia and fg" is the external gravita
tional force while "P'" represent the spray mixture forces
acting on the fiber.
Figure 5 depicts qualitatively a fiber with the acting forces.
fgray
gravitational force
/spray dynarnieforce
s=0 n ra~r
Figure 5: Fiber with external forces (adapted from [10]).
The total aerodynamic force on the fiber equals the integral
of the projection of the stress tensor of the incompressible
NavierStokesequation on the vector normal to the curva
ture, vector n in Figure 5, over the fiber's surface. The
force can be modelled by a typical aerodynamical approach.
comprising the drag coefficient, Bemnoulli's dynamic pres
sure and the cross sectional area of the fiber. We refer to [10,
14] for detailed informations about the modeling of the drag
coefficient and the approach.
Similar to the splitting of the velocity in a mean and a fluc
tuating part as a result of the RANS approach, the external
aerodynamic force f"'" is separated into a deterministic
mean force of the spray and a stochastic turbulent force of
the air only. The first results from the mean spray velocity
and the latter from fluctuating gas velocity: ~ ~
farr(r) = fspray 7)+ fair(r)' (19)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
4. Coupling between Particles and Fibers
In the previous chapter the fluid phases and their coupling
as well as the solid fiber phase and its interaction with a
single fluid phase (air) were introduced. In this chapter a
newly developed approach for the coupling of the particles
with fibers is presented. As already mentioned the apparent
approach would be to detect individual collisions between
the parcels and the fibers. This at first glance apparent ap
proach is however difficult to implement as the parcel con
cept makes it impossible to extract the exact position and
momentum for each particle. Furthermore, even if only the
parcel would be accounted, for the problem of collision
detection and momentum transfer remains. To overcome
this problem both phases (gas and droplets) are averaged or
"homogenized" to a new fluid, called "Hfluid" within the
CFD environment. The cell specific data (p,,eiv uet, k, e)
of the new fluid are then used in the FIDYST environment
to compute the dynamics of falling fibers,
Two variants of this homogenization were tested:
averaged momentum homogenization
discrete momentum homogenization
For both variants the following definitions are valid:
Cell specific volume fraction of both phases:
The conservation of volume within a cell is defined as:
Vceu, = Vas = I1NiVpi (20)
where n denotes the number of parcels within a cell, N, the
number of particles comprised in the parcel i, V,,the vol
ume of a particle, VGasthe volume of the gas phase and V,,en
the volume of a cell.
This equation can be rewritten as:
Vas =1lNiVpi
1 = + (2 1)
Vcell Vcell
The volume fraction a of the dispersed phase is defined as:
(2 = NgVp (22)
Vceul
and should not exceed 1 (a < 1).
This definition leads to the equation of the volume fraction
for the gas phase:
~s= 1 a (23)
Density of a cell containing a homogenized fluid:
The conservation of mass within a cell is defined as:
mceul = mcas =1 C 1 iPi (24)
where m,,eu, mas and myi refer to the mass of a cell, the
mass of the gas phase and the particle within the cell, re
spectively. By making use of the different phase densities
this leads to:
PceulVceul = PcasVcas PP =I1 NiVp, (25)
with fspray(r)
f arr(r)'
=f(relative sprayfibervelocity)
f(relative airfibervelocity, k, E)
PcasVas PP I=1 NiVpi
pcecu = +
Vceul Vceul
(26)
Notice that a fiberfiberinteraction is not yet implemented.
PP CI11NiVp upi
p (1 ')uL +V
Vceu
p (1 a2) + ppa
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Paper No
Pceul = Pcas~l ppa
(27)
the mean cell velocity for the Hfluid is obtained:
pcasl PP a~ Pa
uceul = P P
(33)
where poem, Pas and p, is the density of the cell, the gas
phase and the particle phase, respectively.
With the definition of equation (22), (23) and (27) we can
present the aforementioned two model variants. The two
models differ in their approach for the momentum conserva
tion,
4.1 Averaged momentum homogenization
The averaged momentum homogenization averages the
velocities of the particles present in a cell to a mean velocity
(see equation (28)). This velocity is applied to all the parti
cle in the cell in order to compute the total particle momen
tum in the cell (see equation (30)). The velocity of the ho
mogenized fluid is deduced from the summation of the
momentum contributed by the different fluid phases to a cell
momentum (see equations (3133)). Figure 6 depicts the
mean velocity of a cell computed with the velocity of the
gas phase and the mean velocity of the particles.
SVelocity of the gas phase
SMean velocity of the cell
SMean velocity of liquid phase
Figure 6: Velocities within a cell for the averaged momen
tum homogenization model.
The equation of the mean velocity of the particles u, in
each cell of the computational domain is:
up = 22~ (28)
u,. is the velocity of a particle within the parcel i and A;
denotes the number of particles within the parcel and n
denotes the number of parcels within a cell.
The momentum equation of the gas phase Ias within a cell
can be simplified by making use of equation (20) to:
IGas G masu PGasVGasu PGas ( a Vceullu (29)
The momentum equation of the dispersed phase I, with
use of the mean particle velocity is:
IP = 1): zmpi Up = pp(pr V,)UP = pp(a2Vceu,)U, (30)
The cell momentum can be written as the sum of the
phases' momentum contributions in the following form.
Iceul = IGas + P (1
This definition yields to the equation:
PceulVcellceul = PGas ( a VceluN + PP aVceul)UP
(32)
After the insertion of equation (27) for the mean density,
4.2 Discrete momentum homogenization
The discrete momentum homogenization accumulates the
momentum of every particle in a parcel (see equation (34))
and adds it to the gas phase momentum in order to deter
mine the required cell velocity (see equations (3536)).
Figure 7 depicts this idea of the collection of momentum
contributions to a homogenized cell specific momentum.
This allows the derivation of a cell specific mean velocity.
Velociyof fra sohs
 Velocityofparticles
Figure 7: Velocities within a cell for the discrete momen
tum homogenization model.
The momentum equation of the gas phase Ias within a
cell is equal to equation (29).
The momentum equation of the disperse phase I, accumu
lates the momentum contributions of each particle within a
parcel :
IP = ,C1 Ng Vi upi (34)
The cell momentum can be written as the sum of the phases'
momentum contributions corresponding to equation (31):
PceulVcellceul = PL[(1 a)Vceu]HuL + PP =I1 NiVpi upi
(35)
After the insertion of equation (27) the mean cell velocity
for the Hfluid is obtained:
(36)
ucel =
4.3 Coupling Method
The data features obtained through homogenization are
applied to the fibers using a quadratic Taylor drag approach
for high Reynolds numbers. This empirical motivated ap
proach is used within FIDTST for the airfiber coupling [10].
The spray dynamic force fspray of equation (19) is splitted
in 8 HOrmal f spray and a tangential component:
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Paper No
fspray Pceuld, uceusN (s, t)
d,
(37)
fspray 2 Pceuld, 5.4
Y IuceusN (s, t)I ucr II(s, t)
d,
(38)
uceusN and uceur, are the normal and tangential components
of the cell specific Hfluid velocity uceu, computed with
either the averaged or the discrete model.
In addition, a turbulent airfiber coupling is applied, ac
cording to equation (19). In this equation the fluctuating
stochastic component fai'l depend on the global Gaussian
fluctuation of the gas phase velocity field u and not yet of
the spray velocity field uceu,. Detail on the turbulence mod
eling approach can be found in [10].
5. Description of Simulation Case
Two sequential simulations were conducted, a CFD compu
tation of the fluid flow in the domain and a fiber dynamics
simulation to record the influence of the spray on the be
haviour of the fibers.
5.1.Grid
A container (lm x 0,5m x 3m) was modeled and a hexahe
dral grid (100 x 30 x 200 cells) was drafted in GA.4MBIT. The
outcome was a mesh with a cell volume of 510 m3 and a
total number of 6106 cells. The domain interior was initially
filled with air and five of the six boundaries were defined as
pressure outlet. Only the lower boundary in the XZplane
was defined as a noslip wall to represent the substrate with
wall friction for the gas phase and trapping for the discrete
droplets. The trapping condition means that the parcel tra
jectory computation stops when it gets in contact with the
wall. Gravity was activated in negative ydirection. Figure 8
shows the domain with a symbolically drawn spray cone
and moving direction of the injector representing the spray
nozzle.
Figure 8: Symbolically depiction of the domain with the
spray cone and the moving direction of the injector.
5.2. Gas Phase
The modeling of the moving gas phase inlet in the moving
injector was realized by a determination of the cell contain
ing the injector position in the current time step and a sub
sequent injection of air from that cell into the domain. The
transient computation was executed with a pressure based
solver and the following parameters:
 a time step is set 210's in order to meet a Courant number
of 0.1
 an air mass flow of 0.0225kg/s and an air injection veloc
ity of 60m/s in the cell representing the air injector
 an inj ector velocity of 0.6m/s in zdirection.
5.3. Liquid Phase
The PUR particles were injected in the domain by the "air
blast atomizer model" of Fluent. The atomizer coincided
with the injector and performed the same movement. For the
computation of the transient particle traj ectories
 an inj section timestep of 1 10 s
 a liquid density of 1100kg/m3
 a liquid viscosity of 1kg/(m.s)
 a droplet surface tension of 0.4N/m
 a parcel number per inj section timestep of 50
 a parcel relative velocity vector of [0 ; 10 ; 0]m/s
 a spray nozzle diameter of 710 m
 a spray sheet constant of 12
 a spray ligament constant of 0.5
 and an atomizer dispersion angle of 60
were used.
5.4. Solid Fiber Phase
The mesh and data from the CFD simulation were imported
into with the FIDTST model for the fiber dynamic simula
tion. To complete the modeling three geometric elements
were added:
 a substrate for the fiberwall contact
 a fiber inj ector and a spray j et nozzle for the visualization.
Figure 9 depicts the model with the different geometric
elements and the moving direction vector.
Z
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
struction schemes. These tests showed that the value of the
difference depends strongly on the computing scheme. This
value could be minimized but could not yet be eradicated.
Further analysis of this problem will be conducted in order
to increase the accuracy of the global modeling.
2 B:I:;H+0 1
2otuso B Ye+0 ly ng lue lns ~me132P0 1nr1.21
244e+013,dp bss~, ustsy
L5Bue+10: eito ftecl ae otu lto
Contours of Yelocity Magnitud e [m/s) [Time1.342 Be01 ) L E T B3[ d p p ns, s ru t2ndl
Figure 10: Depiction of the coel based contour plot of the
averalgedl hoogncized sr cn fluid scaled between 0 a d maxi
mumices density giue h atcl est sapr
enl ig tth einis. .eijeto efr hepri
2.1e+ 0 ceeae adtu itiutdi h oe
2.0e *0iaino h cn satiue t h oino h
1ozl. 928+0 le ihte ihs elct r lctdi
1.Be+0tsra ftegspaeweetegsvlct n
1.Bs te+0 mnu rnfrt heprilsi ih h
1.5Bice+0aeesvr bten431~4 n .10
1.47he *0i.W ee gi othsfgr gi ntenx
: :grph
1.3 +
1.0 e+
I :":~ .1
zeoct magiu
Paper No
Movingdirection
Liquid spray nozzle Fiberjet nozzle
a~~ J usrt
Figure 9: Elements of the fiber dynamics computation.
The laterally injected glass fibers are characterized by
 a length of 1.510 'm
 a diameter of 1. 7810" m
 a density of 2103 kg/m3
 an elasticity module of 500N/m2
 an inj section velocity of 3m/s
 an injection angle of 600 with respect to the spray axis
 a moving velocity of 0.6m/s of the fiber inj ector
The transient fiber dynamics computation were conducted
using a time step of 110's.
6. Results and Discussion
As the two approaches were tested the authors distinguish
between the data features common to both of the homoge
nized models and the data features discerning them from
another
Common data features concern the resulting air jet, the
resulting spray and the homogenized density. Assigned
results are presented in chapter 6.1.
The comparative data features are the homogenized veloci
ties of each approach as well as the fiber jet deflection.
These are discussed in chapter 6.2.
6.1 Common Data Analysis
The gas phase:
The data used in FIDYST to compute the dynamic drag
force on the fibers must be node based data. The CFD com
putation stores its data in a cell based manner and interpo
lates them on the cell nodes for visualisation or further
analysis. In order to check the accuracy of the interpolation
two alternative representations are shown and compared,
Figure (10) and (11) show the contours of the air injected
into the domain as a cell based and node based depiction,
respectively. In the cell based depiction the maximum value
in the injection cell is approximately twice as high as the
node based value. However, the node based values in the
actual spray cone possess the same velocity level as in the
cell based depiction. As the fibers are injected into the spray
cone in some distance to the injection cell the observed
differences have not effects on the fiber dynamics.
Another anomaly in the results is the high difference be
tween the predefined air injection velocity of 60m/s in the
injection cell and the plotted cell maximum velocity of
44.4m/s in the cell based velocity depiction of Figure (10).
Some parametrical simulations have been conducted com
bining different discretization schemes in time and space
together with different pressurevelocity and gradient recon
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
air jet itself in Figure (1 1).
The maximum values between the discrete and the averaged
Hfluid do not differ significantly from another with
8.61m/s and 9.89m/s, respectively. But the discrete Hfluid
shows a more wide spread distribution. It can be concluded
that the domain of influence with almost the same velocity
is larger than in the averaged Hfluid.
Paper No
The homogenized density:
Figure (13) shows the homogenized density distribution in
the domain with a density peak of 286kg/m3 at the injection
cell. This can be explained by the high concentration of
particles in that cell at every injection time step (see Figure
(12)). The injected fluid density of 1100kg/m3 is opposed to
the air density of 1.225kg/m3 within a cell,
Because of the iniquitous depiction a rescale of the data is
necessary and another distribution is presented in Figure
(14). It shows that a density between 0 and 15kg/m3 prevails
in the zone of interest which is well distributed in the spray
11gure 10: Depiction or the velocity contour
discrete homogenized fluid.
11gure 14r: Depiction of thle density contour plot or thle
averaged homogenized fluid.
regure lrj: vepiction or mne velocity contour plot or mne
averaged homogenized fluid.
The fiber let deflection:
A video analysis of the spray molding manufacturing proc
ess has highlighted the deflection of fiber flow when it
encounters a spray cone. This jet deflection is used in this
paper as a first validation criterion for the developed ap
proaches.
Figure (17) and (18) depict individual fiber trajectories
through the spray represented by the discrete and the aver
aged Hfluid, respectively. The spray is here illustrated by
the velocity field of the Hfluid. In both figures the fiber
flow is deflected when it enters the spray cone which a good
quality agreement with the real processing.
rlgure le: vepiculon or me aensnry comrour plor or me
averaged homogenized fluid scaled between 0 and
20[kg/m ].
6.2 Comparative Data Analysis
The velocity distribution of the Hfluid:
Figure (15) and (16) show the velocity contour plot of the
node values of the averaged and discrete Hfluid, respec
tively, after a couple of time steps.
The Hfluid velocity is in both model inverse proportional to
its density (see equation (33) and (36)). Thus the areas with
high densities muss have low velocities and the other way
around. This explains the pearshaped high velocity regions
of both approaches in Figure (15) and (16), respectively,
where those areas coincide with the low density area of
Figure (13). This differs from the tubeshaped maximum
velocity region of the air jet itself with a maximum velocity
value at the injection (see Figure (ll)).
Both approaches show a drop in the maximum velocity of
over 50% compared to the maximum node velocity of the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Paper No
fiber cross section
dimensionless constant
drag coefficient
drag coefficient of a sphere
dimensionless constant
dimensionless constant
dimensionless constant
particle, parcel mean diameter (m)
particle, parcel diameter (m)
liquid ligament diameter (m)
fiber's elasticity modulus (Nm 2)
spray dynamic force (N)
meall Spray dynamic force (N)
normal projection of spray dyn. force (N)
18Hgenti81 prOjeCliOn of spray dyn. force (N)
fluctuating aerodynamic force (N)
gravitational force (N)
external force vector
specific virtual mass force (ms 2)
specific pressure gradient force (ms )
drag factor (s )
gravitational acceleration vector (ms 2)
momentum of inertia (m )
cell momentum vector (kgms )
gaS momentum vector (kgms )
dispersed phase momentum vector (kgms)
turbulent kinetic energy (m s )
cell mass (kg)
gas phase mass (kg)
mass of a single particle i (kg)
spread parameter, number of parcels in a cell
vector normal to fiber curvature
number of particles comprised in the parcel
Ohnesorge number
static pressure (Nm )
Reynolds number
fiber position vector (m)
undisturbed droplet radius (m)
random number
fiber arc length (m)
mass source term
time (s)
modified normal force ()
integration time scale (s)
gas velocity vector (ms )
gas velocity fluctuation (ms )
particle velocity vector (ms )
mean particle velocity vector in a cell (ms ')
airdroplet relative velocity (ms )
mean cell velocity (ms )
normal proj section of mean velocity (ms )
tangential proj section of mean velocity (ms ')
particle mean velocity vector (ms )
velocity vector of a particle within a parcel i
cell volume (m )
gas phase volume (m )
volume of a single particle i (m )
droplet distortion (m)
A,
Cb
Ca
Cd~sphere
Cj
C,
do
d,
E,
spray
fspray
spray
N
(spray
T
fair'
fgare"
F
Fz
F2
FD
g
I
Ice,,
IGas
Ip
k
mceen
m
mp
n
N i
Oh
p
Re
r
r
r*
s
t
T
To
u
I,
u,
ucert
uceur
up,
V i
Vp,
11gure 17: Depiction or the 11ber tiow entermg the spray
sillltratP 1 hTw two va~ietyT of two TlierTtP us~l~inia
r figure 15: vepiction of the ~10er nlow entermg me spray
illustrated by the velocity of the averaged Hfluid.
7. Conclusion and further work
The paper describes two efficient approach variations to
homogenize two fluid phases, a gas phase and a droplet
phase, in order to model the polyurethane fiberreinforced
composite spray moulding manufacturing process. This
homogenization avoids the computation of the probability of
individual dropletfiber collision,
The fiber flow is deflected when it enters the region of the
Hfluid which is the first visual validation criteria for any
approach. A further criterion is the evaluation of the fiber
orientation and fiber density on the substrate which should
be in good accordance with experimental data. This analysis
is part of the future work in this proj ect.
However, some collision related events like the mass growth
of a wetted fiber due to the collision with a droplet or the
fiberfiberinteractions are not implemented in the model at
the present time. In addition, a more accurate momentum
transfer which considers the particle density distribution in a
cell of the computational mesh is potentially more accurate
because this would allow an elaborated allocation of mo
mentum transfer between the spray and the fibers.
A coupling of the drag force stochastic term with the ho
mogenized fluid is also planned.
Nomenclature
Paper No
We Weber number
Greek letters
a volume fraction of the dispersed phase
turbulent dissipation rate ()
p, pas gaS density (kgm')
Pcert cell density (kgm ')
pF density of the fiber material (kgm )
pp liquid density (kgm')
Stress tensor
Sdroplet viscosity (I\i n!i
vas gaS phase viscosity ((l\is n, !!
o droplet surface tension (Nm )
f random number
e 7 eddy life time (s)
Subsripts
Sindice
Acknowledgements
The authors would like to thank the German Federal Mimis
try of Education and Research (BMBF) for supporting fi
nancially this project and also the Hennecke GmbH for the
technical support and insight information.
References
[1] Diffo, P., Wulf, P.: Simulation des PUR Langfaser
spriihrerfahrens. in: "PUR 2009", VDIBerichte /
VDITagungsbande Kunststofftechnik, Band Nr. 4301,
Ntimrberg, 2009.
[2] Stiesch, G.: Modeling Engine Spray and Combustion
Processes (Heat and Mass Transfer). Springer, Berlin, 2003.
[3] Fritsching, U.: Spray Simulation Modeling and Nu
merical Simulation of Sprayforming metals. Cambridge
University Press, Cambridge, 2004.
[4] Crowe, C.T., Sommerfeld, M., Tsuji, Y : Multiphase
flows with Droplets and Particles. CRC Press, New York,
1998.
[5] ANSYS CFXSolver Theory and ANSYS CFXSolver
Modeling. Manuals of ANSYS CFX 10.0, ANSYS Inc.,
Waterloo, Ontario, 2005.
[6] ANSYS FLUENT User Guide. Manual of FLUENT 6.3,
FLUENT Inc., Lebanon, 2006.
[7] Reitz, R.D.: Mechanisms of Atomization Processes in
HighPressure Vaporizing Sprays. Atomisation Spray Tech
nology, 3, 1987.
[8] O'Rourke. P.J.: Collective Drop Effects on Vaporizing
Liquid Sprays. PhD thesis, Princeton University, Princeton,
New Jersev, 1981.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
[9] Hietel, D.: FYDIST Fiber Dynamics Simulation Tool,
Innovation durch Simulationsunterstuitzung. Vortrag, 23.
Hofer Vliesstofftage, Hof, 2005.
[10] Marheineke, N.: Turbulent Fibers On the Motion of
Long, Flexible Fibers in Turbulent Flows, PhD thesis, Tech
nische Universitiit Kaiserslautemn, 2005.
[ll] Ferziger, J.H., Peric, M.: Computational Methods for
Fluid Dynamics. Springer, Berlin, 2008.
[12] Lefebvre, A.H.: Atomization and sprays. Hemisphere
NY, 1987.
[13] Lavemnia, J.E, Wu, Yu: Spray Atomization and Deposi
tion. John Wiley & Sons, NewYork, 1996.
[14] Panda, S.: The Dynamics of Viscous Fibers, PhD thesis,
Technische Universitit Kaiserslautern, 2006.
[15] Shih, T.H., Liou, W W., Shabbir, A., Yang, Z., & Zhu,
J.: A New kE EddyViscosity Model for High Reynolds
Number Turbulent Flows Model Development and Valida
tion. ComputersFluids, 24(3):227238, 1995.
[16] Launderand, B.E., Spalding, D.B.: Lectures in Mathe
matical Models of Turbulence. Academic Press, London,
England, 1972.
[17] Schmidt, D.P., Corradini, M.L., Rutland, C.J.: A
TwoDimensional ,NonEquilibrium Model of Flashing
Nozzle Flow. In 3rd ASME/JSME Joint Fluids Engineering
Conference,1999
[18] Liu, A.B., Mather, D., & Reitz, R.D.: Modeling the
Effects of Drop Drag and Breakup on Fuel Sprays. SAE
Technical Paper 930072, SAE, 1993.
[19] Lamb, H.: Hydrodynamics, Sixth Edition. Dover Pub
lications, NewYork, 1945.
[20] Taylor, G.I.: The Shape and Acceleration of a Drop in a
High Speed Air Stream. Technical report, In the Scientific
Papers of G.I. Taylor, ed., G.K. Batchelor, 1963.
[21] O'Rourke, P.J. & Amsden, A.A.: The TAB Method for
Numerical Calculation of Spray DropletBreakup. SAE
Technical Paper 872089, SAE, 1987.
[22] O'Rourke, P.J.: Collective Drop Effects on Vaporizing
Liquid Sprays. PhD thesis, Princeton University, Princeton,
NewJersev, 1981.
[23] Jeffrey, G.B.: The Motion of Ellipsoidal Particles Im
mersed in a Viscous Fluid. Proc. R. Soc. A 102, 161179,
1922.
[24] Folgar F., Tucker III, C.L.: Orientation Behaviour of
Fibers in Concentrated Suspensions. Journal of Reinforced
Plastic Composite 3, S. 98119, 1984.
[25] Advani, S.G., Tucker III, C.L.: A Tensor Description of
Fiber Orientation in Short Fiber Composites. Proceeding
Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
43rd Annual Technical Conference, SPE 1985, S. 11131118,
1985.
[26] Advani, S.G., Tucker III, C.L.: The Use of Tensors to
Describe and Predict Fiber Orientation in Short Fiber Com
posites. Journal of Rheology 31, S. 751784, 1987.
[27] Advani, S.G., Tucker III, C.L.: Closure Approximations
for Three Dimensional Structure Tensors. Journal of Rheol
ogy 34, S. 367386, 1990.
[28] Yamamoto, S.: Matsuoka, T.: Dynamic simulation of
fiber suspensions in shear flow. J. Chem. Phys., Vol. 102, 5,
February 1, 1995.
[29] Rave et al.: Simulation of the fiber spinning process.
ITWMBericht Nr. 26, Kaiserslautern, 2001.
[30] ITWM Jahresbericht 2004, S. 22, Kaiserslautern 2005.
[31] Paschedag, A.R.: CFD in der Verfahrenstechnik, S.79.
WileyVCH, Weinheim, 2004.
