CFD-Modeling of fluid domains with embedded monoliths with emphasis on automotive converters

Journal Pre-proof CFD-Modeling of Fluid Domains with Embedded Monoliths with Emphasis on Automotive Converters Matthias Hettel, Eric Daymo, Tobias Schmidt, Olaf Deutschmann

PII:

S0255-2701(19)30991-2

DOI:

https://doi.org/10.1016/j.cep.2019.107728

Reference:

CEP 107728

To appear in:

Chemical Engineering and Processing - Process Intensification

Received Date:

9 August 2019

Revised Date:

9 October 2019

Accepted Date:

6 November 2019

Please cite this article as: Hettel M, Daymo E, Schmidt T, Deutschmann O, CFD-Modeling of Fluid Domains with Embedded Monoliths with Emphasis on Automotive Converters, Chemical Engineering and Processing - Process Intensification (2019), doi: https://doi.org/10.1016/j.cep.2019.107728

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

CFD-Modeling of Fluid Domains with Embedded Monoliths with Emphasis on Automotive Converters

Matthias Hettel 1, Eric Daymo2, Tobias Schmidt1, Olaf Deutschmann1 1

Institute for Chemical Technology and Polymer Chemistry, Karlsruhe Institute of Technology (KIT), Kaiserstr. 12, 76128 Karlsruhe, Germany Tonkomo, LLC, Gilbert, Arizona, 85297 United States

ro of

2

lP

re

-p

Corresponding Author: Matthias Hettel, Institute for Chemical Technology and Polymer Chemistry, Karlsruhe Institute of Technology (KIT), 76128 Karlsruhe, Germany, Tel.: (+49) 721 608 44269, Fax: (+49) 721 608 44805, E-mail: [email protected]

Jo

ur

na

Graphical Abstract

Highlights -1-

  

New approach for calculating flow domains with embedded monoliths Domain coupling allows exclusion of monolith and usage of different models in domains Approach works perfectly and leads to valuable results RANS modeling leads to good results but transition behind monolith cannot be predicted LES yields all details of flow and captures transition from laminar to turbulent flow

ro of

 

ABSTRACT

ur

na

lP

re

-p

A new approach for calculating flow domains with embedded monoliths is presented, focusing on a system suitable for an automotive converter. Using appropriate boundary conditions, the monolith itself can be excluded from the CFD computational domain, leaving two mapped computational domains upstream and downstream of the monolith. The resulting method enables more detailed results of the downstream flow profile than afforded by the commonly used porous body approach for simulating a monolith, but without the complexities of a full 3D calculation of the entire domain, including the monolith. The present approach was validated with experimental flow data from the literature collected with a prismatic (planar) monolith with approximately 4500 channels. Sensitivity studies were performed to check the influence of the downstream turbulence model, inlet turbulence boundary conditions, and spatial discretization schemes. RANS turbulence models (k- SST and k-) generally predict the experimentally measured flow profile downstream of the monolith although the transition to turbulence cannot be reproduced correctly. Currently, LES is the best approach for adequately describing the characteristics of flow downstream of monoliths.

Jo

Keywords

CFD; OpenFOAM; Catalytic Converter; Monolith Catalyst; Honeycomb Catalyst; Automotive Converter

1 Introduction The application of Computational fluid dynamics (CFD) simulations helps to analyze the interdependency between mass and heat transport and chemical reactions in reactors for chemical and materials synthesis, in particular in the presence of catalytic surfaces. In many -2-

technical processes, e.g. high-temperature catalysis, catalytic combustion, reforming and aftergas treatment, catalytically coated honeycomb monoliths are used. However, even the prediction of flow fields without additional chemical processes is a challenge, particularly for turbulent flows. The focus of the work presented here is on the calculation of the flow field typical of catalyst assemblies used in automotive applications (Fig. 1-1) rather than chemical conversion or mass and heat transport that occurs within the monolith.

Jo

ur

na

lP

re

-p

ro of

The prediction of the flow field is important for several reasons. Firstly, it is highly desirable that the velocity of flow is the same in all channels in order to achieve uniform chemical conversion and utilization of the catalyst (Guojiang & Song, 2005; Holmgren et al., 1997; Ibrahim et al., 2018; S.-J. Jeong, 2014; Karvounis & Assanis, 1993; Shi-jin et al., 2000; Zhang & Tennison, 2008). Flow maldistribution leads to an inhomogeneous ageing of the catalyst (Martin et al., 1998) and influences the light-off characteristics (Chakravarthy et al., 2003; Guojiang & Song, 2005; Su et al., 2013; Zygourakis, 1989). Secondly, probes which measure the air fuel ratio (lambda value) or the OSC (Oxygen Storage Capacity) of the catalyst are used in three-way converters as part of the control system for the motor and the exhaust gas aftertreatment system. The probes should measure data which are representative for the overall status of the system. Therefore, the optimal positioning of the probes in front, inside or downstream of a monolith, or even in between two monoliths, is essential (Hwang et al., 1995; Weltens et al., 1993).

Fig. 1-1:

Schematic representation of the flow in a typical catalyst assembly

Reactors that are based on monolithic structures are comprised of many equally sized small channels. To resolve these structures with CFD methods in detail, grids with high spatial resolution are needed, which lead to large numbers of grid cells, requiring computational times and storage capacities which are not viable. To reduce the computer resources needed for modeling of structured catalytic reactors several strategies with various levels of -3-

complexity and accuracy have been developed. These strategies will be briefly explored, emphasizing CFD of automotive converters.

ro of

The most common approach to limit the number of grid cells is to model the monolith as a continuous media, which is partially fluid and partially solid (Benjamin S. F. et al., 2001; Hayes et al., 2012; Ibrahim et al., 2018; S.-J. Jeong, 2014; Om Ariara Guhan et al., 2016; Su et al., 2013). With this approach, the whole region of the porous body is discretized with a numerical grid that does not reproduce the real geometry of the channels and the solid structure between them. The volumetric ratio between fluid and solid structure defines the porosity. Using volume averaging leads to a loss of information at the micro-level but provides insight into the macroscopic effects (e.g. flow maldistribution inside the monolith). The strongly anisotropic properties of a monolith have to be deduced with additional models. It is common practice to express the pressure loss in the channels with the Hagen-Poiseuille equation. The flow resistance in the transverse direction is set to a large value to prevent flow in that direction. Sometimes CFD models of this type are referred to as Volume Averaged Navier Stokes (VANS) approaches.

lP

re

-p

In most technical applications (e. g. automotive converters) the flow inside the channels of monoliths is laminar because of small channel diameters (typically in the order of one millimeter). The flow distribution inside the monolith is mainly influenced from the flow domain upstream of a monolith (Guojiang & Song, 2005; Hayes et al., 2012; Mat Yamin et al., 2013). As the profile of the axial velocity (into direction of the channels) downstream of the monolith reflects the distribution of the volume fluxes among the different channels inside the monolith, the comparison of measured and calculated profiles downstream of the monolith is of large interest (Fig. 1-1).

Jo

ur

na

Typically, the flow up- and downstream of a monolithic reactor is turbulent. For the CFD calculation of these domains mainly RANS (Reynolds Averaged Navier Stokes) models are applied. Important aspects of turbulence modeling for such applications are near-wall turbulence, flow separation and laminar to turbulent transition. Various turbulence models have been utilized in order to best simulate these flow features. Most literature studies of flow up- and downstream of a monolith employ the well-established k- turbulence model. (Benjamin et al., 1996) show that the usage of the k- RNG (ReNormalization Group) or the Norris-Reynolds (NR) two-layer wall model yields better results than the standard k- model. Based on this finding, the non-linear k- turbulence model coupled to the NR two-layer wall model was used in (Benjamin S. F. et al., 2001; Benjamin et al., 2004). Other authors apply the algebraic turbulence model (Chakravarthy et al., 2003), the Low Reynolds k- model (Tsinoglou et al., 2004), the 2f turbulence model (Porter et al., 2014; Porter et al., 2016), the k- model (Cornejo, Nikrityuk, et al., 2018a; Ibrahim et al., 2018) or the k- SST (Shear Stress transport) model (Cornejo, Nikrityuk, et al., 2018a). Comparisons of results for different turbulence models can be found in (Clarkson, 1995; Cornejo, Nikrityuk, et al., 2018a; S.-J. Jeong, 2014), with the last source favoring the k- SST model. The k- SST model may be somewhat preferred in the recent literature over other RANS approaches -4-

because it combines the benefits from both the k- model (description of far-field) and the k model (description of flow close to the wall), without a large computational penalty. Additional details on the turbulence models can be found in section 3.4.

ro of

In general, the position of flow separation in a diffusor is hard to predict using RANS turbulence models. Since the opening angle of the diffusor is normally large the separation point is well defined (Fig. 1-1). Additionally, wall effects play a minor role. Therefore, the influence of the turbulence model on the flow field in front of a monolith is weak. Nevertheless, small geometrical details, such as a weldseam near the detachment, can play an important role (Hayes et al., 2012; Holmgren et al., 1997). The velocity distribution inside the monolith also depends on the flow profile in the tube upstream of the diffusor. A plug flow profile leads to a more homogeneous distribution than a profile exhibiting a boundary layer at the tube wall (Hayes et al., 2012; Lai et al., 1992; Ozhan et al., 2014).

re

-p

Special interest must be paid to the details of the flow within the first several millimeters of the monolith channel length. Here, turbulence decays quickly (Cornejo, Nikrityuk, et al., 2018a) leading to a fully developed laminar flow profile. When the inlet of the channels is partially misaligned from the direction of the bulk flow, flow separation can occur inside the channels (Fig. 1-1 lower left side). Therefore, some authors add a correction function for the pressure loss in channels where the inflow shows an angle of attack,  (Benjamin S. F. et al., 2001; Benjamin et al., 1996; Chakravarthy et al., 2003; Tsinoglou et al., 2004). Corrections of these types improve the results for simulations of short monoliths.

Jo

ur

na

lP

Since most automotive converters are axisymmetric, the porous body strategy is often implemented in two-dimensions (axial and radial), which is the approach taken by most of the authors referred to above. This approach is feasible as long as the inflow and the outflow are symmetrical. The influence of asymmetrical geometry and the optimization of the flow distribution inside a monolith are investigated in the works of (S.-J. Jeong, 2014; S. Jeong & Kim, 1997; Om Ariara Guhan et al., 2016; Ramadan et al., 2007; Shijin et al., 2000; Will & Bennett, 1992). (Om Ariara Guhan et al., 2016) show that a combination of diffusor angle and a misalignment between the diffusor and monolith can lead to a strong improvement of the homogeneity. The homogeneity of the flow in the monolith can also be enhanced using flow distributors (Guojiang & Song, 2005). The important automotive application of two monoliths in series is addressed rarely in the literature, e. g. in (Lai et al., 1992; Shi-jin et al., 2000; Su et al., 2013; Weltens et al., 1993; Zhang & Tennison, 2008). When using the porous body approach together with turbulence models, special interest must be paid to the handling of turbulence inside and downstream of the monolith. Downstream of the monolith the flow is transitional and changes from laminar to turbulent conditions (Fig. 1-1 lower right side), which can hardly be captured using the RANS approach. Further, (Cornejo, Nikrityuk, et al., 2018a; Ibrahim et al., 2018) show that the values for turbulence -5-

inside the monolith have to be dampened. In contrast, a source term is generally needed to initiate turbulence downstream of the monolith.

ro of

The reason that turbulence initiation is required downstream of the monolith is that the velocity field shows no microscale (on the order of channel diameter) gradients of the axial velocity perpendicular to the main flow direction. This leads to an underprediction of turbulence production. Moreover, (Cornejo, Nikrityuk, et al., 2018b) show that even when some channels are fully discretized, capturing the details of individual jets coming out of the channels, the values of turbulence are now too large for the standard k- model. However, the usage of the k- SST model shows a similar level of turbulent kinetic energy as a LES (Large Eddy Simulation) but with a different spatial distribution. The validation of turbulence modeling downstream of a monolith lacks in the availability of experimental data. Only few measurements of turbulence values have been published up to now. Only in the paper (Eggenschwiler et al., 2009) two-dimensional distributions are given. In (Ibrahim et al., 2018) measurements along a line 5 mm downstream of a monolith are presented.

lP

re

-p

Another possibility is to model only one channel of the whole monolith in 1D/2D/3D. This approach is used for calculations including chemistry and is only successful if one channel is a good representative for all other channels of the monolith. If the conditions inside the monolith are inhomogeneous it is far better to calculate a small number of channels and assume each of them as a representative for a group of channels (Chakravarthy et al., 2003; Deur M. et al., 2002; Hettel et al., 2015; Martin et al., 1998; Tischer et al., 2001). Another approach is to exclude the channels of the catalyst from the CFD-domain and to solve them using 1D/2D grids using a second solver (Hettel et al., 2015; Jan Štěpánek et al., 2011; J. Štěpánek et al., 2012; Sui et al., 2016).

Jo

ur

na

If the geometry of the converter is prismatic (planar), it is possible to model only a representative row of channels in detail (Porter et al., 2014). Here, only the first millimeters of the channels are resolved to capture the effects of flow separation. The remaining part of each channel, which exhibits a developed flow profile, is shrunk to a length of 1 mm (Porter et al., 2016). The pressure loss of this truncated channel is set to the same value as that of the fulllength channel. (Ozhan et al., 2014) replaced the catalytic converter by an interface where the pressure drop as a function of the velocity, according to the Hagen–Poiseuille equation, is applied. Works where all the channels of a real converter are resolved by the numerical grid include restrictions in dimensionality, complexity of models or in the number of channels considered. Either the number of channels resolved is large (some thousands) and there is no chemistry involved, or the number of channels is small (up to some hundreds), which allows the use of computationally expensive models such as LES or detailed chemistry. The 2D calculation of a quarter of a converter with 85 channels using global chemistry is shown by (Agrawal et al., 2012). (Mei et al., 2006) presented the modeling of one sixth of a reactor with ca. 240 channels for the catalytic combustion of CH4 applying a one-step -6-

chemistry mechanism. An example of the calculation of only hydrodynamics in 7539 channels using the Lattice Boltzmann method can be found in (Bertrand et al., 2012). We previously discussed the 3D modeling of a quarter of a CPOX (Catalytic Partial Oxidation) reformer including two monoliths with ca. 300 channels (Hettel et al., 2018). Detailed chemistry was applied and radiation effects between the two solids were considered. The results are compared to detailed measured data. (Choudary & Mazumder, 2014) show the modeling of catalytic combustion of methane in a 13-channel monolith, and (Kumar & Mazumder, 2010) simulated a quarter of a CPOX reformer with 293 channels. Both authors apply detailed chemistry but without comparison to experimental data.

lP

re

-p

ro of

In the present work, a new approach is developed to calculate the flow distribution in catalytic converters by modifying the multi-region solver chtMultiRegionFoam. In addition to modifications to support region coupling of pressure and velocity, the solver was also modified to utilize Local Time Stepping (LTS), whereby the time step is adjusted for each cell in the domain to be as high as possible while still maintaining stability (Espinoza D. E. et al., 2015). As a result of these upgrades to the multi-region solver, along with custom boundary conditions, the monolith is completely excluded from the calculation domain. Similar to porous body approaches the pressure loss of the monolith channels are calculated via the Hagen-Poiseuille equation. However, the channel inlets and outlets are discretized in the CFD domains upstream and downstream of the monolith. Therefore, the flow field near the monolith can be calculated with much more detail than with a porous body approach, with which the microscopic details of the flow cannot be resolved. Furthermore, the presented approach allows the usage of different turbulence models in the flow domains upstream and downstream of the monolith.

na

The functionality of this approach is tested for a case without chemistry from literature (Porter et al., 2014; Porter et al., 2016). The monolith considered consists of approximately 4600 channels.

Jo

ur

2 Measurements from Literature 2.1 Experimental Setup Data from literature (Porter et al., 2014; Porter et al., 2016) were used for the validation of our new approach. The authors investigated a system consisting of a diffusor (domain 1), a monolith and a downstream section (domain 2). Fig. 2-1 shows a schematic diagram of the isothermal flow rig used for the steady state flow measurements, including the dimensions. The system is not axis-symmetric but prismatic. The dimension into z-direction is 96 mm.

-7-

Symmetry plane (xy-plane) of investigated system. The dimension into z-direction is 96 mm. Every flow domain has its own coordinate system. Data was measured in domain 2 at x = 40 mm (indicated with a purple line)

ro of

Fig. 2-1:

2.2 Experimental Data

lP

re

-p

Flow is from left to right. The flow at the inlet originates from an upstream section which is comprised of a flow straightener and a contraction nozzle. Domain 1 (green) consists of a short inlet section and the diffusor with a total opening angle of approximately 60°. It is followed by an unwashcoated cordierite monolith (orange) with the length of 100 mm. The height of the quadratic channels of the monolith is Hchannel=1.12 mm, the thickness of the walls separating the channels is Lwall=0.15 mm. The nominal cell density is 62 cells/cm2 (=400 cpsi or cells per square inch) and the porosity is 0.77. The monolith is followed by an outlet sleeve (domain 2 = blue) to prevent any environmental influences.

ur

na

The experiments of (Porter et al., 2014; Porter et al., 2016) were performed with air at isothermal conditions. The temperature is not given in the papers, so we assumed a value of 298 K. Profiles of the axial velocity were measured 40 mm downstream of the monolith applying hot-wire anemometry (HWA). The authors used a TSI IFA 300 HWA-system including platinum-plated tungsten wires (Dantec 55 P11) with a diameter of 5 m.

Jo

The three Reynolds numbers, Re, applied in the experiments were chosen to be representative of typical passenger sized vehicles and are given in Table 2-1. They are based on a uniform velocity uinlet and the hydraulic diameter of the inlet (left boundary in Fig. 2-1). The hydraulic diameter is calculated to be Dhydraulic=4·cross sectional area / wetted perimeter = 38.4 mm. Also given are the values of gas hourly space velocity (GHSV). Names for the cases listed in the first column will be used in the discussion of the results. The profiles of the axial velocity ux used for the validation of our numerical approach are plotted in Fig. 2-2. In Fig. 2-3 the values are normalized with the inlet velocity uinlet. The HWA measurements were taken across the full width of the outlet sleeve. The profiles were -8-

performed on either side of the symmetry plane (xy-plane) of the system (into positive and negative y-direction) and are plotted versus the modulus of the y coordinate. A variability of around 5% between the two measurements is observable. As stated before, the profiles at the measured position (at x = 40 mm downstream of the monolith) can be interpreted as a good representation of the flow distribution inside the monolith.

ro of

The profiles in Fig. 2-2 and Fig. 2-3 can be explained based on the typical flow characteristics shown in Fig. 1-1. As the flow separates at the beginning of the diffusor a confined jet is produced which traverses the diffusor. This leads to the maximum of velocity in the center region of the monolith. The jet decelerates in front of the monolith and spreads into ydirection. Part of this flow is fed into the recirculation zone. Subsequently, the flow entering the monolith is maldistributed leading to the highest velocity in the channels near the center of the monolith. The flow into y-direction in front of the monolith stagnates towards the diffusor wall. This causes a rise in static pressure and encourages flow through the outermost channels and results in secondary velocity peaks near the outer boundary of the system.

Jo

ur

na

lP

re

-p

It is expected that the flow pattern will be independent of the Re, especially for large values of Re. The plot of the normalized profiles in Fig. 2-3 shows that the Reynolds numbers are too small to yield similarity of the flow field. As the Reynolds number increases the maxima and the minima of the profiles become more pronounced.

Fig. 2-2:

Measured velocity profiles 40 mm downstream of the monolith (non-normalized)

Fig. 2-3:

-9-

Measured velocity profiles 40 mm downstream of the monolith (normalized)

3 Numerical Setup 3.1 CFD-Code The modeling was performed applying the CFD-code OpenFOAM® (OpenFOAM, 2017). The open-source-software offers complete freedom to customize and extend its existing functionality to the user. We used version 6 as a basis for the development of a new solver.

-p

ro of

The present application focuses on the isothermal flow field. Because only the pressure loss of the channels is needed to model the flow-field up- and downstream of the monolith, the monolith is completely removed from the calculation domain by applying the HagenPoiseuille approach to account for the pressure loss in the monolith channels. The flow regions up- and downstream of the monolith are coupled with appropriate boundary conditions (section 3.3). However, the channel inlets and outlets are considered in the CFD grids. Therefore, the flow in the near field up- and downstream of the monolith can be calculated with much more detail than with a porous body approach, where the microscopic details of the flow cannot be resolved. The developed solver allows the calculation of multizonal domains divided by monoliths, including the flow domains up- and downstream and between monoliths.

Jo

ur

na

lP

re

The consideration of detailed chemistry in the channels in the future can be made feasible by applying the approaches outlined here to the computational tool DUO (Benzinger et al., 2019; Hettel et al., 2018; Hettel et al., 2015), which stands for the coupling of the two computer codes OpenFOAM and DETCHEMTM. DETCHEMTM (Deutschmann et al., 2014) is a commercial software package, being rather inexpensive for academic purposes. It is specifically designed for the modeling and simulation of reacting flows based on elementary step mechanisms, particularly for heterogeneous systems such as catalysis, materials synthesis, and fuel cells. With DUO the flow inside the channels of a coated monolith can be calculated by coupling with stand-alone DETCHEMTM executables which handle the detailed chemistry, including additional features such as washcoat models. As only the pressure loss in the channels is needed for the calculation of the flow field up- and downstream of the monolith, the pressure loss can easily be extracted from the calculation of the channels and used for the calculation of the coupled domains.

3.2 Computational Domains and Grids The investigated geometry (Fig. 2-1) is symmetrical with respect to the xy-plane and the xzplane. Therefore the calculation domain comprises only one quarter of the converter system (Fig. 3-1), with the region upstream of the monolith named domain 1, and the region downstream of the monolith called domain 2. Domain 2 was extended in axial direction, relative to the experimental setup, because the measured data was extracted at an axial - 10 -

Geometry of the two calculation domains (red). The coordinate system depicts the directions x, y and z. In each of the two flow domains the origin is positioned in the inlet plane on the symmetry axis of the converter (Fig. 2-1).

na

Fig. 3-1:

lP

re

-p

ro of

position of 40 mm, which is at the outlet plane of the assembly to the atmosphere. Since we do not want the outlet boundary conditions applied in the calculations to influence the results at x = 40 mm, we extended the calculation domain about 20 mm versus the real assembly. Therefore, the outlet boundary plane was positioned at x = 60 mm in the calculations so that the plane for data comparison does not lie on the outlet patch. In doing so, we added a portion of the surrounding atmosphere to the calculation domain. The boundary conditions at the four side planes of the extension were chosen to be symmetrical. The whole monolith consists of approximately 4500 channels. As mentioned before, the monolith is excluded from the calculation, but the inlets and outlets of the monolith are included in domains 1 and 2 respectively. Therefore, with quarter symmetry, domain 1 has about 1100 outlets and domain 2 has about 1100 inlets (i.e., 30 x 36 channels at the inlet and outlet of the monolith).

Jo

ur

Since the geometrical shapes of the flow domains are simple it was possible to generate grids using hexahedral cells with a relatively narrow distribution of cell sizes. The cell counts (ncell) of the grid pairs (one grid for domain 1 and one for domain 2) used are given in Table 3-1. The cell count increases from set coarse to fine. In set very fine the cell count of domain 1 is the same as that for set medium (5 mill. cells), but domain 2 was discretized with around 75 million cells. This latter set was used for LES calculations. To get a feeling of the resolution, a characteristic length scale of the cells, lcell, for each grid is also provided in Table 3-1. This characteristic length can be derived by first dividing the volume of the corresponding flow domain by the number of cells to get a characteristic cell volume. Then, the cubed root of this value is lcell. To get an impression of the grid topologies, renderings of the grid pair medium (Table 3-1) can be found in Fig. 4-3 for domain 1 and in Fig. 4-4 for domain 2.

- 11 -

3.3 Domain Coupling Fig. 3-2 shows the setting of the boundary values (marked with a dot) of a CFD calculation in the usual case of one flow domain including a monolith (upper part), and in the case of two coupled flow domains without a monolith (lower part). In a CFD calculation for a catalytic device consisting of one single flow domain (monolith included), the boundary values have to be set at the inlet and outlet. Normally, at the inlet the volume / mass flux or the velocity uin,0 is set whereas at the outlet the static pressure pout,0 is held constant.

Jo

ur

na

lP

re

-p

ro of

In the present approach, the monolith is completely removed from the flow domain, resulting in two distinct flow domains upstream (domain 1) and downstream (domain 2) of the monolith. To explain the coupling of these domains a part of the original domain (Fig. 3-2, upper picture) is cut out (red rectangle) and plotted enlarged in the lower picture of Fig. 3-2. The region stretches across two channels of the original geometry. These two channels are no longer existent in the grid. The inlets to the original monolith channels are incorporated as outlets in the grid of domain 1. Consequentially, the outlets from the original monolith channels are considered as inlets in the grid of domain 2.

Fig. 3-2:

Principle of domain coupling. Upper picture: setup of boundary conditions for standard CFD including a monolith (one flow domain). Lower picture: exchange of boundary conditions between two coupled domains in the present approach. The boundary conditions which are marked with blue dots have to be set. The boundary conditions which are marked with red dots are part of the solution. As discussed in the text, mass flow is mapped between domain 1 and domain 2, but mass flow is converted to velocity with knowledge of the channel cross-sectional area and fluid density at each monolith outlet (inlet to domain 2). - 12 -

ro of

The two domains are solved sequentially by the solver. Therefore, the boundary conditions which are marked with blue dots have to be set. The boundary conditions which are marked with red dots are part of the solution. The conditions at the inlet of domain 1 (velocity uin,0) and at the outlet of domain 2 (pressure pout,0) are fixed for the whole calculation process. For the first iteration the boundary values at the outlet of domain 1 (n pressure values pout,i for n channels) and the inlet values of domain 2 (n velocity values uin,i for n channels) have to be estimated. The solution of domain 1 yields the n velocity (and mass flux) values at the outlets uout,i which can be used to calculate inlet values uin,i for the next iteration of domain 2. Specifically, the mass flux calculated from the velocity profile at the outlet patch into the ith channel (i.e., the domain 1 outlet for the ith channel) is supplied as input to the velocity inlet boundary condition for the ith channel of domain 2. Then, with knowledge of the channel cross-sectional area and average density, the average velocity of ith inlet (uin avg,i) can be calculated by: Ṁ

uin avg,i = A∙ρ out,i ,

(3-1)

in avg,i

-p

where Ṁout,i is the mass flow rate at the ith outlet of domain 1 (into the monolith), A is the cross-sectional area of a monolith channel and ρin avg,i is the average density at the inlet to the ith channel.

lP

re

The application of equation 3-1 will yield a constant inlet velocity at the outlet of the monolith, but a velocity profile is expected due to fully developed laminar flow within the monolith channels. The fully developed velocity profile in a rectangular channel was thoroughly discussed by (Wörner, 2010), from which the velocity profile approximation of (Purday, 1949) was selected for this work: b+1 a+1 b

a

na

uin,i = uin max,i (1 − Y a )(1 − Z b ) =

uin avg,i (1 − Y a )(1 − Z b ) . (3-2)

ur

uin max,i is the maximum of the velocity profile at the ith inlet to domain 2, Y is defined as the relative inlet face y-coordinate value divided by the channel half-height, and Z is defined as the relative inlet face z-coordinate value divided by the channel half-width. Additional details on the calculation of a and b are found in (Wörner, 2010), but for the present channel

Jo

geometry, a = b = 2.2, and the quantity

b+1 a+1 b

a

equals 2.116.

Next, to conserve mass, a factor is calculated on the fly to ensure the applied velocity profile results in the exact mass flow into each inlet i to domain 2. This factor, called ∅Factor,i , is calculated by numerically integrating the mass flow of each face cell j, of each inlet i to domain 2, and then dividing the calculated mass flux at each inlet i after the application of the velocity profile by the mass flowrate out of domain 1 at each outlet i. ∅Factor,i =

faces ∑# uin i,j ∙ Ai,j ∙ ρi,j j

Ṁout,i

(3-3)

, - 13 -

Finally, the inlet velocity at each inlet face is modified by ∅Factor,i so that the velocity uin i,j at a face j of an inlet i to domain 2 is: uin i,j = ∅Factor,i ∙ 2.116 ∙ uin avg,i (1 − Yj a )(1 − Zj b ) .

(3-4)

Next, the solution of domain 2 allows for the calculation the n pressure values pout,i in domain 1. The values of the static pressure at the outlets of domain 1, pout,i, are calculated by adding the pressure loss ∆𝑝𝑙𝑜𝑠𝑠,𝑖 of the associated channel to the domain 2 inlet pressure of the same channel (at the monolith outlet), pin,i: 𝑝𝑜𝑢𝑡,𝑖 = 𝑝𝑖𝑛,𝑖 + ∆𝑝𝑙𝑜𝑠𝑠,𝑖 .

(3-5)

32

ro of

The calculation of the pressure loss is described by the Hagen-Poiseuille equation: 𝐿

∆𝑝𝑙𝑜𝑠𝑠,𝑖 = 𝑅𝑒 ∙ 𝜌𝑜𝑢𝑡 𝑎𝑣𝑔,𝑖 ∙ 𝑢𝑜𝑢𝑡 𝑎𝑣𝑔,𝑖 2 ∙ 𝐷 .

(3-6)

Here, Re is the Reynolds number, 𝜌𝑎𝑣𝑔 is the average density at the ith outlet patch of domain

lP

re

-p

1, 𝑢𝑜𝑢𝑡 𝑎𝑣𝑔,𝑖 is the average velocity at the ith outlet patch, L is the channel length (100 mm), and D is the hydraulic diameter (1.12 mm). Because the monolith is long, entrance effects are negligible and were not included in the calculation of pressure drop. If necessary entrance effects can be considered using the same general approach by including additional pressure loss terms to equation (3-6) (Benjamin S. F. et al., 2001; Benjamin et al., 1996; Chakravarthy et al., 2003; Tsinoglou et al., 2004).

na

3.4 Turbulence Modeling



ur

The basic equations solved in CFD for modeling turbulent flows are the averaged equation for conservation of mass and the three averaged equations for conservation of momentum (averaged Navier-Stokes equations):



   ui  0. t x j

Jo

(3-7)



  ui u p    ui  uj i     ij  .  l  t x j x j x j  x j 

(3-8)

 l is the laminar viscosity. In the RANS (Reynolds Averaged Navies Stokes) context the

overbar above physical quantities means that they are time-averaged. In the LES (Large Eddy Simulation) context the overbar means that the quantities are filtered or averaged in space.

- 14 -

Due to averaging the equation system is not closed anymore. The matrix of stresses ij contains new unknowns. A turbulence model is an approach to provide these unknowns. The presented approach of domain coupling (section 3.3) allows the usage of different turbulence models in the two flow domains upstream and downstream of the monolith. In flow domain 1 (upstream) the two-equation k- SST model was used in all calculations. In flow domain 2 (downstream) the k- SST model, the k- and the LES approach were evaluated. All three models will be described briefly. For all calculations the assumption of isothermal flow was applied. RANS Modeling

ro of

Within the framework of the RANS approach, the six elements of the matrix ij represent the turbulent stresses:  2    k ij .  3

(3-9)

-p

 u u ij  t   i  j  x j xi 

k is the turbulent kinetic energy and ij the Kronecker symbol.

re

In the k- model (Launder & Spalding, 1974; Rodi, 1993), the turbulent viscosity t is

t   C 

lP

defined as: k2 . 

(3-10)



x j

t   k    2   l    t  S    ,  x j  k  x j 



Jo





ur

  ui  k

na

Two additional transport equations for k and  have to be solved. For steady state conditions and an incompressible medium these are:

  ui   x j

 x j

 t   l   

    2 2 .   C    S  C   1 t  2 k k  x j 

(3-11)

(3-12)

S is the strain rate tensor, the product t  S 2 is the production term. The standard formulation implies the following constants: C1  1.44 , C2  1.92 , C  0.09 , k  1.0 and   1.3 . The k- SST (Shear Stress Transport) model (Menter, 1994) is a combination of the k- model and the k- model, and, as mentioned in the introduction, it generally inherits the - 15 -

advantages of both models: The description of far-field from k- model and the prediction of flow close to wall from the k- model. Therein, the turbulent viscosity t is defined as: t  

k 

1  1 SF  max  * , 2    a1  

. (3-13)

The model also includes a limiter to the production of turbulence. The two additional transport equations for the turbulent kinetic energy k and the turbulence dissipation  are:

x j





  ui  x j



 x j

 t   l  k 

 k  2 *   t S   f* k  ,   x  j 

(3-14)

ro of



  ui k

t     k    * 2 2   l     S   f   2 1  F1    x j    x j  x j x j

(3-15)

-p

Given its complexity and extension, the various constants and functions of the SST transport equations are not described in this paper. They are discussed in (Menter, 1994).

re

LES Modeling

lP

Within the framework of the LES approach, the filtered Navier Stokes equations solve the flow up to a cut-off length scale limit,  f , defined by the characteristic cell length. Below the

na

cut-off size a sub-grid (SGS) scale model is employed. In the present work the sub-grid scale was modeled by the Dynamic Kinetic Energy model (Huang & Li, 2010; Kim, 1995; Yoshizawa & Horiuti, 1985). Here, the matrix ij represents the SGS stresses: 2 ij   kSGS ij  2Ck k1/2  f Sij . 3

ur

(3-16)

Jo

These stresses include the interaction between the resolved flow and the vortices which cannot be resolved from the numerical grid. An additional transport equation for the sub-grid scale turbulent kinetic energy, kSGS , derived to account for the historic effects due to production, dissipation and diffusion, has to be solved: 3/2  kSGS  kSGS kSGS    SGS  kSGS   uj  ij Sij  C     t x j f x j  k x j

  . 

(3-17)

Sij is the resolved strain rate tensor. The sub-grid scale viscosity is defined as:

- 16 -

1/2 SGS  Ck  kSGS f .

(3-18)

For reasons of clarity the overbars over time- or space averaged variables will be not used in the discussion of the results of the RANS and LES calculations.

3.5 Boundary Conditions It is common practice to define the intensity of the turbulent fluctuations of the three velocity components u, v, w in terms of root mean square quantities, that is,

uRMS  u2 , vRMS  v2 ,wRMS  w2 , .

ro of

(3-19)

The degree of turbulence is the relative intensity and usually expressed as a percentage. It is defined as: TU  uRMS u .

-p

(3-20)



 k 2 t

 m2   s3  ,  

 C  k 2

lP

 m2  3 2 k   uRMS  s2  , 2  

re

whereas 𝑢̅ is the time average of the velocity into x-direction. For isotropic turbulence, at equilibrium the turbulent fluctuations are the same (uRMS = vRMS = wRMS ). With this, the turbulent kinetic energy, k , and the specific turbulence dissipation,  , can be deduced by:



t

 m2   s3  .  

(3-21)

na

μt is the turbulent dynamic viscosity and C  0.09 . We estimated that

Domain 1

ur

μt ⁄μl = 10 using the values at room temperature for the dynamic (laminar) viscosity (μl = 17.6 ∙ 10−6 kg/(s ∙ m)) and density (ϱ = 1.17 kg/m3 ) of air.

Jo

As stated before, only the k- SST-model was used for the calculation of the flow domain 1. The boundary conditions of the inlet of domain 1 are summarized in Table 3-2. The uniform inlet velocity uinlet for the three cases under investigation can be calculated from the Reynolds numbers Re given in Table 2-1. Boundary values for the turbulent kinetic energy k and the specific dissipation rate ε for a degree of turbulence of TU=1% were obtained using the relations (3-21). It was verified, that the degree of turbulence at the inlet of domain 1 shows no influence on the results in the tested range 0.1% > TU > 5%. The static pressure at the outlets of domain 1 were set according to the description in section 3.3. - 17 -

Domain 2 The velocity distributions at the inlets of domain 2 were set according to the description in section 3.3. For the turbulence values at the inlets of domain 2, two sets of parameters were utilized (Table 3-3). The set noTurbulence was used to approximate laminar conditions at inlet of domain 2, while the set smallTurbulence was used to introduce weak turbulence at the inlet of domain 2, assuming a degree of turbulence of 1%. A reference velocity, uref (the average velocity in the cross sectional area of domain 2 downstream the monolith) was used for calculating the turbulence values (Table 3-3).

ro of

The static pressure at the outlet of domain 2 was fixed to pout = 1 bar. In both calculation domains the xy-planes and xz-planes were defined to be symmetric. The extension of the original system by 20 mm in axial direction (section 3.2) was enclosed using symmetric conditions. All other boundaries were defined as walls. For the near wall treatment of the flow the standard conditions for the turbulence models have been used (e. g. logarithmic wall law for k- model).

-p

The symmetry, inlet and outlet boundary locations are appropriately marked in Fig. 4-1 for domain 1 and in Fig. 4-2 for domain 2.

re

3.6 Numerics

lP

RANS Calculations

ur

na

For all calculations, OpenFOAM’s PIMPLE algorithm (merged PISO-SIMPLE), which is a combination of PISO (Pressure-Implicit with Splitting of Operators) and SIMPLE (SemiImplicit Method for Pressure Linked Equations) was applied. All results were produced using 2nd order discretization in space except those in section 4.1.4, where the influence of applying 1st order discretization (hybrid scheme) is shown. First order discretization in time was used for RANS calculations, employing the LTS method in both domains 1 and 2 to accelerate the development of the steady-state solution.

Jo

LES Calculations

The LES calculations for the case Re42000 were performed with the PIMPLE solver. 2nd order discretization in space and time (backward scheme) was used. To guarantee a Courant number below 0.5 the timestep t=5·10-6 s was used (Table 3-4). Additionally, the characteristic time scale tchar for one pass through the domain 2 (length L = 60 mm) and the corresponding number of time steps nchar needed to simulate this physical time is given. Averaged results for LES calculations imply physical simulation times of ca. 3 characteristic time scales tchar. The calculations were performed using 500 processor cores, with a total simulation time of about 1 week to reach ca. 3 tchar.

- 18 -

4 Results The presented approach of domain coupling (section 3.3) allows the usage of different turbulence models in the two flow domains, upstream and downstream of the monolith. The focus of the investigations was on the flow field of domain 2 (downstream of the monolith), because measured data were available only downstream of the monolith. In order to limit the influence of too many parameters, for flow domain 1 (upstream) only the two-equation k- SST model was applied in all calculations. In the flow domain 2 (downstream) the k- SST model, the k- and the LES approach were used.

ro of

4.1 RANS Modeling 4.1.1 Flowfield Upstream of the Monolith

Jo

ur

na

lP

re

-p

Fig. 4-1 shows the axial velocity ux in domain 1 for the case Re30000 (Table 2-1) with the grid pair medium (Table 3-1). The flow field shows the characteristics which are discussed in section 1 (Fig. 1-1). A recirculation zone cane be seen at the upper region of the diffusor. The flow near the symmetry axis of the system (x axis at y = 0 and z = 0) reaches the monolith in horizontal direction. Streaklines which start at larger y-positions at the inlet reach the outlets of the domain 1 (inlets to monolith) with a component in the y-direction. The larger the yposition of the streaklines at the entrance, the larger is the angle of attack at the outlets. The enlargement of the region in the blue rectangle, which is plotted in Fig. 4-3, shows that the flow at the inlet of the channels is not equally distributed and the axial component is larger on the upper side of each outlet.

- 19 -

ro of

Axial velocity ux in domain 1 together with boundary conditions and some streaklines for the case Re30000. The k- SST model was applied, using the grid pair medium. The enlargement of the region in the blue rectangle is shown in Fig. 4-3.

re

-p

Fig. 4-1:

4.1.2 Flowfield Downstream of the Monolith

Jo

ur

na

lP

Fig. 4-2 shows the axial velocity ux in domain 2 for the case Re30000 (Table 2-1) with the grid pair medium (Table 3-1). The streaklines indicate that the flow downstream of the domain 2 inlets (outlets of monolith channels) is equally distributed with respect to the y direction. The inlets near the symmetry axis of the domain (z-axis at x = 0 and y = 0) exhibit the largest velocities. The developed laminar profiles of the jets can be indicated for some millimeters and then smear out. It is noted that the velocity profile at each inlet is symmetrical due to the special handling of the boundary conditions described in section 3.3. This can also be seen in the enlargement of the region in the blue rectangle which is shown in Fig. 4-4.

- 20 -

ro of

-p

Axial velocity in domain 2 together with boundary conditions and some streaklines for the case Re30000. The k- SST model was applied, using grid pair medium. The purple line indicates the position where the experimental measurements were done (Fig. 2-1). The enlargement of the region in the blue rectangle is shown in Fig. 4-4.

Jo

ur

na

lP

re

Fig. 4-2:

Fig. 4-3:

Axial velocity ux in domain 1: enlargement of the region in the blue rectangle of Fig. 4-1, showing details of the grid.

Fig. 4-4:

- 21 -

Axial velocity ux in domain 2: enlargement of the region in the blue rectangle of Fig. 4-2, showing details of the grid.

Normalized mass fluxes Ṁout,i,norm leaving the monolith channels (yz-plane at the monolith exit). Walls are the planes y = const. = 39 mm and z = const. = 48 mm. Symmetry is at y = const. = 0 mm and z = const. = 0 mm. The k- SST model was applied in both domains, using the grid pair medium.

na

lP

Fig. 4-5:

re

-p

ro of

Fig. 4-5 depicts the distribution of normalized mass fluxes Ṁout,i,norm leaving the monolith channels (entering domain 2). Every square of the overlaid mesh represents one channel i. The distribution exhibits characteristics previously discussed with respect to Fig. 1-1, where the mass flux peaks at the symmetry axis (y = 0 mm) and near the wall (y = 39 mm).

Jo

ur

This characteristic is still present in the plane x = const. = 40 mm downstream of the monolith. Fig. 4-6 shows the distribution of normalized axial velocity ux at this axial position. Here, the distribution of velocity is homogeneously distributed in the z-direction, neglecting the region near the wall, where the velocity decreases to zero. The assembly investigated by (Porter et al., 2014; Porter et al., 2016) was designed to maintain a two-dimensional flow with respect to the z-direction (Fig. 2-1). Based on the results of the calculation it can be stated that this intention is fulfilled with this converter system.

- 22 -

ro of

Normalized axial velocity ux,norm in the yz-plane at x = 40 mm downstream of the monolith. Walls are the planes y = const. = 39 mm and z = const. = 48 mm. Symmetry is at y = const. = 0 mm and z = const. = 0 mm. Measurements are taken along the y axis at the z =0 (purple line in Fig. 2-1). The k- SST model was applied in both domains, using the grid pair medium.

re

-p

Fig. 4-6:

lP

In the calculations shown up to now we used RANS turbulence modeling in both domains. However, from the measurements it is unknown if the flow downstream of the monolith is laminar or turbulent at the measuring position 40 mm downstream of the monolith. The flow in the channels is laminar (Table 4-1), proven by the Reynolds number of the channels, given by Rechannel:

na

Rechannel   uchannel  Hchannel / l .

(4-1)

Jo

ur

u̅channel is the average velocity in the channels for a plug flow profile, Hchannel=1.12 mm is the height of the quadratic channels, and μl is the laminar dynamic viscosity. u̅channel can be calculated from the average velocity downstream of the monolith, uref (Table 4-1) divided by the porosity, which is 0.77 (section 2.1). Since the Reynolds numbers are below the critical value of 2300 the flow leaves the channels under laminar conditions. In the paper of (Cornejo, Nikrityuk, et al., 2018b) it was found that the state of the flow downstream the monolith is dependent upon the characteristic Reynolds number Rewall (Table 4-1): Rewall   uchannel  2  Lwall / l ,

(4-2)

where Lwall=0.15 mm is the thickness of the wall between the channels. The authors claim that the flow downstream of the monolith is laminar if Rewall < 100. This statement can be - 23 -

interpreted in a way that the flow remains laminar downstream of the monolith up to a distance of minimum 20 mm. This is length of the flow domain which was used in the calculations in (Cornejo, Nikrityuk, et al., 2018b). Based on these findings, simulations have been performed where the flow in domain 2 was treated to be laminar (no turbulence model enabled).

Jo

ur

na

lP

re

-p

ro of

Fig. 4-7 shows the velocity profiles against the measurements. It can be seen that the single jets leaving the monolith can be identified clearly at the measuring position 40 mm downstream of the monolith. The calculation reflects the trend of the experimental profiles if one considers an “average” value of the minima and maxima of the wiggles. Clearly, laminar calculation for case Re30000 fits better with the measurements than the RANS calculation. This is because the turbulent viscosity seems to be overpredicted in the RANS calculations leading to smoother profiles. However, in the experiments, no wiggles can be identified which proves the fact that the flow is turbulent downstream of the monolith, at least at the measuring position. Therefore, the usage of RANS modeling seems to be justified so far. The results achieved with LES, discussed in section 4.2, reveal the characteristics of the flow in greater detail.

Fig. 4-7:

Comparison of calculated and measured velocity profiles downstream 40 mm of the monolith (purple line in Fig. 2-1) for grid pair medium. The k- SST model was applied in domain 1. Compared are the usage of k- SST model and laminar flow in domain 2.

- 24 -

4.1.3 Influence of Grid

na

Influence of grid pairs (Table 3-1). Comparison of calculated and measured profiles of axial velocity ux in domain 2 at x = 40 mm downstream of the monolith (purple line in Fig. 2-1). The k- SST model was applied in both domains.

ur

Fig. 4-8:

lP

re

-p

ro of

Three pairs of grids (Table 3-1: coarse, medium, fine) were also compared. The comparison of calculated axial velocity ux against measurement for the case Re42000, shown in Fig. 4-8, indicates that the solutions with the grid pairs medium and fine seem to be slightly better than the solutions obtained with the grid pair coarse. The gradient of the velocity into y-direction near |y|=20 mm is slightly steeper and the minimum at |y|=30 mm somewhat smaller which fits better with the measurements. However, the influence of the grid is very weak. We used the grid pair medium for most of the calculations.

4.1.4 Influence of Spatial Discretization Scheme

Jo

Two classes of spatial discretization schemes were also compared, 1st order (hybrid scheme) and 2nd order (linear scheme), with each set used in both domains. Fig. 4-9 shows that the influence of discretisation scheme is strong. As expected, the use of 1st order discretisation smooths the profiles. As the gradients of the profiles into y-direction increase with the Reynolds number, the difference between 1st order and 2nd order increases with average velocity. Thus, the usage of 1st order discretization schemes should be avoided except for achieving a rough solution. In all other calculations shown 2nd order discretization schemes were used.

- 25 -

ro of

-p

Influence of discretization scheme for grid pair fine. Comparison of calculated and measured profiles of axial velocity ux in domain 2 at x = 40 mm downstream of the monolith (purple line in Fig. 2-1). The k- SST model was applied in both domains.

re

Fig. 4-9:

lP

4.1.5 Influence of Turbulence Model and Boundary Conditions

na

Whereas in domain 1 the k- SST model was always applied, in domain 2 the influence of the turbulence model (k- versus k- SST) was investigated. Additionally, the turbulence parameters at the 1100 inlets of domain 2 (which represent the outlets of the 11000 channels of the monolith) have been varied using the two parameter sets noTurbulence and smallTurbulence (Table 3-3).

Jo

ur

Fig. 4-10 shows the comparison of the calculated profiles of axial velocity to the experiments for case Re42000. The influence of turbulence model on the velocity profile can clearly be seen near the wall. Here, the k- SST model predicts a steeper increase than the k- model, which agrees better with the experiments. However, we did not prepare the grid such that the first calculation cell was in the correct dimensionless distance from the wall, which is necessary for properly applying the logarithmic wall law associated with the application of the k- model. No influence of turbulence boundary conditions can be seen using the k- model. In the case of the k- SST model, the application of the parameter set noTurbulence leads to a slightly better coincidence with the measurements as compared to the usage of parameter set smallTurbulence. Starting with negligible turbulence at the inlets leads to a slightly smaller overall level of the turbulent kinetic energy (Fig. 4-11) and therefore to a slightly smaller exchange of inertia in the flow. Therefore, the profile is a bit steeper (larger maxima at - 26 -

|y|=0.005 mm and smaller minima at |y|=0.027 m) than the profile obtained with the parameter set smallTurbulence. The trend toward steeper profiles with decreasing turbulence level can also be identified in the laminar solutions shown in Fig. 4-7 if one considers an “average” value of the minima and maxima of the wiggles. These solutions represent the cases with smallest exchange of inertia exhibiting the largest maxima and smallest minima for all Reynolds numbers.

lP

re

-p

ro of

Comparisons of the turbulent kinetic energy on the symmetry axis (symmetry axis of the center channel, at z = 0 mm and y = 0 mm) are shown in Fig. 4-11. For all four cases the level of k increases strongly, with a maximum between 5 – 10 mm downstream the monolith. Farther downstream the values decrease strongly, resting at a small value << 1 m2/s2. For both of the applied turbulence models, the usage of boundary conditions smallTurbulence leads to larger maximum of k values which are positioned a bit earlier downstream of the monolith, as compared to the case noTurbulence. For the k- model, the level of k downstream of x = 20 mm is effectively not influenced by the boundary conditions, and downstream of x = 40 mm k is at the same level as the k- SST calculation for the boundary conditions noTurbulence. The calculation using the k- SST model and boundary condition smallTurbulence exhibits the highest level of k at a downstream position of x = 5 mm. The velocity profile calculated at these conditions does not differ much from the profiles obtained with other turbulence models and boundary conditions, but it exhibits the smoothest velocity profile (profile with smallest derivatives of velocity into y-direction). The production of turbulence downstream of the monolith cannot be captured using porous body models of the monolith. Hence, additional source terms have to be introduced (Cornejo, Hayes, et al., 2018). As the main issue of the paper is the coupling and not turbulence modeling, more details of turbulence shall not be discussed here.

Jo

ur

na

The results of the LES calculations shown in section 4.2 prove that the RANS models cannot predict transition in free flow correctly and show an unphysical distribution of turbulence values. It is therefore surprising that RANS models can approximate the experimental velocity profiles so well. This is owed to the fact that the flow in domain 2 is dominated by the flow distribution coming out of the monolith and the absence of any curvature of the streaklines downstream.

- 27 -

ro of

Jo

ur

na

lP

re

-p

Fig. 4-10: Comparison of calculated and measured velocity profiles downstream 40 mm of the monolith (purple line in Fig. 2-1) for case Re42000 and grid pair medium. The k- SST model was applied in domain 1. Parameters were the RANS turbulence model (k- versus k- SST) in domain 2 and the turbulence level (noTurbulence versus smallTurbulence) at the inlets of domain 2.

Fig. 4-11: Turbulent kinetic energy in domain 2 along symmetry axis (center of central channel at y = z = 0 mm) for case Re42000. Parameters were the RANS - 28 -

turbulence model (k- versus k- SST) in domain 2 and the turbulence level (noTurbulence versus smallTurbulence) at the inlets of domain 2.

4.2 LES Modeling For all the results shown in this section the k- SST model was applied in domain 1 and the Dynamic Kinetic Energy LES model was applied in domain 2. We used the grid pair very fine (Table 3-1) with 5 mill. cells in domain 1 and 74 mill. cells in domain 2. The LES calculations were done for case Re42000. 4.2.1 Instantaneous Velocity

Jo

ur

na

lP

re

-p

ro of

Firstly, we focus on the distribution of instantaneous axial velocity ux in domain 2 on the symmetry plane (xy-plane at z = 0). The snapshot shown in Fig. 4-12 proves that the jets which leave the monolith remain laminar up to a distance of 10-20 mm. No unsteady features can be identified upstream of this position. Afterwards, transition occurs and “real” turbulence can be identified ca. 30 mm downstream of the monolith. Due to the unsteady nature of the process these dimensions vary in time. This behavior is completely different from the results of the RANS models which predict the maximum of the turbulent kinetic energy in around the first 10 mm downstream of the monolith (Fig. 4-11).

Fig. 4-12: Instantaneous axial velocity ux in domain 2 on symmetry plane (xy-plane at z = 0 mm) for case Re42000. Flow is from left to right. Symmetry conditions are on the bottom side, the wall is on the top side.

- 29 -

ro of -p re lP na

Jo

ur

Fig. 4-13: Snapshots of instantaneous axial velocity ux in domain 2 on yz-plane at x = 40 mm downstream of the monolith for case Re42000. View is into direction of monolith (negative x direction). Symmetry conditions are at the bottom and right sides, walls are at the left and top sides. The difference in simulated time between the two pictures is 0.02 s. In the upper picture the major part of the cross-sectional area appears to be in the laminar state, whereas the state in the lower picture is slightly turbulent. A cross sectional slice of the field of ux at x = 40 mm is depicted in Fig. 4-13. The state of turbulence is not equally distributed. Regions with nearly laminar conditions and regions with more turbulent conditions can be identified. As time proceeds, the zone near the center of the system (around z = 0 mm) is permanently turbulent whereas in other regions the flow changes between nearly laminar and turbulent conditions. This intermittency is typical for transitional flows. - 30 -

4.2.2 Time Averaged Velocity

lP

re

-p

ro of

The field of time averaged axial velocity u̅x in domain 2 on the symmetry plane (xy-plane at z = 0) is plotted in Fig. 4-14. It shows that the identity of the single jets coming out of the channels vanishes downstream of ca. x = 20 - 25 mm. The process of time averaging was done during ca. 3 typical time scales tchar of the system (Table 3-4). The identification of small structures in the region near the axis at x = 30 mm shows that the time averaging was not long enough to achieve an optimal time averaged flow field. The figure also depicts a disadvantage of the boundary conditions used for LES. The fluctuation is hindered at the symmetry planes. This will be discussed together with Fig. 4-16.

na

Fig. 4-14: Time averaged axial velocity u̅x in domain 2 on symmetry plane (xy-plane at z = 0 mm) for case Re42000. Flow is from left to right. Symmetry conditions are on the bottom side, the wall is on the top side.

Jo

ur

The cross sectional slice of the field of u̅x at x = 40 mm is depicted in Fig. 4-15. The jets in the zone z > 30 mm and y < 10 mm are nearly laminar at this axial position. But as with the instantaneous velocity plot (Fig. 4-13), even in the more turbulent regions, the matrix of jets can be slightly identified. As mentioned, looking at animations of the flow field, the region near the center of the system is permanently turbulent whereas towards the left wall the flow changes between nearly laminar and weak turbulent conditions, thus exhibiting transitional behavior.

- 31 -

ro of

re

-p

Fig. 4-15: Time averaged axial velocity u̅x in domain 2 on yz-plane at x = 40 mm downstream of the monolith for case Re42000. Symmetry conditions are at the bottom and right sides, walls are at the left and top sides. View is into direction of monolith (negative x direction).

Jo

ur

na

lP

The comparison of the calculated time averaged velocity u̅x at x = 40 mm with the experimental measurements is plotted in Fig. 4-16. It can be seen that the coincidence is nearly perfect for values of |y| > 10 mm and better than the RANS results (e. g. Fig. 4-8). Near the symmetry axis (|y|=0 mm) the influence of symmetry boundary conditions can be readily identified. This leads to a distinct maximum directly on the symmetry axis of the system. The “mixture” of the central jet is hindered because no fluctuations are possible. From the computed results extracted for three different (radial) z-positions it can be seen that the two dimensional nature of the flow (compare with Fig. 4-6) is not fulfilled near the symmetry plane. An interesting fact is that the maximum of the measured velocity is not on the symmetry axis but at ca. |y| = 5 mm. This characteristic, whereby the peak value of u̅x is found at ca. |y| = 5 mm, is reflected from the calculations and is visible in Fig. 4-16. Currently the reason for this detail cannot be explained.

- 32 -

ro of -p

lP

re

Fig. 4-16: Comparison of calculated and measured profiles of time averaged velocity u̅x 40 mm downstream of the monolith for case Re42000. Measurements were taken along the line z = 0 mm (purple line in Fig. 2-1), while calculations are presented for radial positions z = 0/5/10 mm. Note that both the calculation and data reflect that the peak u̅x value is not along the centerline, but at ca. |y| = 5 mm.

ur

na

We are aware using symmetry boundary conditions for LES modeling is not perfect. Nevertheless, the calculations show interesting details and prove that the usage of RANS modeling does not capture all of the relevant characteristics of this flow system. The question of how the LES turbulence conditions are for the two smaller Reynolds numbers (Re22000 and Re30000) remains open. It can be expected that the transition to turbulence is farther downstream leading to weaker turbulence at the position of the measurements.

Jo

5 Conclusions

A new approach is presented for calculating flow domains with embedded monoliths. Similar to porous body approaches, the pressure loss of the monolith channels is calculated via the Hagen-Poiseuille approach. However, by coupling the flow domains up- and downstream of the monolith with appropriate boundary conditions, the monolith itself can be excluded from the CFD computational domain. The channel inlets and outlets are discretized in the CFD domains upstream and downstream of the monolith. Therefore, the flow field near the monolith can be calculated with much more detail compared to a porous body approach where the microscopic details of the flow cannot be resolved be resolved due to missing geometric information. Furthermore, the presented approach allows the usage of different turbulence - 33 -

models, thermophysical models and numerical parameters in the flow domains up- and downstream of the monolith. The modeling was performed with the CFD code OpenFOAM®, for which a new solver was developed. It is shown that the presented methodology works very well and can lead to valuable insights. The functionality of this approach is successfully validated with flow data from the literature for a prismatic (planar) converter system, including a monolith, with consists of approximately 4400 channels. In this work, the coupling of the domains was validated for isothermal flow.

re

-p

ro of

The calculation of laminar flow in domain 2 and the different characteristics of the flow field versus experiments proves that the flow was turbulent in measurements at the position 40 mm downstream of the monolith. It is surprising that the RANS modeling leads to a good coincidence of calculated and measured profiles of axial velocity downstream the monolith since the transition from laminar to turbulent conditions is not well captured by these methods. This behavior is caused by the fact that the flow is dominated by the flow distribution coming out of the monolith and the absence of curvature of the streaklines downstream (i.e., poor mixing in the downstream section). The results are not sensitive to the choice of the turbulence model (k- versus k- SST) in domain 2 and the turbulence boundary conditions at the inlets into flow domain 2. These models presume fully turbulent conditions (large Reynolds numbers) which do not exist in this system at the experimentally validated conditions. A strong increase of turbulent kinetic energy directly downstream of the monolith is predicted. The turbulent viscosity seems to be overpredicted in the RANS calculations leading to smoother profiles versus the experiments.

Jo

ur

na

lP

In contrast, the LES calculations show that the transition of the laminar jets leaving the monolith to turbulence takes place about 20 mm downstream of the monolith. The time averaged velocity profiles agree very well with the measurements. Due to the application of symmetry planes, artifacts can be identified in the results. Nevertheless, the LES calculations show interesting details and prove that the usage of RANS modeling is not adequate for flow systems of this type in the range of operational conditions which are typical for converters in vehicles. Currently, LES seems to be the only approach that can adequately deliver details of flow downstream of monoliths. This prediction is essential for determining the exact position of probes in the aftergas treatment system of cars. In future CFD studies of converter systems, it will be important to include the effect of transition from laminar to turbulent flow.

Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

- 34 -

ACKNOWLEDGMENTS

Jo

ur

na

lP

re

-p

ro of

Electronic data tables with the experimental results were graciously supplied by Professor Svetlana Aleksandrova of Convetry University. We appreciate the permission to use data from (Porter et al., 2014; Porter et al., 2016) in this work. This work was performed on the supercomputer ForHLR-I funded by the Ministry of Science, Research and the Arts BadenWürttemberg and by the Federal Ministry of Education and Research.

- 35 -

REFERENCES Agrawal, G., Kaisare, N. S., Pushpavanam, S., & Ramanathan, K. (2012). Modeling the effect of flow mal-distribution on the performance of a catalytic converter. Chemical Engineering Science, 71, 310-320. Benjamin S. F., Haimad N., Roberts C. A., & J., W. (2001). Modelling the flow distribution through automotive catalytic converters. Proc Instn Mech Engrs, 215, 5. Benjamin, S. F., Clarkson, R. J., Haimad, N., & Girgis, N. S. (1996). An Experimental and Predictive Study of the Flow Field in Axisymmetric Automotive Exhaust Catalyst Systems. In: SAE International.

ro of

Benjamin, S. F., Liu, Z., & Roberts, C. A. (2004). Automotive catalyst design for uniform conversion efficiency. Applied Mathematical Modelling, 28, 559-572. Benzinger, W., Daymo, E., Hettel, M., Maier, L., Antinori, C., Pfeifer, P., & Deutschmann, O. (2019). Reverse water gas shift (RWGS) over Ni – Spatially-resolved measurements and simulations. Chemical Engineering Journal, 362, 430-441.

-p

Bertrand, F., Devals, C., Vidal, D., Préval, C. S. d., & Hayes, R. E. (2012). Towards the simulation of the catalytic monolith converter using discrete channel-scale models. Catalysis Today, 188, 80-86.

re

Chakravarthy, V. K., Conklin, J. C., Daw, C. S., & D’Azevedo, E. F. (2003). Multidimensional simulations of cold-start transients in a catalytic converter under steady inflow conditions. Applied Catalysis A: General, 241, 289-306.

lP

Choudary, C., & Mazumder, S. (2014). Direct numerical simulation of catalytic combustion in a multi-channel monolith reactor using personal computers with emerging architectures. Computers & Chemical Engineering, 61, 175-184.

na

Clarkson, R. J. (1995). A Theoretical and Experimental Study of Automotive Catalytic Converters Unpublished PhD Thesis, Coventry University, Coventry UK.

ur

Cornejo, I., Hayes, R. E., & Nikrityuk, P. (2018). A new approach for the modeling of turbulent flows in automotive catalytic converters. Chemical Engineering Research and Design, 140, 308-319.

Jo

Cornejo, I., Nikrityuk, P., & Hayes, R. E. (2018a). Multiscale RANS-based modeling of the turbulence decay inside of an automotive catalytic converter. Chemical Engineering Science, 175, 377-386. Cornejo, I., Nikrityuk, P., & Hayes, R. E. (2018b). Turbulence generation after a monolith in automotive catalytic converters. Chemical Engineering Science, 187, 107-116. Deur M., Jonnavithula S., Dhanapalan S., Schulz K., Raghunathan B., Nakla H., Meeks E., & P., C. C. (2002). Simulation of Engine Exhaust Aftertreatment with CFD using Detailed Chemistry. In Proc. 12th International Multidimensional Engine Modeling User’s Group, . Engine Research Center, Detroit, MI, USA.

- 36 -

Deutschmann, O., Tischer, S., Correa, C., Chatterjee, D., Kleditzsch, S., Janardhanan, V. M., Mladenov, N., Minh, H. D., Karadeniz, H., & Hettel, M. (2014). DETCHEM Software package. In. Karlsruhe, Germany: www.detchem.com. Eggenschwiler, P. D., Tsinoglou, D. N., J., S., C., B., U., V., & M., G. (2009). Ceramic foam substrates for automotive catalyst applications: fluid mechanic analysis. Exp Fluids, 47, 209–222. Espinoza D. E., Scanlon T. J., & E., B. R. (2015). Validation of Tools to Accelerate HighSpeed CFD Simulations Using OpenFOAM. Proceedings of 20th AIAA International Space Planes and Hypersonic Systems and Technologies Conference (Hypersonics 2015), AIAA 2015-3566.

ro of

Guojiang, W., & Song, T. (2005). CFD simulation of the effect of upstream flow distribution on the light-off performance of a catalytic converter. Energy Conversion and Management, 46, 2010-2031. Hayes, R. E., Fadic, A., Mmbaga, J., & Najafi, A. (2012). CFD modelling of the automotive catalytic converter. Catalysis Today, 188, 94-105.

-p

Hettel, M., Daymo, E., & Deutschmann, O. (2018). 3D modeling of a CPOX-reformer including detailed chemistry and radiation effects with DUO. Computers & Chemical Engineering, 109, 166-178.

re

Hettel, M., Diehm, C., Bonart, H., & Deutschmann, O. (2015). Numerical simulation of a structured catalytic methane reformer by DUO: The new computational interface for OpenFOAM® and DETCHEM™. Catalysis Today, 258, 230-240.

lP

Holmgren, A., Grönstedt, T., & Andersson, B. (1997). Improved flow distribution in automotive monolith converters. Reaction Kinetics Catalysis Letters, 60, 8.

na

Huang, S., & Li, Q. S. (2010). A new dynamic one-equation subgrid-scale model for large eddy simulations. Int. J. Numer. Meth. Engng, 81, 835–865.

ur

Hwang, K., Lee, K., Mueller, J., Stuecken, T., Schock, H. J., & Lee, J.-C. (1995). Dynamic Flow Study in a Catalytic Converter Using Laser Doppler Velocimetry and High Speed Flow Visualization. In: SAE International.

Jo

Ibrahim, H. A., Ahmed, W. H., Abdou, S., & Blagojevic, V. (2018). Experimental and numerical investigations of flow through catalytic converters. International Journal of Heat and Mass Transfer, 127, 546-560. Jeong, S.-J. (2014). A full transient three-dimensional study on the effect of pulsating exhaust flow under real running condition on the thermal and chemical behavior of closedcoupled catalyst. Chemical Engineering Science, 117, 18-30. Jeong, S., & Kim, T. (1997). CFD Investigation of the 3-Dimensional Unsteady Flow in the Catalytic Converter. In: SAE International. Karvounis, E., & Assanis, D. N. (1993). The effect of inlet flow distribution on catalytic conversion efficiency. International Journal of Heat and Mass Transfer, 36, 14951504. - 37 -

Kim, W. a. M., S. . (1995). A new dynamic one-equation subgrid-scale model for large eddy simulation. In 33rd Aerospace Sciences Meeting and Exhibit, Reno, NV, AIAA 950356. Kumar, A., & Mazumder, S. (2010). Toward simulation of full-scale monolithic catalytic converters with complex heterogeneous chemistry. Computers & Chemical Engineering, 34, 135-145. Lai, M. C., Lee, T., Kim, J. Y., Cheng, C. Y., Li, P., & Chui, G. (1992). Numerical and experimental characterizations of automotive catalytic converter internal flows. Journal of Fluids and Structures, 6, 451-470. Launder, B. E., & Spalding, D. B. (1974). The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering, 3, 269-289.

ro of

Martin, A. P., Will, N. S., Bordet, A., Cornet, P., Gondoin, C., & Mouton, X. (1998). Effect of Flow Distribution on Emissions Performance of Catalytic Converters. In: SAE International.

-p

Mat Yamin, A. K., Benjamin, S. F., & Roberts, C. A. (2013). Pulsating flow in a planar diffuser upstream of automotive catalyst monoliths. International Journal of Heat and Fluid Flow, 40, 43-53.

re

Mei, H., Li, C., Liu, H., & Ji, S. (2006). Simulation of Catalytic Combustion of Methane in a Monolith Honeycomb Reactor1 1Supported by the National Natural Science Foundation of China (No.20136010 and No.20376005). Chinese Journal of Chemical Engineering, 14, 56-64.

lP

Menter, F. (1994). Two-equation eddy-viscosity models for engineering applications. AIAA J., 32 (8), 1598-1605.

na

Om Ariara Guhan, C. P., Arthanareeswaran, G., Varadarajan, K. N., & Krishnan, S. (2016). Numerical optimization of flow uniformity inside an under body- oval substrate to improve emissions of IC engines. Journal of Computational Design and Engineering, 3, 198-214.

ur

OpenFOAM. (2017). OpenFOAM-The Open Source CFD Toolbox. In: www.openfoam.org. Ozhan, C., Fuster, D., & Da Costa, P. (2014). Multi-scale flow simulation of automotive catalytic converters. Chemical Engineering Science, 116, 161-171.

Jo

Porter, S., Mat Yamin, A. K., Aleksandrova, S., Benjamin, S., Roberts, C. A., & Saul, J. (2014). An Assessment of CFD Applied to Steady Flow in a Planar Diffuser Upstream of an Automotive Catalyst Monolith. SAE International Journal of Engines, 7, 16971704. Porter, S., Saul, J., Aleksandrova, S., Medina, H., & Benjamin, S. (2016). Hybrid flow modelling approach applied to automotive catalysts. Applied Mathematical Modelling, 40, 8435-8445. Purday, H. F. P. (1949). An Introduction to the Mechanics of Viscous Flow. DoverPublications, New York, 16. - 38 -

Ramadan, B. H., Lundberg, P. C., & Richmond, R. P. (2007). Characterization of a Catalytic Converter Internal Flow. In: SAE International. Rodi, W. (1993). Turbulence Models and their Applications in Hydraulics. A State-of-Arts Review, International Association for Hydraulic Research Monograph, Rotterdam, The Netherlands, 3rd edition. Shi-jin, S., Jian-xin, W., Ren-jun, Z., & Jun-rui, C. (2000). Study on Flow Characteristics of Automotive Catalytic Converters with Various Configurations. In: SAE International. Shijin, S., Jianxin, W., & Renjun, Z. (2000). Numerical Simulation and Optimum Design of Automotive Catalytic Converters. In FISISTA World Automotive Conference. Seoul, Korea.

ro of

Štěpánek, J., Kočí, P., Kubíček, M., Plát, F., & Marek, M. (2011). Spatially 3D simulation of a catalytic monolith by coupling of 1D channel model with CFD. In E. N. Pistikopoulos, M. C. Georgiadis & A. C. Kokossis (Eds.), Computer Aided Chemical Engineering (Vol. 29, pp. 141-145): Elsevier.

-p

Štěpánek, J., Kočí, P., Marek, M., & Kubíček, M. (2012). Catalyst simulations based on coupling of 3D CFD tool with effective 1D channel models. Catalysis Today, 188, 8793.

re

Su, Q., Xie, L., Shuai, S., Wang, J., Song, J., & Li, Z. (2013). Optimization of automotive catalytic converter by numerical modeling and simulation with detailed mechanism. Catalysis Today, 216, 292-298.

lP

Sui, R., Prasianakis, N. I., Mantzaras, J., Mallya, N., Theile, J., Lagrange, D., & Friess, M. (2016). An experimental and numerical investigation of the combustion and heat transfer characteristics of hydrogen-fueled catalytic microreactors. Chemical Engineering Science, 141, 214-230.

na

Tischer, S., Correa, C., & Deutschmann, O. (2001). Transient three-dimensional simulations of a catalytic combustion monolith using detailed models for heterogeneous and homogeneous reactions and transport phenomena. Catalysis Today, 69, 57-62.

ur

Tsinoglou, D. N., Koltsakis, G. C., Missirlis, D. K., & Yakinthos, K. J. (2004). Transient modelling of flow distribution in automotive catalytic converters. Applied Mathematical Modelling, 28, 775-794.

Jo

Weltens, H., Bressler, H., Terres, F., Neumaier, H., & Rammoser, D. (1993). Optimisation of Catalytic Converter Gas Flow Distribution by CFD Prediction. In: SAE International. Will, N. S., & Bennett, C. J. (1992). Flow Maldistributions in Automotive Converter Canisters and their Effect on Emission Control. In: SAE International. Wörner, M. (2010). Approximate residence time distribution of fully develop laminar flow in a straight rectangular channel. Chemical Engineering Science, 65, 3499-3507. Yoshizawa, A., & Horiuti, K. (1985). A Statistically-Derived Subgrid-Scale Kinetic Energy Model for the Large-Eddy Simulation of Turbulent Flows. Journal of the Physical Society of Japan, 54 (8), 2834-2839. - 39 -

Zhang, X., & Tennison, P. (2008). Numerical Study of Flow Uniformity and Pressure Loss Through a Catalytic Converter with Two Substrates. In: SAE International.

Jo

ur

na

lP

re

-p

ro of

Zygourakis, K. (1989). Transient operation of monolith catalytic converters: a twodimensional reactor model and the effects of radially nonuniform flow distributions. Chemical Engineering Science, 44, 2075-2086.

- 40 -

Table 2-1: Experimental Parameters Re [-] 22,000 30,000 42,000

uinlet [m/s] 8.79 12.00 16.80

GHSV [1/h] 97,000 150,000 190,000

ro of

Case Re22000 Re30000 Re42000

Table 3-1: Pairs of grids for domains 1 and 2 used in combination. ncell is the cell count divided by 106. lcell is the characteristic length scale of the cells in [mm].

very fine 5 0.27 74 0.12

-p

fine 10 0.22 16 0.20

na

lP

domain 2

ncell lcell ncell lcell

medium 5 0.27 7 0.25

re

domain 1

coarse 2.5 0.35 3 0.33

Jo

Case Re22000 Re30000 Re42000

ur

Table 3-2: Boundary conditions at the inlet of domain 1 uinlet [m/s] 8,79 12,00 16,80

k [m2/s2] 0.0116 0.0216 0.0432

- 41 -

 [1/s] 77 144 282

Table 3-3: Turbulence boundary conditions at the 1100 inlets of domain 2 smallTurbulence k  2 2 [m /s ] [1/s] 0.0011 7.13 0.0020 13.48 0.004 26.20

ε [1/s] 1E+05 1E+05 1E+05

ε [1/s] 0.0092

ro of

uref [m/s] Re22000 2.7 Re30000 3.7 Re42000 5.2 Case

noTurbulence k  2 2 [m /s ] [1/s] 1E-07 1E+05 1E-07 1E+05 1E-07 1E+05

Table 3-4: Typical time values and number of timesteps for LES calculations t [s] 5.0E-06

tchar [s] 0.031

nchar [-] 6200

na

lP

re

-p

Case Re42000

Table 4-1: Reynolds numbers in and downstream of the channels of the monolith uref [m/s] Re22000 2.7 Re30000 3.7 Re42000 5.2

Jo

ur

Case

u̅channel [m/s] 3.5 4.8 6.75

- 42 -

Rechannel [-] 261 358 504

Rewall [-] 49 68 96