CFD Simulation of Particles in
Pipe Flow and Mixing Tank
Aljaz Janic
Link¨oping University IEI - Department of Management and Engineering
Link¨oping University IEI - Department of Management and Engineering Division of Applied Thermodynamics and Fluid Mechanics Master Thesis 2020|LIU-IEI-TEK-A–20/03862-SE
CFD Simulation of Particles in
Pipe Flow and Mixing Tank
Aljaz Janic
Academic supervisor: Marcus Jansson Industrial supervisor: Fredrik Innings Industrial co-supervisor: Dragana Arlov Examiner: Matts Karlsson
Link¨oping University SE-581 83 Link¨oping, Sweden
Abstract
This project aimed to investigate the capability of the STAR CCM+ software when predicting the flow with particles using Lagrangian Particle Tracking and Discrete Element Method. The first part pertained to rectangular channel flow, with ratio between height of the channel and particle diameter (2hD−1p ) of 15. It was found out that simulations of particles in a channel come with many difficulties. Such as, obtaining accurate pressure drop results using DEM when comparing to DNS simulations including particles within a reasonable computational time.
The second part consisted of a simulation of the off-centred mixing tank. As the use of DEM caused numerical issues, another modeling approach was used. Therefore, the Lagrangian Particle Tracking was used. The outcome of the project is the sen-sitivity study of the forces which can be applied to the particles. The finding was that the Shear Lift force and the Virtual Mass force have a negligible contribution in regards to the particles distribution. In addition to this, it was also discovered that the turbulence model has a large effect on the particles in the near-wall region. Choosing an isotropic turbulence model resulted in clustering of the particles near the wall, therefore an anisotorpic turbulence model needed to be used.
Acknowledgements
During the master thesis project, I have encountered numerous issues and prob-lems which would not be resolved without certain people. Firstly I have to thank my industrial supervisors at TetraPak Fredrik Innings and Dragana Arlov for all the help and guidance given during this project. I would also like to express my gratitude to Alfred Jongsma for his expertise combined with valuable experience in simulating particles. I would also like to thank Mani Johannesson from Siemens, for all the support given with the software Star CCM+. I would also like to thank my academic supervisor from Link¨oping University Marcus Jansson for guidance during this project.
Finally, I would like to thank my family and my girlfriend for all the support given during my master’s degree.
Aljaz Janic June 2020, Lund
Nomenclature
Abbreviations and Acronyms
Abbreviation Meaning
CFD Computational Fluid Dynamics
CAD Computer Aided Design
CPU Central Processing Unit
CFL Courant-Friedrichs-Lewy
DEM Discrete Element Method
DNS Direct Numerical Simulation
FF Fanning Friction Factor
KTH KTH Royal Institute of Technology
LiU Link¨oping University
LPT Lagrangian Particle Tracking
LES Large Eddy Simulation
MRF Moving Reference Frame
RANS Reynolds Averaged Navier-Stokes equations
RKE Realizable k −
RBM Rigid Body Motion
RSM Reynolds Stress Model
SKE Standard k −
VOF Volume of Fluid
Latin Symbols
Symbol Description Units
A Area m2 Cµ, CT, C1, C2 Constant [−] cp Heat capacity W m−2K−1 Cd Drag Coefficient [−] Cl Lift Coefficient [−]
Cf s Static Friction Coefficient [−]
D Reynolds Stress Diffusion [P a]
Dp Diameter of the particle [m]
d1 Length of the particle [cm]
d2 Width of the particle [cm]
E Young’s modulus [P a] Fb Body Forces [N ] Fc Contact Force [N ] FCo Coulomb Force [N ] Fd Drag Force [N ] Fg Gravity Force [N ] Fp Pressure Force [N ]
Symbol Description Units
FM RF Moving Reference Frame Force [N ]
Fs Surface Forces [N ]
Fu User Defined Force [N ]
Fvm Virtual Mass Force [N ]
G Shear modulus [P a]
G Buoyancy production [P a]
h Height of the channel [m]
i Internal Energy [J ]
I Identity Tensor [−]
k Turbulent Kinetic Energy m2s−2
Np Power number [−]
Nq Flow number [−]
N up Nusselt number for particle [−]
NRe Reynolds number mixing tanks [−]
n Number of rotations [RP M ]
P Turbulent Production Tensor [P a]
p Pressure [P a]
P r Prandtl number [−]
R Reynolds Stress Tensor [P a]
Re Reynolds number [−]
Sx Source term in momentum equation [−]
SE Source term in energy equation [−]
t Time [s]
Te Large Eddy Timescale [−]
u+ Dimensionless Velocity [−]
u Velocity ms−1
Vp Volume of the particle m3
y+ Dimensionless Wall Distance [−]
y Wall Distance [m]
Greek Symbols
Symbol Description Units
δij Kronecker delta [−]
Dissipation Rate of TKE m2s−3
λ Thermal conductivity W m−2K−1
µ Dynamic viscosity kgm−1s−1
ν Kinematic viscosity P as−1
νp Poisson’s Ratio [P a]
Πij Pressure-velocity gradient tensor [P a]
ρ Density of the fluid kgm−3
ρc Density of the continuous phase kgm−3
ρd Density of the dispersed phase kgm−3
Symbol Description Units
σ, σk Constant [−]
τij Viscous Stress N m−2
τw Shear stress on the wall N m−2
τ1, τ2, τ3 Particle time scale [s]
Contents
1 Introduction 1
1.1 Background . . . 1
1.2 Purpose and aims of the study . . . 2
1.3 Limitations . . . 2
1.4 Previous Work . . . 3
1.4.1 Particles in the pipe and channel flow . . . 3
1.4.2 Particles in mixing tank . . . 4
2 Theory 7 2.1 Fluid Dynamics . . . 7
2.1.1 Pipe Flow . . . 7
2.1.2 Stirring tank . . . 9
2.2 Computational Fluid Dynamics . . . 11
2.2.1 Governing Equations . . . 11
2.2.2 Modeling Turbulence . . . 12
2.2.3 Near wall flow . . . 17
2.2.4 Rotating Motion . . . 18 2.3 Multiphase flow . . . 19 2.3.1 Particle Forces . . . 19 3 Method 27 3.1 Channel flow . . . 27 3.1.1 Geometry . . . 27 3.1.2 Solver Setup . . . 28 3.1.3 Mesh . . . 30 3.2 Mixing tank . . . 32 3.2.1 Geometry . . . 32 3.2.2 Mesh . . . 33 3.2.3 Solver Setup . . . 34 4 Results 37 4.1 Channel flow . . . 37 4.2 Mixing tank . . . 41 5 Discussion 55 5.1 Channel flow with particles . . . 55
5.2 Mixing tank . . . 58
6 Conclusion 61
7 Future Work 63
A Particle Parameters 67
A.1 Particle Size . . . 67
A.1.1 Potato . . . 67
A.1.2 Pepper . . . 68
1
Introduction
1.1
Background
Particles can be found in different industrial applications, from the pharmaceutical industry, the food industry to the automotive industry. Particles pose problems to engineers when trying to predict their behavior or even simulate them. This was the reason for this project, which was done in collaboration with TetraPak. TetraPak is specialized in food processing and packaging services. For a company like TetraPak, performing simulations is really important to reduce the need for costly experiments. In the food industry, foods with particles are becoming increasingly popular, which is where the need for these simulations comes in.
In order to design an efficient processing system that meets customer demand, one needs to know what is happening inside and understand the physics behind it. Prod-uct development goes through several different stages, and one of them is modeling and simulations of the process. Performing experiments is one way of trying to un-derstand the physics, but Computational Fluid Dynamics is also a tool which can help to explain and visualize what is happening inside the machine, pumps, mixing tanks or pipes. In the end, understanding the implications that come with particles when processing food will distinguish you from others and make you stand out from the rest.
As mentioned above CFD is a tool that helped in the past and will also help engineers in the future to get a better understating of industrial fluid flows. It is constantly evolving and adapting to industrial needs, which makes it even more attractive to use. Furthermore, the latest developments in turbulence models make it more ap-plicable to a wider variety of cases. It will never completely eliminate the need for costly experiments, but it can complement them nicely. In addition to this, CFD gives complete results, while the outcome of the experiment is sometimes limited. Furthermore, when equipment is being scaled up, it is almost impossible to perform experiments on such a large scale, which is another reason why CFD simulations are needed.
This thesis focuses on two parts, the first part covers simulations of particles within a pipe, using program software STAR CCM+. Since TetraPak is the company dealing with all kinds of fluids, with or without the particles, they want to better under-stand how the aforementioned software handles simulations with particles. This is a simple geometry that allows to test different configurations and settings to deter-mine how sensitive the program is to a change of mesh size and applied forces. The parameter which will be investigated is a pressure drop when adding particles in the domain, and that will be compared with results from DNS simulations performed by Ardekani et al. at KTH university [1].
The second part of the thesis will focus on mixing tanks, which represent a large part of the TetraPak portfolio. As mentioned above, to design an efficient mixing process one needs to understand what is happening inside a machine. The food in the mixing tank can contain particles, and homogeneous distribution is desired, which is the reason why a CFD simulation is conducted in this study.
1.2
Purpose and aims of the study
Purpose and aim of this project are:
• Simulate particles in channel flow at moderate Reynolds number using a method-ology developed for higher Reynolds number.
• Develop a CFD methodology when having particles in an off-centred agitated tank.
1.3
Limitations
The limitations of this study are written for easier understanding of what could go wrong during the thesis.
• Limited computational power (TetraPak cluster).
• Limited time, from middle of January to beginning of June. • Lack of experimental data.
• Limitations of the software STAR CCM+.
If any of these problems are encountered, support from a university or company will be needed. That said, in general, one must be aware that limitations exist and have a plan of what to do to solve and overcome them.
1.4
Previous Work
Previous work will be split into two parts, the first section will focus on a literature study and previous work of the particles in pipes and channels and the second part will be about mixing tank.
1.4.1
Particles in the pipe and channel flow
Simulating a fluid with particles and how they affect the pressure drop with the rela-tion of particle size and concentrarela-tion, was already investigated by some researchers. Kuerten in 2011 [2] performed a simulation of turbulent channel flow coupled with particles. The result of the simulation was an enhancement of the heat transfer and a small increase in the friction velocity. Sagar et al. [3] performed experiments in a channel flow, the focus was on pressure drop and velocities in the channel. The ex-periments were performed for three different volume fractions (5, 10, and 20%) and also three different particle sizes (2H/dp=40, 16, and 9). The experiments they
com-pleted coincides with the results from DNS simulation in a channel from Uuhlman et al. [4]. The most recent report available was performed by Ardekani et al. [1] at KTH, in which they investigated heat transfer enhancement and pressure drop in a periodic channel flow using DNS simulations for several different Reynolds numbers, concentrations and particle sizes. The results obtained with this study from DNS simulations were used as a validation for channel flow.
TetraPak, which had done work on simulating particles using DEM in both pipes and ducts that served as a background study for the master thesis project. For the flow with particles, they performed simulations for two different domain shapes, namely a circular pipe and a rectangular channel. It was found that simulating fluid in the domain with particles using DEM is always a compromise between particle size and mesh size, as can be seen in Figure 1.
Figure 1: Mesh size related with particle size
Figure 1 shows that if particle spans over several mesh cells this might lead to un-stable solutions. However, mesh resolution should also be fine enough to resolve near-wall gradients and flow features in the middle of the pipe.
Furthermore, the simulations were only performed for a small range of Reynolds numbers from 10000 to 40000. The main parameter of interest was the pressure drop. The pressure drop was also validated with experimental data from different
researches, for the pipe flow with Martin [5] and the channel flow with Sagar et al. [3]. It showed that in low Reynolds numbers (Re = 10000) discrepancies between results kept increasing from 5.6% without the particles to 9.2% for 5% volume frac-tion and 18.9% for 10% volume fracfrac-tion of particles in the domain. On the other hand, high Reynolds numbers (Re = 40000) simulated with STAR CCM+ predict less than 1% difference compared to the experiments in different particle volume fractions.
To conclude the literature study of simulations with particles inside pipes and chan-nels, it can be said that these geometries were thoroughly investigated in the past. The results of these studies can be used for validation and comparison.
1.4.2
Particles in mixing tank
The second task was concerned with bigger mixing tanks and and the accuracy with which simulations can predict general behavior inside the tank.
Simulating a mixing tank with coupled CFD and DEM approach has been done be-fore by Shao et al. [6], who used Fluent coupled with EDEM to simulate a Rushton blade mixing tank. Simulations were compared with experiments, with the em-phasis on velocities, particle distribution, and rotational velocity of the particles. The concentration which they were simulating was 1% volume fraction and small particles. The conclusion was that Fluent coupled with EDEM predicts accurate values of volume fraction around the domain compared with experiments. Blais et al. [7][8] investigated particles using DEM in an axial stirring tank. Simulations were performed with small particles with an axial agitator at different rotational speeds. Simulation results were compared with experimental results both visually and by measuring pressure at the bottom wall. It was concluded that CFD coupled with DEM gives strong results when comparing it visually.
In addition to the previous article, Blais et al. [9] investigated the properties of the particles and how they affect the results. Properties of the particles which were looked into are stiffness and friction factors of the particles at different rotational speeds of the agitator. It was concluded that the stiffness does not play an important role and that results remain the same no matter how stiff the particle is. However, tangential and rolling friction factors affected the volume fraction of the suspended particles at a lower rotational speed. It was also found out, when the agitator is spinning with higher rotational velocity these parameters do not have an immense effect on the final solution. They will only increase or decrease the simulation time needed to achieve a steady-state solution.
In addition to previous studies which were done using mostly Reynolds Averaged Navier Stokes (RANS) methods, Eng et al. [10] did a LES study on how solid parti-cles affect macro instability frequency in a stirred tank. The partiparti-cles were modeled and resolved with the DEM approach. An axial agitator was used with varying volume concentration (0.4%-3%) and particle size (1 mm and 2 mm) in the mixing tank. It was concluded that macro instabilities also known as low-frequency periodic
flow instability can be detected and that the CFD approach can accurately predict the phenomena. These instabilities influence the flow patterns inside a mixing vessel. Simulating particles in a mixing tank can be done as it was shown in several articles, but it is not as widespread of a case as channel flows.
2
Theory
2.1
Fluid Dynamics
Fluid flow is defined with three main governing equations, which are: Conservation of Mass, Conservation of Momentum and Conservation of Energy [11]. These equations are applied to a control volume, which can be either fixed (Euler approach) or moving with the flow (Lagrange approach). The main difference is that the Eulerian approach is observing the flow from a fixed position. The Lagrangian method is following the control volume, which means that the position of a control volume will change with time.
Figure 2: Euler and Lagrange control volume
Figure 2 shows both control volumes and the way they function. The majority of fluid flow simulations utilize the Eulerian approach, while there are some cases, where using the Lagrangian approach is unavoidable.
There are several important numbers in Fluid Dynamics, one of them is called Reynolds number. Reynolds number is an indicator, whether the fluid is in a turbu-lent or laminar regime. It is the ratio between inertial and viscous forces, as shown in the Eq. 1.
Re = uLρ
ν (1)
Where u is velocity, L is characteristic length (diameter), ρ is density and ν is kinematic viscosity. According to Schlichting [12], the transition from the laminar region in pipe flow occurs at Re = 2300 and the flow is fully turbulent when Re is higher than 2900. In between those numbers, there is a transition region, where flow transforms from laminar to turbulent.
2.1.1
Pipe Flow
To better understand what is happening with the fluid when it is moving, one needs to step back and look at the flow in general. Flows in general are mostly flowing over
or inside a geometry, which causes that fluid interacts with a no-slip wall. When this happens, the flow develops a boundary layer on that surface, which can be also seen in Fig. 3.
Figure 3: Boundary layer [13]
When freestream flow interacts with the wall, which is rough i.e. no slip, the bound-ary layer starts to develop. The thickness of the boundbound-ary layer grows with the distance that fluid travels. This phenomenon is caused by friction between layers of the fluid. The layers of the fluid close to the wall will be affected the most by the no-slip wall. In addition to this, the wall will also experience shear stress caused by the liquid. This is the reason for development of the velocity profile. [14]
The velocity profile is fundamental for a basic understanding of how the fluid be-haves and it can be found in all kinds of the engineering applications. For example, the channel flow is a geometry where the velocity profile plays an important role, due to the fact that the ratio between streamwise and spanwise dimension is very large. On the other hand, the mixing tank has some other distinct flow patterns which are in some cases more important than velocity profile. To sum up, one should know what results to expect from the simulations or experiments and which are the most dominant features that need to be resolved or measured.
To help with validation, there are empirical equations which were derived and devel-oped for a specific case based on experiments. These relations are non-dimensionalized, which makes them applicable to a wide variety of cases. Since the channel flow is a common engineering problem, some of the equations could be used to compare simulation results.
Fanning Friction Factor
Fanning friction factor is one of these empirical relations, which is used to evaluate the pressure drop in pipes and channels. It is calculated with the equation seen in the Eq. 2.
f = τ
ρu22 (2)
Where τ stands for wall shear stress, ρ is the density of the fluid and u is the velocity of the fluid. The expression is defined with streamwise wall shear stress divided by dynamic pressure. It is a non-dimensional value, thus it used to different geometries and cases. Figure 4 shows the Moody chart which relates the Reynolds number with
Fanning Friction factor. Moody chart was constructed with the help of experiments, its purpose it that it can be applied to a wide variety of flow conditions.
Figure 4: Moody diagram [15]
It can be seen that the Fanning friction factor changes with different Reynolds number. The Fanning friction factor is decreasing with higher Reynolds number, this holds for laminar and turbulent region. In between there is a transition region where the Fanning friction factor increases. Moody chart is used to determine the pressure drop with relation to Reynolds number and pipe roughness.
2.1.2
Stirring tank
Stirring tanks are widely used in processing industries for mixing in various ways. There are several agitators design which will have an effect on the performance of the agitation done in the tank.
Figure 5: Different design of the agitators [16]
There are two main directions where the agitator will push the flow, either in the radial direction or in the axial direction, as seen in Figure 5. Axial agitator pushes
the flow in downwards direction which will result in one large circulation region. On the other hand, there is a radial agitator that will push the fluid towards the side wall. Two smaller circulation regions will form in tanks with radial agitators. These tanks in the Figure 5 also have baffles, which will decrease splashing in the tank, but increase pressure losses. In addition to this, baffles also increase mixing efficiency. When baffles are not used, then in low-viscous fluid tanks the agitator is positioned off-center to avoid solid-body-rotation. However, for high-viscous fluid tanks without baffles, the agitator can be positioned in the center.
Power Number
Agitators come with characteristics curves, which are related to Reynolds number for mixing tanks. Equations for Power number (Np) and Reynolds number (NRe) in
mixing tank are shown in Eqs. 3 and 4. Np= P n3 d5 ρ (3) NRe= n d2 ρ µ (4)
Where P is power from the agitator, n is rotations per second, d is the diameter of the agitator, ρ is the density of the fluid and µ is dynamic viscosity. Power number is a dimensionless number that specifies the ratio between resistance force and inertia force. Since this number is dimensionless it is used to compare different agitators and it is usually plotted against Reynolds number.
Flow Number
Flow number is a non-dimensional value of mass flow discharge from the impeller divided by rotations and diameter of the blade, Eq. 5.
Nq=
Q
n d3 (5)
Where Q stands for volume flow discharged from the agitator, n is the number of rotations per second and d is the diameter of the agitator. Using volume flow and the diameter of the tank helps to determine whether the particle will sink or float prior to performing the experiment.
2.2
Computational Fluid Dynamics
Computational fluid dynamics is a tool, which was used in the past and it will also be used in the future. With the help of CFD, some flows and phenomena are easier to understand and use in everyday life. However, simulating cases are getting more and more complex which is why CFD as a tool will be put to the test in the future. Continuous development will ensure that the tool adapts to industrial needs.
2.2.1
Governing Equations
The Continuity equation (Eq. 6) states that the amount of incompressible fluid that comes in the control volume is the same as the amount of the fluid that goes out.
∂ui
∂xi
= 0 (6)
Conservation of Linear momentum (Eq. 7) is derived by applying Newton’s second law to a fluid element (Control volume) in all three directions [11].
∂ρui ∂t + ∂ρ ∂xi (uiuj) = − ∂p ∂xi + ∂τij ∂xj + Sx (7)
The left-hand side of the Momentum Equation (Eq. 7) are forces that come from the shear stress of the fluid. While on the right-hand side are the forces acting on a body, i.e. pressure forces, gravity, and source of external forces. There is an unknown viscous stress term τij, which for Newtonian fluids is proportional to the
time rate of strain. Viscous stresses are caused by molecular diffusion across the boundary around the control volume.
τxx= 2µ ∂u ∂x + λ divu τxy = µ ∂u ∂y + ∂v ∂x (8)
For a Newtonian fluid, they are determined with the Eq. 8. Their formulation is important for the near-wall velocity profile, due to the friction between layers of the flow in the boundary layer. Some of the simulations also require energy equation, because they can not be considered as isothermal, which is the reason to introduce the energy equation (Eq. 9).
ρDE
Dt = −div(pui) + div(k grad T ) + ∂ ∂xj
(uiτij) + SE (9)
The left side of this equation is the change of energy on the fluid particle within time. The right side is the sum of the work done by a fluid particle together with the net rate of the heat added to the fluid and a source term [11]. The specific energy of the fluid is calculated using the Eq. 10.
E = i +1 2 u
2+ v2+ w2
(10) Where i is internal energy and the rest is kinetic energy, there could also be a term including gravitational potential energy.
2.2.2
Modeling Turbulence
Fluid flows in engineering applications are mostly turbulent, which requires that the whirls and vortices in the flow are resolved. There is always a question to what extend the flow needs to be resolved, which will result in a more or less ac-curate solution. The method which resolves all of the vortices and whirls is called Direct Numerical Simulation (DNS), this method also gives the most accurate re-sults. However this high level of accuracy is rarely wanted or possible, which is the reason for Large Eddy Simulation (LES), where larger vortices are resolved, but smaller are modeled. This method reduces the computational cost of the simula-tions. For engineering applications, Reynolds-Averaged Navier-Stokes equations (or RANS equations) are used, because they provide a satisfactory ratio of accuracy and computational cost.
Reynolds-Averaged Navier-Stokes
Reynolds-averaged Navier-Stokes equations are derived from the momentum equa-tion (Eq. 7) when Reynolds decomposiequa-tion is applied to those equaequa-tions. Reynolds Decomposition (Eq. 11) splits a variable into time-averaged part and fluctuating part [11].
φ(t) = φ + φ(t)0 (11)
Derived momentum equation with Reynolds decomposition is shown in the Eq. 12. This equation is also called Reynolds Averaged Navier Stokes equation [11].
ρuj ∂ui xj = ρfi+ ∂ ∂xj −p + µ ∂ui ∂xj +∂uj ∂xi − ρui0u j0 (12) The left side of the equation represents the change in the mean momentum of the fluid. On the right side, the first term is the mean body force, the second is stress from the mean pressure field and the third term are the viscous stresses. The last term is called Reynolds stress, which comes from the fluctuating velocity field and it is nonlinear. It also requires an additional equation to be solved, which is the reason for different turbulence models. Turbulence models can have several equa-tions which close the system, one of the most widespread turbulence models are k − ε and k − ω. Both of these rely on the analogy between the molecular gradient-diffusion process and turbulent motion, by using Boussinesq approximation. On the other hand, there are Reynolds Stress Transport models, which have the poten-tial to predict complex flows more accurately than eddy viscosity models because the transport equations for the Reynolds stresses naturally account for the effects of turbulence anisotropy, streamline curvature, swirl rotation and high strain rates [17]. Realizable k −
Realizable k − ε is a two-equation model, which is used for resolving the turbulent flow. This turbulence model was developed from Standard k − ε (SKE). Realiz-able k − ε turbulence model uses the Boussinesq hypothesis, which states that the Reynolds stresses are proportional to the mean rate of deformation (Eqs. 13 and 14).
− ρu0iu0j = µt ∂u0i ∂x0j + ∂u0j ∂x0i ! −2 3ρkδij (13) k = 1 2 u02+ v02+ w02 (14) One of the improvements from the SKE is that it has a new formulation for the turbulent viscosity, which means that the coefficient Cµ is expressed as a function
of mean flow and turbulence properties, rather than assumed to be constant as in the standard model (Eq. 15) [17].
Te =
k ε µt= ρfµCµTe
(15)
There Te is large eddy timescale, Cµ is a model constant, fµ is damping function
and ρ is density. Large eddy timescale Te suggests that the dissipation rate of
turbulent kinetic energy is constant at all scales. The second difference is that it has a new transport equation for the dissipation rate ε. The equation for kinetic energy production and dissipation can be seen in the Eqs. 16 and 17 [17].
∂ ∂t(ρk) + ∂ ∂xj (ρkuj) = ∂ ∂xi µ + µt σk ∂k ∂xj + Pk− ρε + Sk (16) ∂ (ρε) ∂t + ∂ (ρεuj) ∂xj = ∂ ∂xj µ + µt σε ∂ε ∂xj + 1 Te Cε1Pε−Cε2f2ρ ε Te + ε0 T0 +Sε (17)
In both equations, there are five main terms and all of them represent a contribution from a different cause. The first one is the Rate of change of k or ε, the second one is the transport of k or ε by convection. On the right-hand side, there is a transport of k or ε by diffusion, rate of production of k or ε and the last one which is negative is the rate of destruction of k or ε. The coefficients which are used in those equations can be seen in Table 4.
Table 4: Model coefficients [17]
Coefficient Value CT 1 Cε1 1.44 Cε2 1.9 Cµ 0.09 σε 1.2 σk 1
Reynolds Stress turbulence
Elliptic Blending model is the Reynolds Stress Transport model, which is also known as second-moment closure problems. The transport equation for fluctuating velocity is obtained by subtracting the RANS equation (Eq. 12) from the Navier-Stokes equation. [18] Equation is later multiplied with fluctuating velocity, after simplifying and averaging the equation which was derived can be seen in the Eq. 18.
ui ∂uj0uk0 ∂xi | {z } C = −ui0uk0 ∂uj ∂xi − ui0u j0 ∂uk ∂xi | {z } P − 2ν∂uk 0 ∂xi ∂uj0 ∂xi | {z } DS +p 0 ρ ∂uj0 ∂xk +∂uk 0 ∂xj | {z } PSC + ∂ ∂xi −ui0uj0uk0+ ν ∂ ∂xi (uj0uk0) − p0 ρ (δijuk 0+ δ ikuj0) | {z } DF (18)
Equation 18 can be split into several parts, begin with, the first term marked with C, which stands for convective transport. The first term on the right-hand side stands for Production of Reynolds Stress by mean velocity gradients. This term is one of the few terms which do not require any kind of modeling since it is defined purely with flow quantities. The second term marked with DS stands for the dissipation rate of Reynolds stresses by viscous action as the smallest scales of turbulence. The smallest scales are assumed to be isotropic, which requires some modeling. However, this approximation is not applicable near the walls of surfaces, where the dissipa-tion tensor becomes anisotropic. The third term, named P SC, gives a correladissipa-tion between fluctuating pressure and fluctuating strain rate. Its task is to redistribute the kinetic energy between the stress components. It plays an important role when determining how anisotropic the Reynolds Stresses actually are. The last term DF is a combination of several diffusion terms, which redistribute stresses accordingly. The software uses a slightly revised transport equation, which can be seen in the Eq. 19 [17].
∂
∂t(ρR) + ∇ · (ρRv) = ∇ · D + P + G − 2
3IγM + φ − ρε + SR (19) Where v is the mean velocity, R is Reynolds Stress Tensor, D is Reynolds Stress Dif-fusion, P is Turbulent Production, G is the Buoyancy production and I is identity tensor. Symbol γM stand for Dilatation Dissipation, ε stands for Turbulent
Dissi-pation rate tensor, φ is pressure strain tensor and the last term SR is user specified
source.
To begin with, a term marked D, is Reynolds Stress Diffusion, which is modeled by using an isotropic form of turbulent diffusion, seen in Eq. 20.
D = µ + µt σk ∇R (20)
Where µ is dynamic viscosity and σk is a model coefficient. This formulation of the
diffusion is simpler and considers that diffusion is isotropic. Turbulent viscosity µt
µt=ρCν k2 ε k =1 2tr (R)) (21)
Where tr(R) is the trace of the Reynolds Stress tensor and Cµ is model
coeffi-cient. The next term is Turbulent production P, which is obtained directly, without recourse to modeling.
P = −ρ R · ∇vT + ∇v · R
(22) The Eq. 23 is calculated from mean flow velocity, Reynolds stresses, and density, which is the reason it doesn’t have to be modeled. The next term is called Buoyancy production marked with G.
G = −β µt P rt
∇T g + g∇T
(23) Where β is coefficient of thermal expansion, g is gravitational vector, P rtis turbulent
Prandtl number and T is average temperature. This expression is derived for fluids with constant density.
∂ (ρε) ∂t +∇ (ρεv) = ∇ µ + µt σε ∇ε +ε k 1 2Cε1(tr(P) + Cε3tr(G)) − Cε2ρε (24) The isotropic turbulent dissipation rate is derived from a transport similarly as it was done with k − ε turbulence model. The last term is called Dilatation Dissipation Rate γM.
γM =
ρCMkε
c2 (25)
Where CM is model coefficient and c is the speed of sound. The model coefficients
which are used in the Elliptic blending model can be seen in Table 5.
Table 5: Model coefficients for Reynolds stress Model [17]
Coefficient Value σε 1.15 σk 1 CM 2 CS 0.21 Cε1 1.44 Cε2 1.83 Cµ 0.07
The turbulence model which was used with the Reynolds Stress model was the Ellip-tic Blending model. It is a low Reynolds number model based on an inhomogeneous near-wall formulation of the quadratic pressure-strain term. It uses a blending func-tion to blend the viscous sublayer with the log-layer formulafunc-tion of the quadratic pressure strain and dissipation. The transport equation can be seen in the Eq. 26.
φ − ε = 1 − α3
φw− εw + α3φh− εh
1 =α − L2∇2α (26)
The blending parameter for the elliptic blending model is α which switches from the near-wall layer to the quasi-linear version in the outer region. The pressure strain term formulation for the near-wall region can be seen in the Eq. 27. Once the blending function detects that the distance from the wall is too large, it will switch to formulation seen in the Eq. 28.
φw = −5ρ k R · N + N · R −1 2R : N(N + I) (27) φf = − [C1ρε + C1str(P + G)] A + C3− C3s √ A : A ρkS +C4ρk S · A + A · S − 2 3A : SI + C5ρk W · A + A · WT (28)
The terms in the equations marked with P, G, A, S, W have all been introduced before. Two additional variables are I, which is identity tensor and N, which is computed from wall normal vector. In addition to those matrices, some additional Model coefficients were added. Similar thing happens when deriving the equations for the dissipation rate, seen in the Eq. 29 for near wall region and in Eq. 30 for freestream.
φw = Rε
k (29)
φw = 2
3εI (30)
The coefficients which are added for Elliptic blending model can be seen in the Table 6.
Table 6: Model coefficients for Elliptic Blending Reynolds stress Model [17]
Coefficient Value A1 0.115 C1 1.7 C1s 0.9 C3 0.8 C3s 0.65 C4 0.625 C5 0.2
This turbulence model provides more accurate results since it uses fewer approxi-mations and modeling than a standard two-equation turbulence model. It naturally accounts for the effects of turbulence anisotropy, streamlines curvature, swirls ro-tation, and high strain rates [17]. However this does not come without drawbacks, RSM models need to solve eleven equations compared to two equation turbulence models, which are only solving 6 equations.
2.2.3
Near wall flow
A boundary layer forms on all of the non-slip surfaces which fluid interacts with. The fluid is slowed down because of the shear stress which is caused by friction between the layers of the fluid. Resolving the boundary layer is important because it affects the velocity profiles and its accuracy. Solvers split the near-wall region into several regions. The first region is the viscous sublayer (y+< 5), where the flow is laminar, due to stronger viscous forces. The second region is the buffer region (5 < y+< 30) and the third region is the log-law region (30 < y+< 500) where inertial forces start to prevail over viscous forces. Everything further away is considered as freestream. These regions can be seen in Figure 6. To understand these regions, non-dimensional values of wall distance and velocity need to be introduced. To begin with, friction velocity (Eq. 31), which will be later used to normalize the velocity.
uτ =
r τw
ρ (31)
Where τw is wall shear stress in streamwise direction and ρ is the density. The
equation for u+ can be seen below (Eq. 32).
u+= u uτ
(32)
The last parameter needed, is non-dimensional wall distance, which is marked as y+. This value represents the wall distance and is later used to split the boundary layer into different regions.
y+= y uτ
ν (33)
As aforementioned, different regions in the boundary layer are modeled in a different way. Linear law is used close to the wall and the relation between u+ and y+ can be seen in the Eq. 34.
u+= y+ (34)
However, further away from the wall this relation can no longer be used because it is just not accurate enough. Which is the reason that Log-Law is introduced (Eq. 35).
u+= 1 κln(y
+) + C
1 (35)
Where κ is von Karman’s constant with the value of 0.41, while C1 is also a constant
with the value of 5. With the help of these relations, Figure 6 can be constructed and different regions visualized. The red line shows the DNS results, and the blue lines are just two different laws, namely linear and logarithmic law.
Figure 6: Wall law [19]
2.2.4
Rotating Motion
Rotating motion can be simulated with different approaches. The most basic ap-proach is called Moving Reference Frame, where the domain will remain stationary. Force will be applied to the fluid in the domain, which will mimic the behavior of the spinning domain. This method has a low computational cost compared with others. It also comes with few drawbacks and limited applicability.
A step up from the MRF approach is Rigid Body motion, where the rotating domain actually spins. It will require more computational power because the solver needs to rotate the domain, which will result in mesh misalignment. One should be careful that the spinning domain does not propagate too fast in time, which might result in solver instabilities.
2.3
Multiphase flow
Multiphase flow is a system where different phases (gas, liquid, or solid) coexist and interact with each other. These flows are complicated for solvers since there is an interface dividing one phase from another phase. The software uses two different ap-proaches when it comes to multiphase flows. The first one is the Eulerian multiphase model and the second one is the Lagrange multiphase model. The only difference between these two models is that the Eulerian multiphase grid is solving equations on a fixed grid, while the Lagrangian model is tracking particles. Both models will be presented, distinct features of each model will be pointed out. The last method which is used in DNS codes when resolving particles is called Immersed Boundary Method, which will not be investigated any further in this thesis.
Eulerian Multiphase model
The Eulerian Multiphase model is used for simulations that have several Eulerian phases. There are different ways of how the program treats phases. The first one is the Eulerian Multiphase Mixture, which treats the mass, momentum, and energy as a mixture of quantities. This means that it uses weighted physical properties to solve the transport equation. This model is computationally efficient, due to the fact that it treats several phases as one.
The second model which can be used is called Multiphase Segregated flow. This method solves the conservation equations for every phase separately. Furthermore, the phases are expected to mix on length scales smaller than length scales which are resolved.
The last model in Eulerian multiphase flow is Volume Of Fluid (VOF), which should be used to simulate two unmixable phases. This model can handle free surface, which is the reason, to utilize it for sloshing flow in a water tank, where the free surface always remains smooth. However, this model assumes that the velocity, pressure, and temperature fields are shared within the phases.
Lagrange Multiphase model
Phase continua is always expressed in Eulerian form. On the other hand Lagrangian Multiphase model allows to model and solve arbitrary number of particles of a disperse phase. In Lagrangian model, the particles can be combined together in a parcel which are followed through the continua, as seen in the Fig. 2. These parcels are also known as Lagrangian phases.
2.3.1
Particle Forces
Lagrangian Multiphase model consists out of different forces which can be applied to the particles when they are coupled with the flow. Introducing less forces on the particles will result in less accurate results, since some physics is left out. Fur-thermore, particles can also be coupled to the flow in different ways, which will be presented in this chapter.
Particles in the flow can subjected to the flow behaviour in several different ways. The most basic flow and particle coupling is that mass-less particles are introduced in the domain, which means that the particles will follow the flow and the flow will not be affected by those particles. The next step is to introduce one-way coupling which makes sure that particles follow the flow but the fluid does not see the parti-cles. In other words the Eulerian phase can exchange momentum, mass and energy with Lagrangian phase, but not vice versa. One step further from this is two-way coupling, this introduces fluid-particle and particle-fluid interaction, in other words the Eulerian and the Lagrangian phase can exchange momentum, mass and energy with each other.
With the DEM method, the four-way coupling is introduced. This means that fluid and particles interact with each other. In addition to this interaction, particles can also interact with each other. This is done with a contact force which is described below.
Momentum Equation
The equation of conservation of linear momentum of a particle of mass mp is
pre-sented in the Eq. 36. This equation was derived from Newton’s second law, which states that the sum of forces on the particles will be equal to mass times acceleration.
mp
dvp
dt = Fs+ Fb (36)
Where vp is the instantaneous particle velocity, Fs are surface forces acting on the
particle, and Fb is the resultant of the body forces. Both of the forces can be split
into different terms, which are applied to the particle (Eq. 37) [17].
Fs = Fd+ Fp+ F l + Fvm
Fb = Fg+ FM RF + Fc+ FCo
(37)
Surface forces can be split up in different components, the first one FD is the Drag
force, the second one Fp is Pressure gradient force, the third one Fl is Lift Force
and the last one Fvm is Virtual mass force. A similar thing can be done with body
forces, were the first part Fg stands for gravity, the second contributor FM RF is the
force that comes from moving reference frame. The fourth term Fc is contact force
between particles and the last one FCo is Coulomb force.
Drag Force
Drag force on the particle is calculated using the Eq. 38. Fd=
1
2CdρAp|vs| vs (38)
Where Cdis the drag coefficient of the particle, ρ is the density, Ap is the projected
area of the particle, and vs is the particle slip velocity. There are several different
DEM particles is Di Felice [17]. This model introduces extra term which accounts for the presence of other particles around a particle. The formulation of this method is given by Eq. 39.
ζ = 3.7 − 0.65 exph−0.5 (1.5 − log [Rep])2i Cd= 0.63 + 4.80 piRep !2 2−ζ (39)
Where Rep is Reynolds number from the particle, i is the void faction around a
particle and 2−ζ is the contribution from the other particles around it, which is the reason it is recommended to use with DEM (4-way coupling). While Schiller-Naumann correlation is more appropriate to be used when using only the Lagrangian Particle Tracking method [17].
Cd= 24 Rep 1 + 0.15Re0.687 p Rep ≤ 103 0.44 Rep > 103 (40)
The correlation seen in the Eq. 40 is mostly dependant on particle Reynolds number Rep, and the main difference between Schiller-Naumann and Di-Felice is that latter
takes into account particles around it. Schiller-Naumann correlation is based on Stokes and Newton’s law when it comes to the calculation of the drag coefficient, which is important when assessing the settling velocity [17].
Pressure Gradient Force
Pressure gradient force is defined as it can be seen in the Eq. 41.
Fp = −Vp∇pstatic (41)
Where Vp is the volume of the particle and ∇pstaticis the gradient of the static
pres-sure in the continuous phase. Prespres-sure term also accounts for buoyancy term, which is in some cases important force when particles have different densities compared to the fluid.
Virtual Mass Force
Virtual mass force is defined in Eq. 42. Fvm= CvmρVp Dv Dt − dvp dt (42) Where Cvm is the virtual mass coefficient and the default value is 0.5 for a sphere in
a uniform, inviscid, and incompressible fluid, which is later altered for different con-ditions. The operator inside the brackets is the material derivative. When a body is accelerated through a fluid phase, this body does work. In other words, virtual mass term accounts for the form drag due to acceleration. This work is related to the virtual mass effect. According to Sommerfeld [20] Virtual Mass force become
insignificant when ρc/ρd∼ 10−3.
Spin Lift Force
Lift force can act on particles, there are two ways of producing lift force, the first one is due to the spinning of the particle and the second one is due to the shear force of the fluid on the particle. Lift force due to rotation can be seen in the Eq. 43 and 44. FLR = ρπ 8 D 2 pCLR|vs| Ωp× vs |Ωp| (43) CLR= 0.45 + Rer Rep − 0.45 e−0.5684Re0.4R Re0.3p (44)
Equation 43 is a classic equation for lift force, however, instead of the velocity, there is an angular velocity Ωp, which takes into account the rotations of the particle. Lift
Coefficient is calculated using the Eq. 44, which is obtained with Sommerfeld em-pirical relation [17]. These two equations calculate the lift force which is generated with rotation (Magnus effect). Lift force due to the rotation of the particles can only be accounted for when using DEM particles. Lift force does not necessarily mean that the force is acting in the upward direction, it depends on the rotation direction of the particle.
Shear Lift Force
Lift force which is generated from shear force from fluid to particle Eqs. 45 and 46 is also known as Saffman Lift force [20].
FLS = CLS ρπ 8 D 3 p(vS× Ωp) (45) CLS = 4.1126 Re0.5 s f (Rep, ReS) (46)
To begin with lift coefficient (Eq. 46), it is calculated with Sommerfeld relation which is a function of Rep (Reynolds number of particle) and ReS (Reynolds
num-ber of shear flow). In equation for lift shear force (Eq. 45) there is Ωp which stands
for the curl of the flow and vS is slip velocity.
Contact Force
The first body force which will be described is called Contact Force. It is one of the most important forces in DEM modeling and it is added only in the DEM method. It is calculated with the Eq. 47.
Fc=
X
contacts
Fcm (47)
Where Fcm is a contact model force calculated using Hertz Mindlin [17]. Forces
between particles are split into two components normal and tangential (Eq. 48) as seen in the Figure 7.
Fcm= nFn+ tFt (48)
Figure 7: Contact model for particles [17]
The particle interactions are considered normal damper spring models, in addition to this, there is also friction in the tangential model, which models particle rolling against each other and the boundaries. The equation for normal force and its coef-ficients is shown in the Eq. 49.
Fn= −Kndn− Nnvn Kn= 4 3EeqpdnReq Nn=p5KnMeqNn damp (49)
Where coefficient Knis normal spring stiffness and Nnis the normal damping factor.
When calculating tangential forces, there is first a condition which checks whether the particles slip or not (Eq. 50).
|Ktdt| < |Kndn| Cf s
Ft= −Ktdt− Ntvt
(50) In this equation Ktstands for tangential spring stiffness, dtis overlap of the particles,
Kn is normal spring stiffness and Nn is normal damping factor and Cf s is static
friction coefficient. If this condition holds then tangential force is calculated with Eq. 50. Otherwise, the Eq. 51 is used to model tangential forces in a particle-particle and particle-wall interaction.
Ft= |Kndn| Cf sdt |dt| Kt= 8GeqpdtReq Nt=p5KtMeqNt damp (51)
Where coefficient Kt is tangential spring stiffness and Nt is tangential damping
factor. There are several other coefficients: Nt damp stands for tangential damping
coefficient, Geq is equivalent shear stress module, Req is equivalent radius of two
particles, Eeq is equivalent Young’s modulus and Meq equivalent particle mass.
Force of Gravity
The last important force which is needed is gravitational force (Eq. 52).
Fg = mpg (52)
Where mp is mass of the particle and g is gravitational acceleration vector. In some
cases, this will be the dominant force that will be acting on particles. Moving Reference Frame and Coloumb Force
The last two contributions come from Columb force FCo and Moving Reference
Frame force FM RF. The Columb force is a force that comes from electrical charges,
which results in an additional force acting on the particle. The Moving Reference Frame force is a force that comes from the rotating domain. Every particle that passes through the MRF domain gets a force applied to it, which results in this part FM RF to be added to the momentum equations for particles. Moving reference
frame was never applied to particles in this project. Turbulent Dispersion
Particles in the flow field are affected by many different forces and fluctuating veloc-ity components. The phenomena are similar to Brownian motion which just adds random motion, but when particles are subjected to turbulent eddies the magnitude and the direction vary. If particles are small enough to follow the instantaneous fluid motion, particle dispersion is equal to that of a fluid element. The equation which is used to calculate the eddy velocity scale can be seen in Eq. 53. [20]
ue= lt τt r 2 3 (53)
The turbulence model provides the length and time scales of the turbulence. In addition to the eddy velocity scale, the eddy time scale needs to be added, which measures the lifetime of the eddy. In other words, the maximum interval over which a fluctuating velocity v0 is applied to the particle. It is defined with the diffusion of a scalar, depending on the selected turbulence model, Eq. 54.
τe=
2µt
ρu2 e
(54) Particles can only interact with the eddy for less than the eddy time-scale when they have non-zero slip velocity. It can be said that turbulent dispersion defines how the particle will interact with eddies.
Figure 8 depicts the movement of the particle when it is moving through the flow field. The particle is larger and it does not follow the instantaneous fluid motion.
Figure 8: Particle motion in turbulent velocity field
Additionally, Figure 8 also shows that particles will change the direction with a slight delay since it has its own momentum.
Particle Time Scale
Particle time step depends on three different particle time-scales, which define the timestep that should be used when resolving particle movement. The first one is Rayleigh wave propagation time τ1, which is defined with the Eq. 55.
τ1 = π
R VRayleigh
(55) The equation above limits the time-step for DEM particles, that the force acting on the particle is only affected by its close neighbors. Where R stands for half of the characteristic dimension and VRayleigh is Rayleigh wave velocity, which depends on
material properties. Maximum timestep must be no more than timescale multiplied with τ1. The second time-scale is Impact duration (Eq. 56).
τ2= 2.94 5√2πρ 4 1 − ν2 E !25 R 5 √ νimpact (56) The relation limits the moving particles when colliding, if the timestep would be too big that would result in instabilities because the particles overlap too much. The maximum timestep of the second particle cannot be larger than 0.1τ2. This
timescale is defined mainly with material properties. The last timescale is called Particle transit time which is purely geometric (Eq. 57).
τ3=
R νparticle
(57) This timescale limits the movement of the particles within one timestep, which pre-vents missing collision with particles and walls. Timescale must be lower than 0.1τ3
In the end, user-selected timescale must fulfill all of this aforementioned time scales, if it does not then software automatically selects the highest timescale required.
3
Method
3.1
Channel flow
Channel flow is one of the fundamental engineering flows, which are well known and have experimental data for comparing the solutions. The geometry for channel flow used was the same as the one the study performed by KTH [1]. The same study will be used to compare results obtained from simulations.
3.1.1
Geometry
The geometry which was used for the channel flow simulations can be seen in Figure 9.
Figure 9: Geometry used for channel flow
The geometry is the channel where both edges have the same length (h = 0.015 m) and the length of the domain is 0.5 m, also seen in Table 8. The domain length was determined by simulating the same domain which was extended by 4 times. Domain study was important to cut down computational time and to get a better understanding of what is happening in the domain if it is extended. Wall shear stress was monitored and when the wall shear stress did not change for more than 2% when elongating the channel, the solution was considered as independent from the domain length, these simulations were performed without particles. The domain study can also be seen in Figure 10, where results are compared with empirical relations for Fanning Friction Factor. The relatively low error (less than 2%,) between empirical and simulation results indicated the comparability of the two methods. The final domain was split into two sections, each measuring 0.25m. Only the second section was evaluated in terms of Fanning friction factor in the results.
0 0.5 1 1.5 2 Distance (m) 0.009 0.0095 0.01 0.0105 0.011
Fanning Friction Factor [-]
Domain study-FF
Simulations Moody
Figure 10: Fanning Friction Factor in relation to domain length
3.1.2
Solver Setup
Solver setup was divided into two parts, the first part consists of boundary conditions and settings which were used for the Eulerian phase. Lagrangian phase parameters will follow.
Eulerian phase
The turbulence model was modeled using Realizable k-ε. Boundary conditions pre-scribed on the domain were as follows in Table 7, one face was mass flow inlet, the opposite face was an outlet and everything else was a no-slip wall. The reason behind choosing specified split ratios at the outlet is that it is easier for the solver to handle simulations with particles. Boundary condition outlet with specified split ratios calculates mass flow out of the domain to satisfy the continuity equation. The boundary condition on the inlet was uniform flow velocity, so the velocity profile needs time to develop. The time step used in this study was 0.005 s, the selected timestep was taken from internal reports provided by TetraPak.
Table 7: Boundary conditions
Boundary Condition Type Value Temperature
Inlet Mass Flow Inlet 0.087218 [kg s−1] 300 K
Wall No-slip - 400K
Outlet Specified Split Ratio 1.0
-Mass flow, velocity, and the domain size had to be calculated from data provided by KTH [1], that way it can be ensured that obtained results can be directly comparable with simulations performed with STAR CCM+. Reynolds number was specified, using this number the rest of the parameters needed can be calculated. Calculated parameters, which were used in the set up can be seen in the Table 8.
Table 8: Parameters for Channel flow Coefficient Value Re2h 5600 P r 7 2h Dparticle 15 Dparticle 0.002 m h 0.015 m µ 0.0010383 P a s ρ 1000 kg m−3 Lagrangian phase
Momentum equations are solved for Lagrangian particles in the flow, which means that some parameters need to be set. Some of the important parameters can be seen in Table 9. The shape of the particles is a sphere, which is the same shape used in the study performed by Ardekani et al. [1].
Table 9: Parameters for Particles
Coefficient Value Dparticle 0.002 m Volume Fraction 10 % ρ 1000 kgm−3 ν 0.0 E 74000.0 P a ρ 1000 kg m−3
Normal Restitution Coefficient 0.0 Tangential Restitution Coefficient 0.0 Coefficient of Rolling Resistance 0.001
Cf s 0.01
Forces that were included in the calculation in channel flow were: Drag, Shear Lift, Spin Lift, Virtual mass, Pressure gradient, Gravity, and Contact force. Several dif-ferent models can be used for Drag and Lift coefficient calculation in software, in this case, Di Felice drag coefficient was selected. The reason for selecting Di-Felice is that it takes into account particles close by. Besides, to drag and lift force, rota-tional drag torque, spin lift force, and shear lift force were added and for all of them, Sommerfield formulation was used. Most of the forces have only one formulation which can be used, thus the options are limited. Shear lift force is the only force that has a second option, which is called the Saffman lift force, which is meant to be used in lower particle Reynolds numbers.
Monitoring parameters were very important since this is a transient simulation, there have to be enough inner iterations so the solution converges in inner timestep. Monitoring parameters were Wall Shear Stress in a streamwise direction on two dif-ferent wall sections. A very important parameter to monitor particle convergence is the momentum source of the particles. Particles are moving and rotating in every
timestep differently which requires strict convergence criteria on their momentum generation.
For particle injection in the domain, the injector must be defined, which is a surface from which particles will be injected in the domain. This surface is the face defined as inlet, with a specified mass flow of the particles inside the domain. This helps to control the number of particles in the domain, so for example, if volume fraction of particles is desired to be 10% this means that 90% of the total mass flow will be for the Eulerian phase and the rest will be prescribed to Lagrangian phase. This method can be used if the particles have the same density as the fluid.
The simulation is considered as four-way coupled, which means that particle-particle interactions are computed. As seen in Table 9, all of the restitution coefficients were set to 0, which means that the particles do not bounce, because all of the kinetic energy during the collision is lost.
3.1.3
Mesh
A mesh that was generated for the channel flow consisted out of hex elements and prism layers near the wall to capture the velocity profile accurately. Mesh indepen-dence study was performed firstly for the simulations without the particles, where only the base size in the bulk was varied. The inflation layer remained the same during the grid sensitivity study, if the results would differ from empirical relation significantly, one should consider refining the near-wall region. Transition to the bulk mesh is extreme in coarse mesh, but it gets better with refining the mesh.
Figure 11: Meshes which were used for mesh independence
Mesh M1 is the coarsest and M5 is the finest mesh tested, the mesh was refined up to a point where the results of the Fanning friction factor did not change for more than 1%. Furthermore, the solution was also compared with the empirical relation for Fanning Friction factor in internal flows, the results from simulations deviate for 1% from factor calculated using the Moody chart. The number of elements in each mesh and size of the elements in the bulk can be seen in Table 10. Additionally, there is also a row which contains particle diameter and cell size ratio, which will be important for the future work.
Simulations with particles are very sensitive to mesh size because if a particle is big-ger than the cell it means that the cell will be mostly occupied by the Lagrangian phase. However, STAR CCM+ limits the volume fraction, which is the reason for a thorough investigation of what mesh should be used when simulating flow with
Table 10: Mesh sizing and number of elements Mesh M1 M2 M3 M4 M5 Base size [m] 0.00325 0.0025 0.002 0.0015 0.00125 Number of elements 19712 40800 72000 126920 192000 Dp Cell size 0.62 0.80 1.00 1.33 1.60
particles, and what is the best compromise between computation time and accuracy. The first step in this investigation was to run the particle simulation on the same mesh and observe the behavior.
The software offers an option to smooth the result over several cells, so the afore-mentioned problems do not appear. Smoothing means that the effect of the particle in the cell is smeared over several cells instead of one single cell. This can be useful in the regions where the cells are small i.e. near refinement regions and near the wall. Distributing the particle load on the Eulerian phase might not be physical, but it helps solver to handle particles.
3.2
Mixing tank
3.2.1
Geometry
The mixing tank has an agitator inside which mixes the product inside the tank. The 3D model of the geometry must be accurate since this will affect the results. The geometry for both mixers was provided by Alfa-Laval, which can be seen in Figure 12.
(a) Isometric view (b) Top view
Figure 12: Geometry used for simulations, diameter of the impeller is 400 mm
The shown geometry of the agitator in Fig. 12, has a diameter of 400 mm and was used for the off-centered tank simulation. The off-centered tank has an impeller shaft positioned 0.16 m from the center of the tank. This is common practice to improve mixing efficiency.
(a) Top view (b) Side view
Figure 13: Geometry used for simulations, diameter of the impeller is 400mm
Figure 13 shows geometry which was used for the simulations, the diameter of the tank is 1 m with the height of 0.639 m. Monitoring points that were spread out through the domain can be seen in Figure 13, they were used to monitor velocity in
the domain. In addition to velocity, the power number of the agitator was monitored and compared with manufacturer value.
3.2.2
Mesh
The mesh was constructed using polyhedral elements, since polyhedral mesh cell faces are more likely to align with the flow direction to reduce numerical diffusion. The different numbers of inflation layers were used to capture the velocity profile on the blades accurately which resulted in a value of y+ below 5 around the whole domain.
Figure 14: Four different meshes which were constructed for stirring tank case
Power number was a parameter of interest and also the main indicator for mesh independence. Power number gives a general overview of the flow in the stirring tank. For a detailed view velocity points, which were spread around the mixing tank, were also compared between meshes. Four different meshes were compared, varying only the base size. The Fine mesh was constructed with the help of Siemens support, because Medium mesh was unstable when the particles were injected.
Table 11: Mesh size compared with particle size
Mesh Num of cells Base [m] Phase 1 Dp1/Cell Phase 2 Dp2/Cell
Extra Coarse 165k 0.025 0.012 0.48 0.01 0.4
Coarse 300k 0.014 0.012 0.86 0.01 0.71
Medium 600k 0.01 0.012 1.2 0.01 1
Fine 1800k 0.007 0.012 1.71 0.01 1.43
Table 11 shows four different meshes which were used, Fine mesh has refined regions near the agitator, due to convergence issues with previous meshes.
3.2.3
Solver Setup
The mixing tank consists out of two domains, inner and outer. The motion of the do-main can be simulated in various ways, the first option is Moving Reference Frame, which is also the simplest. The second option is Rigid Body Motion, which is more accurate but requires more computational time because the domain spins. For this case, RBM method was utilized.
Table 12: Parameters rotational domain
Coefficient Value
Tip speed 2.1 m s−1 Impeller Diameter 0.4 m
n 100 RPM
Table 12 contains the values which were used for the rotating domain in this sim-ulation setup. The simulated liquid is water and this simsim-ulation is considered as isothermal. Boundary conditions on the tank walls were a no-slip walls, except the top wall was an exception. In reality, the top surface is in contact with air, which can be defined as a free surface, therefore the top surface is a slip wall. This is an assumption that the agitator is spinning with low speed, which results in a flat top surface.
Lagrangian phase
The particle properties used in experiments tried to be recreated and used also for simulations. In experiments potatoes and peppers were used as particles, and later on referred to as Phase 1 for potatoes and Phase 2 for peppers. In experiments the particles were mainly shaped like cubes, however, this is not possible to recreate in software, which is the reason that the particles will be modeled as spheres.
Table 13: Parameters for Particles in a mixing tank
Coefficient Phase 1 Phase 2 Dparticle 0.0012 m 0.001
Mass Fraction 0.5 % 0.55 %
ρ [kg m−3] 1020 993
Determining some parameters from the Table 13 was hard, due to the fact that it is hard to measure each individual parameter. Density was one of the most impor-tant parameters, since the particles in the tank will be affected by buoyant forces. The density of both phases was calculated with sinking and rising velocity from measurements and later compared with articles [21][22]. The rest of the parameters were approximated by using the experience gained from channel flow. Particle size was calculated from the samples provided from the experiments, see Appendix A. Poisson’s ratio and Young’s modulus were determined from previous simulations of channel flow, where it was discovered that lowering Young’s modulus, did not have an effect on results. Normal and tangential restitution coefficients of the wall were both 1, which means that during the collision the kinetic energy is conserved. Which
will result in rebound particle-wall interaction.
The first step was to simulate low concentration with Lagrangian Particle Tracking (LPT). Using LPT in Star CCM+, one can select which forces to apply to the parti-cles. The simulation plan was to investigate every force in isolation from other forces, so the force effect could be seen and evaluated. Turbulent dispersion and gravity were always utilized. Before simulating a rotating domain and particle interaction with the agitator, a test run was needed with stationary domain. In theory, all of the particles should remain where they are if the density of the particles and water is the same. The simulation plan is the following:
• Test Run: Fluid with zero velocity, Pressure gradient force, Drag force, Shear lift force, Virtual mass force, particles have the same density as water
• Simulation 1: Pressure gradient force, particles have the same density as water • Simulation 2: Pressure gradient force, Drag force, particles have the same
density as water
• Simulation 3: Pressure gradient force, Drag force, Shear Lift force, particles have the same density as water
• Simulation 4: Pressure gradient force, Virtual mass force, particles have the same density as water
• Simulation 5: Pressure gradient force, Drag force, particles have different den-sities compared with water
4
Results
4.1
Channel flow
Mesh independence study is important in CFD methodology because mesh also induces error. Which is the reason for the mesh study presented in Fig. 15. The Moody relation which is used as a reference value was read of Fig. 4.
Mesh independence study without particles
M1 M2 M3 M4 M5 0 0.002 0.004 0.006 0.008 0.01 0.012
Fanning Friction Factor [-]
Numerical simulations Moody relation
Figure 15: Mesh study for flow without particles
Five meshes were tested and as it can be seen after Medium mesh (M3), the results do not change for more than 1%. The same five meshes were used for simulations with particles and the results can be seen in Figure 16.
Mesh independence study with particles
M1 M2 M3 M4 M5 0 0.005 0.01 0.015 0.02 0.025
Fanning Friction Factor [-]
Numerical simulations DNS simulations
Figure 16: Mesh study for flow with particles
As seen from Figure 16, there is no asymptotic behavior, values of Fanning Friction factor keep increasing, regardless of the mesh size. Percentage difference can also be seen in Table 14, which also contains the ratio between cell and particle size.
