A One-Dimensional Model for Particulate Deposition and Hydrocarbon Condensation in Exhaust Gas Recirculation Coolers

Abstract

Cooled exhaust gas recirculation (EGR) is widely used in diesel engines to control engine out NO_x (oxides of nitrogen) emissions. A portion of the exhaust gases is recirculated into the intake manifold of the engine after cooling it through a heat exchanger. EGR cooler heat exchangers, however, tend to lose efficiency and have increased pressure drop as deposit forms on the heat exchanger surface. This adversely affects the combustion process, engine durability, and emissions. In this study, a 1-D model was developed to simulate soot deposition, soot removal, and condensation of several hydrocarbon (HC) species in a circular tube with turbulent gas flow at constant wall temperature. The circular tube, which makes up the computational domain in the model, represents a single channel from any EGR cooler geometry. The model takes into account soot particle deposition due to thermophoresis, diffusion, turbulent impaction, and gravitational drift. However, thermophoresis was found to be the most dominant deposition mechanism for boundary conditions at which EGR coolers typically operate. Soot removal was modeled by considering a force balance between the drag and van der Waals forces. A lognormal distribution of particles was assumed at the tube inlet. The evolution of the particle distribution in the bulk flow along the tube as well as the mass distribution in the deposit layer on the tube walls is predicted by the model. Condensation of six HC species between C₁₅-C₂₄ alkanes was also modeled. Predictions made by the model are in reasonably good agreement with experimental data obtained from a laboratory reactor under the same boundary conditions. There are several assumptions and simplifications built into the model, which can be refined further to improve it.

Copyright 2012 American Association for Aerosol Research

NOMENCLATURE

A	=	Constant used to calculate the Cunningham slip correction factor
Ao	=	Overall heat transfer surface area of the tube (m2)
Acr	=	Cross-sectional area of the EGR cooler tube (m2)
AH	=	Hamaker's constant (J)
B	=	Constant used to calculate the Cunningham slip correction factor
C	=	Constant used to calculate the Cunningham slip correction factor
Cc	=	Cunningham slip correction factor
Cm	=	Constant used to calculate the thermophoretic coefficient
Co	=	Soot mass concentration at the tube inlet (kg/m3)
Cs	=	Constant used to calculate the thermophoretic coefficient
Ct	=	Constant used to calculate the thermophoretic coefficient
Cpg	=	Specific heat of the gas (J/kgK)
Do	=	Tube diameter (m)
D	=	Transient tube diameter (m)
Dp	=	Particle diffusion coefficient (m2/s)
	=	Mass diffusivity of HC species in air (m2/s)
dp avg	=	Average diameter of a particle bin (m)
f	=	Darcy friction factor
FD	=	Particle drag force (N)
FV	=	Van der Waals force (N)
Fb	=	Buoyancy force (N)
Fg	=	Weight of the particle (N)
Fl	=	Lift force (N)
g	=	Gravitational acceleration (9.81 m/s2)
h	=	Convective heat transfer coefficient (W/m2K)
Ho	=	Distance between particles (m)
jg,i	=	HC species mass condensation flux (kg/m2s)
kb	=	Boltzmann's constant (1.381 × 10−23 J/K)
kp	=	Thermal conductivity of soot particle (W/mK)
k deposit	=	Thermal conductivity of deposit layer (W/mK)
K	=	Proportionality constant
Kn	=	Knudsen number
K th	=	Thermophoretic coefficient
Kg	=	Mass transfer coefficient (m/s)
L	=	Tube length (m)
m	=	Deposit mass (kg)
mg	=	Mass flow rate through the tube (kg/s)
MWg	=	Molecular weight of the gas (28.89 g/mol)
MW HC,i	=	Molecular weight of HC species (g/mol)
No	=	Particle number concentration at tube inlet (#/m3)
Np	=	Number concentration in each particle bin in cell (x) and at time step (t) (#/m3)
Nu	=	Nusselt number
Pi ,g	=	Partial pressure of the HC species in the gas (kPa)
Pr	=	Prandtl number
Q	=	Rate of heat transfer (W)
r	=	Radial coordinate
ReD	=	Reynolds number
R th	=	Thermal resistance (m2K/W)
Rf	=	Fouling resistance (m2K/W)
S	=	Sticking probability
Scg,i	=	Schmidt number for HC gas species
t	=	Time step (s)
Tm	=	Mean gas temperature (K)
To	=	Gas temperature at tube inlet (K)
T wall	=	Temperature of the tube wall (K)
T surface	=	Temperature at the surface of the deposit (K)
u*	=	Friction velocity (m/s)
Um	=	Mean gas velocity (m/s)
V	=	Gas velocity at the inlet to the EGR cooler (m/s)
Vd	=	Velocity due to Eddy diffusion (m/s)
Vg	=	Gravitational drift velocity (m/s)
Vt	=	Velocity due to turbulent impaction (m/s)
V th	=	Thermophoretic drift velocity (m/s)
x	=	Distance along the tube or computational cell location (m)
yi ,gas	=	Mole fraction of the HC species in the gas stream
yi ,int	=	Mole fraction of the HC species at the gas/liquid interface

Greek symbols

αg	=	Thermal diffusivity (m2/s)
δ	=	Thickness of the viscous sub-layer (m)
∇T	=	Temperature gradient (K)
ϵ	=	Effectiveness
λ	=	Mean free path (m)
μg	=	Dynamic gas viscosity (kg/ms)
ν g	=	Kinematic gas viscosity (m2/s)
ρ deposit	=	Density of the deposit layer (kg/m3)
ρg	=	Gas density (kg/m3)
ρp	=	Particle density (kg/m3)
τ p	=	Particle relaxation time (s)
τ s	=	Wall shear stress (N/m2)
φ	=	Thickness of the deposit layer (m)
ψ	=	Strength of deposit

Abbreviations

EGR	=	Exhaust Gas Recirculation
LMTD	=	Logarithmic mean temperature difference (K)
NOx	=	Oxides of nitrogen
PM	=	Particulate matter
SOF	=	Soluble organic fraction

1. INTRODUCTION

The widely used methodology for reducing NO_x emissions from diesel engines is to recirculate a portion of the exhaust gases into the combustion chamber after cooling them in a heat exchanger known as an EGR cooler. Diesel exhaust also contains particulate matter, which consists of a solid fraction, a soluble organic fraction (SOF) and sulfates (Witze 2001). The solid fraction is mostly carbon, which consists of molecules with C/H ratios ranging from 4 to 11 (Mehta et al. 2001). The solid fraction is commonly referred to as soot and it makes up the core of the particulate. Morphologically, it consists of nearly spherical primary particles that cluster to form chain like aggregates (Mehta et al. 2001; Witze 2001). Deposition of particulate matter on the walls of the EGR cooler leads to formation of a fouling layer (Hoard et al. 2008). The fouling of EGR coolers results in significant deterioration of the cooler effectiveness and increased pressure drop across the cooler. Since the pressure differential between the exhaust and intake manifolds is the driving mechanism for EGR flow, an increase in pressure drop across the EGR cooler affects control of the desired EGR rate and results in decreased fuel efficiency due to increased pumping work (Hoard et al. 2008). Condensation of hydrocarbons (HC) in diesel exhaust occurs in the EGR cooler as well under certain boundary conditions and is a function of the vapor pressure and concentration of the HC species.

Main transport mechanisms of soot particles in EGR coolers that lead to deposition are thermophoresis, eddy diffusion, turbulent impaction and gravitational drift (Abarham et al. 2010a). It has been extensively reported in the literature that thermophoresis is the dominant mechanism for deposition of submicron particles under non-isothermal conditions (Epstein 1997; He and Ahmadi 1998; Hoard et al. 2008; Abarham et al. 2010a). Thermophoresis is particle motion in the presence of a temperature gradient. Collisions with more energetic gas molecules on the hot gas side versus the cold gas molecules near the tube wall leads to net motion of particles toward the wall.

He and Ahmadi (1998) studied the effect of thermophoretic force on particle deposition and dispersion in laminar and turbulent flows. They also compared the predictions of the Brock–Talbot and modified Cha–McCoy–Wood (MCMW) equations for variation of the thermophoretic velocity as a function of particle diameter. They found that the predictions of the Brock–Talbot and MCMW model equations agree with each other for particle diameters between 50 and 500 nm. For particles below 50 nm, the MCMW model predicts that the thermophoretic velocity decreases as the particle diameter decreases, which is the correct trend, while the Brock–Talbot equation, predicts that the thermophoretic velocity continues to increase. The Brock–Talbot model, predicts slightly higher thermophoretic velocities than the MCMW model for particle diameters above 500 nm.

Housiadas and Drossinos (2005) developed 1-D and 2-D models to investigate thermophoretic deposition in laminar and turbulent pipe flows. They found that in turbulent flow 1-D and 2-D model predictions provided the same degree of accuracy, unless there were large differences in the bulk gas to wall temperatures, in which case the 2-D approach offered significant improvement.

Abarham et al. (2009a, 2009b) developed a 1-D model for thermophoretic soot deposition and HC condensation in turbulent pipe flow. They also did a comprehensive review of the various soot deposition and removal mechanisms in EGR coolers (Abarham et al. 2010a). Predictions made by their model compared well to experimental data obtained from a controlled EGR cooler fouling test performed with sample tubes in the exhaust of a diesel engine (Abarham et al. 2009a; 2009b). Abarham et al. (2010b) also developed an analytical solution for thermophoretic deposition of submicron particles in transition and turbulent flow regimes. They did a parametric study to investigate the effect of boundary conditions on the total particulate mass deposited. It was shown that higher gas inlet temperature and inlet particulate concentration result in greater particulate mass deposited, while higher gas inlet pressure reduces it, when all other boundary conditions are kept constant.

Kern and Seaton (1959) were the first to offer a formulation of the removal term for heat exchanger deposits. They defined an average removal rate, which is assumed to be caused by shearing action of the stream at the surface. Rather than being removed particle by particle, deposit is considered to be sheared off in chunks (Loo and Bridgwater 1985). They postulated that the average removal rate over the whole tube would be roughly proportional to the shear stress and thickness of the deposit. Yung et al. (1989a) did visualization experiments and concluded that fluid drag force was the dominant re-entrainment force. They showed that if the drag force is larger than the adhesion force for small particles, removal occurs. Müller-Steinhagen et al. (1988) and Grillot and Icart (1997) have shown that when the gas velocity is increased particulate fouling is generally reduced. It was found that particulate fouling can be avoided if the gas velocity is above a critical flow velocity (Abd-Elhady et al. 2004, 2010). The critical flow velocity is defined as the bulk velocity above which rolling of the deposited particles occurs.

Our approach in this study was to develop a 1-D model for soot particle deposition, particle removal and condensation of six HC species between C₁₅-C₂₄ alkanes in a circular pipe at constant wall temperature in the turbulent flow regime. There are several differences between the present work compared to previous models, such as those developed by Housiadas and Drossinos (2005) and Abarham et al. (2009a, 2009b). The present model uses a lognormal distribution of soot particles as opposed to single diameter particles. Evolution of the particle distribution along the channel was calculated in the model by grouping the particles into different bins, where each particle bin was represented by an average particle diameter. Particle deposition by thermophoresis, eddy diffusion, turbulent impaction and gravitational drift was calculated as a function of the average particle diameter for each particle bin as opposed to using the bulk particle mass concentration. The model can thus be used to investigate the influence of different particle distributions on mass deposition in the tube or to compare the relative importance of various deposition mechanisms as a function of particle diameter for different boundary conditions. The model can also be used to predict particle deposition in turbulent pipe flows encountered in other applications, where the range of particle diameters might be different and, thus, the relative contribution of different deposition mechanisms.

Another significant difference is that the present model includes a deposit removal mechanism unlike previous models. Lance et al. (2010) have found evidence of erosion and longitudinal grooves in the deposit layer parallel to the direction of exhaust gas flow in field-aged EGR cooler samples obtained from several diesel engine manufacturers. Removal of particles was thus modeled by considering a force balance between the fluid drag and van der Waals forces. Exclusion of a particle removal mechanism will influence the growth of the fouling layer and might lead to an overestimation of fouling for a given exposure time and boundary conditions. In addition, a sticking probability was defined to estimate the fraction of particles approaching the cooler wall that actually adhere to the surface. Since particle deposition and removal equations are solved for each particle bin (or average particle diameter), the evolution of the inlet particle distribution in the bulk flow along the tube as well as the mass distribution in the deposit layer on the tube wall can be tracked in our model.

Predictions made by the model were compared to experimental data obtained from a laboratory reactor under the same boundary conditions. The experimental setup has been designed and built for controlled fouling of EGR coolers. Sub-micron soot particles are produced by a soot generator, which are then mixed with air and HCs. This mixture is subsequently heated in a tubular furnace, to simulate the conditions of exhaust gasses from diesel engines. The gas-particle mixture is cooled in an EGR cooler, and fouling of the cooler is monitored.

2. MODEL DEVELOPMENT

The model simulates particle deposition, removal and HC condensation for turbulent gas flow through a circular pipe at constant wall temperature. The computational domain of the model consists of a circular pipe of diameter (D_o ) and length (L). The pipe was divided into several cells as shown in Figure 1. The number of cells was chosen so that the model predictions were grid independent. Figure 1 also shows a single computational cell. Fluid flow and heat transfer were solved for each cell (x) at a given time step (t). It was assumed that at each time step steady-state was reached quickly. Hence, a steady-state solution was obtained at each time step and transient effects in the equations were neglected. The time-steps in the model were fixed at 1 s. Particle deposition and removal equations were solved for each computational cell and at each time step. HC condensation was solved for each cell and was coupled to the particle deposition. Particle deposition and HC condensation leads to a reduction in the tube diameter. The tube diameter calculated at a particular time step [D(x,t)] was used as input for the next time step (t + 1). The deposit surface temperature, i.e., the temperature at the interface of gas and deposit was determined from an energy balance.

FIG. 1 Computational domain of the model and a single computational cell.

The gas temperature at the tube inlet [T_o (x = 0, t)], soot mass concentration [C_o (x = 0, t)] and mole fractions of the six HC species used in the model [y _HC(i)o (x = 0, t)] were either constant or could vary as a function of time, i.e., the boundary conditions could vary at each time step. This allowed for input of experimental data directly into the model; for example, boundary conditions from an EGR cooler fouling schedule. Due to its high specific heat, engine coolant circulated through EGR coolers typically encounters less than 2°C increase in its temperature. Hence, the tube wall was assumed to be at a constant temperature [T _wall(t)] at any given time step. Other inputs to the model were the average gas pressure in the tube [P _avg(t)], the total number concentration of soot particles at the tube inlet [N_o (x = 0, t)] and the mass flow rate (m_g ) through the tube.

Some of the assumptions made in the model were:

•	Variations in the radial direction were neglected.
•	Fully developed turbulent flow was assumed in the tube and entrance effects were neglected.
•	Properties of air were used for bulk EGR gas properties since it made up more than 95% of the flowing gas mixture in experimental measurements.
•	It was assumed that the properties of the deposit formed by soot deposition and HC condensation, namely thermal conductivity and density, were constant and did not change as a function of temperature, deposit thickness or time, i.e., the deposit properties did not change at each time step.

2.1. Governing Equations for Bulk Fluid Flow and Heat Transfer

Properties of air including density, viscosity, thermal conductivity, and specific heat were used for EGR gas and were calculated at each time step (t) using ideal gas assumption at an average temperature of the gas inlet and wall temperature as:

It was assumed that the properties did not vary as a function of the axial distance along the tube. The governing equation for the mean gas temperature [T_m (x,t)] as a function of axial distance (x) at time step (t) was obtained from energy conservation on a single computational cell (Housiadas and Drossinos 2005; Incropera et al. 2007; Abarham et al. 2009a).

Solving for T_m (x,t) and replacing T _wall(t) by T _surface(x,t), where T _surface(x,t) is the temperature of the surface of the deposit layer at the interface between the gas and the deposit:

The mean velocity of the bulk gas U_m (x,t) was calculated from the mass flow rate and the cross-sectional area of the tube based on the transient tube diameter D(x,t). Reynolds number (Re _D) is given by:

Nusselt number (Nu) is calculated from the Gnielinski correlation, which is valid over a large Reynolds number range, 3000 ≤ Re _D ≤ 5 × 10⁶ (Incropera et al. 2007). Reynolds numbers in cooler passages are around 500 for low EGR rates to above 3500 for high EGR rates. The flow in EGR coolers can thus be laminar or turbulent (Hoard et al. 2008). However, only fully developed turbulent flow is considered in this model. The Reynolds numbers for the test cases described later in the paper vary from 7800–29000.

The Moody (or darcy) friction factor (f) is calculated using the correlation developed by Petukhov (Incropera et al. 2007):

The wall shear stress (τ _s ) and the friction velocity (u*) are defined as (Incropera et al. 2007):

Different expressions have been suggested for the thickness of the viscous sublayer in the turbulent boundary layer by Housiadas and Drossinos (2005) and Byers and Calvert (1969). In our model the thickness of the viscous sublayer (δ) is defined as (Yung et al. 1989b; White 2003; Abarham et al. 2009a):

FIG. 2 Particle distribution at pipe inlet. (a) Fractional distribution; (b) Cumulative distribution.

2.2. Particle Deposition

It was assumed that the soot particles at the inlet of the tube at any time step are lognormally distributed with a mean diameter of 130 nm and a standard deviation of 1.4 as shown in Figure 2a and b. Lognormal distribution is a commonly reported distribution for particles in diesel exhaust (Mehta et al. 2001; Harris and Maricq 2002). The values of mean diameter and standard deviation were chosen to match the particle distribution in the experimental test cases described in detail later in the paper. It was assumed that the particles were spherical with a constant density (ρ_p ) of 1.8 g/cc (Vishwanathan and Reitz 2008; Lance et al. 2009). The total number of particles between 1–1000 nm at the channel inlet was divided into 50 bins to reduce computational requirements. If sufficient computational power is available the particle deposition and removal equations could be solved for each particle diameter. Since the particles found in EGR coolers are typically < 200 nm, a higher resolution of particle bins could be used up to 200 nm, while a coarser resolution could be used from 200 to 1000 nm. As a starting point the width of each bin in the model was set to 20 nm and each bin was represented by an average diameter (dp _avg). The fractional distribution of the particle bins as shown in Figure 3 was determined from the cumulative particle distribution at the channel inlet shown in Figure 2b and was calculated as:

FIG. 3 Particle bin distribution at pipe inlet.

Number concentration in each bin N_p (dp _avg) at the channel inlet was calculated from the fractional distribution of particle bins and the total number of particles [(N _o(t)] at the inlet of the tube. It was assumed that the change in the number concentration within each particle bin along the tube occurred only due to particle deposition on the tube wall or removal of particles from the deposit layer, i.e., there was no coagulation of particles or interaction between the different particle bins. The number concentration in each particle bin for the next computational cell (x + 1) was calculated based on the current cell (x) at any given time step (t).

In every computational cell and at every time step particle deposition and removal equations were solved for each of the 50 particle bins using the average diameter (dp _avg). It was assumed that all the particles in each bin behaved identically and their motion was governed by the average diameter (dp _avg). The model takes into account soot particle deposition in turbulent pipe flow at constant wall temperature by the following mechanisms:

•	Thermophoresis
•	Diffusion
•	Turbulent impaction
•	Gravitational drift

Thermophoresis is particle motion due to thermal gradients. Collisions with more energetic gas molecules on the hot gas side versus the cold gas molecules near the tube wall leads to net motion of particles toward the wall. Two commonly used correlations in the literature for the drift velocity due to thermophoresis are the Brock–Talbot and modified MCMW correlations (Housiadas and Drossinos 2005; Hoard et al. 2008; Abarham et al. 2010a). In our model we have used the Brock–Talbot correlation given by:

The thermophoretic drift velocity (V _th) was calculated for each particle bin, represented by an average diameter, for every computational cell (x) at time step (t). The thermophoretic coefficient (K _th) for each particle bin was calculated as:

where the Cunningham slip correction factor (C _c) was given by:

where the Knudsen number (Kn) for each particle bin was calculated from the mean free path of gas molecules (λ).

It was assumed that the thermal conductivity of soot particles (k_p ) was 0.5 W/mK and the thermophoretic constants A, B, C, C_s, C_m , and C_t were assumed to be 1.2, 0.4, 1.1, 1.14, 1.17, and 2.18 respectively (Lance et al. 2009; Abarham et al. 2010a).

The temperature gradient for a computational cell (x) at time step (t) was given by (Housiadas and Drossinos 2005):

Particle motion can also occur due to concentration gradients. The particle deposition velocity due to diffusion (V_d ) for each particle bin in every computational cell (x) at time step (t) was defined based on Fick's law as:

Particle deposition could also occur inertially if the particle was unable to follow the flow depending on the particle relaxation time. Smaller diameter particles have short relaxation times and are able to follow the flow while larger particles are unable to do so. The drift velocity due to inertial impaction in turbulent flow (V_t ) for each particle bin in every computational cell (x) at time step (t) is defined as (Abarham et al. 2010a):

where the particle relaxation time (τ _p ) was calculated from the Cunningham slip correction factor (C_c ).

Finally the gravitational drift velocity (V_g ) for each particle bin in every computational cell (x) at time step (t) is defined as (Abarham et al. 2010a):

It was found from the above analysis and has been reported by other researchers, Abarham et al. (2010a) and Hoard et al. (2008), that thermophoresis was the dominant deposition mechanism for the range of particles typically found in EGR coolers. The drift velocities due to other deposition mechanisms were retained in the model for completeness and to keep the model versatile enough to predict particle deposition in turbulent pipe flows encountered in other applications, where the range of particle diameters might be different.

2.3. Particle Removal

There are several forces acting on a spherical particle attached to the surface as shown in Figure 4 (Abd-Elhady et al. 2010). Removal of deposit mass was calculated in the model by considering a simple force balance between only the drag and van der Waals forces acting on each deposited particle bin. The lift force and particle weight were found to be not comparable to the drag and van der Waals forces. It was assumed that there was no interaction between the deposited particle bins in the fouling layer. The drag force (F_D ) due to fluid flow acting on a particle bin in a cell (x) at time step (t) is given by (Yung et al. 1989b):

The van der Waals force (F_V ) acting on the particle bin in a cell (x) at time step (t) is given by:

FIG. 4 Different forces acting on a particle resting on a flat surface and exposed to a shear flow. Point O is the fulcrum of rolling. (Abd-Elhady et al. 2010).

Where A _H is the Hamaker's constant, which was assumed to be 10⁻²⁰ J (Abarham et al. 2010a) and H_o is the separation between particles, which was assumed to be one-fiftieth of the particle diameter (dp _avg).

Not all particles that migrate to the tube walls adhere to the surface. A sticking probability (S) was defined to determine the fraction of particles approaching the wall that actually adhere to the surface (Epstein 1997). The sticking probability for a particle bin in a cell (x) at time step (t) was calculated as: If

Else

In reality, the particle agglomerates attached to the surface may not be spherical and can have a range of non-spherical geometries. The assumption of sphericity for the particles will, therefore, limit the accuracy of the forces acting on the particle as well as the sticking probability calculated in the present model. This can be improved further by accounting for the morphology of diesel particulate matter.

2.4. Particle Conservation

The mass deposited by a single particle bin in a computational cell (x) at time step (t) is given by:

So, the total mass deposited by all particle bins in the cell (x) at time step (t) is obtained by a summation over all the particle bins:

The widely accepted removal term for heat exchanger scale or deposit is (Loo and Bridgwater 1985; Bohnet 1987; Royer et al. 2010):

where τ _s is the wall shear stress, ψ is the strength of scale factor, and K is a proportionality constant. In our model, the mass removed from each cell (x) at time step (t) for a single particle bin was defined as:

where (mass_deposit) is the total accumulated mass in the deposit for a single particle bin (dp _avg) in cell (x) up to the previous time step (t − 1). The accumulated mass of each particle bin in the deposit layer is known from particle mass conservation applied to each computational cell at every time step. The strength of scale factor (ψ) was assumed to be 1. The criterion for removal of accumulated mass for each particle bin in the deposit layer in cell (x) at time step (t) was defined as:

The total deposit mass removed from cell (x) at time step (t) for all particle bins is given by:

The net soot mass accumulated in cell (x) at time step (t) is given by:

This approach of solving particle deposition and removal equations for each particle bin enabled the model to predict the evolution of the particle distribution in the gas along the channel as well as the mass distribution in the deposit layer.

2.5. HC Condensation

The EGR cooler feed gas in a diesel engine consists of several HC species. Heavier HC species are more likely to condense at the coolant temperatures (30–90°C) typically encountered in diesel engine EGR coolers. The model took into account condensation of six HC species as listed in Table 1. Sluder et al. (2008) found that the gas chromatograph/mass spectrometer (GC/MS) chromatograms for deposits in surrogate tubes fouled by exhaust from a diesel engine had a peak approximately around C₁₉–C₂₀, both of which are present in diesel fuel in significant quantities. They also found that the tube cooled by coolant at 40°C temperature retained HCs as light as C₁₅. Species like C₂₅ and beyond are present in diesel fuel at insignificant levels and are typical of the lowest end of the lubricating oil HC distribution (Sluder et al. 2008). The mass condensation flux (j_{g ,i} ) from the gas stream to the surface is defined as (Hoard et al. 2008):

TABLE 1 HC species included in the model

Molecular formula	Name
C15H32	Pentacosane
C17H36	Heptacosane
C19H40	Nonadecane
C20H42	Eicosane
C22H46	Docosane
C24H50	Tetracosane

where (y_{i ,gas}) is the mole fraction of the HC species in the gas stream determined from experimental measurements and (y_{i ,int}) is the mole fraction of the HC species at the gas/liquid interface. The mass transfer coefficient (K_g ) is defined as:

Since the heat transfer coefficient (h) is related to the Nusselt number (Nu) the mass condensation flux can be rewritten as:

where, the Schmidt number for each HC species (Sc_{g ,i} ) was calculated using the diffusivity of each HC species in air (). The mass diffusivity of each HC species in air () at standard conditions was obtained from Yaws (2010) and was corrected for temperature and pressure using the average gas temperature and pressure in the EGR cooler assuming ideal gas behavior as shown by Incropera et al. (2007).

The mole fraction of the HC species at the gas/liquid interface (y_{i ,int}) is given by:

The partial pressure of each HC species was calculated using the respective Antoine coefficients for each species similar to the approach of Abarham et al. (2009b):

Antoine coefficients for each HC species were obtained from Yaws et al. (2009). The mass condensed for a HC species in a single computational cell at time step (t) is:

The total condensed HC mass for all 6 species is:

The mole fraction for a particular HC species at the inlet of the next computational cell (x + 1) was calculated based on molar conservation applied to the current cell (x). It was assumed that all the condensed HCs diffuse through the deposit and migrate toward the cold wall as reported by Lance et al. (2010) from their optical microscopy measurements on deposits from field fouled EGR coolers. The cold wall inhibits the evaporation of the condensed HCs. Hence, evaporation of condensed HC species from the deposit was not included in the model.

2.6. Calculation of Various Parameters

The total mass deposited in cell (x) at time step (t) is the summation of the net soot mass and the condensed HC mass:

The local thickness of the deposit layer (φ) is given by:

where (ρ _deposit) is the density of the deposit layer. In our model, the density of the deposit layer was assumed to be 35 kg/m³ as reported by Lance et al. (2009). It was assumed that the density of the deposit layer was constant throughout the tube and did not change with deposit thickness or temperature. The diameter of the tube in cell (x) at the next time step (t + 1) is calculated from the current time step (t):

There is convective heat transfer from the bulk gas to the deposit layer and there is conduction heat transfer across the deposit to the tube wall. The temperature at the surface of the deposit layer, i.e., at the interface between the deposit layer and the gas was calculated from a heat transfer balance by solving the equation:

where the local heat transfer coefficient [h(x,t)] is calculated from the local Nusselt number [N_u (x,t)]. The thermal conductivity of the deposit layer (k _deposit) was assumed to be 0.05 W/mK (Lance et al. 2009). It was assumed that the thermal conductivity of the deposit layer was constant throughout the tube.

Some of the heat exchanger performance metrics were also calculated in the model. The effectiveness (ϵ) was defined as:

The thermal resistance (R _th) was defined as:

where A _o is the surface area for heat transfer, LMTD is the logarithmic mean temperature difference and Q is the heat transfer from the tube. The fouling layer thermal resistance (R_f ) is defined as:

3. MODEL EVALUATION

The model was used to calculate fouling of a single tube exposed to exhaust gases for 6 h. The boundary conditions for a test case used to evaluate the model are given in Table 2, while the concentrations of the six HC species are given in Table 3. The tube geometry and the inlet particle size distribution were chosen to match the experimental data described in the next section. The HC concentration values do not represent actual measurements, but were chosen to test the model. Figure 5 shows the deposit layer thickness along the tube at 3 different exposure times. Maximum deposit thickness occurs at the tube inlet due to the highest temperature gradient and thereby thermophoretic effect.

TABLE 2 Boundary conditions of test case for model evaluation

Parameter
Tube inner diameter	9.2 mm
Tube length	0.21 m
Soot mass concentration	Co: 10 mg/m^3
Soot number concentration	No: ∼10^12 #/m^3
Inlet gas temperature	673 K (400°C)
Coolant (wall) temperature	308 K (35°C)
Flow rate	5.2 kg/hr
Inlet gas velocity	20 m/s
Average pressure	150 kPa
Exposure time	6 h

TABLE 3 HC species concentrations

HC species	Concentration (ppm)
C15H32	10
C17H36	5
C19H40	5
C20H42	5
C22H46	5
C24H50	2

FIG. 5 Variation of deposit layer thickness along the tube length at different exposure times.

The deposit surface temperature, i.e., temperature at the interface of the deposit layer and bulk gas also decreases along the length of the tube, while it increases as a function of time as shown in Figure 6. Due to the thicker deposit layer at the tube inlet, heat transfer through the deposit is inhibited due to its low thermal conductivity leading to the highest deposit surface temperature at the tube inlet. As the deposit layer grows with time, the deposit surface temperature increases from the initial wall or coolant temperature at t = 0 due to inhibition of heat transfer across the deposit.

FIG. 6 Variation of the surface temperature of the deposit layer along the tube length at different exposure times.

Condensation of HC species was coupled to the soot deposition model. The condensed mass for the six HC species (C₁₅–C₂₄) deposited in the tube as a function of time is shown in Figure 7. As mentioned earlier, it is assumed that all the condensed HCs diffuse through the deposit and migrate toward the cold wall (Lance et al. 2010). The cold tube wall inhibits the evaporation of the condensed HCs. Hence evaporation of condensed HC species from the deposit is not included in the model. As the deposit surface temperature increases no further condensation occurs and therefore HC species, C₁₇–C₂₂, reach an asymptote between 0–1.5 h of exposure time, while C₂₄ reaches an asymptote after about 3 h. It was observed that there is very little condensation of C₁₅ at these boundary conditions. This indicates that lighter HC species in the exhaust will not condense under the operating conditions typically encountered in EGR coolers and thereby do not contribute to fouling. It was also observed that the HC condensation was more pronounced toward the tube exit due to the lower deposit surface temperatures, which is in contrast to particle deposition.

FIG. 7 Condensed mass of different HC species in the tube as a function of time.

The net soot mass deposited in the tube is shown in Figure 8 as a function of time. The net soot mass deposited in the tube is about 13% of the total exposure over 6 h. Soot deposition due to thermophoresis continues to occur as long as a temperature gradient exists between the bulk gas and the deposit layer. The total condensed HC mass, which is the summation of the condensed mass for the six HC species considered in the model, is also shown.

FIG. 8 Net soot mass and total HC condensed mass deposited in the tube as a function of time over 6 h.

The total mass deposited in the tube is the summation of total deposited soot mass and condensed HC mass. The effectiveness of the single tube heat exchanger degrades over time due to reduction in heat transfer because of fouling as shown in Figure 9.

FIG. 9 Effectiveness of the single tube heat exchanger as a function of time.

Using the particle bin approach, it was also possible to track the evolution of the particle distribution in the bulk gas through the channel, since it was assumed that change in the number concentration of each particle bin occurred only due to particle deposition or re-entrainment. Figures 10a and b show the distribution of the particle bins in the bulk gas at the tube inlet (x = 0), middle (x = L/2), and tube outlet (x = L) at t = ½ h and t = 6 h. Both figures can be used to visualize changes in the particle distribution in the bulk gas as it moves from channel inlet to outlet at a particular time step. Couple of observations can be made from the figures. First, at any given time step, as the feed gas moves along the channel there is a net decrease in the number concentration of different particle bins as particles are deposited on the channel wall and re-entrained into the gas flow due to shear removal. At both time steps (at t = ½ h and t = 6 h) the evolution of the particle distribution in the bulk flow is captured by the model as the feed gas moves from channel inlet to outlet. Second, at t = 6 h, thermophoretic particle transport, which is the dominant soot deposition mechanism, is reduced compared to t = ½ h due to the increase in the deposit surface temperature and consequent reduction in the temperature gradient between the bulk gas and the deposit. Also, due to the reduction in tube diameter, the bulk gas velocity increases leading to higher shear force on the deposit layer and therefore higher mass removal from the deposit layer.

FIG. 10 Evolution of the particle bin distribution in the gas along the tube. (a) t = ½ h; (b) t = 6 h.

Similarly, the growth of the mass distribution in the deposit layer could also be calculated since it was assumed that there was no interaction between the deposited particle bins on the channel wall. The accumulated mass of each particle bin in the deposit layer is known at every time step from particle mass conservation applied to each computational cell. Figures 11a and b show the evolution of the mass distribution in the deposit layer as a function of time at the tube inlet and outlet. The total mass deposited at the tube outlet was lower compared to tube inlet due to the decreased temperature gradient.

FIG. 11 Evolution of the deposit mass distribution as a function of time at (a) Tube Inlet; (b) Tube Outlet.

This approach can thus be used to investigate the influence of different particle distributions on EGR cooler fouling as well as boundary conditions favorable for removal of particles of different diameters.

4. COMPARISON TO EXPERIMENTAL DATA

4.1. Experimental Setup

To validate the model, theoretical results were compared with the experimental data obtained from a laboratory reactor. The experimental laboratory reactor consisted of a soot generator, gas/particle flow heater, testing section for EGR coolers and an exhaust system as shown in Figure 12. Details of the test setup and measurement procedures are fully described in Abd-Elhady et al. (2010). The soot particles are formed in the soot generator due to the incomplete combustion of ethylene (C₂H₄) and their size distribution is controlled by the stoichiometry of combustion. The mass flow rate of ethylene and air are measured and controlled using flow controllers for each gas. The average diameter of the soot particles was 130 nm with a standard deviation of 55 nm, which is similar to the particle size distributions reported in the literature for diesel engine exhaust by Wentzel et al. (2003); Park et al. (2004); Kittelson et al. (2006); and Lapuerta et al. (2007). The mass concentration of soot particles in the feed gas to the EGR cooler was 100 mg/m³. Soot particles produced by the soot generator were mixed with external air in the mixing tube as shown in Figure 12. The mass flow rate of the injected air was controlled and measured using a mass flow controller. The mixture of air and soot particles was heated up to 400°C in a tubular furnace with 8 kW maximum power, to simulate the conditions of exhaust gasses from diesel engines. The hot mixture of gas and soot particles was then passed to the EGR cooler, which in the present investigation was a shell and tube type heat exchanger consisting of 1–3 corrugated circular tubes. It should be pointed out that the tubes used in the experimental setup were corrugated to improve the rate of heat transfer, while the circular tube used in the model was assumed to have a smooth surface. An independent water circuit was used to maintain the inlet temperature of the cooling water to the EGR cooler at a constant level of 80°C, to simulate the typical coolant temperatures in diesel engines (Zheng et al. 2004). The mass flow rate of the cooling water was controlled by a flow controller and was measured with an accuracy of ± 1%. The gas velocity through the EGR cooler was varied by changing the flow rate of the injected air, via a mass flow controller. The inlet gas velocity (V) through the EGR cooler was calculated from the measured mass flow rate (m_g ) of the flowing gasses and the cross-sectional area (A _cr) of the cooler tubes as follows:

where ρ_g was taken as the density of air, which makes up more than 95% of the flowing gas mixture. The gas velocities reported in this paper are the gas velocities at the inlet to the EGR cooler. The performance of the EGR cooler was monitored through the respective inlet and exit gas/water temperatures, which were measured using calibrated thermocouples with an accuracy of ±0.1°C. The gas pressure was measured before and after the EGR cooler with an accuracy of ±0.4%.

FIG. 12 Schematic flow diagram of the experimental setup (Abd-Elhady et al. 2010).

4.2. Results and Discussion

Abd-Elhady et al. (2004, 2010) found that particulate fouling can be avoided if the gas velocity is above a critical flow velocity. The critical flow velocity is defined as the bulk velocity above which rolling of the deposited particles occurs. Abd-Elhady et al. (2010) used the laboratory reactor described above to experimentally investigate the influence of gas flow velocity on particulate fouling of EGR coolers, by operating the EGR cooler at gas velocities below and above the critical flow velocity. Their results were used to verify the model. In the experimental study three test cases were run, where everything was held constant except for the inlet gas velocity, which was varied by varying the amount of injected air. No HC species were introduced into the feed gas to the EGR cooler for these test cases. Hence, fouling of the cooler in these experiments is only due to soot. The boundary conditions for all three test cases are given in Table 4, which were used as inputs to the model.

TABLE 4 Boundary conditions for experiments

Parameter
Tube inner diameter	9.2 mm
Tube length	0.21 m
Soot mass concentration	C o : 100 mg/m^3
Soot number concentration	N o : ∼10^13 #/m^3
Inlet gas temperature	673 K (400°C)
Coolant (wall) temperature	353 K (80°C)
Flow rate	4.94, 11.5 and 19.75 kg/h
Inlet gas velocity	30, 70, and 120 m/s
Average pressure	101 kPa
Exposure time	6 h

Since the tubes were exposed to the particle laden feed gas for about 6 h in the experimental measurements, the model was also used to calculate soot deposition in the tube for the same time duration. Table 5 shows a comparison between the deposit thickness at the pipe inlet calculated by the model after 6 h of exposure and experimental measurements. The thickness of the deposit layer is highest at the tube inlet and decreases along the length of the tube. For the 120 m/s case the fouling layer thickness was very small and difficult to quantify in the experimental data.

Figures 13a–c show a comparison between the thermal resistance of the fouling layer calculated by the model and experimental data for the three inlet gas velocities. Overall there is good agreement between the model and experimental measurements. The model is able to predict the decrease in the fouling layer thermal resistance with increase in gas velocity as well as the differences in the initial rate of fouling between the three gas velocities.

FIG. 13 Comparison of effect of gas velocity on fouling layer thermal resistance between model and experiment for (a) 30 m/s; (b) 70 m/s; and (c) 120 m/s.

The shear force due to the gas flow is not sufficient to remove much of the deposited particles at 30 m/s inlet gas velocity, leading to a more or less continuous increase in the fouling layer, as shown in Figure 13a, in both experimental measurements and model predictions.

TABLE 5 Comparison of deposit layer thickness

Gas velocity (m/s)	Experiment (mm)	Model (mm)
30	1.0	0.87
70	0.42	0.33
120	–	0.13

Increasing the gas velocity further to 70 m/s, it can be observed that the fouling layer thermal resistance decreases with the flow velocity as shown in Figure 13b. The fouling process can be strongly reduced, but not completely prevented if the gas velocity is increased to 120 m/s. With increase in gas velocity the shear force on the deposit layer increases leading to increased removal of particle mass. An asymptotic value of fouling layer thermal resistance is reached for the 120 m/s case in both the model and experimental measurements as shown in Figure 13c, indicating an equilibrium between particle deposition and removal. The marginal discontinuities in the fouling layer thermal resistance versus time data from experimental measurements for gas velocity of 120 m/s are due to the successive breakage of the fouling layer and re-growth during the experimental run. It should be noted that a velocity of 120 m/s is not practical in current diesel EGR coolers.

There are some differences between the model and experiment especially for the 70 m/s gas velocity and in predicting the initial rate of fouling at 120 m/s. These differences could be due to corrugated tube geometry used in the experiment as opposed to a smooth tube used in the model. Among the assumptions and simplifications made in the model, the assumption of a constant deposit layer density and thermal conductivity along the length of the tube as well as the values used for these two properties could be an important source of the differences observed. Further work is required to understand the variation of these properties with temperature and layer thickness and implement them in the model.

5. SENSITIVITY ANALYIS

There are several assumptions and simplifications made in the model. The assumption of a constant deposit layer density and thermal conductivity along the length of the tube as well as the values used for these two properties are important parameters that affect the calculations made by the model. In reality, the deposit properties change with temperature and layer thickness. In the present model the values selected for these two properties were based on the findings of Lance et al. (2009).

5.1. Thermal Conductivity of the Deposit Layer

The test case described in Table 2 was used to examine the sensitivity of the total deposit mass in the tube, as calculated by the model, to the thermal conductivity of the deposit layer. Figure 14 shows the total mass gain as a function of time calculated by the model for different values of deposit layer thermal conductivity, while the density was held constant at 35 kg/m³. Table 6 shows the variation in the total mass gain calculated by the model, with the deposit layer conductivity of 0.05 W/mK used as the base case. The total deposit mass calculated by the model shows a 10–21% variation with a 40–100% change in the value of the deposit layer thermal conductivity used and a more or less linear relationship with the deposit thermal conductivity. As the deposit layer thermal conductivity is increased, the surface temperature of the deposit decreases, leading to a higher temperature gradient between the bulk gas and the deposit layer and, therefore, increased particle transport due to thermophoresis. The thermal conductivity of the deposit is mainly dependent on its porosity. Further experimental measurements with deposits of varying porosity are required to confirm the influence of deposit layer thermal conductivity on fouling.

TABLE 6 Effect of deposit layer thermal conductivity

Deposit thermal conductivity (W/mK)	Total mass gain (mg)	% Difference
0.03	42.2	12.2
0.05	48.1	–
0.07	52.6	9.5
0.1	58.6	21.9

FIG. 14 Effect of deposit thermal conductivity on total deposit mass calculated by the model as a function of time.

5.2. Deposit Layer Density

The same test case was used to examine the sensitivity of the total mass gain in the tube, as calculated by the model, to the value used for the deposit layer density. Figure 15 shows the total mass gain calculated by the model with variation in the deposit layer density, while the thermal conductivity was held constant at 0.05 W/mK. Table 7 shows the variation in the total mass gain calculated by the model, with the deposit layer density of 35 kg/m³ used as the base case. The total deposit mass calculated by the model shows a 12–16% variation with more than 40% change in the value of the deposit layer density used and a more or less linear relationship with the deposit density. As the deposit layer density is increased, the thickness of the deposit layer decreases leading to lower temperatures at the deposit surface due to increased heat transfer and, therefore, increased mass gain due to thermophoresis.

TABLE 7 Effect of deposit layer density

Deposit density (kg/m^3)	Total mass gain (mg)	% Difference
20	40.2	16.5
35	48.1	–
50	54.1	12.6

FIG. 15 Effect of deposit density on total deposit mass calculated by the model as a function of time.

The total mass gain in the tube calculated by the model is influenced by the thermal conductivity and density of the deposit layer. To make the quantitative predictions of the model more robust, further work is required to implement the variation of deposit layer properties. Further experimental measurements of deposit properties from fouled EGR coolers will help in developing a model for the variation in fouling layer properties.

5.3. Particle Bin Distribution

The particle distribution at the tube inlet was divided into several bins to reduce the run time of the model. Ideally the particle deposition and removal equations could be solved for each particle diameter in the distribution. The base strategy used in the model was to divide the particles between 1–1000 nm into 50 bins of 20 nm width. A different particle bin distribution was used to examine its influence on the model calculations. Since the particles found in EGR cooler feed gas are typically < 200 nm, the particles between 1–400 nm were divided into 20 bins of 20 nm width, while particles between 400–1000 nm were divided into 6 bins of 100 nm width. Each particle bin was represented by a corresponding average diameter. The test case described in Table 2 was used to examine the influence of this staggered particle bin distribution on model predictions. The difference in the total deposit mass over 6 h calculated by the model for the base and staggered particle bin distributions was less than 0.15%. This was true for other parameters calculated by the model as well.

6. CONCLUSIONS

A 1-D model was developed for soot particle deposition, soot removal, and HC condensation in turbulent pipe flow at constant wall temperature. The model took into account soot deposition due to thermophoresis, diffusion, turbulent impaction, and gravitational drift; although thermophoresis was found to be the most dominant deposition mechanism. Soot removal was modeled by considering a force balance between the drag force and van der Waals force. The model was capable of analyzing steady-state as well as transient boundary conditions from an EGR cooler fouling schedule. The model was able to predict growth of the fouling layer as well as several heat exchanger performance parameters under different boundary conditions. For the test case used to evaluate the model at boundary conditions similar to those observed in diesel engines, the model predicted a 21% drop in thermal effectiveness in just 6 h of exposure to soot laden gas flow. This degradation of cooler effectiveness is the main problem of EGR coolers. The total predicted mass gain was about 48 mg, which is about 13% of the total exposure over 6 h. The model could also predict the evolution of the inlet particle distribution in the bulk gas along the tube as well as the mass distribution in deposit layer on the tube walls, which can be used to investigate the influence of particle distribution in the bulk gas on fouling.

Predictions made by the model for deposit layer thickness and fouling layer thermal resistance were in reasonably good agreement with experimental data obtained from a laboratory reactor under the same boundary conditions. The final value of the fouling layer thermal resistance predicted by the model was within 10% of experimental measurements for the three inlet gas velocities. The model was also able to resolve differences in fouling behavior for three different inlet gas velocities.

The assumption of a constant deposit layer thermal conductivity and density as well as the values used for these properties influence the quantitative predictions of the model. The total deposit mass calculated by the model showed a 12–16% variation with a ±40% change in the value of the deposit layer density, while it varied by 10% with a 40% change in the value of the deposit layer thermal conductivity. Further investigations are required to include the variation of deposit layer properties with temperature and thickness in the model.

Acknowledgments

The authors wish to thank Mr. T. Zornek and Dr. Abd-Elhady of the University of Stuttgart for the experimental data and insights.