Studentarbeten – Mekanik och maritima vetenskaper (M2) – Projektarbeten 2022:05 Flow uniformity characterization in catalytic converters under turbulent inlet conditions TME131 Project in Applied Mechanics 2022 Department of Mechanics and Maritime Sciences C H A L M E R S U N I V E R S I T Y O F T E C H N O L O G Y Göteborg, Sweden 2022 Flow uniformity characterization in catalytic converters under turbulent inlet conditions TME131 Project in Applied Mechanics 2022 © A X E L L A R S S O N , J A C O B L A R S S O N , A R V I N D M U R A L I , E H S A N P E Y V A N D I , S U N I L R A N G A S W A M Y , O S K A R T Y L É N , 2022 Studentarbeten – Mekanik och maritima vetenskaper (M2) – Projektarbeten 2022:05 Department of Mechanics and Maritime Sciences Chalmers University of Technology SE-412 96 Göteborg Sweden Telephone +46 (0)31 772 1000 Preface The work in the present report was carried out as a part of the course TME131 Project in Applied Mechanics, which is a mandatory course within the Applied Mechanics Masters programme at Chalmers. The course was carried out during spring semester 2022. The project was supervised by doctoral student Pratheeba Chanda Nagarajan and Associate Professor Jonas Sjöblom at the Division of Combustion and Propulsion Systems and Professor Henrik Ström at the Division of Fluid Dynamics. The simulations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at Chalmers Centre for Computational Science and Engineering (C3SE) partially funded by the Swedish Research Council through grant agreement no. 2018- 05973. Abstract In the sectors utilizing the combustion of fossil fuels, catalytic converters are commonly employed to reduce dangerous gases generated by combustion engines. For greater ef- ficiency, it is critical to understand the exhaust flow homogeneity inside the converter. This work attempts to evaluate the added value of using Detached Eddy Simulation (DES) simulation over Reynolds-Averaged Navier-Stokes (RANS) simulation for analysis of tur- bulent flow inside catalytic converters. The meshing and simulation of the domain was carried out using ANSYS FLUENT and ANSYS Workbench. Initially a RANS simula- tion with k-ω SST turbulence model was carried out for a specific case to understand the behavior of the time averaged flow field with modeled turbulence and its effect on uni- formity index. The results from the RANS simulation were, together with definitions of turbulent length scales, used to develop a mesh for DES simulation. In this process it was found that in order to resolve a sufficient amount of turbulence upstream of the monolith inlet, the DES simulation required a considerably finer grid than the RANS simulation. The DES simulation with k-ω SST model was used to simulate multiple retention times, attempting to achieve a quasi steady flow with resolved velocity fluctuations. However due to constraints on computational power, simulating until statistical convergence was not possible. The time averaged quantities were extracted from the DES simulation in order to make a fair comparison with RANS results. The uniformity indices of RANS and DES were then compared throughout the monolith of the catalytic converter. A difference of ∼1-1.5% in uniformity index was found from this comparison. Further more it was found that the resolved turbulence in the DES generates fluctuations of ∼3-6% in the uniformity index based on instantaneous quantities. The work concludes that for the exact case simulated, the possibility of increased accuracy when predicting uniformity index with DES instead of RANS does not outweigh the increased computational cost. However for different flow conditions, with larger mean velocity magnitude and less uniform flow, the measured difference between RANS and DES could make a larger difference for prediction of catalytic conversion efficiency. Hence in these cases it is possible that DES simulations are worth their computational cost. Whether this is the case or not would have to be examined through further research before coming to any definitive conclusions. Keywords: Catalytic converters, Flow uniformity, Detached-eddy simulation, Reynolds averaged navier stokes simulations TME131 LIST OF ACRONYMS List of Acronyms CFD Computational Fluid Dynamics DDES Delayed Detached Eddy Simulations DES Detached Eddy Simulation HPC High-Performance Computing LES Large Eddy Simulation RANS Reynolds-Averaged Navier-Stokes SAS Scale Adaptive Simulation SST Shear Stress Transportation UI Uniformity Index URANS Unsteady Reynolds-Averaged Navier-Stokes vi TME131 CONTENTS Contents List of Acronyms vi 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.3 Delimitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.4 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Theory 2 2.1 RANS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.2 LES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.3 DES and DDES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.3.1 Evaluation of spatial discretization . . . . . . . . . . . . . . . . . . 5 2.3.2 Time step and time extent . . . . . . . . . . . . . . . . . . . . . . . 6 2.3.3 Inlet condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.4 Scale Adaptive Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.5 Uniformity Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3 Methodology 8 3.1 Geometry and Boundary conditions . . . . . . . . . . . . . . . . . . . . . . 8 3.2 Simulated conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.3 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.4 Spatial and temporal discretization . . . . . . . . . . . . . . . . . . . . . . 10 3.4.1 RANS mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.4.2 DES mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.4.3 Time stepping and simulation time . . . . . . . . . . . . . . . . . . 12 4 Results 12 4.1 RANS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 4.2 DES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4.2.1 Solution convergence . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.2.2 Mesh evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 vii TME131 CONTENTS 4.3 Uniformity index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.4 Computational Cost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 5 Discussion 20 5.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 5.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 References I viii TME131 1 INTRODUCTION 1 Introduction 1.1 Background In many industrial and non-industrial applications of combustion engines, emissions arise from the chemical reaction of hydrocarbon fuel with air. The common emissions from vehicles include carbon dioxide, carbon monoxide, unburned hydrocarbons, and oxides of nitrogen. These are potential green house and poisonous gases and should preferably be chemically converted to less harmful species before being discharged into atmosphere. These emissions are hazardous to people and environment, as they cause health impacts in the former and affect climate change in the latter. Emissions are also governed by emission legislation. To treat emissions, catalytic converters are employed. These are monoliths with multiple parallel channels, whose walls are coated with catalytic material. The active catalytic component in these devices is usually platinum group metals. The performance of the catalytic converter determines the concentration of emissions in the exhaust pipe. This performance, particularly in terms of catalytic conversion of emissions and degradation rate, is influenced by the radial distribution of the flow through the converter [1]. To improve the performance of the catalytic converter, the information about the flow distribution and temperature distribution are hence vital. To obtain this information at an early stage in the design process, reducing prototyping and testing costs, simulations have to be used. However due to the shapes of the common geometries of catalytic converters and the nature of the concept of uniformity, it is of large interest that the simulation results properly account for the turbulence within the domain. In order to truly capture the turbulent dynamics of the flow, models that resolve turbu- lence need to be used instead of Reynolds-Averaged Navier-Stokes (RANS) models which are a popular choice. The drawback of turbulence resolving models compared to RANS is the significantly larger demand for computational resources. This project looks at the choice between Detached Eddy Simulation (DES) and RANS simulations. This includes the effect of model choice on the resulting flow uniformity, but also whether or not the difference between the models is significant enough to justify the increased computational cost of DES. 1.2 Objectives The main objective of the project is to investigate how resolved turbulent velocity fluc- tuations affect flow uniformity at turbulent inlet conditions and to understand the role of turbulence in attenuating or amplifying non-uniformities using DES turbulence model. This objective is broken down into three smaller components listed below: 1. Determine the differences in uniformity index of mean flow between DES and RANS. 2. Determine the magnitudes of the fluctuations in uniformity index induced by tur- bulence for specified steady state conditions. 3. Evaluate the cost and utility of using DES instead of RANS. 1 of (22) TME131 2 THEORY 1.3 Delimitations As a result of the deadlines within the project, some delimitations were set to demarcate its extent: • Only a single muffler model with catalytic converter is analysed for this project. • The project only considers the RANS and DES solver and the k−ω SST turbulence model. The effect of other turbulence models are not studied for this Muffler- Catalytic converter domain. • Comparison will only be made between DES and RANS results obtained within the project. Due to the difficulty of validating experimental results for complex geometries, no comparison will be made with experimental results. • Due to the computational cost of DES-simulations, the difference between RANS and DES will only be evaluated for a single set of inlet conditions. 1.4 Limitations In addition to the delimitations set internally, some external limitations regarding com- putational power were present. The allocation of computational power of computers and cluster to the project, along with the time available for the project, effected the size of the simulation files that could be handled and the time of simulations. Computers with 32GB RAM, 4 core 3.50GHz processor and 25GB of storage space was available in order to prepare and process simulations. HPC cluster Vera at Chalmers University of Tech- nology was available to run simulations on with a maximum of 256 processor cores. The project received an allocation of 8000 core-hours every month on the cluster. 2 Theory The following chapter covers theory such as underlying equations of fluid flow, and dif- ferent methods to solve them. It also covers flow uniformity and previous work on its importance for catalytic performance. Basic understanding of these topics is of impor- tance for understanding of the methodology and results presented in following chapters. 2.1 RANS The RANS equations are based on the assumption that all variables can be separated into mean (time averaged) and fluctuating components ϕ̄i and ϕ′ i. This assumption gives the Reynolds decomposition of the velocity vector vi = v̄i + v′i (2.1) Inserting the decomposition into the the steady state Navier-Stokes equations yield the continuity equation: ∂vi ∂xi = 0 (2.2) 2 of (22) TME131 2 THEORY and the momentum equations [2]: ρ ∂v̄iv̄j ∂xj = − ∂p̄ ∂xi + ρ ∂ ∂xj ( µ ∂v̄i ∂xj − ρv′iv ′ j ) . (2.3) The Reynolds stresses u′ iu ′ j that have appeared in the last term of the momentum equations have to be modeled to close the equation system. When using the k−ω model, as is done in this project, one uses the Boussinesq assumption v′iv ′ j = −2νts̄ij + 2 3 δijk (2.4) where νt, s̄ij and k are the turbulent viscosity, strain rate tensor and turbulent kinetic energy. To obtain k and νt, transport equations are solved for k and the dissipation ε [2]. 2.2 LES In Large Eddy Simulation (LES), rather than using Reynolds time averaging, we use volume averaging. With the volume average of variable ϕ denoted as ϕ̄ the incompressible Navier-Stokes equations become [2]:∂v̄i ∂t + ∂ ∂xj (v̄iv̄i) = −1 ρ + ν ∂2v̄i ∂xj∂xj − ∂τij ∂xj ∂v̄i ∂xi = 0 (2.5) where τij denote subgrid stresses, representing the turbulence which length scale is small enough to be filtered out by the spatial filter. This turbulence has to be modeled using a turbulence model, and this is done in a similar way as when solving the Reynolds stresses in equation 2.3. In this project, the spatial filtering is achieved using the spatial discretization of the CFD- solver as visualized in figure 2.1. Effectively this means that the computational mesh used in ANSYS Fluent is the filter controlling the amount of resolved turbulence. ? (a) Non-resolved turbulent eddy (b) Resolved turbulent eddy Figure 2.1: Schematic of the effect of cell size on the filtering of turbulent eddies. 2.3 DES and DDES The DES method uses a combination of LES and RANS. The method aims to resolve large eddies with LES but uses RANS to resolve boundary layers. This removes the need 3 of (22) TME131 2 THEORY of mesh refinements near walls that are fine enough to resolve the small scale turbulence of the boundary layer. The idea behind using DES is to combine the best characteristics of RANS with the best characteristics of LES; RANS tends to be able to predict attached flows very well with a low computational cost while LES has a high computational cost for these flows. For separated flows on the other hand, RANS tends to give a poor representation of the physical flow, whereas LES gives good flow prediction without the need for the extreme mesh resolution that would be needed to resolve boundary layers. The k−ω SST modeled was also chosen for the DES simulation. The modeled transport equation for turbulent kinetic energy, k, used in k − ω SST DES modeled reads [2]: ∂k ∂t + ∂ ∂xj (v̄jk) = ∂ ∂xj [( ν + νt σk ) ∂k ∂xj ] + P k − β∗kω (2.6) In order for DES to be able to switch the turbulent length scale from a RANS length scale (∝ k1/2/ω) to a LES length scale (∝ ∆) when the grid is sufficiently fine, the dissipation term in k-equation (last term in equation 2.6) needs to be modified. This can be done by introducing a new blending function FDES, according to (β∗kω −→ β∗kωFDES ) that reads: FDES = max { Lt CDES∆ , 1 } (2.7) where Lt = k1/2 β∗ω , ∆ = max{∆x1,∆x2,∆x3} and β∗, CDES are constants equal to 0.09 and 0.61 respectively. Through this modification, one can make sure that the correct scales of turbulence are being modeled instead of being resolved. It can be seen in equation 2.7 that through ∆, the grid determines the location where the switch between RANS and LES occurs, therefore, an appropriate grid is desirable. In order to ensure that RANS is used near walls it is required that the grid spacing in the wall tangential direction is larger than the boundary layer thickness δ [2]. Otherwise it may happen that FDES term switches to LES in the boundary layer if the tangential grid size (i.e. ∆x and ∆z) are too small. This means that the flow in the boundary layer is treated in LES mode which results in a poorly resolved LES as the mesh is coarse near the wall. For compensating this problem, DDES (delayed DES) is introduced. DDES makes sure that DES stays in RANS mode in the boundary layer by going from FDES to FDDES which is defined as: FDDES = max { Lt CDES∆ (1− Fs), 1 } (2.8) where the blending function Fs is taken as F1 or F2 F1 = tanh(ξ4), ξ = min max { √ k β∗ωy , 500ν y2ω } , 4σω2k CDωy2  (2.9) and 4 of (22) TME131 2 THEORY F2 = tanh(ξ2), ξ = max { 2 √ k β∗ωy , 500ν y2ω } (2.10) For the details and a deeper derivation of different variables in equations 2.9 and 2.10, the reader is referred to literature [2]. 2.3.1 Evaluation of spatial discretization To evaluate the DES mesh resolution, one can look at a number of different properties. These range from length scales and measurements relating to the turbulent energy spec- trum to simple comparisons between properties of the resolved turbulence and properties of the modeled sub grid trubulence. Relative comparisons A method to evaluate how well the mesh is resolving the tur- bulence is by comparing the amount of modeled and resolved turbulent kinetic energy. By statistical time averaging of the DES flow field we can find the resolved turbulent kinetic energy kres = 1 2 v′iv ′ i (2.11) and compare it to the modeled turbulent kinetic energy of equation 2.4 as kratio = kres kres + k (2.12) A common threshold used in evaluation is that kratio should be greater than 80 % [3]. Although the method is popular its accuracy has been questioned [4], hence it should preferably be used cautiously and in conjunction with other evaluation methods. Taylor microscale Apart from the uncertainty surrounding its accuracy, a large draw- back of the methods using the ratio of turbulent kinetic energy is that it requires a DES solution to base conclusions on. For this project, it would be preferable to make a good guess of the required DES mesh size based only on RANS solutions. This can be done as in equation 2.13 by estimating the Taylor microscale λRM i (RM denoting that values can be obtained with RANS methods) using the turbulent kinetic energy. λRM i = √ 15 ν ε ⟨v′i⟩rms (2.13) By using the definition of k from equation 2.11 an equation for λ that uses k and ϵ from RANS can be derived: λRM = √ 10ν k ε (2.14) The aim is to have cell lengths equal or close to λ. This approach places the cutoff of the turbulence resolution on the low energy side of the inertial sub range. It has, although for simpler flows than the one in the recirculation region of the catalytic inlet, been shown to provide good results [5]. 5 of (22) TME131 2 THEORY Using the Taylor microscale as a target for cell size also ensures that the cells are several times smaller than the integral length scale, which is the upper limit of the inertial sub range and is defined as: l0 = k3/2 ε . (2.15) 2.3.2 Time step and time extent The size of the time steps is of importance in order to resolve the fluctuations in the flow. An estimate of sufficient time steps can be based on turbulent time scales or CFL (Courant-Friedrichs-Lewy) number, equation 2.16, given that the mesh is sufficiently re- fined to resolve the turbulence. A CFL < 1 will resolve all turbulence that is possible given the mesh. CFL = |v|∆t ∆ , ∆ denotes cell length (2.16) For implicit time discretization the solver is unconditionally stable independent of time step which alows for higher CFL number but with the risk that the solution is incorrect and turbulence might not be resolved. In a comparison of different CFL numbers on a 2D-hill case resolved with DES it has been shown that an increased time step to CFL = 5 did not result in any major differences for time averaged properties compared to lower CFL numbers [6]. When solving with a transient solver as in DES, efforts have to be made to ensure that a solution is obtained for a long enough physical time to capture the main characteristics of the flow. In this project the concept of retention time, the average time it takes for a fluid particle to move through the entire domain, is used. The retention time tr is defined by the inlet volumetric flow rate and domain volume V as tr = V (vA)inlet (2.17) where v is inlet normal velocity and A is area. Given the fact that equation 2.17 gives the average flow through time. Simulations should be run for at least a few retention times to ensure that a majority of the flow passes through the domain at least once. 2.3.3 Inlet condition With homogeneous inlet condition there is a risk that there are not any disturbances to the flow that cause turbulence and that the flow will stay laminar. The flow that comes to the catalyst is already turbulent and it would be desirable to consider that in specification of the inlet. An effective way of obtaining a developed turbulent flow is to induce perturbations on the inlet. Turbulence generated on the inlet is specified with turbulent intensity and length scale according to recommendations from the Ansys Fluent user guide [7] by assuming developed pipe flow. The turbulent intensity is estimated as I = 0.16Re −1/8 D (2.18) 6 of (22) TME131 2 THEORY where ReD is the Reynolds number based on inlet velocity and inlet pipe diameter. The turbulent length scale for developed channel flow is defined as l = 0.07L C 3/4 µ (2.19) where the characteristic length L, as for the Reynolds number, is the inlet pipe diameter. 2.4 Scale Adaptive Simulation Scale Adaptive Simulation (SAS) models solve unsteady flow with both RANS and LES content in a single model. The SAS concept is based on introducing von Kármán length scale into the equation of turbulence scale. The information provided by the von Kármán length scale allows SAS models to dynamically adjust to resolved structures in an Un- steady Reynolds Averaged Navier Stokes (URANS) simulation, which provides LES-like behaviour in unsteady regions of the flow domain. In addition, the SAS model provides capabilities of standard RANS in the flow regions which are stable [8]. 2.5 Uniformity Index The uniformity index indicates how varied the values for a variable ϕ is over a surface. An index of 1 indicates the highest uniformity where the value is the same for every single cell. As shown by equations 2.20-2.23 [9] the index can be evaluated as area- or mass-weighted where the first captures the quantity variation and the latter captures the flux variation. The area-weighted uniformity index γa is given by γa = 1− Σn i=1[(|ϕi − ϕa|)Ai] 2|ϕa|Σn i=1Ai (2.20) where i is the mesh grid area index for a surface consisting of n such areas and ϕa is the average value for the studied variable over the whole surface. It is calculated as ϕa = Σn i=1ϕiAi Σn i=1Ai . (2.21) The mass-weighted uniformity index differs from the area weighted equation in that it incorporates flux terms: γm = 1− Σn i=1[(|ϕi − ϕm|)(|ρiv⃗iAi|)] 2|ϕ̄m|Σn i=1(|ρiv⃗iAi|) (2.22) where ϕm is the average flux over the whole surface: ϕm = Σn i=1[ϕi(|ρiv⃗iAi|)] Σn i=1[|ρiv⃗iAi|] . (2.23) Previous work [1][10] has shown that the flow distribution through the monolith of a catalytic converter plays a crucial role in the catalytic performance and degradation rate. In addition to the fact that uniformity impacts performance, Guojiang [10] has shown that the impact of changes to uniformity differs at different conditions. As flow velocity increases and uniformity decreases, smaller differences in γa have a larger impact on catalytic conversion efficiency. 7 of (22) TME131 3 METHODOLOGY 3 Methodology This chapter describes the methods used in order to obtain comparable results from RANS and DES solutions. The main strategy of the process was to first obtain RANS results for the geometry and conditions in question, and then use these results to make sure that the DES results would be as reliable as possible. 3.1 Geometry and Boundary conditions The simulations will be carried out with a model that consists of in inlet pipe followed by a circulation area and a monolith modeled as a porous medium. These regions, together with the velocity inlet at which the conditions were changed and the pressure outlet can be seen in figure 3.1. Figure 3.1: Geometry used for simulations. The velocity inlet is marked in blue and the pressure outlet is marked in red, the green part is the monolith. The inlet of the monolith is located at z=0.02355m and the outlet is located at z=0.1735m. The z-axis is the axis parallel to and positive in the direction of the flow inside the mono- lith. The x-axis and y-axis are then the aligned to the width and height of the monolith respectively. The monolith is represented as a porous media instead of modeling it exactly which would generate a very detailed geometry. The porous media is a region with low permeability, but relatively higher in the axial direction compared to the other. This achieves the effect of the monolith with a laminar flow in the axial direction. The cordierite is used as the material for porous media. This material is used to convert pollutants like CO, HC and NOx into less dangerous gas. The walls of the catalyst are coated with layers of insulating mat, glass wool and steel. The insulating mat gives protection against the mechanical shock and vibration, the glasswool is an insulating material and the steel provide protection to the catalytic converter. The walls of domain before and after the monolith are defined as steel with a specific convection coefficient. This wall boundary condition helps to replicate the exact model of catalytic converter used in automotive industry. 8 of (22) TME131 3 METHODOLOGY 3.2 Simulated conditions The boundary condition used represents realistic steady operating conditions in terms of velocity and temperature in a catalytic converter. The inlet condition together with the turbulent properties, calculated as in section 2.3.3, are defined in table 3.1. The turbulent properties are only active for DES. Table 3.1: Boundary condition at inlet T [K] V [m/s] I [%] l [m] 423 25 4 0.01 The surrounding domain and the pressure outlet are set to atmospheric conditions with a temperature of 293K. In addition to to the conditions specified in table 3.1, an additional set of inlet conditions was run using a RANS solver only which are defined in table 3.2. Table 3.2: Additional RANS conditions used to understand the relevance of the magnitude of uniformity changes T [K] V [m/s] 423 75 By examining the effect of inlet conditions on uniformity index, the group was able to get an idea of the importance of the magnitude of the uniformity differences between RANS and DES simulations as well as the value of considering the impact of the chosen inlet conditions on the results. 3.3 Comparison To ensure the fairest possible comparison between results from the time averaged RANS solution and the inherently transient DES solution, time averaged quantities were ex- tracted from the DES solution. This was automated using the built in functionality for statistical time sampling within ANSYS Fluent. Uniformity indices defined by equations 2.20-2.21 were calculated using average quantities as ϕ. The built in definitions for mass- weighted UI use the instantaneous density and velocity to calculate mass flux, causing the indices of average quantities to be effected by instantaneous data. Because of this, use of mass-weighted UI was avoided, and the area-weighted formulation was the primary choice. As it is within the monolith of the catalytic converter that the uniformity index is of interest when evaluating performance, the comparison of UI was made on seven cross sections within the monolith. In addition to analyses of uniformity indices, contours of the examined quantities were also analyzed in order to understand differences between RANS and DES results. This also aids the understanding of how the uniformity is changing in time when evaluating the DES simulations. 9 of (22) TME131 3 METHODOLOGY 3.4 Spatial and temporal discretization For spatial discretization of the domain, a combination of tetrahedral and quadrilateral meshing was used. Tetrahedral cells are used in the domain between the inlet and monolith due to its ability to adapt to changing geometries. Quadrilateral cells are used in and after the monolith where the geometry is simple and the flow is predominantly laminar. All walls are refined with prism layers in order to accurately calculate the wall shear stresses and provide a smooth transition of the cell volume. 3.4.1 RANS mesh The computational mesh used for the RANS simulations was selected through a mesh convergence study. Four meshes with a cell count between 0.7 to 2.2 million were com- pared based on uniformity index of velocity at the monolith outlet. Identical boundary conditions and convergence criteria were maintained. Figure 3.2: Uniformity index of mass weighted velocity at monolith outlet obtained from RANS simulations on meshes of different sizes As seen in figure 3.2, relative convergence of the results seems to have been achieved when passing a cell count of about 1.2 million cells. As a result of this, the 1.2 million cell mesh used in the mesh convergence study was used for the RANS simulations. 3.4.2 DES mesh Using equations 2.14 and 2.15, it could be determined whether or not the mesh was also suitable for DES simulations. The results on the yz-plane in the center of the domain is shown in figure 3.3. 10 of (22) TME131 3 METHODOLOGY (a) l0 ∆ (Integral length scale l0) (b) λ ∆ (Taylor length scale λ) Figure 3.3: Turbulent length scales compared with cell length ∆ from the mesh used for the RANS simulation. As made evident by figure 3.3b, using the RANS mesh of 1.2 million cells for the DES simulation would result in a solution where a majority of the cells in the important re- circulation region after the inlet pipe would not be able to resolve turbulence at the Taylor length scale. In figure 3.3a, one can also see that in critical regions of the domain not even the largest eddies would be resolved properly. To deal with this, the mesh from the RANS simulation was refined based on the integral and Taylor length scales. After iterative refinement, this resulted in the length scale resolution shown in figure 3.4, with a mesh count of about 6 million cells. Mesh refinement is specifically done in the non boundary layer parts of the inlet pipe and flow re-circulation region. To obtain the smallest sufficient cell count, mesh in the other regions is not refined. (a) l0 ∆ (Integral length scale l0) (b) λ ∆ (Taylor length scale λ) Figure 3.4: Turbulent length scales compared with cell length ∆ for the refined mesh. After the refinement the cells in the inlet pipe and circulation region were slightly above or below the size of Taylor micro-scale except at the walls where RANS is supposed to be used. The cells are also at least five times smaller than the integral length in large parts of the re-circulation region and 2-3 times smaller in the beginning of the separated shear layer, allowing for resolution of the largest eddies in the energy spectrum. 11 of (22) TME131 4 RESULTS 3.4.3 Time stepping and simulation time After refining the mesh using turbulent length scales, the time step has to be set to give a good CFL number together with the new mesh and expected velocity magnitudes. Figure 3.5: CFL number based on the mesh defined for DES, the velocity magnitudes from the steady state RANS results and timestep ∆t = 0.0001. In figure 3.5, we see that with a time step of ∆t = 0.0001 the CFL number does not go above 5. Given the previous work on CFL effects on mean quantities [6], this time step is considered fine enough for the task at hand, and it is chosen for the DES simulations. For the transient DES simulations, the steady state RANS results were used as initial conditions. However they were slightly altered to include synthetic turbulent fluctuations generated by ANSYS Fluent. In order to allow for the solution of physically correct turbulent fluctuations, the solver was run for two retention times (see equation 2.17) before initializing the statistical time sampling used for comparison with RANS results. For the specified case, two retention times equate to about 0.25 seconds. 4 Results The results obtained from the process described in the previous section are here presented in the following order; RANS results, DES results, model comparison. 4.1 RANS The time averaged results from RANS are studied by observing the velocity and temper- ature contours at the cross sections where UI is measured. The progression of the velocity and temperature field inside the catalyst can be seen in figures 4.1 and 4.3. 12 of (22) TME131 4 RESULTS (a) Catalyst inlet, z=0.05m. (b) z=0.1m. (c) Catalyst outlet, z=0.1735m. Figure 4.1: Velocity throughout the catalyst. The velocity field is very similar in all cross sections which is expected due to continuity since the flow is predominantly laminar and only travailing in the axial direction trough the porous media. The velocity field inside the monolith is characterized by the high velocity towards the bottom and the low velocity in the center of the upper half. It is intuitive that there will be higher velocity towards the bottom due to momentum of the flow from the inlet that comes from above. 29.16 21.87 14.58 0.00 7.29 ms-1 Figure 4.2: Vector plot of flow velocity in the recirculation region. The color scale of the vectors represent velocity magnitude. The re-circulation and chaotic motion of flow before entering the monolith is the main 13 of (22) TME131 4 RESULTS reason for the non uniform distribution of velocity. The vector contour in figure 4.2 shows the re-circulation happening before the monolith. From the vector contour, it is clear that a majority of the flow enters the region with high velocity and impinges on the wall to recirculate. The smaller concentration of flow in the upper part of the catalytic converter can be explained by the fact that the recirculating flow is blocked by the jet flow entering the re-circulation region from the inlet. The temperature, figure 4.3, is effected more than the velocity throughout the catalyst since the change in temperature is time dependent and the flow at the outlet has been exposed to the far field temperature diffusion for a longer time than the flow at the inlet. (a) Catalyst inlet, z=0.05m. (b) z=0.1m. (c) Catalyst outlet, z=0.1735m. Figure 4.3: Static temperature throughout the catalyst. The temperature decreases in areas close to the wall which has a lower temperature than the flow throughout the domain. The temperature also decreases more in areas with low velocity and thus has more time to cool down when travailing between cross sections. The temperature field hence looks as expected based on the velocity distribution. The RANS simulation gives an insight of the time averaged flow inside the monolith. The effect of instantaneous variation of uniformity index is still unknown. Therefore it is one of the reasons for analysing this project with a DES modelling approach. 14 of (22) TME131 4 RESULTS 4.2 DES The results from DES is, unlike those from RANS, presented with a time history. The instantaneous solution are used in order to evaluate fluctuations and the time averaged statistics allows for later comparison with RANS. Figure 4.4: Area-Wighted uniformity index of mean and instantaneous velocity in the middle of the catalyst (z = 0.1). t = 0 is the time at which statistical sampling was initialized. The solution of the DES consists of fluctuations with different amplitude and frequency. In figure 4.4 there are high frequency fluctuations with an amplitude of ∼ 0.03 in UI and there appears to be lower frequency fluctuations with an amplitude of ∼ 0.06. Contours of velocity inside the catalyst in figure 4.5 are captured at 0.36s and 0.48s with a high value of UI and at 0.40s and 0.55s with a low value of UI after the high amplitude drops that can be seen in figure 4.4. 15 of (22) TME131 4 RESULTS (a) t = 0.36s (b) t = 0.40s (c) t = 0.48s (d) t = 0.55s Figure 4.5: Velocity in the middle of the catalyst (z=0.1m) at different times. The right hand side contours are captured right after the drops in uniformity seen in figure 4.4. The left hand side contours are captured just before the same drops. The flow is distributed with a region of high velocity at the bottom of the catalyst and lower velocity in the rest of the catalyst. In the DES simulation the region of high velocity changes in size and moves slightly right and left with time. The change of the high velocity region correlates with the change of UI, as the region gets more prominent the UI decreases. It can also be seen that a thin curve of slightly higher velocity flow is present in an arc from this high velocity region to the upper part of the monolith. This curve is always located on the same side of the monolith as the high velocity region at the bottom (in figure 4.5 (d) the high velocity region is present on both sides so two such arcs are seen). The rest of the cross section does not change with a clear pattern. The time averaged velocity is studied in figure 4.6 in order to compare with RANS. Figure 4.6: Mean velocity in the middle of the catalyst (z=0.1m) at the last time step. 16 of (22) TME131 4 RESULTS The time averaged result from the DES also has the region of high velocity at the bottom wall. There is also a distinct region in the middle of the upper half where the velocity is low. Less distinct but still noticeable low velocity regions are also found on just below the middle on the right and the left side of the monolith. All these low velocity regions are the mean result of the oscillating thin arcs seen in the instantaneous results in figure 4.5. These are characteristics that could also be identified in the RANS results in figure 4.1b. One interesting difference between the result for DES in figure 4.6 and for RANS in figure 4.1 (b) is that the high velocity region at the bottom has a nearly vertical gradient while the DES gives a more crescent shaped velocity gradient that is also biased towards the left. The uniformity index of temperature, figure 4.7, behaves in a different way compared to the velocity. Figure 4.7: Area-Wighted uniformity index of mean and instantaneous temperature in the middle of the catalyst. t = 0 is the time at which statistical sampling was initialized. The changes of UI are small in absolute values, it is however clear that the changes are caused by the solution changing at a larger time scale than the simulated time rather than resolved fluctuations. (a) t = 0.00s (b) t = 0.70s Figure 4.8: Temperature in the middle of the catalyst (z=0.1m) at first and last time instance of figure 4.7. The temperature distribution in figure 4.8 is visually identical at the first and last time instance of figure 4.7. This is thus the reason of the the low changes in UI. The contour of 17 of (22) TME131 4 RESULTS mean temperature is also visually identical to the instantaneous and therefore not shown. 4.2.1 Solution convergence The simulated flow time was decided based on including several retention times and limited by allocation of cluster access. The time history of flow quantities can be studied in order to analyze whether the solution has statistically converged or not. The quantities of interest in this report are velocity and temperature of which the history in a point roughly in the center of the monolith is presented in figure 4.9a & 4.9b. (a) Velocity history (b) Temperature history Figure 4.9: Time history in point x = y = 0, z = 0.1. The time history of velocity and temperature indicates that the mean statistics has not converged during the simulated time. The fluctuations of velocity are however resolved by the simulated time while the fluctuations of temperature appears to not have been resolved if there are any. The velocity has an upward trend both in figure 4.9a and in figure 4.4 the velocity UI also has the same upward trend which indicates that the solution is not perfectly statistically converged for the velocity. Both 4.9 and 4.7 shows that the temperature has not statistically converged during the simulation time. Due to the small changes of both instantaneous and mean temperature it is difficult to tell if the cause is fluctuations with longer time scale than the simulated time or that the solution is still adapting from incorrect initial conditions. 4.2.2 Mesh evaluation The results of the DES is evaluated with the methods described in section 2.3.1. Figure 4.10 shows that the ratio of resolved turbulent kinetic energy kratio is larger than 80 % everywhere except in the boundary layer of the inlet pipe and parts of the wall in the re-circulation region below the inlet pipe and on the inlet of the monolith. 18 of (22) TME131 4 RESULTS Figure 4.10: Volume rendering of kratio (see equation 2.12). 4.3 Uniformity index By use of the statistical mean quantities obtained from the DES solver, comparison could be made against the uniformity indices from the RANS simulation as done in figure 4.11. Figure 4.11: Area-Weighted uniformity index of mean velocity across the seven cross- sections within the monolith. Comparison between RANS and DES results. The result shows that in all monolith cross-sections, the DES simulations gave a more uniform flow than the RANS simulations. The difference in terms of UI is between ∼0.010 and ∼0.015 in all cross sections. However as one can note in figure 4.4, there could be some oscillation present at timescales larger than the temporal extent of the performed simulation. From the same figure we also see that the regular oscillations (∼0.03) are at least two times bigger in size than the difference between the averaged DES and RANS results (∼0.015) for uniformity index. The difference in UI between RANS and DES can be put in perspective by analysing the difference in UI based on inlet condition. In figure 4.12 the RANS result is compared with an additional RANS simulation with different inlet conditions, table 3.2. 19 of (22) TME131 5 DISCUSSION Figure 4.12: Area-Weighted uniformity index of mean velocity across the seven cross- sections within the monolith. Comparison of RANS results at different inlet velocities. Figure 4.12 shows the effect of increased inlet velocity on the area weighted uniformity index within the monolith. A larger velocity gives a lower uniformity index, for this case at a difference of ∼0.1. It is clear that the inlet conditions have a larger effect on the uniformity index than the choice between DES and RANS. 4.4 Computational Cost As the change between RANS and DES was made, the changes to the mesh presented in section 3.4.2 resulted in a considerably larger equation system for the CFD solver to solve in each iteration. In addition to this, the transience of the DES simulations together with the small time step requirements meant that more iterations had to be run for the DES case. For the RANS simulation, approximately 4000 iterations were run to reach convergence for the case in question. In comparison, the DES simulation required approximately 230 000 iterations, and should have been run further for better statistical convergence. Together with the increased mesh size, this caused a significant increase in core hours used. The RANS simulation required 26 core hours to converge whilst the results obtained with DES required 15000 core hours, an increase of ∼50 000 %. 5 Discussion In section 1.2, three main objectives were set. After analyses of the simulation results obtained within the project they can be evaluated: 1. Determine the differences in uniformity index for mean flow between DES and RANS. From the results presented in figure 4.11, it seems that running DES rather than RANS does result in a difference in obtained flow uniformity through the monolith of the catalytic converter. However the change of UI is not of a highly significant magnitude. As shown by Guojiang [10], changes of this magnitude to a flow distri- bution with this high uniformity will only have very minor effects on the catalytic conversion efficiency of some gases. 20 of (22) TME131 5 DISCUSSION On the other hand, the same work implies that the impact of small uniformity changes are larger as the uniformity goes down and velocity increases. Combining this knowledge with the results presented in figure 4.12, it becomes possible that the magnitude of UI change between RANS and DES could become very significant at higher inlet velocities. 2. Determine the magnitudes of the fluctuations in uniformity index in- duced by turbulence for specified steady state condition. For the simulated case, the magnitude of the regular oscillations of the uniformity index for velocity seems to be of size ∼0.03 with irregular jumps at size ∼0.06. Similarly to the difference of mean flow UI between DES and RANS, this is rather small in comparison with the measured mean uniformity index at 0.90-0.92. How- ever as stated beforehand the impact of small uniformity changes are larger when the uniformity goes down. Therefore even these fluctuating differences (especially the larger irregular ones) could become significant to the instantaneous catalytic conversion efficiency at higher inlet velocities. 3. Evaluate the cost and utility of using DES instead of RANS. As reported in section 4.4, the computational cost of running DES simulations at a sufficiently fine resolution in time and space is significantly larger in comparison to the cost of running steady state RANS simulations. Hence one has to carefully evaluate the possible added value of DES simulations to their analyses of catalytic flow uniformity before going down that road. It is unlikely that sufficient value to make up for the additional cost will be added unless the flow is expected to be highly non-uniform. On top of that, high non- uniformity is likely to be a result of higher inlet and mean velocities within the domain. In turn, this will have an impact on the turbulent length scales and CFL numbers presented in sections 2.3.1 and 2.3.2 which could further increase the com- putational cost of simulations of sufficient resolution. Furthermore the time statistics in the DES solution are not fully converged. This is because the mean values are still changing and there might be lower frequency fluctuations which has not been captured. The left bias for the high velocity region in figure 4.6 is another telling sign that the simulation has not run for a long enough physical time. This is because nothing in the geometry or the inlet velocities justify such a bias and any turbulence bias should have been averaged out to zero if the simulation has reached a statistical quasi-steady state. One question that remains completely unanswered by this project is whether DES is the best alternative to RANS. It is possible that other options for scale resolving are better for the task at hand. One option for this could be the scale adaptive simulations briefly mentioned in section 2.4, as the possibility to transition between RANS and LES in the recirculation region could possibly reduce the computational cost when compared with DES. 5.1 Conclusions From the results we see that the DES results showed a different uniformity index for velocity than the RANS results. However, the difference was minor in comparison with 21 of (22) TME131 5 DISCUSSION the difference a change in inlet velocity would cause in RANS simulations. The difference was also minor in the effect it will have on the catalytic conversion efficiency. Similarly the fluctuations in the velocity uniformity index in the DES simulations were minor in comparison with the mean uniformity index. The cost of DES compared to RANS also turned out to be significantly higher. The computational time to run was several orders of magnitude higher for DES than for RANS and the ratio is estimated to be even worse for higher inlet velocities. It was also unclear if the simulated time was enough for a converged result in the DES simulation. Based on these findings the DES results cannot be said to be a better alternative to RANS for time averaged values for the studied inlet boundary condition, as the cons of computational cost seem to weigh heavier than the slight possible increase in accuracy. However the DES results did show that the time dependent flow field and the fluctuations were differentiable in comparison with the mean RANS results. DES might hence be worth the high computational cost if the effect of these fluctuations are important. This could for example be when simulating different operating conditions or when designing new catalytic converters. For most applications however the cost is prohibitive for what likely will be a marginal increase in accuracy. 5.2 Future work The results presented in this report concerns a single case with one velocity and one temperature inlet boundary value. A better comparison between RANS and DES could be made if more cases with different velocities and temperatures were simulated. Partly due to time constraints and partly due to the increasingly refined meshes needed to resolve the turbulence for increased flow velocities as given by the length scale in equation 2.13 this was not done for this report. For deeper analyses and understanding of the case simulated in this project, the DES simulation should also be run until a statistically converged flow had been reached. To have in depth comparison of the uniformity index between DES and RANS results it would be recommended to run simulations for a number of cases. The results of the present work do however indicate that if more cases are run, maximum value will likely be achieved if their inlet velocity conditions are much higher. This would allow for analysis of cases where any difference in uniformity index would be of larger importance. Another area to look into is the use of other turbulence resolving approaches, such as SAS which performs well where high separations are present in the flow[11]. This is relevant in this case since separation occurs when the pipe flow enters the circulation region before the monolith. 22 of (22) References [1] John S Howitt and TC Sekella. “Flow effects in monolithic honeycomb automotive catalytic converters”. In: SAE Transactions (1974), pp. 1067–1075. [2] Lars Davidson. “Fluid mechanics, turbulent flow and turbulence modeling”. Mar. 2022. [3] Stephen B Pope. “Ten questions concerning the large-eddy simulation of turbulent flows”. In: New journal of Physics 6.1 (2004), p. 35. [4] Lars Davidson. “Large eddy simulations: how to evaluate resolution”. In: Interna- tional Journal of Heat and Fluid Flow 30.5 (2009), pp. 1016–1025. [5] Yacine Addad et al. “Optimal Unstructured Meshing for Large Eddy Simulations”. In:Quality and Reliability of Large-Eddy Simulations. Ed. by Johan Meyers, Bernard J. Geurts, and Pierre Sagaut. Dordrecht: Springer Netherlands, 2008, pp. 93–103. isbn: 978-1-4020-8578-9. doi: 10.1007/978- 1- 4020- 8578- 9_8. url: https: //doi.org/10.1007/978-1-4020-8578-9_8. [6] Michael Kornhaas, Dörte C. Sternel, and Michael Schäfer. “Influence of Time Step Size and Convergence Criteria on Large Eddy Simulations with Implicit Time Dis- cretization”. In: Quality and Reliability of Large-Eddy Simulations. Ed. by Johan Meyers, Bernard J. Geurts, and Pierre Sagaut. Dordrecht: Springer Netherlands, 2008, pp. 119–130. isbn: 978-1-4020-8578-9. doi: 10.1007/978-1-4020-8578- 9_10. url: https://doi.org/10.1007/978-1-4020-8578-9_10. [7] Ansys® Fluent User Guide. Version 2021 R1, January 2021. [8] T. Belamric M. Younsi A. Djerradab and F. Menterd. “Application of the SAS tur- bulence model to predict the unsteady flow field behaviour in a forward centrifugal fan”. In: International Journal of Computational Fluid Dynamics (2008). [9] Ansys® Fluent Theory Guide. Version 2021 R1, January 2021. [10] Wu Guojiang and Tan Song. “CFD simulation of the effect of upstream flow distri- bution on the light-off performance of a catalytic converter”. In: Energy Conversion and Management 46.13-14 (2005), pp. 2010–2031. [11] Weilin Zheng et al. “Comparative assessment of SAS and DES turbulence modeling for massively separated flows”. In: Acta Mechanica Sinica 32.1 (2016), pp. 12–21. I https://doi.org/10.1007/978-1-4020-8578-9_8 https://doi.org/10.1007/978-1-4020-8578-9_8 https://doi.org/10.1007/978-1-4020-8578-9_8 https://doi.org/10.1007/978-1-4020-8578-9_10 https://doi.org/10.1007/978-1-4020-8578-9_10 https://doi.org/10.1007/978-1-4020-8578-9_10 List of Acronyms Introduction Background Objectives Delimitations Limitations Theory RANS LES DES and DDES Evaluation of spatial discretization Time step and time extent Inlet condition Scale Adaptive Simulation Uniformity Index Methodology Geometry and Boundary conditions Simulated conditions Comparison Spatial and temporal discretization RANS mesh DES mesh Time stepping and simulation time Results RANS DES Solution convergence Mesh evaluation Uniformity index Computational Cost Discussion Conclusions Future work References