Developing Computational Fluid Dynamics ( CFD ) Models to Evaluate Available Energy in Exhaust Systems of Diesel Light-Duty Vehicles

Around a third of the energy input in an automotive engine is wasted through the exhaust system. Since numerous technologies to harvest energy from exhaust gases are accessible, it is of great interest to find timeand cost-efficient methods to evaluate available thermal energy under different engine conditions. Computational fluid dynamics (CFD) is becoming a very valuable tool for numerical predictions of exhaust flows. In this work, a methodology to build a simple three-dimensional (3D) model of the exhaust system of automotive internal combustion engines (ICE) was developed. Experimental data of exhaust gas in the most used part of the engine map in passenger diesel vehicles were employed as input for calculations. Sensitivity analyses of different numeric schemes have been conducted in order to attain accurate results. The model built allows for obtaining details on temperature and pressure fields along the exhaust system, and for complementing the experimental results for a better understanding of the flow phenomena and heat transfer through the system for further energy recovery devices.


Introduction
Around a third of the fuel energy in a light-duty diesel engine is wasted through the exhaust system [1].Rising awareness of environmental issues together with the fuel economy have encouraged research on energy recovery [1,2] from exhaust gas.An obvious first step for energy recovery is evaluating the heat source.The nature of exhaust flow changes during engine operation.Hence, there is a need for models that allow evaluating exhaust gas in a cost-and time-efficient manner under different engine conditions.
Computational fluid dynamics (CFD) is becoming a very popular tool for numerical predictions in the automotive industry.CFD models are employed in five major areas: vehicle aerodynamics, thermal management (cooling and climate control), cylinder combustion, engine lubrication and exhaust system performance.The main applications of CFD in exhaust systems are: flow distribution in the front of the monolith, pressure loss through exhaust system, skin temperature prediction, heat loss analyses and emission-related studies [3].
High levels of energy lost through engine exhaust have brought attention to heat transfer in exhaust systems and how to model associated phenomena.
In Konstantinidis et al. [4], the status of knowledge regarding heat transfer phenomena in automotive exhaust systems was summarized, and a transient model covering diverse exhaust piping configurations was presented.Using this model, Kandylas et al. [5] studied heat transfer in automotive exhaust pipes for steady and transient conditions.It is stated that the code is suitable to support a complete and efficient methodology for design optimization.
Not only are exhaust heat transfer and temperature models important for energy recovery, but also for on-board control and diagnosis of the after-treatment system.In the work of Guardiola et al. [6], diesel oxidation catalyst (DOC) inlet temperature was modelled depending on fuel mass flow and engine speed using an engine map look-up table.
Traditionally, one-dimensional (1D) simulations using in-house or commercial codes have been used to study exhaust flow.In Shayler et al. [7], a 1D model of system thermal behavior was developed.The exhaust system is modelled as connected pipe and junction elements with lumped thermal capacities.Heat transfer correlations for quasi-steady and transient conditions were investigated.The computational model supports studies of exhaust and after-treatment system design, the investigation of thermal conditions, and performance characteristics.Model predictions and experimental data were in good agreement.
In Kapparos et al. [8], a 1D heat transfer model in an automotive diesel exhaust is presented.The external natural convective heat transfer coefficient was given by Churchill and Chu's correlation.A sensitivity analysis of the main variables (heat transfer coefficient, external pipe emissivity and others) affecting heat transfer in exhaust pipes was carried out.Fu et al. [9] developed a 1D model of an exhaust pipe (without the catalytic converter).Apart from typical variables studied, effects of pipe geometry and flow regime were investigated.
One-dimensional tools have the advantage of shorter computational time, but the accuracy of the results is limited by their inability to simulate 3D flow effects in regions such as the inlet and outlet cone.Usually strong turbulent flow exists when angled and asymmetrical cones are presented [3].
Three different approaches were used in Fortunato et al. [10] to simulate the exhaust and underbody of the car: 1D, 3D finite elements and 3D complete thermo-fluid-dynamics.Both 3D models allowed to detect the most critical points in terms of higher temperatures In Chuchnowski et al. [11], temperature distribution and heat flow analyses on the exhaust system of a mining diesel engine to determine technical parameters of a heat recuperation system are presented.Three-dimensional numerical simulations show it is possible to develop a suitable design of an energy recovery device.
Simulations in a coupled three-dimensional thermal model including the underhood, underbody and exhaust systems of a vehicle were conducted by Guoquan et al. [12].Pulsating flow and steady flow effects were compared.Pulsating exhaust showed an enhancement of more than a 10% in heat transfer.Influence of this effect was higher at the final components of the exhaust, such as the muffler.
Research has also been focused on finding optimal position for exhaust components.Durat et al. [13] compared experimental and CFD heat transfer results in an exhaust system of a spark-ignition engine.An optimal position of a close-coupled catalytic converter in terms of light-off time was found using the CFD model.Liu [14] developed a 3D exhaust model to study the optimal position of a heat recovery device.Three different simulations were carried out varying the position of a thermoelectric generator along the system.Higher surface temperature and lower backpressure were obtained when the recovery device was placed between catalytic converter and muffler.The influence of materials and insulation on heat transfer in catalytic converters has also been discussed employing a 3D model [15].
In the same way preceding authors have employed CFD simulations successfully to find an optimal position for exhaust components or to assist in their design, CFD can also be used to gain knowledge of the relevant characteristics of the flow in terms of its usage in energy recovery.Previous works [16] state that knowledge of heat transfer within the gas is critical to enhance harvested energy.
In exhaust systems, not only should heat transfer to surroundings be reproduced but also chemical reactions in catalytic processes, since after-treatment devices can modify gas temperature.Most mathematical kinetics models rely on the Langmuir-Hinshelwood expressions derived on pellet-type Pt catalysts for the CO and HC oxidation with modified parameters (activation energies and activity factors), best suited to the case modeled.This approach is widely employed because it provides acceptable accuracy in the most common automobile range of operation.
Measurements on pellet-type platinum-alumina catalysts have been processed [17] to derive kinetic rate expressions for the oxidation reactions of CO and C 3 H 6 under oxygen-rich conditions for a better understanding of platinum oxidation catalysts in automotive emission control systems.Terms accounting for CO, C 3 H 6 and NO inhibition were included.The isothermal reactor approach in a numerical integration-optimization method was used to fit the kinetic parameters.These expressions, especially as written by Oh et al. [18], have been widely used [19][20][21][22].This approach has also been used to build simpler kinetic models over a range of interest [23].
Detailed kinetics for a Pt/Rh three-way catalyst were described in Chatterjee et al. [24].The model consists of 60 elementary reaction steps and one global reaction step, involving 8 gas species and 23 site species.It considers the steps of adsorption of the reactants on the surface, reaction of the adsorbed species, and desorption of the reaction products.This approach includes a mechanism of C 3 H 6 oxidation on Pt/Al 2 O 3 , a mechanism of NO reduction on Pt, and a mechanism for NO reduction and CO oxidation on Rh.This model has been used when more accuracy is required, as in Kumar et al. [25].
Nevertheless, due to the variety of washcoat materials and loadings of the monoliths and to ageing effects of after-treatment devices, the kinetic parameters published have difficulties in fitting to conditions or catalytic converters different from the study cases.To overcome the difficulty of finding kinetic parameters fitted to their application, it would be of interest to have a simple method for researchers to develop their own kinetic constants for CFD exhaust models.
The real physical space of a full-size monolith converter contains thousands of tiny channels and could present an extremely time-consuming problem to solve.There are several approaches to modeling the monolith of a catalytic converter.For simulations of the converter performance, most models use a volume-averaging method and treat the monolith as a continuous porous medium [23][24][25][26].With the porous medium, approach channels are not simulated but a macroscopic understanding of the flow is achieved while keeping the computational requirements affordable.
The scope of this work is to develop a combined theoretical and experimental methodology to build an exhaust model to obtain information about temperature and heat losses along the system.Results were used to provide an insight into those flow characteristics that must be taken into account when harvesting energy, especially but not exclusively for thermoelectric modules.The model does not only intend to reproduce the behavior of the exhaust system at a single operation point but in the whole most-used part of the engine map.This information is needed to evaluate the exhaust gas energy available for recovery under real operating conditions.

Engine
Tests were carried out in a test-bench with a Nissan YD22 four-stroke, turbocharged, four-cylinder diesel engine.The engine bore and stroke are, respectively, 86 mm and 130 mm, and the compression ratio is 16.7:1.The exhaust system (Figure 1) is equipped with a diesel oxidation catalyst (DOC) and a muffler.

Measurement Devices
Piezo resistive pressure sensors and K-type thermocouples were used for measuring pressure and temperature of the gas along the exhaust system (see Figure 2).The exhaust mass flow rate was calculated from the addition of fuel and air mass flow rates.Skin temperatures used in calculations were provided by the IR camera Gobi-384 GigE.
Pollutants are measured using an ENVIRONNEMENT manufacturer equipment.A TOPAZE 32M model was used as NOx analyser, based on the chemiluminescence effect from NO oxidation by ozone (O3).The GRAPHITE 52M gas analyser measures total hydrocarbons (THCs) by flame ionization detection while a MIR 2R gas analyser measures CO and CO2 species, detecting the molecules absorption in the infrared spectrum.Other species compositions are estimated via chemical balances from fuel and air consumption.

Test Plan
The velocity profile imposed by the New European Driving Cycle (NEDC) used for light-duty vehicle certification was translated into engine operating conditions (torque and engine speed), as shown in Figure 3 (black squares), employing longitudinal dynamics equations [27,27].Then, a matrix of nine steady-state modes (see Figure 3, circles) covering the most-used quarter of the engine map under both urban and extra-urban NEDC conditions was selected for testing.

Measurement Devices
Piezo resistive pressure sensors and K-type thermocouples were used for measuring pressure and temperature of the gas along the exhaust system (see Figure 2).The exhaust mass flow rate was calculated from the addition of fuel and air mass flow rates.Skin temperatures used in calculations were provided by the IR camera Gobi-384 GigE.

Measurement Devices
Piezo resistive pressure sensors and K-type thermocouples were used for measuring pressure and temperature of the gas along the exhaust system (see Figure 2).The exhaust mass flow rate was calculated from the addition of fuel and air mass flow rates.Skin temperatures used in calculations were provided by the IR camera Gobi-384 GigE.
Pollutants are measured using an ENVIRONNEMENT manufacturer equipment.A TOPAZE 32M model was used as NOx analyser, based on the chemiluminescence effect from NO oxidation by ozone (O3).The GRAPHITE 52M gas analyser measures total hydrocarbons (THCs) by flame ionization detection while a MIR 2R gas analyser measures CO and CO2 species, detecting the molecules absorption in the infrared spectrum.Other species compositions are estimated via chemical balances from fuel and air consumption.

Test Plan
The velocity profile imposed by the New European Driving Cycle (NEDC) used for light-duty vehicle certification was translated into engine operating conditions (torque and engine speed), as shown in Figure 3 (black squares), employing longitudinal dynamics equations [27,27].Then, a matrix of nine steady-state modes (see Figure 3, circles) covering the most-used quarter of the engine map under both urban and extra-urban NEDC conditions was selected for testing.Pollutants are measured using an ENVIRONNEMENT manufacturer equipment.A TOPAZE 32M model was used as NO x analyser, based on the chemiluminescence effect from NO oxidation by ozone (O 3 ).The GRAPHITE 52M gas analyser measures total hydrocarbons (THCs) by flame ionization detection while a MIR 2R gas analyser measures CO and CO 2 species, detecting the molecules absorption in the infrared spectrum.Other species compositions are estimated via chemical balances from fuel and air consumption.

Test Plan
The velocity profile imposed by the New European Driving Cycle (NEDC) used for light-duty vehicle certification was translated into engine operating conditions (torque and engine speed), as shown in Figure 3 (black squares), employing longitudinal dynamics equations [27,28].Then, a matrix of nine steady-state modes (see Figure 3, circles) covering the most-used quarter of the engine map under both urban and extra-urban NEDC conditions was selected for testing.The classical approach that links pressure drop and velocity in flows through porous media is the Forchheimer equation (Equation ( 1)), that can be derived from the Navier-Stokes equation for one-dimensional, incompressible and steady laminar flow of a Newtonian fluid in a rigid porous medium [29]: The second term in the right can be interpreted as a second-order correction to account for the contribution of inertial forces, but, at sufficiently low velocities, this effect is negligible and Equation (1) can be reduced to Darcy's law (Equation ( 2)) [30]: Darcy's law can be rewritten in the form of Equation (3), relating the average fluid velocity through the pores with pressure drop Δ along a segment of length .
If the equation from the linear fit of the scatter plot of − vs. is (Equation ( 4)): then Darcy's constant can be derived as follows (Equation ( 5)): (b) Heat transfer to surroundings.
Natural convection to pipe surroundings in the test-bench needs to be included.Wall temperatures along the exhaust pipe were measured using thermography.Figure 4 shows, as an example, views of the external surface of the pipe and the DOC.The classical approach that links pressure drop and velocity in flows through porous media is the Forchheimer equation (Equation ( 1)), that can be derived from the Navier-Stokes equation for one-dimensional, incompressible and steady laminar flow of a Newtonian fluid in a rigid porous medium [29]: The second term in the right can be interpreted as a second-order correction to account for the contribution of inertial forces, but, at sufficiently low velocities, this effect is negligible and Equation ( 1) can be reduced to Darcy's law (Equation ( 2)) [30]: Darcy's law can be rewritten in the form of Equation ( 3), relating the average fluid velocity v through the pores with pressure drop ∆p along a segment of length L.
If the equation from the linear fit of the scatter plot of − ∆p µL vs. v is (Equation ( 4)): then Darcy's constant κ can be derived as follows (Equation ( 5)): (b) Heat transfer to surroundings.
Natural convection to pipe surroundings in the test-bench needs to be included.Wall temperatures along the exhaust pipe were measured using thermography.Figure 4    An average skin temperature was used to calculate an effective convection heat transfer coefficient.This skin temperature simplification is considered as reasonable, since calculation of heat transfer rates is not very sensitive in the values of pipe wall temperatures [5].
Total transferred heat from pipe wall to surroundings is (W) (Equation ( 6)): Here, all terms in the right side of Equation ( 6) are grouped in one single equivalent convection term (Equation ( 7)): The effective convection coefficient accounts for convection and radiation heat losses, although radiation heat transfer to the environment is expected to be low, since wall temperature is below 400 °C [4].A similar approach is explained in Kapparos et al. [8].
Based on the measured exhaust gas and pipe wall temperatures, an energy balance for the exhaust gas is employed for the calculation of the heat flux.The resulting heat fluxes are employed in the estimation of a mean gas-to-wall convection coefficient.
Since temperature at both ends of the exhaust pipe are measured, total heat losses (W) can be obtained as Equation ( 8 Equaling Equations ( 7) and ( 8), ℎ can be derived as Equation ( 9): (c) Estimation of kinetic parameters.
Under light-off temperatures (about 200 °C) catalytic processes remain basically inactive.Until activation temperature is reached, reactions are chemically controlled.On the other hand, post-light-off reaction rates are limited mainly by mass transfer and, consequently, conversion efficiency depends on the residence time within the monolith, the surface to volume ratio of the monolith and the mass transfer [31].In operating conditions from test matrix, catalytic reactions are normally within the light-off band.Thus, estimation of kinetic parameters is needed.
Two main reactions occurring in DOCs, as in Voltz et al. [17], are modeled (Equations ( 10) and ( 11)): An average skin temperature was used to calculate an effective convection heat transfer coefficient.This skin temperature simplification is considered as reasonable, since calculation of heat transfer rates is not very sensitive in the values of pipe wall temperatures [5].
Here, all terms in the right side of Equation ( 6) are grouped in one single equivalent convection term (Equation ( 7)): .
The effective convection coefficient accounts for convection and radiation heat losses, although radiation heat transfer to the environment is expected to be low, since wall temperature is below 400 • C [4].A similar approach is explained in Kapparos et al. [8].
Based on the measured exhaust gas and pipe wall temperatures, an energy balance for the exhaust gas is employed for the calculation of the heat flux.The resulting heat fluxes are employed in the estimation of a mean gas-to-wall convection coefficient.
Since temperature at both ends of the exhaust pipe are measured, total heat losses .Q (W) can be obtained as Equation ( 8 Equaling Equations ( 7) and (8), h e f f can be derived as Equation ( 9): (c) Estimation of kinetic parameters.
Under light-off temperatures (about 200 • C) catalytic processes remain basically inactive.Until activation temperature is reached, reactions are chemically controlled.On the other hand, post-light-off reaction rates are limited mainly by mass transfer and, consequently, conversion efficiency depends on the residence time within the monolith, the surface to volume ratio of the monolith and the mass transfer [31].In operating conditions from test matrix, catalytic reactions are normally within the light-off band.Thus, estimation of kinetic parameters is needed.
Two main reactions occurring in DOCs, as in Voltz et al. [17], are modeled (Equations ( 10) and ( 11)): Appl.Sci.2017, 7, 590 Propylene is representative of the easily oxidized hydrocarbons, which constitute about 80% of the total hydrocarbons found in a typical exhaust gas.Other saturated hydrocarbons (typically represented as methane or propane) that are resistant to oxidation usually make up the remaining 20% [17].Since only total hydrocarbon (THC) data is available and fast oxidation hydrocarbons represent the majority of the THC, propylene alone was used as representative of the HC (as is common [32]).
Usually adopted forms of the rates of CO and HC oxidations are as follows (Equations ( 12) and ( 13)): where G is a term accounting for NO, CO and HC inhibition effects on oxidation (Equation ( 14)): Since the intended target is an empirical model to fit experimental data and power-law reactions showed good performance, an inhibition term was not included (Equations ( 15) and ( 16)): Arrhenius kinetic constants k m are defined as (Equation ( 17)): The problem is reduced to obtain the pre-exponential factors A m and activation energies E a,m .For both species, an iterative process varying the pre-exponential factor in CFD simulations is conducted until the deviation of outlet concentration value from the experimental is as much as 0.1%.Once this is achieved, values of k m are obtained.This tuning process is done for the different operating conditions, and then k m of each reaction is obtained as a function of monolith temperature (provided by CFD results).This iterative process needs to be followed just once (to obtain the above-mentioned constants and implement them in the model) and not for every prospective 3D simulation.
The simulations were performed on a 2D axisymmetric DOC model.The 2D model maintained the same characteristic lengths but with a simplified geometry (see Figure 5), allowing fast iterative simulations.Other aspects about set-up of these 2D simulations are the same than those presented below for the 3D calculations.
Appl.Sci.2017, 7, 590 7 of 20 Propylene is representative of the easily oxidized hydrocarbons, which constitute about 80% of the total hydrocarbons found in a typical exhaust gas.Other saturated hydrocarbons (typically represented as methane or propane) that are resistant to oxidation usually make up the remaining 20% [17].Since only total hydrocarbon (THC) data is available and fast oxidation hydrocarbons represent the majority of the THC, propylene alone was used as representative of the HC (as is common [32]).
Usually adopted forms of the rates of CO and HC oxidations are as follows (Equations ( 12) and ( 13)): where G is a term accounting for NO, CO and HC inhibition effects on oxidation (Equation ( 14)): Since the intended target is an empirical model to fit experimental data and power-law reactions showed good performance, an inhibition term was not included (Equations ( 15) and ( 16)): Arrhenius kinetic constants are defined as (Equation ( 17)): The problem is reduced to obtain the pre-exponential factors and activation energies , .For both species, an iterative process varying the pre-exponential factor in CFD simulations is conducted until the deviation of outlet concentration value from the experimental is as much as 0.1%.Once this is achieved, values of are obtained.This tuning process is done for the different operating conditions, and then of each reaction is obtained as a function of monolith temperature (provided by CFD results).This iterative process needs to be followed just once (to obtain the abovementioned constants and implement them in the model) and not for every prospective 3D simulation.
The simulations were performed on a 2D axisymmetric DOC model.The 2D model maintained the same characteristic lengths but with a simplified geometry (see Figure 5), allowing fast iterative simulations.Other aspects about set-up of these 2D simulations are the same than those presented below for the 3D calculations.From results, the kinetic exponential law (Equation ( 17)) can be obtained, since (Equation ( 18)): where C i and D i are constants from the linear fit.Sought kinetic constants are (Equations ( 19) and ( 20)): The numerical solution of the continuity, momentum, energy and species equations was computed using a CFD proprietary code (ANSYS Fluent 16), based on the finite volume method.In this work, a steady state, three-dimensional, viscous, turbulent and incompressible (since the maximum Mach number is below 0.3) flow was assumed.Pressure-velocity coupling was taken care of by the segregated, pressure-based solver, the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm [33].Summary details of spatial discretization are presented in Table 1.The convergence criteria were set to 10 −6 for the thermal energy and chemical residuals and 10 −4 for residuals from mass, momentum, turbulence kinetic energy and turbulence energy dissipation rate.Relevant physical and chemical quantities were monitored to assure convergence.The governing equations solved are the continuity equation (Equation ( 21)), momentum equation (Equation ( 22)), energy equation (Equation ( 23)) and conservation equation for chemical species (Equation ( 24)): ∂ρ ∂t Notice that for steady simulations, time-dependent terms become zero.

(b) Turbulence
A laminar regime was forced in the monolith because of small Reynolds numbers due to low velocity and very small hydraulic diameter inside channels.
In the other parts of the domain, Reynolds-averaged Navier-Stokes (RANS) approach was employed with the realizable k − ε model [34] selected.This model has the same turbulent kinetic energy equation as the standard k − ε model but holds an improved equation for ε.Compared to standard k − ε, it shows better performance for flows involving: planar and round jets (predicts round jet spreading correctly), boundary layers under strong adverse pressure gradients or separation, rotation, recirculation or/and strong streamline curvature [34].
(c) Computational domain Simulation domain comprehends the area of interest for energy recovery.Great energy is lost in the muffler [35] and energy harvesting must be performed before it.After-treatment processes, such as chemical reactions within the DOC, might be affected by a temperature change caused by some sort of energy recovery device [36].Consequently, energy recovery should take place between after-treatment devices and the muffler, as pointed out in [14].The simulation domain includes the DOC and the exhaust pipe (see Figure 6).DOC must be included, since exothermal chemical reactions taking place in it can modify gas temperature.In order to obtain a properly-developed velocity profile at the inlet, an extruded entrance zone is added to the physical domain.
Appl.Sci.2017, 7, 590 9 of 20 A laminar regime was forced in the monolith because of small Reynolds numbers due to low velocity and very small hydraulic diameter inside channels.
In the other parts of the domain, Reynolds-averaged Navier-Stokes (RANS) approach was employed with the realizable − model [34] selected.This model has the same turbulent kinetic energy equation as the standard − model but holds an improved equation for .Compared to standard − , it shows better performance for flows involving: planar and round jets (predicts round jet spreading correctly), boundary layers under strong adverse pressure gradients or separation, rotation, recirculation or/and strong streamline curvature [34].
(c) Computational domain Simulation domain comprehends the area of interest for energy recovery.Great energy is lost in the muffler [35] and energy harvesting must be performed before it.After-treatment processes, such as chemical reactions within the DOC, might be affected by a temperature change caused by some sort of energy recovery device [36].Consequently, energy recovery should take place between after-treatment devices and the muffler, as pointed out in [14].The simulation domain includes the DOC and the exhaust pipe (see Figure 6).DOC must be included, since exothermal chemical reactions taking place in it can modify gas temperature.In order to obtain a properly-developed velocity profile at the inlet, an extruded entrance zone is added to the physical domain.(d) Computational mesh Given the nature of the geometry, a hybrid Tet/Hex mesh (see Figure 7) has been used to define the computational domain.Complex 3D features such as inlet and outlet DOC cones have been meshed by means of tetrahedral elements, and hexahedral meshes have been used for the monolith and the duct.In order to ensure the accuracy of the calculations, a grid independence study was conducted and mesh independency was achieved with a 6.5 × 10 5 elements grid with an average mesh size of 3 × 10 −3 m.(d) Computational mesh Given the nature of the geometry, a hybrid Tet/Hex mesh (see Figure 7) has been used to define the computational domain.Complex 3D features such as inlet and outlet DOC cones have been meshed by means of tetrahedral elements, and hexahedral meshes have been used for the monolith and the duct.In order to ensure the accuracy of the calculations, a grid independence study was conducted and mesh independency was achieved with a 6.5 × 10 5 elements grid with an average mesh size of 3 × 10 −3 m.A laminar regime was forced in the monolith because of small Reynolds numbers due to low velocity and very small hydraulic diameter inside channels.
In the other parts of the domain, Reynolds-averaged Navier-Stokes (RANS) approach was employed with the realizable − model [34] selected.This model has the same turbulent kinetic energy equation as the standard − model but holds an improved equation for .Compared to standard − , it shows better performance for flows involving: planar and round jets (predicts round jet spreading correctly), boundary layers under strong adverse pressure gradients or separation, rotation, recirculation or/and strong streamline curvature [34].
(c) Computational domain Simulation domain comprehends the area of interest for energy recovery.Great energy is lost in the muffler [35] and energy harvesting must be performed before it.After-treatment processes, such as chemical reactions within the DOC, might be affected by a temperature change caused by some sort of energy recovery device [36].Consequently, energy recovery should take place between after-treatment devices and the muffler, as pointed out in [14].The simulation domain includes the DOC and the exhaust pipe (see Figure 6).DOC must be included, since exothermal chemical reactions taking place in it can modify gas temperature.In order to obtain a properly-developed velocity profile at the inlet, an extruded entrance zone is added to the physical domain.(d) Computational mesh Given the nature of the geometry, a hybrid Tet/Hex mesh (see Figure 7) has been used to define the computational domain.Complex 3D features such as inlet and outlet DOC cones have been meshed by means of tetrahedral elements, and hexahedral meshes have been used for the monolith and the duct.In order to ensure the accuracy of the calculations, a grid independence study was conducted and mesh independency was achieved with a 6.5 × 10 5 elements grid with an average mesh size of 3 × 10 −3 m.

(e) Boundary conditions
Table 2 shows boundary condition types selected for the model.A commonly employed porous medium approach was used for the monolith, with a one-directional pressure gradient.Boundary conditions values are presented in Table 3.A total amount of seven species are considered in gas composition (see Table 4).Oxidations of carbon monoxide and propylene, as explained in the previous section, are accounted for.A mesh with a special refinement (three layers) near the walls and other without it were tested.Only a part of the domain was used (DOC and a section of the exhaust pipe) in order to reduce computational requirements.As can be seen in Table 5, refinement along the boundary layer adds a considerable number of extra cells and does not provide significant variations.Boundary layer refinement results differed from the base case by less than 1% and required five times more cells.It can be stated that the model is not boundary layer-sensitive and, hence, boundary layer refinement is omitted.
(g) Wall treatment sensitivity analysis.
Three different wall approaches were tested to study their influence in the results of the model: standard, non-equilibrium and enhancement wall functions (Figure 8).
The standard wall functions proposed by Launder and Spalding [37] have been widely used, but the assumption of logarithmic velocity distribution treatment may not be adequate for complex non-equilibrium flows.To overcome this, non-equilibrium wall functions are based on pressure-gradient sensitized Launder and Spalding's [37] log-law for mean velocity.
Enhanced wall treatment is a near-wall modelling method that combines a two-layer model with enhanced wall functions.A one-equation relationship is used to evaluate the laminar sub-layer with fine mesh and transition to log-low function for the turbulent part of the boundary layer.The restriction that the near-wall mesh must be suitably fine might impose large computational requirements.Boundary layer refinement results differed from the base case by less than 1% and required five times more cells.It can be stated that the model is not boundary layer-sensitive and, hence, boundary layer refinement is omitted.
(g) Wall treatment sensitivity analysis.
Three different wall approaches were tested to study their influence in the results of the model: standard, non-equilibrium and enhancement wall functions (Figure 8).
The standard wall functions proposed by Launder and Spalding [37] have been widely used, but the assumption of logarithmic velocity distribution treatment may not be adequate for complex non-equilibrium flows.To overcome this, non-equilibrium wall functions are based on pressure-gradient sensitized Launder and Spalding's [37] log-law for mean velocity.
Enhanced wall treatment is a near-wall modelling method that combines a two-layer model with enhanced wall functions.A one-equation relationship is used to evaluate the laminar sub-layer with fine mesh and transition to log-low function for the turbulent part of the boundary layer.The restriction that the near-wall mesh must be suitably fine might impose large computational requirements.More accurate wall functions such as non-equilibrium or enhanced showed no more than a 3% discrepancy with the standard approach, while requiring eight times more mesh elements.Consequently, the standard wall functions were selected.

Pressure Loss Coefficient
Different exhaust mass flows and their corresponding pressure drops were measured to obtain accurate predictions.Empirical data (Figure 9) show a clear linear correlation between velocity and pressure drop (being the horizontal axis and − the vertical axis) and Darcy's law can be employed.The resulting Darcy's constant is = 1.695 × 10 −8 m 2 .More accurate wall functions such as non-equilibrium or enhanced showed no more than a 3% discrepancy with the standard approach, while requiring eight times more mesh elements.Consequently, the standard wall functions were selected.

Pressure Loss Coefficient
Different exhaust mass flows and their corresponding pressure drops were measured to obtain accurate predictions.Empirical data (Figure 9) show a clear linear correlation between velocity and pressure drop (being v the horizontal axis and − ∆p µL the vertical axis) and Darcy's law can be employed.The resulting Darcy's constant is κ = 1.695 × 10 −8 m 2 .

Convection Heat Transfer Coefficient
Effective convection coefficients derived in each operation mode are presented in Table 6.

Kinetic Parameters
Complete and zero species consumption at DOC are not useful for estimating kinetic parameters.Therefore, data close to 0 or 100% conversion were excluded from calculations.Given that light-off bands seemed to be narrow, particularly for CO, two sets of experiments needed to be used.Due to intended repeatability of the experiments, similar tests would result in similar results, adding no new information.It was decided to repeat the test matrix under very different external ambient conditions, so that different exhaust temperature results in the range of interest are obtained (see Figure 10).

Convection Heat Transfer Coefficient
Effective convection coefficients derived in each operation mode are presented in Table 6.

Kinetic Parameters
Complete and zero species consumption at DOC are not useful for estimating kinetic parameters.Therefore, data close to 0 or 100% conversion were excluded from calculations.Given that light-off bands seemed to be narrow, particularly for CO, two sets of experiments needed to be used.Due to intended repeatability of the experiments, similar tests would result in similar results, adding no new information.It was decided to repeat the test matrix under very different external ambient conditions, so that different exhaust temperature results in the range of interest are obtained (see Figure 10).Resulting kinetic parameters are shown in Table 7: (m 3 /(mol•s) 2.30 × 10 3   (J/mol) 1.98 × 10 4

Chemical Results
Kinetic model results are shown in Table 8.Measured and modelled conversions of CO and propylene are included.

Temperature and Heat Results
To assess the complete model, experimental tests were conducted in five operating points from the test matrix.Modelled outlet temperature , and heat loss through the exhaust pipe are compared with experimental results in the five calibrating points (see Table 9).Resulting kinetic parameters are shown in Table 7:

Chemical Results
Kinetic model results are shown in Table 8.Measured and modelled conversions of CO and propylene are included.

Temperature and Heat Results
To assess the complete model, experimental tests were conducted in five operating points from the test matrix.Modelled outlet temperature T g, out and heat loss through the exhaust pipe .Q are compared with experimental results in the five calibrating points (see Table 9).The remaining four operating points of the test matrix are used to evaluate the model out of training points.Results are shown in Table 10.As commented in Section 2.2.3, it is suggested that energy recovery processes should take place downstream of after-treatment devices in order not to interfere with their operation.However, they should also occur sufficiently close in order to minimize further thermal energy dissipation to the surroundings.Therefore, the analysis of the flow abandoning the DOC is essential to know the working conditions of energy-harvesting devices.
Low uniformity in inlet flow occurring in most catalysts leads to different flow paths with different residence times.Flow inside the DOC is subject to cooling because external convection and heating because of chemical reactions in the DOC.These different paths (see Figure 11) caused by the inlet cone cause a gradient of temperature in the outlet section of the catalytic converter, as can be seen in Figure 12.

Flow Temperature Distribution at the DOC Outlet
As commented in Section 2.2.3, it is suggested that energy recovery processes should take place downstream of after-treatment devices in order not to interfere with their operation.However, they should also occur sufficiently close in order to minimize further thermal energy dissipation to the surroundings.Therefore, the analysis of the flow abandoning the DOC is essential to know the working conditions of energy-harvesting devices.
Low uniformity in inlet flow occurring in most catalysts leads to different flow paths with different residence times.Flow inside the DOC is subject to cooling because external convection and heating because of chemical reactions in the DOC.These different paths (see Figure 11) caused by the inlet cone cause a gradient of temperature in the outlet section of the catalytic converter, as can be seen in Figure 12.Coefficient of variation (standard deviation divided by mean value) of temperature distribution across the DOC outlet area (distribution shown in Figure 12b) was calculated from results of the CFD model.This statistical coefficient accounts for the dispersion from the average temperature flow.As can be seen in Figure 13, variation in temperature is enhanced in engine conditions with low mass flows, since flow thermal inertia is also lower.Engine modes in which catalytic reactions are active (see black dots in Figure 13) present more uneven distribution than those in which they are inactive (see white dots in Figure 13).Coefficient of variation (standard deviation divided by mean value) of temperature distribution across the DOC outlet area (distribution shown in Figure 12b) was calculated from results of the CFD model.This statistical coefficient accounts for the dispersion from the average temperature flow.As can be seen in Figure 13, variation in temperature is enhanced in engine conditions with low mass flows, since flow thermal inertia is also lower.Engine modes in which catalytic reactions are active (see black dots in Figure 13) present more uneven distribution than those in which they are inactive (see white dots in Figure 13).

Temperature Loss along the Exhaust Pipe
Sometimes, because of limitations in available space or ease of installation, the selected place for energy recovery devices moves away from the location with the maximum temperature available (immediately before the after-treatment systems).Average temperature of the exhaust gas in the first 50 cm downstream of the DOC was obtained from the validated model to quantify the loss in temperature when moving the energy recovery device away from the DOC.
Test-bench results with natural convection were compared with external forced convection (as in a moving vehicle with velocity ).Forced convection was simulated as boundary condition in pipe walls.Two different modes (lowest and highest engine power within the test matrix) were  Coefficient of variation (standard deviation divided by mean value) of temperature distribution across the DOC outlet area (distribution shown in Figure 12b) was calculated from results of the CFD model.This statistical coefficient accounts for the dispersion from the average temperature flow.As can be seen in Figure 13, variation in temperature is enhanced in engine conditions with low mass flows, since flow thermal inertia is also lower.Engine modes in which catalytic reactions are active (see black dots in Figure 13) present more uneven distribution than those in which they are inactive (see white dots in Figure 13).

Temperature Loss along the Exhaust Pipe
Sometimes, because of limitations in available space or ease of installation, the selected place for energy recovery devices moves away from the location with the maximum temperature available (immediately before the after-treatment systems).Average temperature of the exhaust gas in the first 50 cm downstream of the DOC was obtained from the validated model to quantify the loss in temperature when moving the energy recovery device away from the DOC.
Test-bench results with natural convection were compared with external forced convection (as in a moving vehicle with velocity ).Forced convection was simulated as boundary condition in pipe walls.Two different modes (lowest and highest engine power within the test matrix) were

Temperature Loss along the Exhaust Pipe
Sometimes, because of limitations in available space or ease of installation, the selected place for energy recovery devices moves away from the location with the maximum temperature available (immediately before the after-treatment systems).Average temperature of the exhaust gas in the first 50 cm downstream of the DOC was obtained from the validated model to quantify the loss in temperature when moving the energy recovery device away from the DOC.
Test-bench results with natural convection were compared with external forced convection (as in a moving vehicle with velocity v ∞ ).Forced convection was simulated as boundary condition in pipe walls.Two different modes (lowest and highest engine power within the test matrix) were studied for each external condition.Forced convection coefficients were calculated with external parallel flow correlations [38] taking into account velocities of the vehicle at those engine conditions.Notice that natural convection conditions can also be relevant in moving vehicles since the part of the exhaust pipe adjacent to the DOC could not be in contact with car underbody air.
Average temperature decreases (∆T/L) for external natural convection (test-bench) range from 0.9 • C/dm in the lower power mode to 1.2 • C/dm in the higher power mode.For forced convection, decrease in temperature ranges from 1.2 • C/dm to 2.6 • C/dm (see Table 11).

Accuracy of the Model
As can be seen, maximum error in chemical conversion was 11.58% and 11.81% for the CO model and the C 3 H 6 model, respectively.Notice also that at relative low temperature modes with very low conversion rates, the model shows a 0% conversion.This is due to exclusion of near-to-zero points while developing the kinetic model, in order to better predict the whole operating range.Similarly, the CO kinetic model tends to show near total conversion for all high conversion modes, since these points were excluded when developing it for the same reason as above.
As can be seen, the maximum error at calibrating points was 4% in outlet temperature and 4.4% in heat losses.Out of training conditions, the maximum error in outlet was 5.9% in outlet temperature and 5.4% in heat losses.
It is difficult to find published kinetic parameters that fit the performance of a particular DOC, due mainly to the different washcoat materials and loadings of the monoliths and also to ageing effects of the DOC.A method to develop tailored kinetic parameters in the area of interest using numerical simulations was presented.The results show that although error in chemical conversion predictions reaches almost 12%, the error in the overall performance of the thermal model is lower.This is due to the low concentration of reactant species in the exhaust gas (see Table 4): small errors in conversion do not involve big changes in total generated enthalpy of chemical reactions in the DOC.It can be concluded that the model shows good agreement with empirical data.
Regarding pressure drop, the linear correlation fits the data appropriately but for high velocities, the deviation can be larger (as seen in Figure 9).
Complex numerical schemes and boundary layer refinement were proven not necessary to obtain reliable results.Standard wall-functions showed good results in comparison with other, more demanding approaches.

Energy Recovery Considerations for Exhaust Systems
DOC outlet flow analysis showed that modules in thermoelectric generator devices placed after catalytic converters may have different performances not only because of the internal design of the devices but also because of the nature of the exhaust flow.It was found that the smaller the flow, the more uneven the flow temperature distribution is.
This can affect recovery devices being placed in the exhaust system.For instance, modules in a thermoelectric generator expected to have the same behaviour (for instance, upper and lower part thermoelectric modules in a flat-shaped, symmetrical design) may show a different performance due to exhaust flow structures.
Awareness of this effect be used to allow for electrical connections between modules depending on differences in voltage output from modules caused by variations in temperature.Even various types of modules can be selected in each zone to optimize the power output result for the different ranges of temperature.
Loss in gas temperature with distance from DOC outlet was quantified.This will help to assess the trade-off between placing the recovery device near after-treatment systems (higher temperatures) or not (easiness in installation).Notice that temperature falls fast in a relatively short distance, reducing the amount of achievable harvested power.

Conclusions
Experimental exhaust data obtained for most common light-duty diesel engines was used to construct a 3D model.A complete methodology to build 3D CFD exhaust models is presented, including how to estimate main parameters needed and numerical approaches.A 2D axi-symmetrical complementary model was shown to be useful to overcome the necessity of chemical parameters fitted to a specific model of catalytic converter.
A computational model incorporating these findings will lead to accurate and faster numerical calculations of analysis of recoverable energy in exhaust systems.This will allow more designs to be explored in a given time.
It has been pointed out that the particular nature of the exhaust flow leaving catalytic converters has to be taken into account for energy-harvesting calculations, since the end of after-treatment systems is the indicated position to place energy recovery devices.
Furthermore, temperature losses caused by placing the recovery device distant from the DOC were evaluated.This information can be used in the evaluation process of the position of a recovery device in an automotive exhaust system.Since temperature falls promptly, it is advised to place the device as close as possible.
In works regarding simulations of recovery devices, if not enough experimental information is available, it is encouraged to include upstream elements of the exhaust that could alter the flow significantly in order to better predict their performance when installed in vehicles.

Figure 1 .
Figure 1.View of the engine test bench including (a) the exhaust system pipe and (b) the diesel oxidation catalyst (DOC).

Figure 2 .
Figure 2. Sketch of the exhaust with measurement points involved in the development of the model.The exhaust system is 4 m long and pipe has a diameter of 5 cm (figure not to scale).

Figure 1 .
Figure 1.View of the engine test bench including (a) the exhaust system pipe and (b) the diesel oxidation catalyst (DOC).

Figure 1 .
Figure 1.View of the engine test bench including (a) the exhaust system pipe and (b) the diesel oxidation catalyst (DOC).

Figure 2 .
Figure 2. Sketch of the exhaust with measurement points involved in the development of the model.The exhaust system is 4 m long and pipe has a diameter of 5 cm (figure not to scale).

Figure 2 .
Figure 2. Sketch of the exhaust with measurement points involved in the development of the model.The exhaust system is 4 m long and pipe has a diameter of 5 cm (figure not to scale).

Figure 3 .
Figure 3. Test matrix.Circles represent selected engine operating points.
shows, as an example, views of the external surface of the pipe and the DOC.

Figure 4 .
Figure 4. Example of infrared images of the exhaust line in (a) the exhaust pipe and (b) the DOC.Temperatures in the scale are in °C.Squares seen in (a) are wall temperature measurement points.

Figure 4 .
Figure 4. Example of infrared images of the exhaust line in (a) the exhaust pipe and (b) the DOC.Temperatures in the scale are in • C. Squares seen in (a) are wall temperature measurement points.

Figure 5 .
Figure 5. Geometry of the 2D axisymmetric DOC model for the tuning process of kinetic parameters.Darker shade of grey represents the monolith.

Figure 6 .
Figure 6.Three-dimensional geometry model of the exhaust system.

Figure 6 .
Figure 6.Three-dimensional geometry model of the exhaust system.

Figure 6 .
Figure 6.Three-dimensional geometry model of the exhaust system.

Figure 7 .
Figure 7. Mesh detail.Tetrahedral (DOC inlet and outlet cones) and hexahedral zones (monolith and pipe) can be distinguished.

Figure 8 .
Figure 8. Results of test runs with different wall functions.No great disagreement between the three approaches was obtained.

Figure 8 .
Figure 8. Results of test runs with different wall functions.No great disagreement between the three approaches was obtained.

Figure 11 .
Figure 11.Streamlines and inlet velocity distribution in the monolith.An example of usual bad flow uniformity caused by inlet cones in automotive catalysts can be seen.

Figure 11 .
Figure 11.Streamlines and inlet velocity distribution in the monolith.An example of usual bad flow uniformity caused by inlet cones in automotive catalysts can be seen.

Figure 12 .
Figure 12.Temperature distribution for the 2400 rpm-110 Nm mode in (a) monolith, outlet cone and exhaust pipe and (b) cross-section of DOC outlet (before exhaust pipe).

Figure 13 .
Figure 13.Plot for coefficient of variation in temperature distribution at DOC outlet.White dots represent engine modes in which chemical processes were not active.

Figure 12 .
Figure 12.Temperature distribution for the 2400 rpm-110 Nm mode in (a) monolith, outlet cone and exhaust pipe and (b) cross-section of DOC outlet (before exhaust pipe).

Figure 12 .
Figure 12.Temperature distribution for the 2400 rpm-110 Nm mode in (a) monolith, outlet cone and exhaust pipe and (b) cross-section of DOC outlet (before exhaust pipe).

Figure 13 .
Figure 13.Plot for coefficient of variation in temperature distribution at DOC outlet.White dots represent engine modes in which chemical processes were not active.

Figure 13 .
Figure 13.Plot for coefficient of variation in temperature distribution at DOC outlet.White dots represent engine modes in which chemical processes were not active.

Table 1 .
Summary of the spatial discretization.

Table 2 .
Boundary conditions included in the model.

Table 3 .
Measured temperatures, mass flow values and outlet pressure.These values have been employed as boundary conditions for the model.

Table 4 .
Modeled species and volumetric fraction range in gas inlet composition within the test matrix.

Table 5 .
Boundary layer sensitivity analysis results.

Table 5 .
Boundary layer sensitivity analysis results.

Table 6 .
Experimentally-derived effective convection heat transfer coefficients for each operating mode in test matrix.

Table 6 .
Experimentally-derived effective convection heat transfer coefficients for each operating mode in test matrix.

Table 7 .
Obtained kinetic parameters for CO and C3H6 oxidations.

Table 8 .
Experimental versus chemical model conversion results.

Table 7 .
Obtained kinetic parameters for CO and C 3 H 6 oxidations.

Table 8 .
Experimental versus chemical model conversion results.

Table 9 .
Experimental versus model thermal results in calibration engine modes.

Table 10 .
Experimental versus model thermal results out of calibration engine modes.

Table 9 .
Experimental versus model thermal results in calibration engine modes.The remaining four operating points of the test matrix are used to evaluate the model out of training points.Results are shown in Table10.

Table 10 .
Experimental versus model thermal results out of calibration engine modes.

Table 11 .
Cooling in the exhaust gas for the part of the pipe adjacent to the DOC (first 50 cm).