DEGREE PROJECT IN CHEMICAL ENGINEERING AND TECHNOLOGY, SECOND LEVEL,
STOCKHOLM, SWEDEN 2018
19/3-2018
KTH ROYAL INSITUTE OF TECHNOLOGY KTH CHEMICAL SCIENCE AND ENGINEERING
CFD for mixing efficiency in
commercial and industrial
advanced air oxidation
JOHAN BERGSTRÖM
i
DEGREE PROJECT
Master of Science in
Chemical Engineering and Technology
Title: CFD for mixing efficiency in commercial and
industrial advanced air oxidation
Keywords: Computational fluid dynamics, Mass transfer,
Chemical engineering, Advanced oxidation processes,
Transport Phenomena, Simulation, COMSOL
Work place: Ozonetech
Supervisor at the
work place: Francesco Montecchio
Supervisor at KTH Matthäus U. Bäbler
Student Johan Bergström
Date: 2018-03-19
Examiner Matthäus U. Bäbler
ii
Abstract
Ozone treatment of pollutants in air is a relatively young technology with limited literature available. To the authors knowledge no literature discussing simulations of commercial ozone treatment without UV-lamps in gas phase has been published up to this point.
The purpose of this project was to identify issues and propose recommendations related to the distribution of ozone in industrial ducting and the injection of ozone into commercial ducting. The injection rate of the ozone mixture is small relative to the flow of the treated air stream. In such cases ozone is easily swept away and confined to a limited section of the ducts, affecting overall efficiency. The injection inlet and ducting were simulated together in both 2D and 3D environments using COMSOL Multi-physics with CFD, CAD and transport of diluted species modules.
Improving mixing in industrial ducting was simple in comparison to the commercial ducting where the flow is complex. For the commercial application simulations
showed that the mixing efficiency varies greatly between injection positions. Based on CFD data in the commercial ducting two models for injection point analysis were developed.
2D and 3D simulations showed different result for injections inside the main duct, the 3D case could properly simulate rotating flows inside the main duct which makes certain injection points in the main duct more effective that predicted in 2D.
This master thesis project was done in cooperation between Royal Institute of Technology and Ozonetech in Sweden.
i
Table of content
1 Introduction ... 1
2 Theoretical background ... 2
2.1 Computational Fluid Dynamics (CFD) ... 2
2.1.1 Direct Numerical simulation (DNS) ... 4
2.1.2 Turbulence models ... 4
2.1.3 Wall functions ... 8
2.2 Ozone properties and reactions ... 9
2.2.1 Ozone reactions in VOC treatment ... 9
2.3 Mass transfer – Basic overview ... 11
2.3.1 Turbulent diffusion ... 12
2.4 COMSOL Multi-physics ... 13
3 Materials and methods ... 14
3.1 Industrial ducting... 14
3.1.1 2D Industrial ducting ... 14
3.1.2 Mesh refinement study ... 15
3.1.3 Injector placement ... 17
3.2 3D Industrial ducting ...19
3.3 Commercial ducting ...19
3.3.1 2D Commercial ducting ... 20
3.3.2 3D Commercial ducting ... 20
3.3.3 Inflow distribution in commercial ducting ... 21
3.4 Error estimation ... 22
3.4.1 CFD based tool for general injection point recommendations ... 22
4 Results – 2D simulations ... 22
4.1 2D - Industrial ducting ... 22
4.1.1 2D Injector placement - Industrial ducting ... 23
4.2 2D Commercial ducting ... 26
5 Results - 3D simulations ... 29
5.1 3D Industrial ducting ... 29
5.2 3D Commercial ducting ... 29
5.3 Injection point prediction based on CFD data ... 32
5.4 Simulation of inlet distribution ... 36
6 Discussion ... 39
ii
6.1 Industrial air duct ... 39
6.2 Commercial ducting ... 39
7 Conclusion ... 41
8 References ... 42
Appendix I - Thermal Decomposition ... 45
1
1 I NTRODUCTION
In the early days of chemical industries waste flows were left untreated to a large extent, following the advancements in environmental science, legal and economic incentives have increased the need for pollution control. Today industries have regulations imposed on them. Furthermore, unregulated odour emitted from both industrial and commercial applications also needs to be treated based on complaints from the public.
To comply with environmental regulations a wide variety of emission controlling technologies are competing. In the case of VOC emissions, thermal incineration is the most straightforward way to completely neutralize organic compounds in a gas and solid phase. Unfortunately, incineration is not always suitable due to economical, technical or safety reasons. Ozone treatment provides a comparatively gentle oxidation at lower temperatures. Unlike incineration, ozone treatment does not oxidize all carbon chains to carbon dioxide or carbon monoxide. Ozone is useful in such cases where it is sufficient to break down toxic compounds and VOCs partially to remove odour or toxicity. The implementation of ozone air treatment is simple, an ozone generator needs to be installed and the ozone can then be injected into already present tubing or ventilation.
For this thesis project Royal Institute of Technology cooperates with Ozonetech who supplies the kind of ozone treatment introduced above. Injecting small volumetric flows of ozone gas mixtures into a much larger industrial gas stream is expected to have the typical scale up issues of poor mixing and potentially dead zones affecting the efficiency. Up to this point Ozonetech has relied on their experience to ensure satisfactory results for their installations, the general goal is to contribute to the understanding of the systems and provide recommendations based on simulations.
The aim of the project is to implement Computational Fluid dynamic (CFD)
simulations coupled with turbulent mass transfer in 2D and 3D environments. The simulations are done in COMSOL Multi-physics to investigate the distribution of injected ozone at different points after injection. To improve mixing problems different placements are simulated and compared. Conclusions about the current issues and solutions are made both injection in industrial sized ducting and odour treatment in commercial applications.
Parts of the report has been removed, figures and text has been removed or altered due to containing sensitive information.
2
2 T HEORETICAL BACKGROUND
2.1 C
OMPUTATIONALF
LUIDD
YNAMICS(CFD)
The movement of gases and liquid has great importance in many fields, points of interest could be how a particle gets carried by a flow, how a solid object affects mixing or perhaps how liquid infiltrates small pores. Due to the difficulty of physically observing and quantifying these processes, computational models have been developed to observe these in a digital environment instead, using methods referred to as Computational Fluid Dynamics (CFD).
According to classical fluid mechanics, flows exhibit different behaviour depending on factors such as speed, viscosity and density of the flow. The different flow types are laminar, turbulent and laminar-turbulent transition flows. Whereas the laminar flow has smooth flow profile without backflow, while turbulent models has a more chaotic flow with vortices [1].
The fundamental equations in CFD are based on the principles that mass is conserved, the acceleration of an object is proportional to force according to Newton’s second law and that energy is conserved. The implication of these principles when deriving governing equations is that when defining a fluid flow in a finite volume element, the mass in and out of each one of these are the same and that there are forces exhibited on the volume element that provides momentum [2].
There are two ways of viewing the volume element, either it is locked in place with flow in and out, or it is seen as a moving volume element where the inflow is not balanced against the outflow out the element. These two forms are referred to as being in conservation form or non-conversion form respectively. Conservative forms are more suited for dealing with systems that have shocks, chemical reactions or other sharp interfaces [3].
The “go to equation” in CFD is the famous Navier-Stokes equation which describes the momentum of a viscous flow. Combining the original Navier-Stokes equation with Newton’s shear stress equations for Newtonian fluids yields what is called the complete Navier-Stokes equation [2], seen below is the momentum equation in z direction.
∂(𝜌𝑤)
∂𝑡 +∂(𝜌𝑢𝑤)
∂𝑥 +∂(𝜌𝑣𝑤)
∂𝑦 +∂(𝜌𝑤2)
∂𝑧
= −∂𝑝
∂𝑧+ ∂
∂𝑥[𝜇 (∂𝑢
∂𝑧+∂𝑤
∂𝑥)] + ∂
∂𝑦[𝜇 (∂𝑤
∂𝑦 +∂𝑣
∂𝑧)] + ∂
∂𝑧(𝜆 ∗ ∇ · 𝑉⃗ + 2𝜇∂𝑤
∂𝑧) + 𝜌𝑓𝑧
The symbols are the following, ∂ – partial derivative sign, 𝜌 – fluid density, 𝜇 - fluid viscosity, 𝜆 –bulk viscosity coefficient, ∇ – vector operator taking partial derivative in x,y and z direction , 𝑉⃗ – vector velocity field in the cartesian space, u – velocity in x direction, 𝑣 – velocity in y direction, 𝑤 – velocity in z direction, t – time, 𝜌𝑓𝑧 body force per unit mass acting on the fluid element in z direction.
To solve the complete Navier-Stokes equation, other equations are required as well, namely equations for viscous flow, continuity, energy and momentum. Together these
3 equations form a system, which is solved numerically [2]. The equations are given below:
∂𝜌
∂𝑡 + ∇(𝜌𝑉⃗ ) = 0
The continuity equation provides mass flux 𝜌𝑉⃗
∂(𝜌𝑢)
∂𝑡 + ∇(𝜌𝑢𝑉⃗ ) = −∂𝑝
∂𝑥+∂𝜏𝑥𝑥
∂𝑥 +∂𝜏𝑦𝑥
∂𝑦 +∂𝜏𝑧𝑥
∂𝑧 + 𝜌𝑓𝑥 x component flux of momentum 𝜌𝑢𝑉⃗
∂(𝜌𝑣)
𝑑𝑡 + ∇(𝜌𝑣𝑉⃗ ) = −∂𝑝
∂𝑦+∂𝜏𝑥𝑦
∂𝑥 +∂𝜏𝑦𝑦
∂𝑦 +∂𝜏𝑧𝑦
∂𝑧 + 𝜌𝑓𝑦 y component flux of momentum 𝜌𝑣𝑉⃗
∂(𝜌𝑤)
𝑑𝑡 + ∇(𝜌𝑤𝑉⃗ ) = −∂𝑝
∂𝑧+∂𝜏𝑥𝑧
∂𝑥 +∂𝜏𝑦𝑧
∂𝑦 +∂𝜏𝑧𝑧
∂𝑧 + 𝜌𝑓𝑧 z component flux of momentum 𝜌𝑤𝑉⃗
For Newtonian fluids, the shear stresses (𝜏) are given as below
∂𝜏𝑥𝑥 = 𝜆∇ · 𝑉⃗ + 2𝜇∂𝑢
∂𝑥
∂𝜏𝑦𝑦 = 𝜆∇ · 𝑉⃗ + 2𝜇∂𝑣
∂𝑦
∂𝜏𝑧𝑧 = 𝜆∇ · 𝑉⃗ + 2𝜇∂𝑤
∂𝑧 𝜏𝑥𝑦= 𝜏𝑦𝑥 = 𝜇 (∂𝑣
∂𝑥+∂𝑢
∂𝑦) 𝜏𝑥𝑧 = 𝜏𝑧𝑥 = 𝜇 (∂𝑢
∂𝑧+∂𝑤
∂𝑥) 𝜏𝑦𝑧= 𝜏𝑦𝑥 = 𝜇 (∂𝑤
∂𝑦 +∂𝑣
∂𝑧)
Stokes made the hypothesis that the bulk viscosity coefficient 𝜆 is the following:
𝜆 = −2 3𝜇
For ideal gases the pressure term is described as follows:
𝑝 = 𝜌𝑅𝑇
Different variations on the complete Navier-Stokes equations are made to involve different physics and phenomena. For example, if mass transfer is to be added then additional continuity equations in addition to an additional term for the energy transport due to the diffusion of species. [2]
4 2.1.1 Direct Numerical simulation (DNS)
To directly solve the Navier-Stokes equations provides the most accurate results in turbulent flow simulations. Direct numerical simulations (DNS) numerically provides the instantaneous and three-dimensional flow fields without trying to reduce
computational cost by averaging flow fields and having a lower resolution. Although DNS provides accurate solutions it also uses the most memory and computational power of all the CFD methods presented in this report. As a result, DNS mostly find usage in smaller simulations and not in engineering designs [4].
With the use of DNS, vortices larger than the grid size can be followed as it moves in time and space, this data can then be used to produce time average flow fields [3].
The resolution of the DNS simulation determines its accuracy, the use of a low resolution does not capture the phenomena smaller than the grid size which in turn could cause the kinetic energy to be transferred to larger vortices by the simulation when it is dissipated by smaller vortices. Which is why DNS require such a fine grid size and consequently large computational power and memory [3].
2.1.2 Turbulence models
Turbulent models should model turbulent flows, which in practice means that it can predict turbulent stress and vortices for these conditions. In contrast to DNS these turbulence models use approximations to model flows smaller than the grid size and provides simplifications which reduce the computational cost. The choice of a turbulence model is important for both computational power and accuracy of the solution. There is unfortunately no universal turbulence model that works for all turbulent flows [3].
Information about how the most common turbulence models work, their advantages and disadvantages are presented under the following headings.
2.1.2.1 Reynolds-averaged numerical simulation (RANS)
Reynold-averaging deals with the large-scale dynamics of a flow field. The mathematical definition is made so that the averaging operator 𝑢̅𝑘 which averages a flow field 𝑢𝑘 that has variations 𝑢′𝑘 = 𝑢𝑖− 𝑢̅𝑖 , the operation should satisfy the following conditions:
𝑢′̅
𝑖 = 0, 𝑢′̅̅̅̅̅̅ = 0,𝑖𝑢̅𝑗 𝑢̿𝑖 = 𝑢̅𝑖
The ensemble average that satisfies the stated conditions are called Reynold average.
It removes fluctuating components from the flow field variables.
By using the Reynold average, the continuity equation below:
∂𝑢𝑘
∂𝑥𝑘 = 0 Becomes
∂𝑢̅𝑘
∂𝑥𝑘 = 0
5 Another part of the Reynold average is Reynolds stress equation, which gets the average stress by employing the average velocity operator to stress equations. Taking the difference between the momentum equation and its Reynolds averaged equation followed by re-expressing the velocity fluctuation derivative terms yields the expression below:
𝐷̅𝜏𝑖𝑗
𝐷̅𝑡 = 𝑃𝑖𝑗 + Π𝑖𝑗 − 𝜀𝑖𝑗 +∂𝐽𝑗𝑖𝑘
∂𝑥𝑘
Where 𝐷̅ is the velocity gradient tensor based on the average velocity, 𝜏𝑖𝑗 is the stress, 𝑃𝑖𝑗 the pressure gradients, Π𝑖𝑗 is the pressure-strain correlation tensor, 𝜀𝑖𝑗 is the dissipation, 𝐽𝑗𝑖𝑘 is the diffusive flux. Π𝑖𝑗, 𝜀𝑖𝑗 and 𝐽𝑗𝑖𝑘 are defined as follows:
𝑃𝑖𝑗 = −𝜏𝑖𝑘 ∂𝑢̅𝑗
∂𝑥𝑘− 𝜏𝑖𝑘 ∂𝑢̅𝑖
∂𝑥𝑘 Π𝑖𝑗 =𝑝′
𝜌 (∂𝑢′𝑗
∂𝑥𝑗 +∂𝑢′𝑗
∂𝑥𝑖)
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜀𝑖𝑗 = 2𝑣∂𝑢′𝑖
∂𝑥𝑘
∂𝑢′𝑗
∂𝑥𝑘
̅̅̅̅̅̅̅̅̅̅̅
𝐽𝑖𝑗𝑘 = 𝐽(𝑇)𝑖𝑗𝑘+ 𝐽(𝑃)𝑖𝑗𝑘 + 𝐽(𝑉)𝑖𝑗𝑘 𝐽(𝑇)𝑖𝑗𝑘 = −𝑢′̅̅̅̅̅̅̅̅̅̅̅ , - Velocity fluctuation 𝑖𝑢′𝑗𝑢′𝑘
𝐽(𝑃)𝑖𝑗𝑘 = −1
𝜌(𝑝̅̅̅̅̅̅ ∂′𝑢′𝑖 𝑗𝑘+ 𝑝′𝑢′̅̅̅̅̅̅ ∂𝑗 𝑗𝑘 – Pressure fluctuation 𝐽(𝑉)𝑖𝑗𝑘 = 𝑣∂𝜏𝑖𝑗
∂𝑥𝑘 - Viscosity fluctuation
Due to nonlinearity in the Navier-Stokes equations it becomes necessary to write truncate higher order 𝜏𝑖𝑗 correlations and represent them with lower order- correlations. Expressing 𝜏𝑖𝑗 for higher order correlations results in different turbulence models [3].
𝜀𝑖𝑗, Π𝑖𝑗 and the turbulent diffusion flux 𝐽(𝑇)𝑖𝑗𝑘+ 𝐽(𝑃)𝑖𝑗𝑘 require some form of modeling since they are not expressed with 𝑢̅𝑖, 𝑝̅ and 𝑢′̅̅̅̅̅̅, such a model is referred to as the 𝑖𝑢̅𝑗 Reynolds stress model (RSM) [3].
2.1.2.2 Eddy Viscosity model
In 1877 Boussinesq introduced an equation called the Boussinesq approximation which utilizes a dynamic viscosity coefficient 𝜇𝑇 [3]as seen below;
𝜌𝜏𝑖𝑗 = 2
3∂𝑖𝑗𝜌𝑘 − 2𝜇𝑇𝐷̅𝑖𝑗
The equation in turn is used to derive the eddy-viscosity models Reynold averaged momentum equation below [3]:
6 𝐷̅𝑢̅𝑖
𝐷̅𝑡 = −1 𝜌
∂𝑃̅
∂𝑥𝑗+ ∂
∂𝑥𝑗[2(𝑣 + 𝑣𝑇)𝐷̅𝑖𝑗]
Where 𝑣𝑇 = 𝜇𝑇/𝜌 is the kinematic eddy viscosity. The pressure variable is expressed as follows, with 𝑝̅ as the dynamic pressure:
𝑃̅ = 𝑝̅ +2 3𝜌𝑘
The unknown stress tensor from Navier-Stokes equations has been into a single scalar by relating the Reynolds stress and the average velocity gradient. The eddy viscosity can be modelled with equations of different orders and appearance as seen in the following turbulence models [3].
Two equation models such as the k- ε can model the eddy viscosity using either length(l), time(t) and velocity (q) as follows:
𝑣𝑇 =𝑙2
𝑡 = 𝑙𝑞 = 𝑞2𝑡 2.1.2.3 Standard k-ε model
The eddy viscosity 𝑣𝑇 can be related to the turbulent kinematic energy k and the energy dissipation rate ε as follows.
𝑣
𝑇= 𝐶
𝜇𝑘2𝜀
Where 𝐶𝜇 is a non-dimensional constant, 𝑘
𝜀 can be interpreted as the characteristic time scale for turbulence to be sustained. The k-ε model is based on these scaling arguments and chooses 𝑞 = √𝑘, 𝑡 = k ε⁄ and 𝑙 = 𝑘2/3
Following some steps, the turbulent kinetic energy and the energy dissipation are determined by solving the following:
𝐷̅𝑘
𝐷̅𝑡 = 𝑃𝑘− 𝜀 + ∂
∂𝑥𝑗[(
𝑣𝑇
𝜎𝑘+ 𝑣)∂𝑘
∂𝑥𝑗] 𝐷̅𝜀
𝐷̅𝑡=(𝐶𝜀1𝑃𝑘− 𝐶𝜀2𝜀)𝜀 𝑘+ ∂
∂𝑥𝑗[(
𝑣𝑇
𝜎𝜀+ 𝑣) ∂𝜀
∂𝑥𝑗] For the Standard k-ε model, the constants are the following.
𝐶𝜇 = 0.09, 𝜎𝑘= 1.0, 𝜎𝜀= 1.3,𝐶𝜀1= 1.44, 𝐶𝜀2= 1.92
The k-ε model is widely accepted for many applications, the computational cost is lower by reducing the number of grid points near walls by using wall law. The eddy viscosity is always positive which provides stability to numerical calculations [3]. It is due to the k-ε model’s high convergence rate that it is sometimes recommended to first run a simulation using the k-ε just to get good initial conditions for alternative models that are less stable such as the k-ω model. [5].
7 It is recommended to use the k-ε model when the flows are stationary that have
average streamlines that are relatively straight and when the flow have low levels of acceleration [2]. The k-ε model and eddy-viscosity models in general does not predict the flow behaviour correctly for example when a flow hits a wall or the flow is inside a rotating channel [3] [5].
2.1.2.3.1 Low Reynold number k-ε model
The Low Reynold number k-ε model is based on the k-ε model but does not require wall functions and can solve the flow everywhere. For problems where turbulent quantities are important near the wall regions it is necessary to use non-slip condition at the walls and solve the flow in the wall boundary [3].
The model uses dampening function fµ and the correction functions f1 and f2 for the production and dissipation terms in the ε-equations. Many low Reynolds numbers models has been proposed, but the general form can be written as follows [3];
𝑣𝑇 = 𝐶𝜇𝑓𝜇𝑘2 𝜀 𝐷̅𝑘
𝐷̅𝑡 = 𝑃𝑘− 𝜀 + 𝐷 + ∂
∂𝑥𝑗[(𝑣𝑇
𝜎𝑘 + 𝑣)∂𝑘
∂𝑥𝑗] 𝐷̅𝜀
𝐷̅𝑡 = (𝐶𝜀1𝑓1𝑃𝑘− 𝐶𝜀2𝑓2𝜀)𝜀
𝑘+ 𝐸 + ∂
∂𝑥𝑗[(𝑣𝑇
𝜎𝑘 + 𝑣) ∂𝜀
∂𝑥𝑗]
Low Reynold number k-ε models needs a dense mesh everywhere where its turbulence dampening effect is used. It can model lift and drag forces as well as heat fluxes more accurately than the normal k-ε model [5].
2.1.2.3.2 k-ω model
The k-ω model is based on modelled transport equations for the turbulence kinetic energy k and the specific dissipation rate ω. The model incorporates low-Reynolds- number effects that is applicable to wall bounded flows and free shear flows [4].
There are several variations of the k-ω model [3], the standard form by Wilcox for a constant density is given below:
𝑣𝑇 = 𝑘 𝜔 𝐷̅𝑘
𝐷̅𝑡 = 𝑃𝑘− 𝛽∗𝑘𝜔 + ∂
∂𝑥𝑗[(
𝑣𝑇
𝜎𝑘+ 𝑣)∂𝑘
∂𝑥𝑗] 𝐷̅𝜔
𝐷̅𝑡 = 𝛼𝜔
𝑘𝑃𝑘− 𝛽∗𝑘𝜔2+ ∂
∂𝑥𝑗[(
𝑣𝑇
𝜎𝑘+ 𝑣) ∂𝑘
∂𝑥𝑗] 𝛼 =5
9, 𝛽 = 3
40, 𝛽∗ = 0.09, 𝜎𝑘∗ = 2, 𝜎𝑤 = 2 𝜀 = 𝛽∗𝜔𝑘
8 The k-ω model performs well in predicting boundary layer flows with adverse pressure gradient and separation. It is reported that the model has difficulty predicting spatially developing free shear flows [3].
Since the model is more non-linear that the k-ε model it is harder for numerical solvers to converge and sensitivity to the initial guesses of the solution is increased. Initial guesses could potentially be provided by first running a k-ε model simulation or similar and then use the results as initial guesses for the k-ω model. The k-ω model is useful in internal flows, flows with strong curvature, separated flows and jets [5].
2.1.2.4 Shear stress Transport (SST) model
To deal with the problem of modelling adverse pressure gradients and separation with the k-ω model, the shear stress transport (SST) model was proposed. Which combines the strengths of the k-ε and k-ω models [3]. The model takes advantage of the k-ω model’s accurate formulation near wall regions and the k-ε model’s free stream independence in the free stream. With the addition of three refinements;
namely dampened cross-diffusion derivative term, the definition of the turbulent viscosity affected wall region and by having different model constants [4]. The expression is seen below;
𝐷̅𝜔 𝐷̅𝑡 = 𝛾𝜔
𝑘𝑃𝑘− 𝛽∗𝑘𝜔2+ ∂
∂𝑥𝑗[(
𝑣𝑇
𝜎𝑘+ 𝑣)∂𝜔
∂𝑥𝑗] + 2(1 − 𝐹1)𝜎𝜔21 𝜔
∂𝑘
∂𝑥𝑗
∂𝜔
∂𝑥𝑗
The expression is very similar to that of the k-ω model, with the addition of 𝛾 and a new set of terms to the right. The new terms blend the k-ε and k-ω models, the
function F1 takes values between zero and 1. When F1 = 0, the k-ε part of the model is used, and while F1 = 1 the k-ω model is used. The constants are determined with a weighted average between the constants of the k-ω and k-ε model as a function of F1
[3].
2.1.2.5 Large-eddy simulation (LES)
In Large-eddy simulations large eddies are directly computed while small eddies are modelled, direct numerical solutions in turn require a finer mesh. Therefore, LES requires less memory (RAM) and CPU time than DNS, but more than the RANS models. For LES computations high-performance computing is necessary, especially for industrial applications. LES is believed to be more accurate and reliable than RANS for flows in which large-scale unsteadiness is dominant such as the flow over bluff bodies [4].
Filtered Navier-Stokes equations are commonly used for LES with the filtering width set to be close to the size of the mesh, the filtered parts are then modelled [3].
Close to a solid wall a suitable turbulence model or wall function is needed to eliminate the need for an extremely fine mesh in its vicinity. Examples of these are logarithmic law of the wall, such as the Werner and Wengle wall functions [4].
2.1.3 Wall functions
Flows in contact with solid bodies such as walls are not always simulated accurately by models, or the modelling requires a very fine mesh at these boundaries that increases the computational cost. When the boundary next to the wall is simulated non-slip
9 boundary condition is imposed on the wall surface, the condition states that the velocity at the closest point on the wall has a velocity zero [3].
The logarithmic wall law is a simple wall function that uses the expression below to approximate the velocity tangential to the wall, in a boundary layer starting at a point 𝑦𝑝+ = 30 − 200 grid points normal to the wall.
𝑢̅𝑝+ = 1
Κ𝑙𝑛𝑦𝑝++ 𝐵 Where 𝑦+ ≡ yu𝑇/𝑣 and 𝑢+ ≡ u
u𝑇, with uT being the shear velocity, For a flat plate K = 0.4 and B=5.5. The starting value is taken from the point where the wall function is first imposed in the border to the bulk and is readily solved [3].
2.2 O
ZONE PROPERTIES AND REACTIONSIn the industry ozone’s (O3) oxidative properties is utilized for disinfection, odour reduction and VOC removal, among others. There are different methods of producing ozone, but common to the commercial methods is that a source of oxygen, such as air or pure oxygen gas is used. The oxygen molecules are then excited with corona discharge or UV light (185 nm) which then combines to ozone, this method produce relatively dilute ozone mixtures. Ozone is considered an unstable molecule with an explosion risk at gas concentrations over 30% [6].
Ozone can react directly or indirectly, the direct reaction usually involves an electrophilic attack on unsaturated bonds, while the indirect path involves thermal or UV initiated decomposition to free oxygen atoms that in turn forms hydroxyl radicals starting a chain reaction with radicals [7]. The reaction system involving radicals becomes extremely complicated which becomes apparent by the huge number of kinetic studies compiled in different review papers [8] [9].
2.2.1 Ozone reactions in VOC treatment
The reactions of ozone has been extensively studied [8] [9], most often with the intent of improving the understanding the role of ozone, VOC, moisture and NOx among others gases in the atmosphere. Using this kinetic data, models used to simulate the dynamics in the atmosphere has been made with varying results [10] [11] [12].
Unfortunately, there are still uncertainties about the behaviour of different gas mixtures for a large range of concentrations, pressure and temperatures [8] [9]. It is therefore difficult to make entirely convincing kinetic model involving all reactions in engineering applications.
Searching for exact mechanisms and kinetic data involving organic compounds quickly becomes complicated when considering the radical reaction chains. Reactions take different pathways depending on the location of ozone addition or radical propagation [13] [14] [15].
If the organic compound contains unsaturated bonds such as in alkenes, the ozone molecule can undergo the reactions seen in Figure 1:
10
Figure 1: Addition reaction of an alkene with ozone, picture taken from [15], resulting in different alkoxy radicals
The alkoxy radicals can then continue to react in many ways, reactions include isomerization, radical propagation to other organic or inorganic molecules and termination of the radical propagation [16] [17] [18] [19].
In the cases when Ozone does not directly react with an organic molecule, the thermal decomposition plays a role in the production of free oxygen which in turn forms radicals [9].
In absence of UV, the thermal decomposition of ozone is the following [20] [21]
A compilation of kinetic paRAMeters and reactions used to describe these reactions can be seen in Appendix I - Thermal Decomposition.
The free oxygen atom reacts with moisture to produce hydroxyl radicals [9]
𝑂(1𝐷) + 𝐻2𝑂 → 2𝑂𝐻′
The hydroxyl radical then initiates many reactions, most of which are presented in [22] [8].
𝑂3 ⇄ 𝑂 + 𝑂2 𝑂 + 𝑂3 → 2𝑂2
11
2.3 M
ASS TRANSFER– B
ASIC OVERVIEWIn the simplest case of mass transfer, a molecule is said to diffuse from a region of high concentration to a region of low concentration as described by Fick’s law seen below
𝐽 = −𝐷
𝑧(𝐶1− 𝐶2)
Where J is the molecular flux in units [mol/m2*s], D is the molecules diffusivity in units [m2/s], C stands for the concentration 0f the molecule [mol/m3] in area 1 and 2 respectively and z is the length over which diffusion takes place. The generalized form for a diffusion over length z can be described with the differential [23]:
𝐽 = −𝐷𝑑𝐶 𝑑𝑧 Where J is the flux in unit (mol/m2*s).
To estimate the diffusion coefficient of species 1 in species 2, at a total pressure P in bar and temperature T in Kelvin, the Chapman-Enskog formula below can be used with the following equations [24].
𝐷12= 0.266 ∗ 10−6𝑇3/2 𝑃𝑀𝐴𝑉𝐸1/2𝜎122 Ω MAVE is the average molecular weight as follows
𝑀𝐴𝑉𝐸1/2 = √2𝑀1𝑀2 𝑀1+𝑀2
𝜎12is the characteristic length involved with the two gases 𝜎122 = (0.5(𝜎1+ 𝜎2))2 Ω is the collision integral
Ω = 1.06036
(𝑇∗)0.1561+ 0.193
𝑒0.47635𝑇∗+ 1.76474 𝑒3.89411𝑇∗ T* is then described as follows
𝑇∗= 𝑇
√ε1ε2 𝑘
Constant values for a few gas species are given in Table 1 below.
12
Table 1: Compilation of the calculated mass transfer coefficients for different gas species in air, using the Chapman-Enskog formula at 1 atm total pressure and 295 K. Diffusion coefficient calculated for diffusion in air
Species ε/k σ D12*10-5
[m2/s] (296 K)
Ref Coefficient
values.
Air 78.60 3.711 - [25]
Water 809.10 2.641 2.436 [25]
Oxygen 106.70 3.467 2.051 [25]
Ozone* 233.86 3.709 1.616 [26]
Carbon
monoxide 91.70 3.941 2.027 [25]
Carbon dioxide 195.20 3.941 1.576 [25]
Methane 148.60 3.758 2.219 [25]
Ethane 215.70 4.443 1.494 [25]
Benzene 412.30 5.349 9.414 [25]
2.3.1 Turbulent diffusion
In turbulent flows mixing occurs at macro- and micro scale. Macro scale refer to larger revolving flows (vortices) produced by instability in the average flow, which is characterized by its velocity U and length L. Micro scale vortices are generated by the dissipation of larger vortices which transfers kinetic energy to smaller vortices, this process is known as the energy cascade. The energy cascade describes the process which vortices break up into smaller ones which in turn break up until the smaller vortices dissipates into heat through friction as viscosity becomes an increasing factor [27]. As a measure of energy transferred between vortices the rate of energy
dissipation Ɛ [energy/time] is an important variable for mixing characterization [28]
as it relates to vortex speed and size.
Figure 2: Visualization of the energy cascade process. Larger vortices break down into smaller vortices which finally levels out from friction losses
13 Different from molecular diffusion, mass transport in a turbulent flow is much faster.
The turbulent diffusivity of a passive compound in stationary and homogenous turbulent flow can be described by:
𝜕𝑋̅
𝜕𝑡 = (𝐾 + 𝐷)𝜕2𝑋̅
𝜕𝑥𝑗2
Where 𝜕𝑋̅ is the concentration change space over a change in time 𝜕𝑡, where K is the turbulent diffusivity m2/s, D is the molecular diffusivity m2/s. xj is the direction. K is proportional to the velocity U and length L of vortices. But for practical reasons it is difficult to predict the varying and chaotic nature of real vortices. According to the gradient-diffusion hypothesis the turbulent diffusivity can be described as the ration between the turbulent viscosity and Schmidt number 𝐾 = 𝜈𝑇
𝑆𝐶𝑇 for gases. [29]
2.4 COMSOL M
ULTI-
PHYSICSCOMSOL Multi-physics is a simulation software useful for problems that varies in time and/or space. It mainly uses finite element methods and comes with packages for many types of simulations including fluid dynamics, structural stress and chemical reactions among others. The software can couple several physics, such as fluid flow with diffusion, flow past a membrane for example could involve CFD with mass transfer or reactions and so on. Simulations can be done for 1D, 2D or 3D geometries and commonly as a stationary solution or time depending simulation [30].
For single phase flows COMSOL’s CFD module use Navier Stokes equations, RANS, k- ε, low-Reynolds k- ε, k-ω, SST, v2-f and Spalart-Allmaras turbulence models. All equations in the CFD module can be edited and constants changed to suit modelling needs [31]. By using the Chemical Engineering module diffusing species and reactions can be modelled [32].
COMSOL also has wall treatment in CFD and as of version 5.3 it also includes automatic wall treatment to save computational power. The automatic wall treatment will assign the use of low-Reynolds number k- ε models to accurately simulate the flow behaviour near walls that have a dense mesh, such as in bends and corners. Then for straight walls the logarithmic wall functions can be used. An effect of using this function is that there is no need to use the low-Reynolds number k- ε model at every wall, which requires a dense mesh at these surfaces and subsequently increase computational demand [33].
14
3 M ATERIALS AND METHODS
Several of the CFD models could be used, the k-ε model is widely accepted, memory efficient and robust. But for this project it is interesting to get a clear picture of how the flow around an injector tube of varying lengths submerged in the gas flows. Which makes the low-Reynold number models and especially the SST model interesting. The Large-eddy simulation and DNS approaches are not supported in COMSOL Multi- physics and they would also require an enormous computational effort for the large volumes that needs to be simulated.
In this project many simulations are made using k-ε due to its speed and robustness.
Since COMSOL includes automatic wall functions, these will always be used to reduce the computational effort required.
2D simulations will be used to represent 3D systems in the projects early stages, the much smaller computational effort to solve a 2D system is used to preform many tests at different flow rates and to provide a basic understanding of the system. The conclusions are then used to guide the test cases for the 3D systems.
For the used computers, 16GB of RAM was available. When dealing with 3D
simulations the mesh had to stay below 2.8 million mesh elements for COMSOL to solve the problems on RAM memory for a fully coupled problem when using default iterative solvers. If a pure CFD problem is solved the mesh could be around 3 million elements instead.
3.1 I
NDUSTRIAL DUCTINGThe tests on the industrial ducting are based on the dimensions of a typical industrial sized duct. The details for the model case is seen in Table 2.
Table 2: Shows the dimensions of the industrial ducting and working conditions
PaRAMeter Value
Duct size [mm] 700 x 800 Total air inlet rate [m3/h] 20 000 Duct inlet air speed [m/s] 10 3.1.1 2D Industrial ducting
The industrial duct was represented with the geometry seen in Figure 8 below. The geometry consists of two rectangles merged together with 5 cm radius rounded corners to improve the numerical properties of the system.
In the 2D industrial duct, 4 different cut sections of the duct as seen in Figure 3 were set as the point of extracting the standard deviation of the ozone concentration.
15
Figure 3: Cross sections used to evaluate the standard deviation of the ozone concentration after injection (red colour). Positioned at 1 meter after entrance, 5 cm before the bend (195 cm), 5 cm from the bend and at the outlet
3.1.2 Mesh refinement study
A mesh refinement study was done on the 2D industrial duct by building the mesh in COMSOL with different degrees of resolution. The k-ε model was used for all
simulations to do mesh refinement and all simulations were done with the injector tube at the top of the duct as shown previously in Figure 8
For a quadratic mesh, the expressions in Table 3 was used to control the mesh. The expressions were made from curve fittings in excel based on values from COMSOLs predefined fluid dynamics size settings for “normal”, “fine”, ”finer”, “extra fine” and
“extremely fine” as represented by “refine” values between 1-5. Values above 5 are extrapolations outside the COMSOL predefined settings based on the curve fittings.
Table 3: Shows the expressions put into different fields in the COMSOL mesh (size) settings, where refine is the variable to refine the mesh
PaRAMeter
name Equation Meaning
Curvature -4e-16*refine^3-0.0036*refine^2-0.0036*refine+0.31 Curvature setting value Maxel_size (refine<6)*(-
0.0149*refine+0.0832)+(refine>5)*0.038*exp(- 0.663*refine)
Maximum element size Maxelgrow (refine<7)*(-0.025*refine+1.177)+(refine>6)*1.0015 Maximum
element growth rate minelsize (refine<4)*(-4*10^(-5)*refine^3+0.0005*refine^2-
0.0029*refine+0.0054)+(refine>3)*(0.0127*exp(- 1.111*refine))
Minimum element size Using refine values between 1-5.3 produced the meshes seen in Figure 4 below.
16
Figure 4: Shows a collection of mesh using the refining factor 1,2,3,4,5, 5.1, 5.2 and 5.3 with 1 being the coarsest and 5.3 the finest
As seen from Figure 4, the meshing tool gradually made the mesh elements smaller.
But for refine = 5 mesh near the bend becomes much coarser than the rest of the mesh. Furthermore, for refine 5.2 and 5.3 areas seen in black along the middle of the duct started using trigonal elements instead when it cannot adapt quadratic
elements. In Figure 5 below, results from using the different meshes are shown.
17
Figure 5: Comparison on the result of identical cases using different mesh settings. Standard deviation of the ozone concentration on y-axis and number of mesh elements with refine value on x-axis
Figure 5 show that there is no large difference when “refine” varies between 1-4, but when refine = 5 there is a noticeable change in the result for each iteration. In Figure 6 below the result is presented using percentage units based on deviation from the average from the standard deviation between all tests.
Figure 6: Comparing how much the standard deviation of ozone concentration varies compared to the average value for the standard deviation of ozone concentration.
Figure 6 show that for lower refine value between 1-4, the mesh has a small influence on the solution, for example the deviation from the average only changed by 0.8 % between 3 and 4. When the total mesh elements becomes order of 5 the answer varies more with refined mesh. 0.6 % between 5 and 5.1 and 7.1% between 5.2 and 5.3.
3.1.3 Injector placement
To investigate how ozone distribution varies with injector position the position of a single injector was varied between different heights as seen in Figure 7.
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04
2133 2985 4301 8088 40462 39712 63698 147710
1 2 3 4 5 5.1 5.2 5.3
Standard deviation [mol/m^3]
# Mesh elements (top) and refine value (bottom)
Mesh refinement influence
11.3
9.4 10.2 11.0
8.4 9.0
14.7
21.8
0.0 5.0 10.0 15.0 20.0 25.0
2133 2985 4301 8088 40462 39712 63698 147710
1 2 3 4 5 5.1 5.2 5.3
Absolute deviation from average [%]
# Mesh elements (top) and refine value (bottom)
Mesh refinement influence
18
Figure 7: Qualitative picture, shows how all injector placements are varied from top to bottom with each simulation
Two different geometries are tested as seen in Figure 8 and Figure 9 below
Figure 8: 2D representation of the model system, bulk inflow from the left with the injector represented as a block as seen inside the red borders
19
Figure 9: Straight duct 2D representation, inlet to the left, outlet to the far right, injector is unchanged from the bent duct
3.2 3D I
NDUSTRIAL DUCTINGIn the industrial duct the injection is done 30 cm after the inlet in all cases using the geometry seen in Figure 10 below. The geometry is 2.3 meter long.
Figure 10: Industrial ducting with inflow from the left, outlet to the right
3.3 C
OMMERCIAL DUCTINGFor the commercial ducting simulations, the model case has the following dimensions, seen in Table 4.
Table 4: Initial design paRAMeters for the commercial ducting simulations
PaRAMeter Value
Duct size [mm] 700 x 300
Inlet tubes diameter [mm] 315 Total inlet air flow [L/s] 792-1000
Number of air inlets 3
Air inlet speed [m/s] 3.4
For the final 2D and 3D simulations inlet flow rate of 1000 L/s was used instead, changing the Air inlet speed to 4.3 m/s in each inlet tube.
20 3.3.1 2D Commercial ducting
The 2D ducting was solved with a custom mesh using quadratic elements and COMSOLs extremely fine mesh for fluid dynamics. Tests of different paRAMeter values were done for the cases listed in Table 5
Table 5: Tested values if geometric values in 2D commercial duct simulations
Total flow [L/s]
Inlet tube
width [mm] Main duct
height [mm] Inlet width / duct height ratio
Comment
750 &
1000 315 500 0.63 Preliminary
study 750 &
1000
315 600 0.53 Preliminary
study 750 &
1000
315 700 0.45 Preliminary
study
1000 260 500 0.52 Main study
The points tested in the main study were done with results from preliminary studies as guidance. The injection points are iterated using a paRAMetric sweep in COMSOL.
3.3.2 3D Commercial ducting
In the preliminary simulations the geometry seen in Figure 11 below, was used.
Figure 11: Original ducting geometry, 1-meter straight section followed by 0.5 meter laid as a bend and 1 meter horizontally
It was then decided that a straight duct as seen in Figure 12 gives a more general representation and can use finer mesh within a reasonable computational time.
21
Figure 12: Geometry of the simplified ducting, main duct 2.4 m in total
For the straight duct, several areas were tested. Injections was made from the long side, short side and below at different positions and injector length.
Most of the simulations were successful, but some were skipped due to convergence problems or that the points were not necessary to test anymore.
“Case H” was added to test hypothesis formulated during the development of “mixing predictions models” which try to predict the location of good injection zones based on simulations only using CFD physics.
3.3.3 Inflow distribution in commercial ducting
To predict the path and quantify flow distribution between the inlets when connected to a common atmosphere, CFD simulations with boundary conditions as inlet velocity to outlet pressure or inlet pressure to outlet velocity. For 2D simulations the
geometry seen in Figure 13 below, was used.
Figure 13: Inlet flow distribution simulation geometry, inlet from below, filters at the start of the inlets seen as small squares and outlet to the right of the main duct
22 Calculations are made with the inlet condition as 0 Pa (gauge pressure) or a constant inlet velocity of 3.1 m/s with an outlet condition of either -170 Pa (gauge pressure) or a constant outlet velocity of 4.8 m/s. In the simulations inlet and outlet conditions were velocity in with pressure out and pressure in with velocity out. The simulations were made with k- Ɛ to provide initial conditions for the final simulation using k-ω.
The resistance over the filters were expressed as -abs(k*v^2) [N/m3] in y direction.
Where abs take the absolute value, k [N*s2/m5] is a pressure drop constant and v is the velocity [m/s] in y direction.
Simulated values for the pressure drop constant was -10, -15, -20, -25, -30, -40, -45, - 50, -55 and -60. Resulting in 20 unique solutions, with velocities measured at the end of the inlet tubes and pressure measured over the filters and between inlet and outlet.
3.4 E
RROR ESTIMATIONDuring the 3D simulations some simulations were completely solved for a low resolution (coarse) and high resolution (fine) mesh to find out if the used mesh is good enough. A few cases at different mesh resolution are compared in Table 6 below.
Table 6: Showing change in the standard distribution of ozone concentration at the outlet based on used mesh, H cases (commercial ducting) are solved using symmetry, case H1 full is solve without using symmetry
Case H1 H2 H3 HD H1 Full
# mesh elements 316168 266704 342757 447537 370928
# mesh elements 1176203 1158629 1160117 1160839 1504517 Change outlet standard deviation % 32.08 9.56 1.06 2.35 38.14 The table show that the cases have different mesh reliance and that the result difference is much smaller compared to the increase in mesh elements.
For some of the cases, the position of the injector allows the duct to be cut in half along a symmetry axis. By comparing the results between a simulation using
symmetry to reduce computational effort with a full geometry simulation, numerical errors can be estimated.
3.4.1 CFD based tool for general injection point recommendations Two tools were developed to provide recommendations based on a single CFD simulation.
4 R ESULTS – 2D SIMULATIONS
In the beginning of the project, 2D simulations were carried out for both cases. The industrial case was used to compare results between different models and mesh settings. After which different placements of the injection nozzle were tested in 2D environments to guide development of following 3D simulations.
4.1 2D - I
NDUSTRIAL DUCTINGA flow profile for the ducting is seen in Figure 14 below, the air enters at 10 m/s and separate into different velocity zones near the bend.
23
Figure 14: Velocity magnitude with legend in units [m/s] as read by the colour legend, red arrows show local direction of the flow with black stream lines
Figure 14 show that air velocities can be as high as 18 m/s and as low as 1 m/s near the bend.
4.1.1 2D Injector placement - Industrial ducting
To estimate a good injection depth different positions were tested in a 2D case. The injector position along the height of the duct was varied and ozone distribution was compared between the cases. For the ducting with a bend the change of injection position affects the ozone spread as seen in Figure 15 below.
Figure 15: Shows the standard deviation before the bend and at the outlet of the duct with a bend. Lower standard deviation equals to more uniform ozone distribution
Based on Figure 15 determining which position that gives the best ozone spread depends on which cross section that is analysed. Looking at the ozone distribution
24 before the bend show that the position 630 mm from the top gave the most even distribution while doing the same analysis at the outlet show that the position 630 mm from the top gave the worst distribution in terms of mixing. From the outlet point of view positions between 105 to 560 mm from the top produce similar results where the 560 mm case resulted in the best distribution. Figure 16 below show the concentration profile at different injection depths.
Figure 16: Shows ozone concentration profile with red as high and blue as low concentration respectively, for positions 70, 105, 560 and 630 mm from the top respectively
Concentration profiles seen in Figure 16 enters the low velocity zones at the top corner of the duct and just after the inner bend as seen previously in Figure 14. The used geometry is too small to show where the flow becomes uniform after the bend, which was also seen previous in Figure 14.
Based on results above it was assumed that the optimal position depends on the distance from the injection point to the bend. A straight ducting was simulated as well to investigate the best injection spot in a case where the injector is far away from any bends. For a straight 2-meter duct the position from the top affects the
distribution as seen in Figure 17.
25
Figure 17: Shows the relative difference between the standard deviation of the ozone injection, efficiency values rounded to closest whole value. Calculated as 100*standard deviation (lowest) /standard deviation of current
point
Figure 17 shows that when the injection spot is more than 70 mm away from the top or bottom wall, the ozone spread is near constant. With the best case being 525 mm from the top. Figure 18 below, shows the concentration profile for injection 70, 105 385 and 630 mm from the top
Figure 18: Shows the concentration profile from the injection point at 70, 105 385 and 630 mm from the top respectively, red represent high concentrations and blue is low concentrations of ozone
26 Figure 18 shows that in a 2-meter distance the result is not sensitive to position if the concentration profile does not touch the walls with more than just the edge of the concentration profile.
4.2 2D C
OMMERCIAL DUCTINGThe 2D commercial ducting represent a 3D geometry that has cylindrical inlets connecting to a rectangular main duct. To get an average main duct velocity as expected in the 3D case, the height of the 2D main duct was set to 500 mm as the average between height and width of the 3D duct, then the width of the width of the 2D inlet tubes were set to 150 mm to maintain a correct area ratio between inlets and main duct.
Figure 19 below shows the final 2D representation as described above.
Figure 19: Shows 2D representation of the commercial ducting. Inlet flow from the bottom of the 3 rectangles at the bottom and outlet at the right side of the main duct
Different positions were simulated inside the inlet cylinder and first half of the main duct. Before the 2D geometry was adapted and realistic velocities were set different inflow velocities and inlet width to main duct ratio were used. Velocity profiles from preliminary studies as shown in Figure 20, Figure 21 and Figure 22 below.
27
Figure 20: Velocity magnitude [m/s] shown in colour legend and red flow arrows pointing in the flow direction.
Main duct height is 500 mm and width of inlets 315 mm each
Figure 21: Velocity magnitude [m/s] read using the colour legend and red flow arrows pointing in the flow direction. Main duct height is 700 mm and width of inlets 315 mm each
Figure 22: Velocity magnitude [m/s] read using the colour legend and red flow arrows pointing in the flow direction. Main duct height is 800 mm and width of inlets 315 mm each
In the figures above, it is seen that the change in main duct height to inlet width ratio affects the velocity profile and size of low velocity zones. Although, overall the profiles
28 follow the same shape and position. Injecting in the low velocity zone at the far end opposite the outlet give concentration profiles as seen in Figure 23.
Figure 23: Ozone concentration profile for 500 and 800 mm main duct height with the injection point in the main duct low velocity zone to the left
Figure 23 show that injection inside the low velocity zone opposite the outlet makes the ozone take a path along the top of the duct.
For the case with correct inlet width to main duct height ratio a series of simulations were made. The injection points were graded according to the standard deviation of ozone at the outlet boundary, then a figure was made based on the result that maps out injection zones based on expected mixing efficiency.
The concentration profile of the best case can be seen in Figure 24 below.
Figure 24: Concentration profile of the best case
In Figure 24 the ozone follows the bulk flow from the last duct with the higher ozone concentration entering the middle. In contrast Figure 25 below show the
concentration profile of a less efficient injection spot.
29
Figure 25: Concentration profile for a less efficient injection point
5 R ESULTS - 3D SIMULATIONS
3D simulations produce different results than that of the 2D simulations when the flow has an additional length dimension to flow, leading to revolving flows not seen in the 2D simulations.
5.1 3D I
NDUSTRIAL DUCTINGFor the industrial duct, flow velocity is high, but with simple changes the mixing efficiency was improved drastically.
5.2 3D C
OMMERCIAL DUCTINGSolving the preliminary 3D system using a mesh with “normal” size for fluid dynamics provided Figure 26 and Figure 27 below.
Figure 26: Preliminary commercial ducting. Slices show the velocity according to the colour legend and stream lines originating from the inlets represent flow paths also coloured according to legend
30 The figure shows that the flow entering the main duct undergoes a revolving motion and goes towards the top or along the sides.
In the preliminary case the injection tube was used. Resulting concentration profile can be seen in Figure 27.
Figure 27: Shows the concentration profile as colours, red is high concentration, blue is low concentration
Figure 27 show that ozone gradually spreads to the other side of the duct in a
revolving motion in the first part of the duct. The ozone distribution after the end of the straight section is decent and gets well mixed at the outlet. A new smaller
geometry was made to enable the use of a finer mesh and focus simulation efforts on quick mixing.
The simulations in the straight duct geometry showed that mixing efficiency can vary up to 60 % between cases and that very good mixing can be achieved in a short distance. An example of very good mixing is seen in Figure 28 below.
Figure 28: Best case, concentration profile in slices
31 Figure 28 show that the injected flow has spread out evenly to both sides after the 3 inlets, the same case with a 2.5 cm longer injector is seen in Figure 29 below.
Figure 29: 2.5 cm longer injector compared to the best case, concentration profile in slices
Figure 29 show that unlike the best case in most of the ozone stays on the right side of the duct. Using a 5 cm shorter injector compared to the best case results in the
concentration profile seen in Figure 30 below.
Figure 30: 5 cm shorter injector than best case, slices show the concentration profile of ozone
Figure 30 shows an injection where the ozone stays on the right side in a tight rotation along the right side of the duct.
The second best injection spot gets the ozone distribution seen in Figure 31 below.
32
Figure 31: Second best case, slices show the concentration profile of ozone
In Figure 31 injected ozone spreads out to both sides after injection along the top of the duct, the streamlines separate shortly after the inlet. Figure 28 and Figure 31 both get very good ozone distributions but when looking at the first slice it’s clear the ozone takes different paths over the third inlet.
5.3 I
NJECTION POINT PREDICTION BASED ONCFD
DATAIn a symmetrical duct geometry with a symmetric inlet flow the flow should in theory be mirrored along the middle of the duct. In the commercial geometry the flow should split symmetrically along the middle and other places along the duct length.
By looking at streamlines, attempts can be made to explain the behaviour of injection flows seen earlier in the 3D commercial.
Figure 32 below show the stream lines exiting the first inlet tube and based on the appearance and destination of the streamline red lines has been added to
approximate where the flow change direction and destination.
33
Figure 32: Stream lines seen from above the first inlet. Red lines are added to approximate potential borders where the stream lines change direction and/or destination
The figure shows that the flow splits along the middle (vertical line) and the horizontal lines split flow taking different pathways such as going into the corner, doing a large rotation with slight back flow and rotation towards the outlet without backflow.
Figure 33 below shows the second duct seen from the side, analysed in the same manner.
Figure 33: Stream lines from second inlet duct, side view. Red lines are added to approximate potential borders where the stream lines change direction or destination
The figure shows that there are no large back mixing flows, the stream lines to the right splits into the rotating flow and above the rotation. Figure 34 show the second inlet from above.
34
Figure 34: Second inlet with streamlines, top view. Red lines are added to approximate potential borders where the stream lines change direction or destination
The figure show that the flow splits along the middle (vertical line), the tilted lines are attempts at separating zones that do clockwise and counter clockwise turns.
Figure 35 show the second inlet as seen from the outlet.
Figure 35: Second inlet with streamlines, seen from the outlet. Red lines are added to approximate potential borders where stream lines shift direction or destination
The figure show that the flow separates along the middle (middle vertical line), and the vertical lines next to it separates fields where the streamlines take a path along the top or straight into the rotating bulk flow. Horizontal line separates areas where the streamlines enter the small inner rotation the bottom. Figure 36 below show the second inlet looking towards the outlet.
