# Development of a Particulate Filter Model for the Prediction of Backpressure: Improved Momentum Balance and Entrance and Exit Effect Equations.

1. INTRODUCTIONParticulate Matter (PM) [1,2,3,4] emitted in the exhaust of internal combustion engines is detrimental to the human health and the environment. Studies have shown that exposure to PM may cause health problems such as premature death, asthma, lung cancer and other pulmonary and cardiovascular problems [3,4,5]. Diesel PM is also a major source of carbon black, which has a high global warming potential (albeit a short lived one) [6]. PM influences atmospheric visibility and can lead to the soiling of buildings [7]. Diesel Particulate Filters (DPF) [4,8] have been extensively and successfully used to remove PM from the exhaust of diesel engines for around 20 years [8,9,10]. More recently, with the introduction of legislated limits for particulate number, Gasoline Particulate Filters (GPF) are being introduced to control ultrafine particulate emissions from gasoline direct injection engines [4,11,12],

As this work applies equally to DPFs and GPFs, the term Particulate Filter (PF) will be used to refer to both.

Computer simulation has become an important tool for designing aftertreatment systems [13,14], enabling a wide range of parameters and scenarios to be investigated more quickly and cost effectively than a purely experimental, "trial and error" approach. Simulation can also be used gain insights into the emissions control system which are not readily obtained by measurement.

This paper concerns the development of a PF model for backpressure prediction. Most one-dimensional PF models are based on the pioneering work of Bissett [15], Since that work numerous models for PFs with various improvements have been published, e.g. [16,17, 18,19,20,21,22,23,24,25], An extensive review of the development of PF models has been written by Koltsakis et al. [26].

Probably the most significant advance in PF modelling in recent years has been made by Bissett, Kostoglou and Konstandopoulos [27,28]. Previously, PF models have assumed that the gas flow in the channels is not significantly influenced by the flow through the porous walls. These workers have investigated this, publishing correlations which enable the effect of the gas flow out of or into channels through the porous walls to be incorporated into the usual one-dimensional PF model framework.

This work describes the development and application of a one-dimensional model for a PF. The model makes two innovations. Firstly, the term for momentum convection in the gas momentum balance equations includes the loss or gain of axial momentum in the direction perpendicular to the channels; neglecting this results in the momentum convection term being too large. Secondly, equations for the pressure change due to the abrupt contraction at the PF entrance and for abrupt expansion at the exit are derived which take into account the fact that the velocity profile across the channels is not flat.

The rest of the paper is structured as follows: The experimental section (section 2) describes the collection of data for model calibration and for model input for drive cycle simulations. The development of the model is described in section 3 and model calibration in section 4. The following section (section 5) describes application of the model to investigate the relative magnitude of the different contributions to backpressure, viz. across-wall losses (Darcy and Forchheimer), along-channel losses (viscous and inertial) and entrance and exit effects, and how these vary with flow rate, PF length, temperature, etc. Finally, (section 6) the effect of the model modifications proposed by Bissett, Kostoglou and Konstandopoulos [27,28] to account for the impact of wall flow is investigated.

2. EXPERIMENTAL

2.1 Cold Flow Experiments

Cold flow measurements were made on uncoated, 4.5 and 6.0 inch (114.3 and 152.4 mm) long, 4.66 in (118.4 mm) diameter cordierite PFs supplied by the same manufacturer. The cell density was 300 cpsi (465 000 cells per [m.sup.2]) and the wall thickness 8 thou (0.20 mm). The parts were of the same type/material, i.e. with the same porosity, mean pore diameter, etc.

Backpressure measurements were made with a SuperFlow[R] SF1020 ProBench flow bench [29]. The key design points of this equipment for this work are shown in the schematic in Fig. 1 (top). During testing, air is drawn into the top of the part from the room. On leaving the part, this air expands into a large plenum, which has a cross-sectional area approximately 60 times that of the PFs tested in this work. Since, the inlet of the PF is at atmospheric pressure, it is clear that the pressure within the PF is sub-atmospheric. This contrasts with a PF on a vehicle (Fig. 1, bottom) where the exit is at (or close to) atmospheric pressure and the pressure within the PF is above atmospheric pressure.

2.2 Hot Flow Data

As input for transient simulations data for a 1.6 L gasoline passenger car measured on a chassis dynamometer over a class 3 Worldwide harmonised Light vehicles Test Cycle (WLTC) [30,31] was used. The test was run according to standard procedures. This data file was used purely for input to the model (gas temperature and mass flow); no measurements of PF backpressure or outlet temperature were available for comparison with simulation results. A temperature measured 25 mm in front of the catalyst position was used as input to the PF model.

3. DEVELOPMENT OF MODEL

This section describes the derivation of a one-dimensional model for a PF. The model is based on the one by Bissett [15], but with various improvements, such as allowing for the constriction of the inlet channel by the soot cake and compressible flow. The model makes two innovations: Firstly, the term for momentum convection in the gas momentum balance equations includes the loss (inlet channel) or gain (outlet channel) of axial momentum in the direction perpendicular to the channels; neglecting this results in the momentum convection term being too large. Secondly, equations for the pressure change due to the abrupt contraction at the PF entrance and for abrupt expansion at the exit are derived which take into account that fact that the velocity profile across the channels is not flat.

While this paper only covers simulations of symmetric filters, where the inlet and outlet channels have the same dimensions, for completeness the model will be derived (with the exception of the across-wall pressure drop) for the case where the inlet and outlet channels have different dimensions. Similarly, while this paper does not include soot-loaded PFs, the model derivation will include the presence of soot in the channels.

The model makes the following assumptions:

* Uniform flow distribution at the entrance to the PF.

* Negligible radial temperature gradient across the catalyst.

* Reactant gases are dilute so that changes in the amount of substance (and hence gas volume) with reaction can be ignored, i.e. this is a trace system.

* The channels are square in cross section.

* The flow profile in the channels is unaffected by the flow to/from the porous walls of the channels.

* The flow profile in unaffected by the presence of PM due to the low concentration of PM in the exhaust gas.

* The gas velocity, pressure and temperature profiles along the channels are at steady-state as the gas residence time is short compared to the transient behaviour of the PF.

* Negligible pressure gradients across the width of the channels.

* The pressure drop across the soot cake and PF wall can be described by Darcy's law with the Forchheimer extension; any effect of slip flow is neglected.

* Adiabatic behaviour, i.e. no heat losses to the surroundings.

* Axial transport of heat in the solid phase by conduction.

* Heat is transported along the channels by convection; axial conduction through the gas in the channels is negligible in comparison.

* Conduction of heat through the gas in the channels to the channel walls in the direction perpendicular to the channels is significant and can be described by heat transfer terms.

* Negligible heat transport by radiation since those portions of the long thin channels that have significant view factors with each other are expected to have similar temperatures.

* Uniform temperature across the PF wall and soot cake in the direction perpendicular to the channels at each axial location, i.e. the PF wall and soot cake temperature varies axially, but not in the through-wall direction.

* The gas passing through the soot cake and PF wall has the same temperature as the soot cake and PF wall due to efficient heat transfer.

* Soot is only located only as a cake; soot in the wall is not considered.

* Ideal gas behaviour within the PF, i.e. the gas within the PF is compressible.

* The gas is incompressible during contraction as the gas enters the PF and during expansion as the gas leaves the PF.

* Heat capacities and viscosity are temperature dependant.

The focus of this paper is on backpressure prediction; simulation of reactions occurring over a catalytic coating on the PF and soot accumulation and removal by oxidation will not be discussed in this work.

3.1 Dimensions / Geometry

The model is one-dimensional, describing a single inlet channel-outlet channel pair. Figures 2 and 3 show cross-sectional and longitudinal schematics of the PF with the dimensions marked.

The channel width of the of the soot-loaded inlet channel, [d.sub.i], is related to the soot-free width of the inlet channel, d, and the soot loading, S (mass of soot per unit volume of PF), by:

[d.sub.i] = [square root of [d.sup.2] - 2S / [rho]cell [rho]s] (1)

where [[rho].sub.s] is the packing density of the soot cake and [[rho].sub.cell] is the cell density, i.e. the number of channels per unit cross-sectional area of the PF. For segmented PFs, [[rho].sub.cell] should be for the whole filter (including the cement between the segments), not just a segment.

The open cross-sectional area of the soot-free PF in the region away from the plugs (i.e. the cross-sectional area of the empty channels) divided by the cross-sectional area of PF, s, is given by:

[??] = [([d.sup.2] + [d.sup.2.sub.0])[[rho].sub.cell]/2] (2)

The open area fraction of the front and rear faces of the PF, [[epsilon].sub.i] and [[epsilon].sub.o], are given by:

[[??].sub.i] = [d.sup.2.sub.i][|.sub.z=0][rho]cell / 2 (3)

[[??].sub.0] = [d.sup.2.sub.0][rho]cell / 2 (4)

where [d.sup.2.sub.i|z=0] indicates the value of [d.sup.2.sub.i] at z=0, i.e. the front of the PF. Since ideal gas behaviour is assumed, gas density is given by,

[rho] = [[PM]/[RT]] (5)

where M is the molecular mass of the gas and R is the molar gas constant.

3.2 Gas Phase Balance Equations

When deriving the gas balance equations, it is important to remember that the axial velocity of the gas in the channels is not uniform across the width of the channels, but rather varies across the channels, being fastest along the centre of the channels and slowest at the channel walls. Thus, the variables [V.sub.i] and [V.sub.o] represent the average gas velocity in the inlet and outlet channels (respectively) at a given axial location. The average velocity is given by the volumetric flow divided by the cross-sectional area of the channel. The mass flux of the gas in the channels, [[phi].sub.i] and [[phi].sub.o], will also be referred to; again these quantities refer to averages across the channel (mean mass flux = mass flow in channel divided by the cross-sectional area of the channel). In the derivation, it is necessary to refer to the variation in gas velocity across the channel. In this case the variable u will be used. Note that [V.sub.i] = <[u.sub.i]> and [V.sub.o] = <[u.sub.o]>, where the angle brackets indicate the average over the cross section of the channels.

As with gas velocity, the temperature of the gas in the channels will also vary across the width of the channels, e.g. because of heat transfer with the walls of the PF and the velocity distribution across the channels. Thus, the variables [T.sub.i] and [T.sub.o] refer to the mixing cup or bulk temperature of the gas in the inlet and outlet channels.

Variations in pressure across the width of the channels will be ignored. Thus, [P.sub.i] and [P.sub.o] refer to the pressure of the gas in the inlet and outlet channels at a given axial location; the same value applies across the whole width of the channels.

3.2.1 Gas Mass Balance

The mass balance for the inlet channels is given by:

[[partial derivative][d.sup.2.sub.i][[rho].sub.i][V.sub.i]/[[partial derivative]z]] = -4d[[rho].sub.g,W][V.sub.W] (6)

where [V.sub.W] is the superficial through-wall velocity and [[rho].sub.g,W] the density of the gas as it passes through the soot cake and the PF wall. The term on the left represents the net rate of mass loss per unit length of PF due to gas convection along the inlet channel and the term on the right the rate of mass loss per unit length due to flow through the wall. The "4" arises from the fact that the channels have four sides. Note that since [d.sub.i] varies along the length of the PF, [d.sub.i] must be included in the derivative.

The inlet channel mass balance can also be expressed in terms of the inlet channel and through-wall mass fluxes, [[phi].sub.i] and [[phi].sub.W], where [[phi].sub.i] = [[rho].sub.i] [V.sub.i] and [[phi].sub.W] = [[rho].sub.g,W] [V.sub.W]. Thus,

[[partial derivative][d.sup.2.sub.i][[phi].sub.i]]/[[partial derivative]z] = [-4d[[phi].sub.W]] (7)

Similarly, the outlet channel mass balance is given by,

[[partial derivative][[rho].sub.o][V.sub.o]]/[[partial derivative]z] = [4d[[rho].sub.g,W][V.sub.W]/[d.sup.2.sub.o]] (8)

Or in terms of the mass fluxes,

[[partial derivative][[phi].sub.o]/[[partial derivative]z]] = [4d[[phi].sub.W]/[d.sup.2.sub.o]] (9)

The outlet mass balance equations (Eq. (8) and (9)), are similar to those for the inlet channel balance (Eq. (6) and (7)), except for the change in sign of the wall flow term (gas is entering the channel through the walls rather than leaving) and that since [d.sub.o] is a constant it can be taken out of the derivative. Note that both [V.sub.W] and [[phi].sub.W] are referenced to the dimension d, so that they still appear with d in the outlet channel mass balance.

3.2.2 Gas Momentum Balance

Typically, PF models have used axial momentum balance equations resembling the following:

[[partial derivative][P.sub.i]]/[[partial derivative]z] = [1/[d.sup.2.sub.i]]/[[partial derivative][d.sub.i][[rho].sub.i][V.sup.2.sub.i]]/[[partial derivative]z] - [F[[micro].sub.i][V.sub.i]]/[d.sup.2.sub.i] (10)

[[partial derivative][P.sub.o]]/[[partial derivative]z] = [[partial derivative][[rho].sub.o][V.sup.2.sub.o]]/[[partial derivative]z] - [F[[micro].sub.o][V.sub.o]]/[d.sup.2.sub.o] (11)

where [[micro].sub.i] and [[micro].sub.o] are the viscosity of the gas in the inlet and outlet channels.

Each term in these equations can be interpreted as the net rate of gain in axial momentum per unit volume by different causes. Alternatively, they can be interpreted as the resultant force in the z-direction applied to the gas per unit volume by different sources. From left to right, the terms correspond to i) the force applied to the gas in a volume element by pressure; ii) the rate of momentum gain by convection (or alternatively the force required to accelerate the gas); and iii) the force due gas viscosity, which retards flow. The form of the viscosity term in explained in Appendix A.

The momentum convection terms in Eq. (10) and (11) describe the net rate of gain of z-momentum (i.e. momentum in the axial direction) per unit volume across a volume element due to gas entering and leaving in the z-direction. However, for the case of the inlet channel, this term neglects the loss in z-momentum due to the net flow of gas out of the element in directions perpendicular to the channel axis towards the wall. As a result, the net rate of change in momentum due to convection is overstated and this term is larger than it should be. The same applies for the outlet channel, but with the "signs reversed", i.e. the outlet channel term neglects the gain in z-momentum due to the net flow of gas into the element in directions perpendicular to the channel axis away from the walls and towards the centre of the outlet channel. Again, this results in the net rate of change in momentum due to convection being overstated and a term which is larger than it should be. Therefore, we propose a new momentum convection term, which is much closer to including the loss or gain of axial momentum in the direction perpendicular to the channels. The reason why this has been previously neglected, is that PF models are generally derived assuming that the flow field in the channels is not changed by the gas flow to/from the walls so that gas velocities in the direction perpendicular to the axis are taken to be zero.

In this development, only the momentum convection term will be considered. The momentum convection term from the most general form of the momentum balance equation for fluids is -[nabla] * [rho]uu [32]. The z-component of this is:

- [[partial derivative][rho][u.sub.x][u.sub.z]]/[[partial derivative]x] - [[partial derivative][rho][u.sub.y][u.sub.z]]/[[partial derivative]y] - [[partial derivative][rho][u.sub.z][u.sub.z]]/[[partial derivative]z] (12)

Applying the product rule to each derivative gives:

- [rho][u.sub.x][[partial derivative][u.sub.z]]/[[partial derivative]x] - [rho][u.sub.y][[partial derivative][u.sub.z]]/[[partial derivative]y] - [rho][u.sub.z][[partial derivative][u.sub.z]]/[[partial derivative]z]***

- [u.sub.z]([[partial derivative][rho][u.sub.x]]/[[partial derivative]x] + [[partial derivative][rho][u.sub.y]]/[[partial derivative]y] + [[partial derivative][rho][u.sub.z]]/[[partial derivative]z]) (13)

The mass balance under steady-state conditions is [nabla] * [rho]u = 0 [32]. Therefore, the bracketed term in the above equation is zero, which leaves:

- [rho][u.sub.x][partial derivative][u.sub.z] / [partial derivative]x - [rho][u.sub.y][partial derivative][u.sub.z] / [partial derivative]y] - [rho][u.sub.z][partial derivative][u.sub.z] / [partial derivative]z (14)

Along the centre line of the channel, [u.sub.x][partial derivative][u.sub.z]/[partial derivative]x and [u.sub.y][partial derivative][u.sub.z]/[partial derivative]y are zero on symmetry grounds; away from the centre line these terms will be non-zero. The exact values of these terms will depend on the flow field (i.e. the value of u as a function of position). However, for the inlet channels it seems likely that [u.sub.x][partial derivative][u.sub.z]/[partial derivative]x and [u.sub.y][partial derivative][u.sub.z]/[partial derivative]y will increase in magnitude with proximity to the channel walls; [u.sub.x] and [u.sub.y] will likely increase towards the walls with the gathering flow to the walls, while the typical parabolic-like velocity profile associated with laminar flow (which still broadly occurs for square channels with wall flow, albeit with more some more complicated behaviour close to the walls of the inlet channels [27,28]) will result in [partial derivative][u.sub.z]/[partial derivative]x and [partial derivative][u.sub.z]/[partial derivative]y (mostly) increasing in magnitude from the centre line to the walls.

Nevertheless, we chose to neglect [u.sub.x][partial derivative][u.sub.z]/[partial derivative]x and [u.sub.y][partial derivative][u.sub.z]/[partial derivative]y, making the usual assumption in PF modelling that the flow field in the channel is unaffected by the (relatively slow) flow to/from the walls, which would result in these terms being zero. Neglecting these terms at this stage is conceptually no different from neglecting [partial derivative][rho][u.sub.x][u.sub.z]/[partial derivative]x and [partial derivative][rho][u.sub.y][u.sub.z]/[partial derivative]x in Eq. (12) as is done to give the momentum convection terms in Eq. (10) and (11). This leaves just one term:

-[rho][u.sub.z][partial derivative][u.sub.z]/[[partial derivative]z] (15)

Strictly speaking, this term only refers to the convection of z-momentum at a single set of coordinates. For including in the z-momentum balance equation of a one-dimensional model, this equation needs to be averaged over the channel cross section by integrating over the cross section and dividing by the cross-sectional area of the channel. To do this, it is assumed that the z-component of the velocity ([u.sub.z]) can be expressed as [u.sub.z] = <[partial derivative][u.sub.z]> f(x,y) [27,28].

Substitution of this into Eq. (15) gives:

-[f.sup.2][rho]<[u.sub.z]>[partial derivative]<[u.sub.z]>/[[partial derivative]z] (16)

Note that in this equation only f is a function of x and y, so it is only this term that needs averaging over the cross section. Averaging Eq. (16) over the cross section of the channel gives:

-[alpha][rho]<[u.sub.z]>[partial derivative]<[u.sub.z]> / [[partial derivative]z] (17)

where [alpha]= [[[integral][integral].sub.Channel][f.sup.2]dxdy] / [[[integral][integral].sub.Channel]dxdy] = <[u.sup.2.sub.z]> / [<[u.sub.z]>.sup.2]

The right hand term in this expression is obtained by substituting [u.sub.z] = <[u.sub.z]> f(x,y) for f in the integral. Equation (17) is the final form of the improved momentum convection term. Thus the new momentum balance equations are:

[partial derivative][P.sub.i] / [[partial derivative]z] = -[[alpha][[rho].sub.i][V.sub.i]] / [d.sup.2.sub.i][[partial derivative][d.sup.2.sub.i][V.sub.i]]/[[partial derivative]z] - [F[[mu].sub.i][V.sub.i]]/[d.sup.2.sub.i] (18)

[partial derivative][P.sub.o] / [[partial derivative]z] = -[alpha][[rho].sub.o][V.sub.o] / [[partial derivative][V.sub.o]]/[[partial derivative]z] - [F[[mu].sub.o][V.sub.o]]/[d.sup.2.sub.o] (19)

The [alpha] term in these equations is the momentum flux correction factor which serves to correct for the non-flat velocity profile across the channel (remember <[u.sup.2]> [not equal to] [<u>.sup.2]). As Bissett, Kostoglou and Konstandopoulos [27,28] have pointed out, such a term should also appear in the momentum convection terms of the "old" equations (Eq. (10) and (11)), but this is rarely seen in the literature. For developed laminar flow in a square channel with non-porous walls, a=1.377 [28].

Another way of looking at the "old" and "new" momentum balance equations is that if one starts with the momentum convection term from most general form of the momentum balance equation for fluids, -[nabla] * puu [32], and takes the term for convection of z-momentum in the z-direction, one gets [partial derivative][rho][u.sup.2.sub.z] / [partial derivative]z, the term in the old momentum balance (Eq. (10) and (11)). If, on the other hand, one starts with the momentum convection term from the Navier-Stokes equation, [rho]u * [nabla]u [32], and takes the term for convection of z-momentum in the z-direction, one gets [rho][u.sub.z][partial derivative][u.sub.z] / [partial derivative]z, the term in the new momentum balance (Eq. (18) and (19)). Note that the use of this Navier-Stokes term does not require the assumption of incompressible flow since only steady-state conditions are being considered. (1)

The only other model we have found to consider the convection of z-momentum in the transverse direction was by Wurzenberger and Kutschi [20]. This model includes the effect in the inlet channel momentum balance, but not in the outlet channel balance. Their momentum balance equations like most neglect the velocity profile across the channels. However, although differently expressed, their inlet channel momentum balance is mathematically equivalent to Eq. (18) with [alpha]=1, except that they have [d.sup.2.sub.i] appearing outside of the derivative in the momentum convection term instead of inside.

It is useful to compare the magnitude of the "new" and "old" momentum convection terms. Under isothermal conditions and when the pressure change across a volume element is not large enough to cause a significant change in gas density, [rho] can be considered a constant. Under these conditions, [partial derivative][rho][V.sup.2] / [partial derivative]z (the "old term") is equal to 2[rho]V [partial derivative]V/[partial derivative]z, which is twice the new momentum convection term.

A comparison of predicted backpressure as a function of flow rate for the new (Eq. (18) and (19)) and old (Eq. (10) and (11)) momentum balance equations for both [alpha]=1 (flat velocity profile) and [alpha]=1.377 (realistic velocity profile) is shown in Fig. 4. At low flow rates, the lines coalesce indicating that the momentum convection term is insignificant under these conditions. At higher flows, it is clear that the new momentum convection term results in a lower backpressure prediction than the old one and that allowing for a velocity profile across the channel ([alpha]=1.377) results in a higher prediction than if a flat velocity profile is assumed ([alpha]=1).

Since the gas velocity decreases along the inlet channel, the momentum convection term in both old and new the momentum balance equations (Eq. (10) and (18)) is positive. As a result, at sufficiently high flow rates, these equations will predict (perhaps surprisingly) that pressure will increase along the inlet channel. As the new momentum convection term is smaller than the old one the tendency for this to happen is reduced, but the effect is not prevented.

Since [d.sub.i] varies along the length of the PF, [d.sub.i] must be included in the derivative of the momentum convection term ([partial derivative][d.sup.2.sub.i][V.sub.i]/[partial derivative]z or [partial derivative][d.sup.2.sub.i][[rho].sub.i][V.sub.i]/[partial derivative]z)in the inlet channel momentum balance (both new and old forms, Eq. (10) and (18)), although not everybody does this [17,18,19,23,25]. In contrast, [d.sub.i] should not appear in the [partial derivative]P / [partial derivative]z term. Most workers have written their [partial derivative]P / [partial derivative]z term without [d.sub.i] in the derivative [17,18,19,20,21,23,25] (although there are exceptions [22]), but, in general, the reason for this has not been made clear. The only explanation we have been able to find is that this is required "to balance the forces" [21,26]. We have found that including [d.sub.i] in the pressure derivative results in the prediction unphysical pressure and velocity profiles during soot loading, so this is an important point.

Therefore, it seems worth explaining why [d.sub.i] does not appear in the [partial derivative]P/[partial derivative]z term. Consider Fig. 5, which represents two elements of the inlet channel of length [delta]z across which the width of the channel changes due to variation in the soot loading. The forces acting on the gas in the left hand element in the axial (z) direction are indicated. The net (or resultant) force acting on the gas in the axial direction is given by:

Net force = [A.sub.z][P.sub.z] - [A.sub.z + [delta]z][P.sub.z + [delta]z] - ([A.sub.z] - [A.sub.z + [delta]z]) [P.sub.z + [delta]z] (20)

where [A.sub.z] & [A.sub.z + [delta]z] and [P.sub.z] & [P.sub.z] +[delta]z are the cross-sectional area and the pressure at positions z & z+[delta]z. The first two terms after the equals sign are the forces acting on the gas in the element from the left by the upstream gas and from the right by the downstream gas. The third term is the force on the gas exerted by its "container", in this case the (soot-lined) channel walls. Since [[delta].sub.z][right arrow]0, in calculus terms the channel width changes in a stepwise fashion across the element (Fig. 5), rather than there being a continuous change in cross section.

Simplifying the equation:

Net force = [A.sub.z]([P.sub.z] - [P.sub.z + [delta]z]) (21)

Dividing this by [delta]z and as [delta]z[right arrow]0, gives the derivative -A[partial derivative]P/[partial derivative]z, i.e. the pressure derivative with the cross-sectional area appearing outside of the derivative. Note that the same conclusion would be reached if the channel width changed in the opposite direction to that depicted in Fig. 5. Incidentally, this consideration did not arise in Bissett's original model [15] as he assumed that the inlet channel width did not change with soot loading.

The momentum balance equations (Eq. (18) and (19)) have so far been presented in a form that is not practical for solving with a Boundary Value Problem (BVP) solver as the derivatives [partial derivative][d.sup.2.sub.i][V.sub.i]/[partial derivative]z, and [partial derivative][V.sub.o]/[partial derivative]z are not readily available. For the inlet channel balance, this can be remedied by writing [partial derivative][d.sup.2.sub.i][V.sub.i]/[partial derivative]z as [partial derivative]([d.sup.2.sub.i][[phi].sub.i]/[[rho].sub.i])/[partial derivative]z, substituting in Eq. (5) for the density and applying the product rule to the resulting derivative. Substitution of this result into Eq. (18) gives:

(1 - [alpha] [phi].sub.i][V.sub.i] / [P.sub.i])[partial derivative][P.sub.i] [partial derivative]z = -[V.sub.i][[alpha](1 / [d.sup.2.sub.i][partial derivative][d.sup.2.sub.i][[phi].sub.i] / [partial derivative]z + [[phi].sub.i] / [T.sub.i][partial derivative][T.sub.i] / [partial derivative]z) + F[[mu].sub.i] / [d.sup.2.sub.i]] (22)

This equation contains only derivatives that are known; it defines [partial derivative]P/[partial derivative]z, while the other derivatives are given by Eq. (7) and (26). A similar expression applies for the outlet channel:

(1 - [[alpha][[phi].sub.o][V.sub.o] / [P.sub.o])[partial derivative][P.sub.o] / [partial derivative]z] = -[V.sub.o][[alpha]([partial derivative][[phi].sub.o] / [partial derivative]z + [[phi].sub.o] / [T.sub.o] [partial derivative][T.sub.o] / [partial derivative]z) + F[[mu].sub.o] / [d.sup.2.sub.o]] (23)

3.2.3 Gas Energy Balance

In terms of the specific enthalpy, [H.sub.i], the gas energy balance for the inlet channel is:

[[partial derivative][d.sup.2.sub.i][[phi].sub.i][H.sub.i]]/[[partial derivative]z] = 4[d.sub.i][h.sub.i]([T.sub.w] - [T.sub.i]) - 4d[[phi].sub.w][C.sub.p,g,w][T.sub.w] (24)

where [h.sub.i] is the heat transfer coefficient for heat transfer between the gas in the inlet channel and the PF wall, [T.sub.w] is the temperature of the soot cake and the PF wall and [C.sub.pg,w] is the specific heat capacity at constant pressure of the gas passing through the soot cake and the PF wall.

From left to right, the terms in Eq. (24) correspond to the rate of enthalpy loss per unit length of PF by convection along the inlet channel; the rate of energy gain per unit length due to heat transfer between the gas in the channel and the PF wall (or soot cake for a soot-loaded PF); and the rate of enthalpy loss per unit length due to flow through the wall.

Substituting in d[H.sub.i] = [C.sub.pg,i] d[T.sub.i.sup.2], where [C.sub.pg,i] is the specific heat capacity of the gas in the inlet channel, gives,

[C.sub.pg,i] [[partial derivative][d.sup.2.sub.i][[phi].sub.i][T.sub.i]]/[[partial derivative]z] = 4[d.sub.i][h.sub.i]([T.sub.w] - [T.sub.i]) - 4d[[phi].sub.w][C.sub.pg,w][T.sub.w] (25)

Note that [C.sub.pg,i] appears outside the derivative. Applying the product rule to the derivative and substituting in Eq. (7), gives the final equation for the inlet channel gas energy balance,

[mathematical expression not reproducible] (26)

Analogously, the outlet channel energy balance is:

[mathematical expression not reproducible] (27)

This differs from the inlet channel balance predominantly in the effective change in the sign of the [[phi].sub.w] term as gas is entering the channel through the wall rather than leaving it.

The heat transfer coefficients, [h.sub.i] and [h.sub.o], are calculated from the Nusselt number, [N.sub.u], and the thermal conductivity of the gas in the inlet and outlet channels ([[lambda].sub.i],[[lambda].sub.g]), thus from the definition of [N.sub.u]:

[d.sub.i][h.sub.i] = Nu [[lambda].sub.i] (28)

[d.sub.o][h.sub.o] = Nu [[lambda].sub.o] (29)

In this model, Nu is assumed to be a constant. The effect of Nu being a function of the wall Reynolds number, [Re.sub.W], will be considered in section 6.

In deriving the inlet channel energy balance, it has been assumed that the gas leaving the inlet channels through the walls has the same temperature as the PF wall (Eq. 24)). This recognises that the gas temperature varies across the width of the channel with the gas close to the PF wall (or soot cake, if there is one) having the same temperature as the wall (or soot cake) due to efficient heat transfer with the wall and the low velocity of the gas close to the wall (due to the no slip boundary condition (i.e. gas in the channel adjacent to the wall having zero axial velocity) and the low through-wall velocity resulting from the large filtration area). This approach was used in Bissett's original model [15] and has been used by other workers [24,35]. However, this assumption also results in an apparent paradox in that if [[phi].sub.w] is sufficiently large, Eq. (26) predicts that [T.sub.i] will increase along the inlet channel when [T.sub.i] > [T.sub.w], which clearly contradicts thermodynamics. To avoid this apparent paradox, many workers have developed models in which it is assumed that the gas leaving the inlet channel through the walls has a temperature of [T.sub.i] instead of [T.sub.w] [18,21,22,23,25,26]. More recently, Bissett, Kostoglou and Konstandopoulos [27,28] have resolved this apparent paradox by showing that [h.sub.i] in fact increases with [[phi].sub.w] (more specifically with [Re.sub.w]) at a rate which ensures that [partial derivative]T/[partial derivative]z can never have "the wrong sign". Therefore, there is no reason not to use the original approach where the gas leaving the channel through the PF wall has a temperature of [T.sub.w] and indeed Bissett et al. [27,28] argue that this approach is preferable for a number of reasons. (However, others remain unconvinced [26], stating that the Bissett, Kostoglou and Konstandopoulos' work relies on a fixed gas temperature at the channel wall boundary, the value of which is not known a priori).

As Bissett et al. [27] have pointed out, the terms in the final form of the inlet channel energy balance (Eq. (26)) correspond (from left to right) to i) the enthalpy change in the portion of the gas entering a volume element that continues along the channel, ii) the total heat flux at the walls (due to conduction through the gas [15,28]), and iii) the change in enthalpy of the remaining portion of the gas entering the volume element which leaves through the PF wall. For the latter term, the enthalpy change is associated with the fact the gas enters the element with a mixing cup temperature of [T.sub.i], but leaves at a temperature of [T.sub.w]. With the first term, the enthalpy change is associated with the change in mean gas temperature ([T.sub.i]) across the element. This physical interpretation is possible because of the splitting the derivative in Eq. (25) into terms corresponding to gas passing along the channel and gas exiting through the wall.

3.2.4 Boundary Conditions for Gas Balance Equations

Since there are six differential equations (Eq. (7), (9), (22), (23), (26), (27)), six boundary conditions are required. For the on-vehicle configuration (Fig. 1, bottom), these are:

[d.sup.2.sub.i][[phi].sub.i][|.sub.z=0] = 2[??]/[rho]cell A (30)

[[phi].sub.i][|.sub.z=L] = 0 (31)

[P.sub.0][|.sub.z=0] = 0 (32)

[P.sub.0][|.sub.z=L] = ([P.sub.5] - [P.sub.6]) + [P.sub.atm] (33)

[T.sub.i][|.sub.z=0] = [T.sub.In] (34)

[T.sub.0][|.sub.z=0] = [T.sub.w][|.sub.z=0](35)

where A is the cross-sectional area of the PF, [??] is the total mass flow into the PF and ([P.sub.5] - [P.sub.6]) is the pressure change due to expansion of gas as it leaves the PF; [P.sub.5] and [P.sub.6] refer to pressures at specific locations defined in Fig. 6. The evaluation of ([P.sub.5] - [P.sub.6]) is described in section 3.5.1.

In words, the boundary conditions are that the inlet mass flux and temperature are known (Eq. (30) and (34)); that the gas mass flux at the end of the of the inlet channel and the start of the outlet channel is zero (Eq. (31) and (32)); that the pressure at the end of the outlet channel is equal to the external pressure plus the change in pressure due to expansion as the gas leaves the PF (Eq. (33)); and that since the gas at the start of the outlet channel (if there were any) is stagnant (Eq. (32)) it must have the same temperature as the PF wall at the same axial location (Eq. (35)).

For the cold flow configuration (Fig. 1, top), where atmospheric pressure is at the PF inlet rather than the outlet, Eq. (33), is replaced by:

[P.sub.i][|.sub.z=0] + ([P.sub.1] - [P.sub.2]) = [P.sub.atm] (36)

where ([P.sub.1] - [P.sub.2]) is the pressure loss due to contraction of the gas as it enters the PF (Fig. 6). The evaluation of this is described in section 3.5.2. In words, this boundary condition is that the pressure at the start of the inlet channel plus the pressure change due to gas contraction at the PF entrance is equal to the external pressure.

3.3 Across-Wall Pressure Drop

The pressure drop across the soot cake and PF wall can be described by Darcy's law with the Forchheimer extension. Since gas density and hence through-wall gas velocity (at constant mass flow) vary with pressure, this equation is written in differential form:

- [partial derivative]P / [partial derivative]y = [mu][V.sub.W] / k + [beta][[rho].sub.g,W][V.sub.2.sub.W] (37)

where y is the distance through the soot cake and PF wall, k is the permeability and [beta] is the Forchheimer constant. Both k and [beta] are properties of the porous media through which the gas flows. The first term on the right hand side of the equals sign is Darcy's law and the second term the Forchheimer term.

Darcy's law (3) [37] and the Forchheimer extension (4) [38] are both empirical laws originally developed for the flow of water across beds of sand. However, they have been extensively used for the prediction of pressure drop across all sorts of porous media. Both Darcy's law [39,40] and Darcy's law with the Forchheimer extension [41,42] have subsequently been derived from first principles. Darcy's law describes pressure drop across the porous media due to viscous losses and the Forchheimer extension the pressure drop due to inertial effects.

Equation (37) needs to be solved for the soot cake and PF wall. The soot cake in the inlet channel can be divided into four trapezoidal sections (Fig. 2), each of which sits on a slab shaped section of wall. This geometry, together with the dimensions, is depicted in Fig. 7. In this work we only consider flow through the wall of symmetric PFs, where the wall cross section is rectangular; asymmetric PFs represent an additional level of complexity as the flow path through the PF wall is non-trivial [35]. Since this model allows for the variation in gas density (and hence through-wall velocity) with pressure, it is necessary to rewrite Eq. (37) in terms of the through-wall mass flux. Thus, for the PF wall Eq. (37) becomes:

- [partial derivative]P / [partial derivative]y = [[mu].sub.W][[phi].sub.W] / [k.sub.W][[rho].sub.g,W] + [[beta].sub.W][[phi].sup.2.sub.W] / [[rho].sub.g,W] + [RT.sub.W] / PM[[mu].sub.W][[phi].sub.W] / [k.sub.W] + [[beta].sub.W][[phi].sup.2.sub.W] (38)

In the right hand expression, Eq. (5) has been substituted for [[rho].sub.g,W*] [k.sub.w] and [[beta].sub.w] refer to the permeability and Forchheimer constant of the wall.

Similarly, for the soot cake Eq. (37) becomes:

-[[partial derivative]P]/[[partial derivative]y] = [RT.sub.W]/[PM][[[mu].sub.W][[phi].sub.W]]/[k.sub.S][d/[d.sub.i] + [2y]] + [[beta].sub.S][[phi].sup.2.sub.W][d.sup.2]/[([d.sub.i] + [2y]).sup.2](39)

where [k.sub.s] and [[beta].sub.s] refer to the permeability and Forchheimer constant for the soot cake. Note that [d.sub.i] + 2y is the width of the soot cake at a distance y from the "top" (or channel side) of the soot cake as depicted in Fig. 7. Here it has been assumed that the flow through the soot cake can be modelled by assuming that the gas flow simply spreads out to fill the expanding cross-sectional area of the soot cake. Thus the gas velocity through the soot cake varies with d/([d.sub.i] + 2y). Note also that since the temperature is assumed to be constant through the soot cake and PF wall, [[micro].sub.w] and [T.sub.w] apply to both Eq. (38) and (39).

Equations (38) and (39) are first order differential equations with separable variables, which can be solved to give:

[P.sup.2.sub.i] - [P.sup.2.sub.o] = [RT.sub.w] / M [([2w.sub.w] / [k.sub.w] + d / [k.sub.s]]ln([d / [d.sub.i]))[[mu].sub.w][[phi].sub.w] + *** + (2[[beta].sub.w][w.sub.w] + [[beta].sub.s][d.sup.2](1 / d.sub.i - 1 / d))[[phi].sup.2.sub.w]] (40)

Thus the pressure drop across the soot cake and PF wall, is a quadratic in [[phi].sub.w] This can be solved using the usual equation for the solution of a quadratic to give [[phi].sub.w] as a function of [P.sub.i] and [P.sub.o], which is what is actually required for numerically solving the gas balance equations.

Finally, it is worth noting that the plugs of the PF can be modelled by setting [[phi].sub.W] = 0 in the plugged region. With a plug in the outlet channel, Eq. (7), (9) and (32) mean that [d.sup.2.sub.i][[phi].sub.i] will be constant and [[phi].sub.o] will be zero in the plugged region. Similarly, with an inlet channel plug, [[phi].sub.o] will be constant and [[phi].sub.i] will be zero in the plugged region.

3.4 Solid Phase Energy Balance Equation

The solid phase (PF wall plus soot cake) energy balance is given by:

[mathematical expression not reproducible] (41)

where [[rho].sub.W] is the density of the PF wall; [C.sub.pw] and [C.sub.ps] are the specific heat capacities of the PF wall and the soot cake; [[lambda].sub.w] and [[lambda].sub.s] are the thermal conductivities of the PF wall and soot cake; and Q is the rate of heat generation per unit volume by reaction. The "2" arises because each square channel has four walls for heat transfer, but only half of the channels are inlet or outlet channels.

All the terms in this equation correspond to the rate of gain in energy of the solid phase per unit volume of PF by different causes. From left to right, the terms on the right hand side correspond to i) axial conduction of heat along the PF walls, ii) axial conduction of heat along the soot cake, iii) heat transfer from the gas in the inlet and outlet channels and iv) heat generation by reaction (included here for completeness).

(1-[epsilon]) is the fraction of the cross section occupied by the PF wall and S/[[rho].sub.S] the fraction of the cross section occupied by the soot cake; these give the fraction of the cross section available for conduction.

In principal, Eq. (41) should also have terms for the energy gain due to the gas entering the soot cake (if there were one) or the wall from the inlet channel (equal to +[2[rho].sub.cell] d[[phi].sub.w][C.sub.p,g,W] [T.sub.W]) and for gas leaving the PF wall and entering the outlet channel (equal to -[2[rho].sub.cell] d[[phi].sub.w][C.sub.p,g,W] [T.sub.W] However, because it has been assumed that gas entering the PF wall (or soot cake, if there is one) has a temperature of [T.sub.W] (see section 3.2.3), these terms cancel.

Since the solid energy balance has a second order spatial derivative, two spatial boundary conditions are required. These are that there is no heat flux due to conduction beyond the ends of the PF, i.e.:

[partial derivative][T.sub.w] / [partial derivative]z[|.sub.z=0,L] = 0 (42)

The solid energy balance given above applies for an uncoated PF with a soot cake. However, the thermal mass of a washcoat or stored ash could easily be added in the same way as has been done for soot. Similarly, thermal conductivity along a washcoat or ash layer could be handled in the same way as soot conductivity. The thermal mass of the plugs can be included by adding another term to the bracket in front of [partial derivative][T.sub.W]/[partial derivative]t, which is only non-zero in the plugged region. Opitz et al. [24] report that the thermal mass of the plugs has a significant impact on GPF light-off.

The sold energy balance (Eq. (41)) relies on two assumptions: firstly, that the temperature is uniform across the soot cake and PF wall at any given axial location; and secondly, that the gas entering the soot cake or PF wall from the inlet channel already has the same temperature as the soot cake and PF wall due efficient heat transfer between the gas in the channel (adjacent to the soot cake or wall) and the soot cake or PF wall. Without the latter assumption, the terms for heat convection into and out of the wall would not cancel. Together these assumptions mean that the gas passing through the soot cake and the PF wall has a constant temperature (at a given axial location), which is the same as that of the soot cake and PF wall (at the same axial location). These assumptions have been commonly applied [26]. Calculations by Bissett and Shadman [43] and Depcik and Assanis [21] suggest that gas passing through the soot cake or wall very rapidly takes on the temperature of the soot cake or wall, validating the assumption that the gas passing through the soot cake or wall has the same temperature as the soot cake or the wall. Bissett and Shadman's [43,44] calculations also support the assumption of uniform temperature across the soot cake and wall. On the other hand, simulations by Haralampous and Koltsakis [45] show significant temperature gradients across the soot cake during active regeneration with very high initial soot loadings (> 10 g [L.sup.-1]), suggesting that the assumption of uniform temperature across the soot cake and wall may not hold under extreme conditions. Despite the prediction of a temperature gradient across the soot cake, practically uniform temperature was still predicted across the wall due to its higher thermal conductivity [45]. Of course, whether or not a temperature gradient is predicted across the soot cake depends on the value used for the thermal conductivity of the soot.

3.5 Entrance and Exit Effects

Pressure changes associated with the contraction of the gas flow at the entrance to the PF and expansion of the gas as it exits the PF can have a significant contribution to the total pressure drop across a PF at higher flows (see section 5). Despite this, apart from a couple of CFD studies [19,35,46], entrance and exit effects have received relatively little attention. These pressure changes are generally modelled using expressions of the form [zeta][rho][V.sup.2]/2, where [zeta] is an expansion or contraction loss coefficient ([[zeta].sub.o] or [[zeta].sub.i], respectively). When modelling PFs, these coefficients have either been taken as empirical fitting constants [19,36,47], calculated from CFD [19,35,46], obtained using specified [17] or unspecified [22,48] standard relationships from the literature or given an unspecified value [20,23]; alternatively, these effects have been ignored altogether [15,18,21,25]. Where standard relationships have been used, these are typically for abrupt contraction and expansion in pipes [49], despite the fact that these relationships are only valid at high Re (fully developed turbulent flow), where the velocity profile can be assumed to be flat across the channel/pipe; the flow in the channels of a PF is laminar, so the velocity profile is clearly not flat and hence these relationships do not apply. When using these standard relationships, it is also important to be clear whether they are for pressure loss or mechanical energy loss.

The expressions derived in this section for the pressure changes at the entrance and exit of the PF are based on work by Kays [50] on pressure drops at the entrance and exit of multi-tube heat exchangers; these equations have been validated against experimental data [50]. A square channelled heat exchanger is geometrically similar to a PF at the entrance and exit, so this approach seems reasonable.

The following assumptions will be made when deriving expressions for pressure change during expansion and contraction:

* Incompressible flow, i.e. gas density remains constant during expansion and contraction.

* Inertial effects are assumed to dominate during gas expansion and contraction, so that viscous losses can be neglected, i.e. the flow is assumed to be inviscid.

* No pressure variation across width of channel

* Developed flow is achieved in the inlet channels in a length which is very short compared with the length of the PF.

* The velocity profile across the channels is unaffected by the flow to/from the porous walls of the channels.

The assumption of incompressible inviscid flow is made largely for convenience as it considerably simplifies the problem, enabling analytical expressions to be derived. However, the error due to the assumption of incompressible flow can at least be estimated by taking the difference in predicted entrance or exit pressure change if the density calculated with the upstream pressure is used instead of the downstream pressure (or vice versa) and dividing by the predicted entrance or exit pressure change calculated with the density normally used. For the simulations presented in this work with mass flows typical of a real application for the size of PF used (m [less than or equal to] 0.06 kg [s.sup.-1]), which includes the WLTC simulations, the maximum estimated error in the entrance and exit effect pressure change is no more than 0.7 and 0.3%, respectively, which would seem acceptable. For simulations going to higher mass flow (0.3 kg [s.sup.-1]), the maximum errors were somewhat larger, at 4.3 and 1.8% for the entrance and exit pressure changes (respectively) for the on-vehicle configuration (Fig. 1, bottom) and 5.0% for the entrance pressure change for the cold-flow configuration (Fig. 1, top); as will be shown in the next section, the exit pressure change predicted for the cold-flow configuration is zero.

3.5.1 Gas Expansion at the PF Exit

The situation of the gas as it expands on leaving the PF is shown schematically in Fig. 8. The streamlines are unable to follow the abrupt change in conduit shape as the gas leaves the PF, but take a smooth course instead. An equation for the pressure change can be derived by applying the macroscopic mass and momentum balance equations [49] between "plane 5", located at the outlet channel just before the gas exits and where the flow is unaffected by the forthcoming expansion (marked in Fig. 6 and 8), and "plane 6", where there is fully developed flow after the expansion.

The mass balance for the expansion is:

[A.sub.5][rho][V.sub.5] = [A.sub.6][rho][V.sub.6] (43)

where [V.sub.5] and [V.sub.6] are the mean velocities at locations 5 and 6 (Fig. 6 and 8) and [A.sub.5] and [A.sub.6] are the cross-sectional areas available for flow at the same locations. Note that [A.sub.5] is the total cross-sectional area for all the outlet channels in the PF. Obviously, the left and right sides of this equation correspond to the gas just before the expansion and just after completion of the expansion (re-establishment of developed flow). Since incompressible flow is assumed, the density ([rho]) is the same on both sides of the equation; this cancels, but is included for completeness.

The equation for macroscopic conservation of momentum [49] in the axial direction is:

[A.sub.5][rho][[alpha].sub.5][V.sup.2.sub.5]+[A.sub.5][P.sub.5] + ([A.sub.6]-[A.sub.5])[P.sub.5] = [A.sub.6][rho][[alpha].sub.6][V.sup.2.sub.6]+[A.sub.6][P.sub.6] (44)

The first term on each side of the equals sign is the rate of flow of momentum into (left side) and out of (right side) the region under consideration (i.e. the region between planes 5 and 6 in Fig. 6 and 8). The second term on each side is the force applied to the gas in the region of interest by gas upstream (left) and downstream (right) of the region. The third term on the left is the force exerted on the gas in the region by the solid portion of the rear face of the PF. For subsonic flow, the pressure at this face is equal to the static pressure just before expansion, i.e. [P.sub.5] [50,51]. [[alpha].sub.5] and [[alpha].sub.6] are correction factors to allow for the fact that gas velocity is not constant across the cross section; remember <[u.sup.2]> [not equal to] <[u.sup.2]>. By definition, [[alpha].sub.5] = <[u.sup.2.sub.5]>/[<[u.sub.5]>.sup.2] = <[u.sup.2.sub.5]>/[V.sup.2.sub.5]; [[alpha].sub.6] is similarly defined.

Solution of Eq. (43) and (44) gives an expression for the pressure change due to the expansion at the PF exit:

[P.sub.5] - [P.sub.6] = [[sigma].sub.0]([[alpha].sub.6][[sigma].sub.0] - [[alpha].sub.5]) [[rho].sub.0][V.sup.2.sub.0][|.sub.z=L] (45)

where [[sigma].sub.0] (= [A.sub.5] / [A.sub.6]) is the expansion ratio, i.e. the ratio of the cross-sectional area for flow before and after the expansion. Note that [[rho].sub.0][V.sup.2.sub.0][|.sub.z=L] = [rho][V.sup.2.sub.5]

The value of [[sigma].sub.0] depends on the configuration. For the cold flow configuration (Fig. 1, top), the gas leaving the PF expands into a plenum much larger in cross section than the PF (see section 2.1), so [[sigma].sub.0] is effectively zero. For the on-vehicle configuration (Fig. 1, bottom), [[sigma].sub.0] is the open area fraction of the rear face of the PF ([[epsilon].sub.0]). To summarise:

[mathematical expression not reproducible] (46)

Since [[sigma].sub.0] is zero for the cold flow configuration, the pressure change due to expansion at the PF exit is zero, according to Eq. (45). To evaluate the exit effect pressure change for the on-vehicle configuration, values are required for [[alpha].sub.5] and [[alpha].sub.6]. Calculation of [[alpha].sub.5] and [[alpha].sub.6] requires knowledge of the velocity profile. Since the flow in the channels is laminar we assume [[alpha].sub.5]=1.377, the value for developed laminar flow in a square channel with non-porous walls [28]; for laminar flow this value is independent of the gas velocity. The Reynolds number for the gas in the circular conduit (or rather the canning) after the PF is sufficiently high for the flow to be turbulent. At high Re, a flat velocity profile is expected so [[alpha].sub.6]=1. At intermediate Re, the velocity profile will not be completely flat and [[alpha].sub.6] will be greater than unity and a function of Re. There are correlations for the velocity profile as a function from Re from which [[alpha].sub.6] could be calculated, e.g. Kays [50] gives an equation for [[alpha].sub.6] as a function of Re obtained in this way. However, in this work [[alpha].sub.6] will be assumed to be 1 for simplicity; since simulation results for the on-vehicle configuration are not compared to measured data, the exact value used is not critical.

Since the velocity profile for the laminar flow in the channels is expected to be more pronounced that that for the turbulent flow in the canning following the PF, it is expected that [[alpha].sub.5] will be greater than [[alpha].sub.6]. This means that Eq. (45) predicts that [P.sub.6]>[P.sub.5] (except for [[sigma].sub.0]=0), i.e. that the pressure actually rises as the gas leaves the PF. The same conclusion would be reached for flat velocity profiles ([[alpha].sub.5]=[[alpha].sub.6]=1). This pressure rise is associated with the deceleration of the gas as it leaves the PF. Assuming there is a cone downstream of the PF, this pressure rise will be "lost" as the gas accelerates into the cone, the exact pressure change depending on how well mechanical energy is conserved during this process. Finally, note that while the expansion results in in a pressure rise, mechanical energy is still lost in the process.

It is interesting to see how the exit effect pressure change predicted by Eq. (45) depends on the values of [[alpha].sub.5] and [[alpha].sub.6]. This is shown in Fig. 9, which shows the expansion loss coefficient,

[[zeta].sub.0] ( = ([P.sub.5] - [P.sub.6])/1/2[[rho].sub.0][V.sup.2.sub.0][|.sub.z=L]) as a function of the expansion ratio. [[alpha].sub.5]=1.377, [[alpha].sub.6]=1 corresponds to the values used in the simulations in this work.

Comparison of the lines for [[alpha].sub.5]=1.377, [[alpha].sub.6]=1 with those for [[alpha].sub.5]=1, [[alpha].sub.6]=1, which correspond to a flat velocity profile in the channels, shows that neglecting the velocity profile in the outlet channels results in a considerable under prediction of the exit effect pressure change. [[alpha].sub.6]=1.333 corresponds to laminar flow in the (circular) canning after the PF [50,52]. While laminar flow in the canning is unlikely, comparison of the lines for [[alpha].sub.5]=1.377, [[alpha].sub.6]=1 and [[alpha].sub.5]=1.377, [[alpha].sub.6]=1.333 shows the absolute maximum (but practically unlikely) range of values that [[zeta].sub.0] can have depending on the velocity profile in the canning after the PF. In the range of expansion ratios of relevance for PFs (0.2< [[sigma].sub.0]=[[epsilon].sub.0] < 0.4 for on-vehicle configuration, [[sigma].sub.0]=0 for cold-flow (Eq. (46)), there is not that much difference between these extreme values of [[zeta].sub.0], i.e. in practice [[zeta].sub.0] is not that sensitive to the velocity profile in the canning. Finally, note that if [[alpha].sub.5][not equal to][[alpha].sub.6], [[zeta].sub.0] is predicted to be non-zero when expansion ratio is one, despite there being no change in the area for flow.

3.5.2 Gas Contraction at the PF Entrance

In flow across a sudden contraction [50,52], the fluid stream contracts to a minimum cross section (vena-contracta (VC)) somewhat narrower than the width of the channel, before re-expanding to fill the channels of the PF. The gas accelerates to the vena contracta before decelerating as the gas expands. This flow pattern results as the streamlines are unable to follow the sharp change in shape at the front face of the PF. According to the literature, the fluid acceleration is approximately isentropic up to the vena-contracta point and mechanical energy loss takes place predominantly during the deceleration associated with re-expansion after the vena-contracta point [50,52].

The situation of the gas as it contracts on entering the PF is shown schematically in Fig. 10. The derivation will consider what happens between "plane 1" (marked in Fig. 6 and 10), located upstream of the PF where the flow is unaffected by the forthcoming contraction and "plane 2", where there is fully developed flow after the expansion following the VC. The plane of the vena contracta, i.e. the point where the diameter of the stream is least and the velocity highest, is also important for the derivation.

The mass balance for the contraction is:

[A.sub.1][rho][V.sub.1] = [A.sub.VC][rho][V.sub.VC] = [A.sub.2][rho][V.sub.2] (47)

where the subscript VC indicates values at the vena contracta. As before, the density ([rho]) is the same at all points as incompressible flow has been assumed. Note that [A.sub.2] is the total cross-sectional area for all the inlet channels in the PF at z=0 (or at least a point near the front of the PF after the completion of the contraction). Similarly, [A.sub.VC] is the total of the cross-sectional areas of the VC in all of the inlet channels.

Since the contraction is approximately isentropic up to the vena contracta [50,52], mechanical energy (i.e. the kinetic and potential energy of the gas) is assumed to be conserved between plane 1 and the VC, which gives [49]:

[A.sub.1][[alpha].sub.e1][rho][V.sup.3.sub.1] / 2 + [A.sub.1][P.sub.1][V.sub.1] = [A.sub.vc][rho][V.sup.3.sub.VC] / 2 [A.sub.VC][P.sub.VC][V.sub.VC] (48)

This represents the rate of flow of mechanical energy (i.e. kinetic and potential energy) into the region of interest at plane 1 (left side of the equals sign) and out from the VC (right of the equals). The first term on each side the equals describes the rate of flow of kinetic energy and the second term the rate of flow of potential energy; this potential energy is due to the pressure. [[alpha].sub.e1] is a correction factor to allow for the fact that gas velocity is not constant across the cross section; [[alpha].sub.e1] = <[u.sup.3.sub.1]>/[<[u.sub.1]>.sup.3] Due to the nozzle like flow to the VC, the velocity profile is assumed to be essentially flat at the VC so no correction factor is required [50].

Applying the equation for macroscopic conservation of axial momentum [49] to the expansion of the gas from the VC to plane 2 gives:

[A.sub.VC][rho][V.sup.2.sub.VC] + [A.sub.VC][P.sub.VC] + ([A.sub.2] - [A.sub.VC])[P.sub.VC] = [A.sub.2][rho][[alpha].sub.2][V.sup.2.sub.2] + [A.sub.2][P.sub.2] (49)

Note that this equation is equivalent to Eq. (44) since it also describes momentum conservation during an abrupt expansion. [[alpha].sub.2] is a correction term to allow for the gas velocity at plane 2 not being constant across the cross section; as already stated, no such correction term is required for the velocity at the VC.

Combing Eq. (47), (48) and (49) to eliminate [P.sub.VC] gives an expression for the pressure change due to gas contraction at the PF entrance:

[P.sub.1]- [P.sub.2] = 2[[alpha].sub.2] - [[alpha].sub.e1][[sigma].sup.2.sub.i] - 2 / [C.sub.C] + 1 / [C.sup.2.sub.2] [[rho].sub.i][V.sup.2.sub.i][|.sub.z=0]/2] 50)

where [[sigma].sub.i] (=[A.sub.1] / [A.sub.2]) is the expansion ratio, i.e. the ratio of the cross-sectional area for flow after and before the contraction; and [C.sub.C](=[A.sub.VC] / [A.sub.2]) is the coefficient of contraction, i.e. the cross-sectional area of the VC divided by the area of the orifice through which gas flows to form the VC. Note that [[rho].sub.i][V.sup.2.sub.i][|.sub.z=0] = [rho][V.sup.2.sub.2].

The value of [[sigma].sub.i] depends on the configuration. For the cold flow configuration (Fig. 1, top), gas is effectively entering the PF from an infinite reservoir (see section 2.1), so [[alpha].sub.i] is zero. For the on-vehicle configuration (Fig. 1, bottom), [[alpha].sub.i] is the open area fraction of the front face of the PF ([[epsilon].sub.i]). To summarise:

[mathematical expression not reproducible] (51)

Evaluation of the pressure loss due to contraction (Eq. (50)) requires a value for [C.sub.C]. Calculated values for [C.sub.C] as a function of [[alpha].sub.i] have been taken from the literature [53]. These values were calculated for two-dimensional flow, but the reference indicates that these values give a good approximation to three-dimensional flow. Kays [50] used the same source for [C.sub.C]. These values for [C.sub.C] can be fitted to the following function (Fig. 11):

[C.sub.C] = 0.611-0.517[[sigma].sub.i] / 1-0.90[[sigma].sub.i] (52)

Values of [[alpha].sub.e1] and [[alpha].sub.2] are required for evaluating the pressure change due to contraction. [[alpha].sub.2] has the same value as [[alpha].sub.5], 1.377 [28]. For the cold flow configuration (Fig. 1, top), [[alpha].sub.i] is zero (Eq. (51)), so the value of [[alpha].sub.e1] is unimportant (Eq. (50)). The discussion above on the value of [[alpha].sub.6] (section 3.5.1) applies equally to [[alpha].sub.e1]. Thus, [[alpha].sub.e1] can be calculated using correlations for the velocity profile as a function of Re, e.g. Kays [50] gives an equation for [[alpha].sub.e1] as a function of Re. However, in this work [[alpha].sub.e1] will be assumed to be 1 for simplicity; since simulation results for the on-vehicle configuration are not compared to measured data, the exact value used is not important.

It is interesting to see how the entrance effect pressure loss predicted by Eq. (50) depends on the values of [[alpha].sub.e1] and [[alpha].sub.2]. This is shown in Fig. 12, which shows the contraction loss coefficient, [[zeta].sub.i] ( = ([P.sub.1] - [P.sub.2]) / 1 / 2[[rho].sub.i][V.sup.2.sub.i][|.sub.z=0]) as a function of the contraction ratio. [[alpha].sub.e1]=1, [[alpha].sub.2]=1.377 corresponds to the values used in the simulations in this work.

Comparison of the lines for [[alpha].sub.e1] = 1, [[alpha].sub.2] = 1.377 with those for [[alpha].sub.e1] = 1, [[alpha].sub.2] = 1, which corresponds to a flat velocity profile in the channels, shows that neglecting the velocity profile in the inlet channels results in a considerable under prediction of the entrance effect pressure loss. [[alpha].sub.e1]=2 corresponds to laminar flow in the circular conduit (or rather the canning) in front of the PF [50,52]. While laminar flow in the canning is unlikely, comparison of the lines for [[alpha].sub.e1] = 1, [[alpha].sub.2] = 1.377 and [[alpha].sub.e1] = 2, [[alpha].sub.2] = 1.377 shows the absolute maximum (but practically unlikely) range of values that [zeta] can have depending on the velocity profile in the canning in front of the PF. In the range of expansion ratios of relevance for PFs (0[less than or equal to][[alpha].sub.i]<0.4), there is not that much difference between these extreme values of [[zeta].sub.i], i.e. in practice [[zeta].sub.i], like [[zeta].sub.0], is not that sensitive to the velocity profile in the canning. Note that if [[alpha].sub.e1][not equal to][[alpha].sub.2], [[zeta].sub.i] is predicted to be non-zero when the contraction ratio is one, despite there being no change in the area for flow.

In an attempt to partially validate the equation for the entrance effect pressure drop, cold flow data for two PFs with different lengths, but otherwise the same, has been considered. This is the same data as used for model calibration (section 4). In section 4, the model is calibrated using an entrance effect given by Eq. (50). However, to test the entrance effect equation another parameter estimation was carried out in which [[zeta].sub.i] was optimised (in addition to [k.sub.W] and [[beta].sub.W]). This gave a value of 2.385 for [[zeta].sub.i], which compares favourably with the value of 2.159 predicted by Eq. (50) for the cold flow configuration ([[alpha].sub.i]=0). While encouraging this result falls short of fully validating Eq. (50); further validation is required.

3.5.3 Entrance and Exit Effect Prediction: Comparison with Literature

We have attempted to compare values of [[zeta].sub.i] and [[zeta].sub.0] predicted by the method developed in this work with the literature. However, the range of values quoted in the literature (Table 1) makes it difficult to draw any conclusions about the accuracy or "correctness" of the equations developed here (Eq. (45) and (50)). Even the CFD studies gave contradictory predictions; for example, Konstandopoulos et al. [46] found the exit losses to be approximately twice the entrance losses, while Bouteiller et al. [19] found contraction losses to be greater than exit losses, which is in line with this work (Fig. 9, 12).

3.6 Solution of Equations

The gas balance equations (Eq. (7), (9), (22), (23), (26), (27)), together with equations for the pressure drop across the wall (Eq. (40)), the soot-loaded channel width ([d.sub.i], Eq. (1)) and the gas densities (Eq. (5)), form a BVP. This can be solved with Matlab[R]'s bvp4c BVP solver, which is a finite difference code that implements the three-stage Lobatto Ilia formula [55,56]. The solution has a continuous first derivative and is uniformly fourth-order accurate over the whole spatial domain.

The solid energy balance (Eq. (41)) is a partial differential equation in time (t) and space (z). This was solved by discretising the PF in the axial direction to eliminate the spatial derivatives and hence produce a series of coupled Ordinary Differential Equations (ODEs) in t (method of lines). The discretisation was based on a uniform mesh. The resulting ODEs were solved using Matlab[R]'s ode15s [57,58], a variable time step, variable order (1-5), ODE solver designed for stiff problems.

The physical properties for the model have been taken from the literature (Table B.1). On the subject of physical properties, in case it is of importance to the reader, we note that the correlation given for the viscosity of air in Bissett's original paper [15] actually gives a very poor fit to measured data [59,60], e.g. viscosity is over predicted by approx. 30% at 20[degrees]C.

4. CALIBRATION OF BACKPRESSURE PREDICTION

To be able to predict backpressure across a PF two empirical parameters need to be determined, viz. the wall permeability and Forchheimer constants ([k.sub.W] and [[beta].sub.W] respectively). These were determined by optimising [k.sub.W] and [[beta].sub.W] to minimise the sum of the square of the differences between measured and simulated backpressure using a simplex algorithm (Matlab[R]'s fminsearch [61,62]). Cold flow data for two lengths of PF of the same type and diameter were used. These PFs were uncoated and clean (soot-free). The simulations were run with a realistic plug length (5 mm). Figure 13 compares measured and simulated backpressure both as backpressure versus flow rate (Fig. 13, top) and as a parity plot (Fig. 13, bottom); good agreement between measurement and simulation was achieved.

The values of [k.sub.W] and [[beta].sub.W] obtained are given in Table 2. [k.sub.W] lies within the range suggested by Konstandopoulos et al. [46] of around [10.sup.-13] to [10.sup.-12] [m.sup.2]. The range of values given in the literature for the [[beta].sub.W] are much wider than those given for [k.sub.W]. Thus, the value obtained for [[beta].sub.W] here is an order of magnitude greater than the range of [10.sup.5] to [10.sup.6] [m.sup.-1] suggested by Konstandopoulos et al. [46], but an order of magnitude smaller than the value of 5 x [10.sup.8] [m.sup.-1] given in an earlier work by the same author [16].

Since the model gives a good prediction for PFs with two lengths, this gives some confidence that the model is correctly predicting the contributions to total backpressure of the entrance effect (which is largely unaffected by PF length) and the along-PF terms. (Remember that the exit pressure change is zero for the cold flow configuration (section 3.5.1)). However, it must be noted that the backpressure is not very sensitive to length in the range of lengths investigated (Fig. 13). Although, that is at least partially due to the fact that the entrance effect makes up a large portion of the total backpressure (see section 5). This concept of using measurements on PFs with a range of lengths to separate entrance and exit effect from along-PF contributions to backpressure has been used by other workers [19,47].

Ideally, the model would be validated over a wider range of lengths, particularly shorter ones. However, unfortunately PFs are only commercially available in a limited number of lengths. Some workers have overcome this problem either by shortening PFs by cutting them and re-plugging [47] or even by using silica balls as "synthetic ash" to shorten the effective length of the channels [19].

There are three non-linear (in velocity) terms contributing to the backpressure, viz. the entrance and exit effects, the momentum convection term in the momentum balance equations and Forchheimer across-wall pressure drop. These, together with the effects of compressible flow, result in the plot of backpressure versus mass flow being curved (Fig. 13); without these terms/effects backpressure would be linear in mass flow. Of the three non-linear contributions, only the Forchheimer contribution has an adjustable parameter. Thus, if the momentum convection and entrance and exit effect terms in the model were too large, the predicted backpressure would increase more steeply with mass flow than the measured data even with zero [[beta].sub.W]; it does not (Fig. 13). If, on the other hand, the momentum convection and entrance and exit effect terms in the model were too small compared to reality, then this could be compensated for by increasing [[beta].sub.W]. The fact that the model matches the experimental data (Fig. 13), at least tells us that the momentum convection and entrance and exit effect terms in the model are not too large.

5. RELATIVE CONTRIBUTIONS OF DIFFERENT SOURCES OF PRESSURE DROP

In this section, the effect of flow rate, PF length, temperature, etc. on predicted backpressure will be investigated, as will the relative magnitude of the different contributions to backpressure, viz. across-wall losses (Darcy and Forchheimer), along-channel losses (viscous and inertial) and entrance and exit effects. All simulations were run using the values of [k.sub.W] and [[beta].sub.W] given in Table 2; note that these values were determined for a clean, uncoated PF.

The first thing to realise is that the magnitude of the different contributions to backpressure (excluding entrance and exit effects) depends on the path taken by the gas through the PF. This is illustrated in Fig. 14, which shows how the different contributions vary with the point at which the gas crosses the PF wall. This plot was calculated from the axial profiles of [P.sub.i] and P Thus, gas that crosses the wall at the front of the PF will have zero contribution to backpressure from the inlet channel; conversely, gas that crosses at the end of the PF will have zero contribution from the outlet channel. The curve for the across-wall pressure drop reflects the through-wall mass flux (Fig. 15); clearly, a higher through-wall mass flux means a higher pressure drop across the wall according to Eq. (40). Since the total pressure difference across the PF (for constant mass flow and temperature) is constant irrespective of the path taken by the gas, the curve for the total along-channel backpressure contribution (i.e. inlet + outlet channel) mirrors the curve for the across-wall pressure drop. Note that for the permeability and flow rate used, the along-channel pressure losses are greater than the across-wall loses, except at the very end of the PF.

It is worth explaining how the contributions to backpressure can vary with path taken by the gas, while the total backpressure remains constant. Figure 6 defines a series of locations across the PF at which it useful to define the pressure. Note that while [P.sub.2] and [P.sub.5] are fixed at the inlet channel entrance and the outlet channel exit (respectively), [P.sub.3] and [P.sub.4] can be at any point along the PF provided that they are both at the same axial coordinate. Using the pressures at these locations, the total pressure across the PF (excluding entrance and exit effects) can be expressed in terms of the contributions:

[mathematical expression not reproducible] (53)

As indicated, the bracketed terms are the contributions to backpressure due to the inlet channel, the across-wall pressure drop and the outlet channel, respectively. Since [P.sub.3] and [P.sub.4] cancel, it is clear that the total pressure drop is the same irrespective of the path taken by the gas through the PF.

With one exception, the simulations in this section have been run for the on-vehicle configuration (Fig. 1, bottom), rather than the cold flow configuration (Fig. 1, top) used for calibrating the model, as this makes the results more applicable to real applications. The on-vehicle configuration is somewhat idealised, which is possible as the predictions are not being compared to measured data. The PF sits in a canning of the same internal diameter as the external diameter of the PF with the exit of the canning at atmospheric pressure (Fig. 1, bottom). There are no cones to act as extra sources of pressure change. The canning extends far enough in front of the PF for there to be fully developed flow before the gas contracts to enter the PF and far enough after the PF for fully developed flow to be established after the gas expands on exiting the PF.

5.1 Steady-State Simulations

Since the contributions to backpressure depend on the path taken by the gas through the PF (see above), it is important to be clear how the different contributions to backpressure for the whole filter were calculated. Further, since compressible fluid behaviour is assumed, the different contributions to backpressure are not additive if the effects of gas compressibility are significant. For each study, three sets of simulations were run with different contributions removed, viz.:

1. Darcy's law with Forchheimer extension only

2. Darcy's law with Forchheimer extension and along-channel viscous losses

3. Full model

All simulations were run with entrance and exit effects included. Exit and entrance effect contributions were taken from the full model simulation (#3); clearly, entrance and exit effect pressure change contributions are easier to isolate than contributions to pressure losses occurring within the PF. Simulation #1 was run with [partial derivative]P / [partial derivative]z = [partial derivative]P / [partial derivative]z = 0 instead of Eq. (18) and (19). This result in a flat through-wall velocity profile along the PF. The combined contribution of Darcy's law and the Forchhheimer term can be obtained by subtracting the entrance and exit pressure changes from the total pressure drop from this simulation (#1). The flat through-wall velocity profile means that this result could be obtained directly (and with identical results) from Eq. (40), with [P.sub.0] calculated from the exit effect pressure change from simulation #3 and atmospheric pressure.

Inspection of Eq. (40) shows that the contributions of Darcy's law and the Forchheimer term to the across-wall pressure drop are not additive. This is due to the change in density with pressure across the wall being accounted for (compressible flow). However, the contributions of Darcy's law and the Forchheimer term to [P.sup.2.sub.i] are additive, which means the contribution of Darcy's law to the across-wall pressure drop, ([P.sub.i] - [P.sub.0])[|.sub.Darcy], can be calculated from:

[mathematical expression not reproducible] (54)

where [mathematical expression not reproducible] is calculated from Eq. (40) with [[beta].sub.W]=0 and [mathematical expression not reproducible] is calculated from Eq. (40) with both Darcy and Forchheimer terms in included and with P obtained from the exit effect pressure change from simulation #3. Note that Eq. (54) has used the difference of two squares, [P.sup.2.sub.i] - [P.sup.2.sub.0] = ([P.sub.i] - [P.sub.0])([P.sub.i] + [P.sub.0]). The contribution of the Forchheimer term to the across-wall pressure drop is similarly given by:

[mathematical expression not reproducible] (55)

Note that addition of Eq. (54) and (55) gives Eq. (40), as it should.

Simulation #2 was run with [alpha]=0 in Eq.. (18) and (19). The along-channel viscous contribution to backpressure was taken by subtracting the results of simulation #1 from #2 (each with the entrance and exit effects already subtracted). Note that including the along-channel losses results in the through-wall velocity profile no longer being flat, which in turn will affect the across-wall pressure drop (Darcy's law and Forchheimer term). The calculation method used means that any change in the across-wall pressure drop contribution associated with the change in through-wall velocity profile due to introduction of viscous losses along the channels is counted as being part of the contribution of the along-channel viscous losses. Having said that, if the contribution of the Forchheimer term is negligible and the pressure change across the PF wall is not large enough to cause significant change in the gas density, then Eq. (40) reduces to the across-wall pressure drop being directly proportional to [[phi].sub.W], in which case the average across-wall pressure drop along the PF will be the same irrespective of the through-wall mass flux distribution.

Finally, the contribution of the momentum convection term is taken as being the difference between simulation #3 and #2 (again with the entrance and exit contributions already subtracted from both). As with the along-channel viscous losses, the method of calculation means that the momentum convection term includes any changes in the other contributions caused by the change in the velocity field on including this contribution.

5.1.1 Effect of Including Plugs

Figure 16 compares predicted backpressure as a function of flow rate for simulations with 0, 5 and 10 mm long plugs, the latter being somewhat longer than expected in reality. Increasing the length of the plugs is predicted to increase the backpressure, but only slightly. Since the effect is so small, all the simulations in the rest of this paper, have been run with zero length plugs. (Note that while the plugs have zero length in the simulation, the channels are still plugged.)

Increasing the length of the plugs increases the backpressure since the plugs reduce the effective length of wall for through-wall gas flow thereby increasing the through-wall mass flux and hence the pressure drop across the wall according to Eq. (40). Reducing the wall permeability increases the effect of plug length on backpressure (not shown) as it increases the pressure drop across the wall.

The other thing to note in Fig. 16 is that the predicted backpressure is much lower than in the simulation/data used to calibrate the model (Fig. 13). This is because the model calibration was done with the cold flow configuration (Fig. F top) where the PF entrance is at atmospheric pressure, while Fig. 16 is for the on-vehicle configuration (Fig. F bottom) where the PF exit is at atmospheric pressure. This means that the pressure within the PF is lower for the cold flow case, which means that the gas densities will be lower so a higher gas velocity is required to achieve the same mass flow and hence the backpressure will be greater; all contributions to backpressure increase with gas velocity.

5.1.2 Effect of Flow Rate on Backpressure

The magnitude of the different contributions to backpressure and how they vary with flow rate are shown in Fig. 17. This is expressed both in terms of absolute backpressure (Fig. 17. top) and as the percentage contribution to backpressure (Fig. 17. bottom).

As often stated [17,21,22,36,47], but rarely demonstrated [19], the contribution of the Forchheimer term is tiny, amounting to only 2.5% of the total backpressure even at the very high mass flow of 0.3 kg [s.sup.-1]; under the conditions of this simulation, this term could safely be neglected. As might be expected, the importance of the contributions which are non-linear in velocity (momentum convection, entrance and exit effects (5)) increase with flow rate over the contributions which are linear in velocity (Darcy's law, along-channel viscous losses). Above about 0.14 kg [s.sup.-1] the entrance losses become the largest contribution to backpressure. However, the entrance loss is partially offset by the exit effect, which, as already discussed (section 3.5.1). actually results in a rise in pressure as gas leaves the PF.

The ratio of the entrance effect to the exit effect is fairly consistent across the full flow range, falling from 2.66 to 2.52 as the flow increases. This change is due to the increase in pressure at the PF entrance with increasing flow rate resulting in an increase in the gas density and hence a slightly lower velocity being required to achieve the same mass flow which in turn results in a slightly slower increase in entrance loss than in exit pressure change with increasing mass flow.

For the value of [k.sub.W] used in this simulation (Table 2). the viscous losses along the channel are a bigger contribution to total backpressure than the across-wall (Darcy's law) pressure drop. This will change if the permeability is reduced; in this case the Darcy's law contribution exceeds the viscous losses along the channels somewhere between a value of [k.sub.W] of 3.0 x [10.sup.-13] and 2.0 x [10.sup.-13] [m.sup.2].

The contribution of along-channel viscous losses to backpressure as a function of mass flow is a straight line through the origin. This can be rationalised as follows. Under isothermal conditions, for a symmetric (d=[d.sub.0]), soot-free ([d.sub.i]=d) PF, when the pressure change across the PF is not large enough to cause a significant change in the gas density (incompressible gas behaviour), [V.sub.i] + [V.sub.0] is a constant along the PF (= [V.sub.i][|.sub.z-0] = [V.sub.o][|.sub.z-L] By adding Eq. (18) and (19) to give an expression for [partial derivative]([P.sub.i] + [P.sub.o]) / [partial derivative]z, setting the momentum convection term to zero ([alpha]=0) and solving the resulting equation under these conditions, one obtains an expression for the pressure drop contribution due viscous losses along the channels, viz. F [mu][V.sub.i][|.sub.z=0]L / [d.sup.2]. Thus, under the stated conditions the backpressure contribution due to along-channel viscous losses is linear in flow rate and PF length.

The contribution of Darcy's law is also linear in mass flow. This not immediately expected from Eq. (40). However, if the change in gas density with pressure across the PF wall not significant, i.e. the gas is effectively behaving as if it were incompressible, then the Darcy's law contribution will indeed appear to be linear in mass flow.

The simulation in Fig. 17 was run for the same flow range as the cold flow data used for model calibration (0 to 0.3 kg [s.sup.-1], Fig. 13). However, this flow range is much larger than would likely be encountered by a PF on a vehicle over a drive cycle. Hence, Fig. 18 shows a zoomed in version of Fig. 17 covering the range of mass flows likely to be encountered in a real application. This figure goes up to 0.6 kg [s.sup.-1], which is just below the maximum mass flow reached in the WTLC test used for running transient simulations (section 5.2).

In this more practical mass flow range, the predicted contributions of both the momentum convection and Forchheimer terms are negligible, amounting to about 2 and 1% of the total backpressure at 0.06 kg [s.sup.-1] (Fig. 18). Entrance and exit effects are still significant in this flow range, but are smaller than the Darcy's law contribution, except at high mass flow (>0.055 kg [s.sup.-1]), and are much smaller than the viscous losses along the channels.

5.1.3 Effect of PF Length on Backpressure

The predicted effect of PF length on backpressure is shown in Fig. 19. These simulations were run for 20[degrees]C and a mass flow of 0.06 kg [s.sup.-1] (which, as already stated, is just below the maximum mass flow reached in the WTLC test in section 5.2).

As the length of the PF is reduced (Fig. 19), the through-wall gas velocity increases, resulting in an increase in the across-wall pressure drop. The Darcy's law contribution rises rapidly once the length is reduced below 50 - 60 mm. The Forchheimer contribution starts to rise sharply at a shorter length than the Darcy's law contribution at around 30 - 40 mm, but increases more rapidly with decreasing length since this term has a stronger wall mass flux dependency (Eq. 40).

The contribution of viscous losses along the channels increases linearly with length, as expected from the analysis in the previous section (section 5.1.2).

The entrance and exit effect pressure changes are, as might be expected, little effected by the length. The exit effect shows no change with length, while the entrance effect shows a very slight decrease at the very shortest lengths. This fall in the entrance effect is due to the increase in gas density at the PF entrance. The entrance effect pressure change is proportional to [[rho].sub.i][V.sup.2.sub.i][|.sub.z=0]

(Eq. 50). [rho][V.sup.2] can be written as [rho]VxV; [rho]V is independent of gas density, while V is inversely proportional to density at fixed mass flow. Thus, the increase in pressure at the PF entrance results in an increase in density and hence a (slight) decrease in the entrance effect pressure drop at very short PF lengths. In contrast, the pressure at the PF exit remains at close to 1 atm irrespective of PF length, so the exit pressure change remains constant.

At this mass flow (0.06 kg [s.sup.-1]), the momentum convection term has negligible contribution to backpressure, amounting to no more than 2.3%.

As the PF length is varied, the total backpressure exhibits a shallow minimum at about 100 mm. Reducing the length much below this results in a rapid rise in backpressure due to the dramatic increase in across-wall pressure drop. Increasing the length above the 100 mm results in a slow increase in backpressure, as the rise in the along-channel viscous losses with increasing length is offset by the falling contribution of the Darcy's law across-wall pressure drop.

5.1.4 Effect of Temperature on Backpressure

The effect of temperature is explored in Fig. 20. These simulations were run for a mass flow of 0.06 kg [s.sup.-1]. All contributions to predicted backpressure increase in magnitude with temperature, although the Forchheimer and momentum convection terms remain negligible (no more than 1.0 and 2.2% of total backpressure). Note that the while the exit effect pressure change increases in magnitude with temperature, it actually becomes more negative in value.

As the temperature increases, the pressure drop due to the linear backpressure terms (along-channel viscous losses, Darcy's law) increase more rapidly than the non-linear terms (entrance and exit effects, momentum convection, Forchheimer term) (Fig. 20, top), resulting in the relative contribution to backpressure of the viscous losses along the channels and Darcy's law rising with temperature (Fig. 20, bottom), while the relative contributions of entrance and exit effects fall. This is because the linear pressure drop terms are proportional to [micro]V, while the non-linear ones are proportional to [rho][V.sup.2]. Since pV (equal to the mass flux) is independent of temperature and V increases linearly with temperature (for a fixed mass flow and provided the change in pressure is not so large that it causes a significant change in gas density), both V and [rho][V.sup.2] increase linearly with temperature. However, the gas viscosity also increases with temperature (Table B.1), so the linear (viscous) contributions increase faster with temperature than the non-linear (inertial) terms. The temperature dependency of [rho][V.sup.2] also explains why the absolute values of the entrance and exit effects increase linearly with temperature (Fig. 20, top).

Finally, the ratio of entrance effect to the exit effect decreases very slightly with increasing temperature from 2.62 to 2.49 as the temperature is increased from 20 to 875[degrees]C. This is due to the gas density at the PF entrance increasing slightly faster with temperature than the gas density at the exit due to the increase in pressure at the PF entrance; the pressure at the PF exit remains close to 1 atm.

5.1.5 Effect of Gas Compressibility on Prediction

In Bissett's original model [15], gas densities were calculated at 1 atm, assuming that the change in pressure across the PF was too small to have significant impact on the density. This section considers the significance of this assumption. Figure 21 compares predicted backpressure as a function of flow rate at 20[degrees]C for the on-vehicle configuration for gas densities calculated at 1 atm and the actual pressure (compressible gas behaviour). Since the temperature is constant, the "densities at 1 atm simulation" corresponds to incompressible gas behaviour. For a compressible gas, as the pressure increases the gas density also increases, so the gas velocity for a given mass flow is lower and hence the pressure loss will also be lower. So, in principal, at high flow rates where the backpressure is greater, the predicted backpressure would be expected to be lower if the gas density were calculated at the actual pressure (compressible flow) than if it were calculated at 1 atm. However, while Fig. 21 does show that the predicted backpressure is lower at high flow rates for compressible flow, the effect is extremely slight and could be safely neglected. Note that this lack of dependency of the prediction on the pressure used for the density calculation is support for the validity of the assumption of incompressible flow used when deriving expressions for the entrance and exit effect pressure changes (section 3.5).

On the other hand, the effect of the density calculation is much more pronounced for the cold-flow configuration (Fig. 1, top; Fig. 22) than for the on-vehicle configuration (Fig. 1, bottom; Fig. 21). This is because with the cold-flow configuration the pressure along the PF is decreasing from atmospheric pressure at the inlet, rather than decreasing to atmospheric pressure at the exit, so the change in pressure relative to 1 atm and hence the change in gas density is greater. Note that the predicted backpressure for the cold flow configuration is higher for compressible flow than that for densities calculated at 1 atm (Fig. 22), while the converse is true for the on-vehicle configuration (Fig. 21). This is because the pressure in the PF is below 1 atm in the cold flow configuration, but above 1 atm for the on-vehicle situation.

5.2 Transient Simulations

While steady-state simulations make it easy to see how backpressure varies with flow rate, PF length, temperature, etc., transient simulations are perhaps of more practical interest. Since it is difficult to see the effect of different contributions over a transient cycle, a different approach will be taken; these simulations will look at finding what the simplest possible model able to predict backpressure and outlet temperature is.

These simulations have been run with gasoline input data rather than diesel as this enables a wider temperature range to be covered. Figure. 23 shows the mass flow and inlet temperature over the WLTC used as model input.

5.2.1 Effect of Density Calculation

The effect of the density calculation on the prediction over the WLTC is shown in Fig. 24. This compares predictions for compressible flow (densities calculated at actual temperature and pressure, according to Eq. 5) and for densities calculated at 1 atm and the actual temperature. In the second-by-second plot (Fig. 24, top), there little obvious difference between the two scenarios. The differences are clearer in the parity plot (Fig. 24, bottom). This shows a clear (albeit slight) shift in the points below the diagonal at high backpressures for the compressible flow simulation, as expected (see section 5.1.5). However, the effect is small. (Note that the parity plot shows some scatter about the diagonal. This seems to be evenly distributed about the diagonal and so does not indicate a trend but rather numeric error between the simulations. This plot has approx. 1800 points (one for each second of the WLTC), so it is easy for a relatively few slightly off-diagonal points to cover the majority that lie on or very close to the diagonal.)

5.2.2 Effect of Different Contributions to Backpressure

In this sub-section we have tried to determine the simplest model required for accurate prediction of a transient simulation. Since measured backpressure data was not available, a series of simulations were run with various terms removed to determine which terms had a significant contribution to the prediction and which could be safely neglected.

Figure 25 compares predictions for four versions of the model, viz. i) the full (or base) model; ii) a model with only Darcy's law; iii) a model with Darcy's law and viscous losses along the channels only; and iv) a model with Darcy's law, viscous losses along the channels and entrance and exit effects. From the second-by-second plot (Fig. 25, top), it is clear that a model with just Darcy's law will grossly under predict the backpressure (assuming the same permeability is used). The results for other simulations are quite similar, so it is necessary to look at the parity plot to distinguish between the simulations (Fig. 25, bottom). The points for the model with Darcy's law and viscous losses along the channels only lie above the diagonal at higher backpressures, indicating that this model is predicting a lower backpressure than the base model, while the model with Darcy's law, viscous losses along the channels and entrance and exit effects lies on the diagonal, indicating that this model is giving an equivalent prediction to the base model. From this it can be concluded that for this simulation, the Forchheimer and momentum convection terms can safely be neglected from the model without noticeably impacting the prediction, but the entrance and exit effects cannot be neglected.

The predicted outlet temperatures for the four models are compared in Fig. 26; there are no noticeable differences in the predictions.

6. EFFECT OF THE LATEST PF MODEL DEVELOPMENTS BY BISSETT, KOSTOGLOU AND KONSTANDOPOULOS

Probably the most significant advance in PF modelling in recent years has been made by Bissett, Kostoglou and Konstandopoulos [27,28]. Previously, PF models have assumed that the gas flow in the channels is not significantly influenced by the flow through the porous walls so that values of F and Nu for flow in square channels with non-porous walls can be used. These workers have numerically solved the fully developed laminar flow and heat transfer within long square channels with porous walls (subject to the assumptions of constant wall velocity along the PF, incompressible flow and constant viscosity) to investigate the influence of the wall flow. They suggest the following changes to the PF model:

* Recognition that <[u.sup.2]> [not equal to] [<u>.sup.2], so that a correction factor, [alpha], needs to be applied to the momentum convection term. [alpha] varies with the wall Reynolds number, [Re.sub.W].

* That the friction factor for flow in square channels (and hence F) is no longer a constant, but varies with [Re.sub.W]

* That the Nusselt number is also no longer a constant but a function of [Re.sub.W]

Thus the momentum balance equations (Eq. (18) and (19)) become:

[partial derivative][P.sub.i] / [partial derivative]z = - [[alpha].sub.i]([Re.sub.W])[[rho].sub.i][V.sub.i] / [d.sup.2.sub.i] [partial derivative][d.sup.2.sub.i][V.sub.i] / [partial derivative]z - [F.sub.i](R[e.sub.W])[[mu].sub.i][V.sub.i] / [d.sup.2.sub.i] (56)

[partial derivative][P.sub.0] / [partial derivative]z = -[[alpha].sub.o]([Re.sub.W)[[rho].sub.o][V.sub.o][[partial derivative][V.sub.o] / [partial derivative]z - [F.sub.o]([Re.sub.W])[[mu].sub.o][V.sub.o] / [[d.sup.2.sub.o] (57)

Thus, [alpha] and F are now functions of [Re.sub.W] and are different for the inlet and outlet channels. Here the new form of the momentum balance equations has been modified, but obviously Bissett et al. [27,28] started from the old form (Eq. (10) and (11)). However, the modifications apply equally to both forms; the viscous loss term is the same for both old and new momentum balance equations so [F.sub.i] and [F.sub.o] will be the same for both and [[alpha].sub.i] and [[alpha].sub.0] are calculated directly from the three-dimensional velocity profile (rather than being related to the form of the momentum convection term). [Re.sub.W] is defined as [d[phi].sub.W] / [[mu].sub.i] or d[[phi].sup.W] / [[mu].sub.o], depending on the channel. Note that since [[phi].sub.W] is referenced to the width of the soot-free inlet channel (d), the dimension in [Re.sub.W] is the same irrespective of whether the inlet or outlet channel is being considered or whether the PF is soot loaded or not.

For ease of implementation of these changes into the one-dimensional model framework, these workers have published correlations for [alpha], F and Nu as a function of [Re.sub.W]. It is worth noting that the correlations published in these two papers by the same authors (albeit in a different order) are different, e.g. in some cases the order of the polynomials used is different. Therefore, it is useful to compare the correlations graphically.

The correlations for F as a function of [Re.sub.W] from the two papers are very similar (Fig. 27), although there are small differences. Nu depends on Pr, as well [Re.sub.W]. Since Pr for air is 0.70[+ or -]0.01 between 0 and 725[degrees]C [63], Nu was calculated with Pr=0.7. Values of Nu thus calculated from the correlations in both papers appear identical for Pr=0.7 at least, despite the equations looking very different (Fig. 28). The correlations for [alpha] show greater differences (Fig. 29), although the general trends are very similar.

To see if what sort of effect these changes to the model have on the predictions, some example simulations have been run. These simulations were run using the correlations published in [27] rather than [28]. In the first simulation, the mass flow was increased at constant temperature (Fig. 30). Since the correlations are only valid up to [Re.sub.W]=3, the simulation was only run until the largest value of [Re.sub.W] along the PF exceeded 3. (Note that [Re.sub.W] varies along the length of the PF with the through-wall mass flux (Fig. 15)). Thus the maximum mass flow in this simulation is much smaller than the comparable simulation run in the previous section (Fig. 17). While this limit on the maximum flow rate for the simulation may seem restrictive, it is worth noting that the same work [27,28,64] questions whether one-dimensional models are really valid for approximating the three-dimensional flow in a PF when [Re.sub.W]>3, given that there appears to be no solution for fully developed flow in the inlet channels from the three-dimensional model when [Re.sub.W]>3 [27,28].

The simulations show only the slightest difference in predicted backpressure between the two models at the highest mass flow (Fig. 30) with the [Re.sub.W]-dependant model predicting fractionally lower backpressure. While these simulations were run using the new form of the momentum balance, the same conclusion was reached when the same modifications were made to the old momentum balance equations (not shown).

Moving on to transient simulations, Fig. 31 and 32 compare predicted backpressure and outlet temperature for the two models. In these simulations [Re.sub.W] remains well within the range of validity of the correlations ([Re.sub.W]<3); apart from a single excursion in maximum [Re.sub.W] to 2.4 during the initial PF warm-up, the highest maximum [Re.sub.W] encountered over the cycle was less than 1.6. Note that [Re.sub.W] varies along the length of the PF with the through-wall velocity (Fig. 15) and temperature profile, so the values mentioned here refer to the maximum value of [Re.sub.W] along the PF. When predicted backpressure as a function of time during WLTC is compared (Fig. 31, top), no significant difference between the two models can be seen. A parity plot comparison of the predictions of the two models (Fig. 31, bottom), which makes any differences between the two models more obvious, shows only very small differences in the prediction.

Turning to the outlet temperature prediction of the two models (Fig. 32), the second-by-second plot shows no difference between the two models (Fig. 32, top), while the parity plot (Fig. 32, bottom) shows absolutely no deviation from the diagonal.

Thus for the simulations run here, it is concluded that the improved model in which [alpha], F and Nu are functions of [Re.sub.W] offers no advantage over a model where [alpha], F and Nu are constants and that the extra model complexity is unwarranted.

In a study on GPF warm-up, Opitz et al. [24] also report that [Re.sub.W]-dependant heat (and mass) transfer coefficients are unnecessary. This conclusion was based on the very good agreement found between the predictions of their one-dimensional model and those of a two-dimensional model, which does not rely on heat and mass transfer coefficients.

On the other hand, in their example simulations, Bissett et al. [64,65] report some differences in prediction between "standard" and "[Re.sub.W]-dependant" models, but these were largely in predicted internal properties, such as axial distributions, rather than in the outlet predictions, where they also found only small differences. One could argue that it is only improvements in the predictions of quantities that can be validated against measured data that is important. However, there may be other scenarios where these improvements in the model could be important for the accuracy of the predictions. Nevertheless, while the effects of [Re.sub.W]-dependant [alpha], F and Nu may appear not be important, it is important that possibility that they could be has been investigated.

7. CONCLUSIONS

A one-dimensional model for a PF has been developed based on the approach of Bissett [15]. The model makes two innovations.

The first innovation is that the term for momentum convection in the gas momentum balance equations includes the loss/gain of axial momentum in the direction perpendicular to the channels; neglecting this results in the momentum convection term being too large. Under isothermal conditions the new momentum convection term is half the size of the old term.

The second innovation is that the equations for the pressure change due to the abrupt contraction at the PF entrance and for abrupt expansion at the exit are derived which take into account the fact that the velocity profile across the channels is not flat; other workers have used equations appropriate for high Reynolds numbers which assume flat velocity profiles. Neglecting the velocity profile across the channels results in underestimation of the entrance and exit effect pressure changes.

The model has been calibrated (i.e. values for the wall permeability and Forchheimer constant determined) using cold flow data for PFs of two lengths. The use of more than one length gives some confidence that the model is correctly differentiating between pressure losses occurring along the PF and those occurring at the entrance and exit, although ideally the model would be validated for a larger range of lengths.

A simulation study was carried out to evaluate the importance of the different contributions to backpressure and how these change with flow rate, temperature and length. The main findings were:

* The contribution of the Forchheimer term is negligible, except for unrealistically small PFs. This term can safely be neglected

* The contribution of the momentum convection term is relatively small and this term can probably be neglected for the conditions encountered on a vehicle drive cycle. However, it remains significant for the type of cold flow test used in this work.

* Entrance and exit effect pressure changes are significant and should not be neglected.

* Gas expansion at the exit of the PF results in a pressure rise.

* For the PF used in this study, viscous losses along the channels were a bigger contribution to backpressure than the pressure drop across the PF wall. Of course, this depends on the wall permeability of the PF.

* Viscous contributions to backpressure (Darcy's law and along-channel viscous losses) increase more rapidly with temperature than inertial contributions (momentum convection, entrance and exit effects and Forchheimer term).

* For vehicle drive cycle simulations, whether gas densities were calculated at the actual pressure or at 1 atm had no significant effect on the backpressure prediction.

Finally, including a [Re.sub.W] dependency on the momentum flux correction factor ([alpha]), viscous loss coefficient (F) and Nusselt number (Nu), as described in the literature [27,28], was found to have no noticeable effect on the predictions.

REFERENCES

[1.] DieselNet, "Diesel Particulate Matter," https://www.dieselnet.com/tech/dpm.php, accessed Jan. 2017.

[2.] Niessner, R., "The Many Faces of Soot: Characterization of Soot Nanoparticles Produced by Engines," Angew. Chem. Int. Ed. 53(46):12366-12379, 2014, doi:10.1002/anie.201402812.

[3.] Resitoglu, I.A., Altinisik, K., and Keskin, Ali, "The pollutant emissions from diesel-engine vehicles and exhaust aftertreatment systems," Clean Techn. Environ. Policy 17:15-27, 2015, doi:10.1007/s10098-014-0793-9.

[4.] Guan, B., Zhan, R., Lin, H., and Huang Z., "Review of the state-of-the-art of exhaust particulate filter technology in internal combustion engines," J. Environ. Man. 154:225-258, 2015, doi: 10.1016 j.jenvman.2015.02.027.

[5.] DieselNet, "Health Effects of Diesel Particulates," https://www.dieselnet.com/tech/health_pm.php, accessed Jan. 2017.

[6.] Tollefson, J., "Soot a major contributor to climate change," Nature News: 15th Jan. 2013, doi: 10.1038/nature.2013.12225.

[7.] Kittelson, D.B., "Engines and Nanoparticles: A Review," J. Aerosol Sci., 29(5/6):575-588, 1998, doi:10.1016/S0021-8502(97)10037-4.

[8.] DieselNet, "Diesel Particulate Filters," https://www. dieselnet.com/tech/dpf.php, accessed Jan. 2017.

[9.] Hawker, P., Huthwohl, G., Henn, J., Koch, W. et al., "Effect of a Continuously Regenerating Diesel Particulate Filter on Non-Regulated Emissions and Particle Size Distribution," SAE Technical Paper 980189, 1998, doi:10.4271/980189.

[10.] Johnson, T., "Diesel Engine Emissions and Their Control," Platinum Metals Rev. 52(1):23-37, 2008, doi:10.1595/147106708X248750.

[11.] Morgan, C., "Platinum Group Metal and Washcoat Chemistry Effects on Coated Gasoline Particulate Filter Design," Johnson Matthey Technol. Rev. 59(3):188-192, 2015, doi: 10.1595/205651315X688109^

[12.] Johnson, T., "Vehicular Emissions in Review," SAE Int. J. Engines 9(2):1258-1275, 2016, doi:10.4271/2016-01-0919.

[13.] Chatterjee, D., Schmei[ss]er, V., Frey, M., and Weibel, M., "Perspectives of the Automotive Industry on the Modeling of Exhaust Gas Aftertreatment Catalysts," in "Modeling and Simulation of Heterogeneous Catalytic Reactions: From the Molecular Process to the Technical System, First Edition," (Weinheim, Wiley-VCH Verlag, 2012), Chap. 10, 303-343, doi:10.1002/9783527639878.ch10.

[14.] Ahmadinejad, M., Etheridge, J.E., Watling, T.C., Johansson, A, and John, G., "Computer Simulation of Automotive Emission Control Systems," Johnson Matthey Technol. Rev. 59(2):152-165, 2015, doi:10.1595/205651315X687876.

[15.] Bissett, E.J., "Mathematical Model of the Thermal Regeneration of a Wall-Flow Monolith Diesel Particulate Filter," Chem. Eng. Sci. 39(7-8): 1233-1244, 1984, doi:10.1016/0009-2509(84)85084-8.

[16.] Konstandopoulos, A., Skaperdas, E., Warren, J., and Allansson, R., "Optimized Filter Design and Selection Criteria for Continuously Regenerating Diesel Particulate Traps," SAE Technical Paper 1999-01-0468, 1999, doi:10.4271/1999-01-0468.

[17.] Haralompous, O.A., Kandylas, I.P., Koltsakis, G.C., and Samaras, Z.C., "Diesel particulate filter pressure drop. Part 1: modelling and experimental validation," Int. J. Engine Res. 5(2):149-162, 2004, doi:10.1243/146808704773564550.

[18.] Mohammed, H., Triana, A., Yang, S., and Johnson, J., "An Advanced 1D 2-Layer Catalyzed Diesel Particulate Filter Model to Simulate: Filtration by the Wall and Particulate Cake, Oxidation in the Wall and Particulate Cake by NO2 and O2, and Regeneration by Heat Addition," SAE Technical Paper 2006-01-0467, 2006, doi:10.4271/2006-01-0467.

[19.] Bouteiller, B., Bardon, S., Briot, A., Girot, P. et al., "One Dimensional Backpressure Model for Asymmetrical Cells DPF," SAE Technical Paper 2007-01-0045, 2007, doi:10.4271/2007-01-0045.

[20.] Wurzenberger, J. and Kutschi, S., "Advanced Simulation Technologies for Diesel Particulate Filters, A Fundamental Study on Asymmetric Channel Geometries," SAE Technical Paper 2007-01-1137, 2007, doi:10.4271/2007-01-1137.

[21.] Depcik, C., and Assanis, D., "Simulating Area Conservation and the Gas-Wall Interface for One-Dimensional Based Diesel Particulate Filter Models," ASME J. Eng. Gas Turbines Power 130(6):062807-(1-18), 2008, doi:10.1115/1.2939002.

[22.] Shejbal, M., Marek, M., Kubicek, M., and Koci, P., "Modelling of diesel filters for particulates removal," Chem. Eng. J. 154(1-3):219-230, 2009, doi:10.1016/j.cej.2009.04.056.

[23.] Watling, T.C., Ravenscroft, M.R., and Avery, G., "Development, validation and application of a model for an SCR catalyst coated diesel particulate filter," Catal. Today 188(1):32-41, 2012, doi:10.1016/j. cattod.2012.02.007.

[24.] Opitz, B., Drochner, A., Vogel, H., and Votsmeier, M., "An experimental and simulation study on the cold start behaviour of particulate filters with wall integrated three way catalyst," Appl. Catal. B: Environ. 144:203-215, 2014, doi:10.1016/j.apcatb.2013.06.035.

[25.] Tronconi, E., Nova, I., Marchitti, F., Koltsakis, G., Karamitros, D., Maletic, B., Markert, M., Chatterjee, D., and Hehle, M., "Interactions of N[O.sub.X] Reduction and Soot Oxidation in a DPF with a Cu-Zeolite SCR Coating," Emiss. Control Sci. Technol. 1(2):134-151, 2015, doi:10.1007/s40825-015-0014-y.

[26.] Koltsakis, G., Haralampous, O., Depcik, C., and Ragone, J.C., "Catalyzed diesel particulate filter modeling," Rev. Chem. Eng. 29(1):1-61, 2013, doi:10.1515/revce-2012-0008.

[27.] Bissett, E.J., Kostoglou, M., and Konstandopoulos, A.G., "Frictional and heat transfer characteristics of flow in square porous tubes of wall-flow monoliths," Chem. Eng. Sci. 84:255-265, 2012, doi:10.1016/j.ces.2012.08.012.

[28.] Kostoglou, M., Bissett, E.J., and Konstandopoulos, A.G., "Improved Transfer Coefficients for Wall-Flow Monolithic Catalytic Reactors: Energy and Momentum Transport," Ind. Eng. Chem. Res. 51(40):13062-13072, 2012, doi:10.1021/ie3011098.

[29.] SuperFlow[R], "SF-1020 ProBench", http://www.superflow.com/aspx/prodDetail.aspx?prodid=17&catid=4&navid=10, accessed Jan. 2017.

[30.] DieselNet, "Emission Test Cycles: Worldwide Harmonized Light Vehicles Test Procedure (WLTP)," https://www.dieselnet.com/standards/cycles/wltp.php, accessed Jan. 2017.

[31.] Tutuianu, M., Bonnel, P., Ciuffo, B., Haniu, T., Ichikawa, N., Marotta, A., Pavlovic, J., and Steven, H., "Development of the World-wide harmonized Light duty Test Cycle (WLTC) and a possible pathway for its introduction in the European legislation," Trans. Res. D: Trans. Environ. 40:61-75, 2015, doi:10.1016/j.trd.2015.07.011.

[32.] Bird, R.B., Stewart, W.E., and Lightfoot, E.N., "Transport Phenomena, 2nd Edition," (New York, John Wiley & Sons, 2002), 75-87, ISBN 0-471-41077-2.

[33.] Bird, R.B., Stewart, W.E., and Lightfoot, E.N., "Transport Phenomena, 2nd Edition," (New York, John Wiley & Sons, 2002), 822, ISBN 0-471-41077-2.

[34.] Bird, R.B., Stewart, W.E., and Lightfoot, E.N., "Transport Phenomena, 2nd Edition," (New York, John Wiley & Sons, 2002), 286, ISBN 0-471-41077-2.

[35.] Konstandopoulos, A., Kostoglou, M., Vlachos, N., and Kladopoulou, E., "Progress in Diesel Particulate Filter Simulation," SAE Technical Paper 2005-01-0946, 2005, doi:10.4271/2005-01-0946.

[36.] Konstandopoulos, A., "Flow Resistance Descriptors for Diesel Particulate Filters: Definitions, Measurements and Testing," SAE Technical Paper 2003-01-0846, 2003, doi:10.4271/2003-01-0846.

[37.] Darcy, H., "Les Fontaines Publiques de la Ville de Dijon", (Paris, Victor Dalmont, 1856), 570, https://books.google.co.uk/books?id=DOwbgyt_MzQC&pg=PA1&redir_esc=y#v=onepage&q&f=false, accessed Jan. 2017.

[38.] Forchheimer, P., "Wasserbewegung durch Boden," Zeitschrift des Vereines Deutscher Ingenieure 45:1781-1788, 1901.

[39.] Whitaker, S., "Advances in Theory of Fluid Motion in Porous Media," Ind. Eng. Chem. 61(12):14-28, 1969, doi:10.1021/ie50720a004.

[40.] Whitaker, S., "Flow in Porous Media I: A Theoretical Derivation of Darcy's Law," Trans. Porous Media 1(1):3-25, 1986, doi:10.1007/BF01036523.

[41.] Ahmed, N., and Sunada, D.K., "Nonlinear Flow in Porous Media," J. Hydr. Div. Proc. ASCE 95(HY 6):1847-1856, 1969.

[42.] Whitaker, S., "The Forchheimer Equation: A Theoretical Development," Trans. Porous Media 25(1):27-61, 1996, doi:10.1007/BF00141261.

[43.] Bissett, E.J., and Shadman, F., "Thermal Regeneration of Diesel-Particulate Monolithic Filters," AIChE J. 31(5):753-758, 1985, doi:10.1002/aic.690310508.

[44.] Bissett, E.J., "Thermal Regeneration of Particle Filters with Large Conduction," Math. Model. 6(1):1-18, 1985, doi:10.1016/0270-0255(85)90018-1.

[45.] Haralampous, O., and Koltsakis, G.C, "Intra-layer temperature gradients during regeneration of diesel particulate filters," Chem. Eng. Sci. 57(13):2345-2355, 2002, doi:10.1016/S0009-2509(02)00135-5.

[46.] Konstandopoulos, A., Skaperdas, E., and Masoudi, M., "Inertial Contributions to the Pressure Drop of Diesel Particulate Filters," SAE Technical Paper 2001-01-0909, 2001, doi:10.4271/2001-01-0909.

[47.] Masoudi, M., Heibel, A., and Then, P., "Predicting Pressure Drop of Wall-Flow Diesel Particulate Filters - Theory and Experiment," SAE Technical Paper 2000-01-0184, 2000, doi:10.4271/2000-01-0184.

[48.] Konstandopoulos, A., Kostoglou, M., Skaperdas, E., Papaioannou, E. et al., "Fundamental Studies of Diesel Particulate Filters: Transient Loading, Regeneration and Aging," SAE Technical Paper 2000-01-1016, 2000, doi:10.4271/2000-01-1016.

[49.] Bird, R.B., Stewart, W.E., and Lightfoot, E.N., "Transport Phenomena, 2nd Edition," (New York, John Wiley & Sons, 2002), 197-216, ISBN 0-471-41077-2.

[50.] Kays, W.M., "Loss Coefficients for Abrupt Changes in Flow Cross Section With Low Reynolds Number Flow in Single and Multiple Tube Systems," Trans. ASME 72:1067-1074, 1950.

[51.] Benedict, R.P., Wyler, J.S., Dudek, J.A., and Gleed, A.R., "Generalized Flow Across an Abrupt Enlargement," Trans. ASME J. Eng. Power 98(3):327-334, 1976, doi:10.1115/1.3446171.

[52.] Abdelall, F.F., Hahn, G., Ghiaasiaan, S.M., Abdel-Khalik, S.I., Jeter, S.S., Yoda, M., and Sadowski, D.L., "Pressure drop caused by abrupt flow area changes in small channels," Expt. Therm. Fluid Sci. 29(4):425-434, 2005, doi:10.1016/j.expthermflusci.2004.05.001.

[53.] Rouse, H., "Elementary Mechanics of Fluids," (New York, John Wiley & Sons, 1946), 57, ISBN 047174316X.

[54.] "VDI Heat Atlas, Second Edition," (Berlin, Springer-Verlag, 2010), 1065-1067, doi:10.1007/978-3-540-77877-6.

[55.] Mathworks[R], "bvp4c," http://uk.mathworks.com/help/matlab/ref/bvp4c.html, accessed Jan. 2017.

[56.] Kierzenka, J., and Shampine, L.F., "A BVP Solver Based on Residual Control and the Matlab PSE," ACM TOMS 27(3):299-316, 2001, doi:10.1145/502800.502801.

[57.] Mathworks[R]'s, "ode15s," http://uk.mathworks.com/help/matlab/ref/ode15s.html, accessed Jan. 2017.

[58.] Shampine, L.F., and Reichelt M.W., M.W., "The Matlab ODE Suite," SIAM J. Sci. Comput. 18(1):1-22, 1997, doi:10.1137/S1064827594276424.

[59.] Kadoya, K., Matsunaga, N., and Nagashima, A., "Viscosity and Thermal Conductivity of Dry Air in the Gaseous Phase," J. Phys. Chem. Ref. Data 14(4):947-970, 1985, doi:10.1063/1.555744.

[60.] Perry, R.H., and Green, D.W., "Perry's Chemical Engineers' Handbook, 7th Edition," (New York, McGraw-Hill, 1998), 2-208, ISBN 0-07115982-7.

[61.] Mathworks[R], "Optimizing Nonlinear Functions: fminsearch Algorithm", http://uk.mathworks.com/help/matlab/math/optimizing-onlinear-functions.html, accessed Jan. 2017.

[62.] Lagarias, J.C., Reeds, J.A., Wright, M.H., and Wright, P.E., "Convergence Properties of the Nelder--Mead Simplex Method in Low Dimensions", SIAM J. Optim. 9(1):112-147, 1998, doi:10.1137/S1052623496303470.

[63.] Perry, R.H., and Green, D.W., "Perry's Chemical Engineers' Handbook, 7th Edition," (New York, McGraw-Hill, 1998), 2-325, ISBN 0-07115982-7.

[64.] Bissett, E.J., Kostoglou, M., and Konstandopoulos, "Better Friction Coefficients and Nusselt Numbers for DPFs and When Not to Use Them at All," Cross-Cut Lean Exhaust Emissions Reduction Simulations (CLEERS), Dearborn, April-May 2012, http://www.cleers.org/workshops/workshop2012/presentations/Bissett_CLEERS2012.pdf, accessed Jan. 2017.

[65.] Bissett, E.J., and Wang, W., "On the Implications of Wall Reynolds Number Dependant Nusselt Number and Friction Factor on the Accuracy of Wall-Flow DPF Modeling," 3rd International Symposium on Modeling of Exhaust-Gas After-Treatment MODEGAT III, Bad Herrenalb, Germany, 8th-10th Sept. 2013.

[66.] Bird, R.B., Stewart, W.E., and Lightfoot, E.N., "Transport Phenomena, 2nd Edition," (New York, John Wiley & Sons, 2002), 178-184, ISBN 0-471-41077-2.

[67.] Lemmon, E.W., Jacobsen, R.T., Penoncello, S.G., and Friend, D.G., "Thermodynamic Properties of Air and Mixtures of Nitrogen, Argon, and Oxygen From 60 to 2000 K at Pressures to 2000 MPa," J. Phys. Chem. Ref. Data. 23(3):331-385, 2000, doi:10.1063/1.1285884. '

[68.] Shah, R.K., and London, A.L., "Laminar Flow Forced Convection Heat Transfer and Flow Friction in Straight and Curved Ducts - A Summary of Analytical Solutions," Tech. Rep. No. 75, Dept. Mech. Eng., Stanford Univ., 1971, http://www.dtic.mil/dtic/tr/fulltext/u2/736260.pdf, pp. 117.

[69.] Cybulski, A., and Moulijn, J.A. (Eds.) "Structured Catalysts and Reactors, 2nd Edition," (Boca Raton, CRC Press, 2005), 52, doi:10.1201/9781420028003.

CONTACT INFORMATION

Dr Timothy Watling

Senior Principal Scientist

Johnson Matthey Technology Centre

Blount's Court

Sonning Common

Reading RG4 9NH UK

tim.watling@matthey.com

Tel: +44 (0)118 924 2139

ACKNOWLEDGMENTS

The authors wish to thank Johnson Matthey PLC for permission to publish this paper. Neil Edwards, Johnson Matthey Technology Centre, is thanked for machine translating the Forchheimer reference [38].

ABBREVIATIONS / SYMBOLS

ABBREVIATIONS

atm - Atmosphere, 101325 Pa

BVP - Boundary value problem

CFD - Computational fluid dynamics

cpsi - Cells per square inch

GPF - Gasoline particulate filter

DPF - Diesel particulate filter

ODE - Ordinary differential equation

PF - Particulate filter

PM - Particulate matter

Re - Reynolds number

SI - Systeme international d'unites

VC - Vena contracta

WLTC - Worldwide harmonised light vehicles test cycle

[DELTA]P - Backpressure

SYMBOLS

A - Cross-sectional area of the PF, [m.sup.2]

[A.sub.j] - Cross-sectional area available for flow at location j in Fig. 6, [m.sup.2]

[C.sub.C] - Coefficient of contraction, i.e. the cross-sectional area of the VC divided by the area of the orifice through which gas flows to form the VC, -

[C.sub.pg,i], [C.sub.pg,o] - Specific heat capacity at constant pressure of the gas in the inlet and outlet channels, J [kg.sup.-1] [K.sup.-1]

[C.sub.pg,W] - Specific heat capacity at constant pressure of the gas passing through the soot cake and PF wall, J [kg.sup.-1] [K.sup.-1]

[C.sub.pS], [C.sub.pW] - Specific heat capacity at constant pressure of the soot and PF wall, J [kg.sup.-1] [K.sup.-1]

d - Width of soot-free inlet channel, m

[d.sub.i] - Width of soot-loaded inlet channel, m

[d.sub.0] - Width of outlet channel, m

F - Viscous loss coefficient, -

f - Transverse dependence of axial velocity, -

[H.sub.i] - Specific enthalpy of the gas in the inlet channel, J [kg.sup.-1]

[h.sub.i], [h.sub.o] - Heat transfer coefficients for heat transfer between the gas in the inlet and outlet channels and the PF wall, W [m.sup.-2] [K.sup.-1]

[k.sub.S], [k.sub.W] - Permeability of the soot cake and of the PF wall, [m.sup.2]

L - PF length, m

[L.sub.Plug] - Length of the plugs at the ends of the channels, m

[??] - Gas mass flow into the PF, kg [s.sup.-1]

M - Molar mass of gas, kg [mol.sup.-1]

Nu - Nusselt number, [d.sub.i][h.sub.i] / [[lambda].sub.i] or [d.sub.o][h.sub.o] / [[lambda].sub.o], -

[P.sub.atm] - External (atmospheric) pressure, Pa

[P.sub.i], [P.sub.o] - Gas pressure in the inlet and outlet channels, Pa

[P.sub.j] - Pressure at location j in Fig. 6, Pa

Pr - Prandtl number, [micro] [C.sub.p] / [lambda], -

Q - Rate of heat generation per unit volume by reaction, W [m.sup.-3]

R - Molar (or ideal) gas constant, J [mol.sup.-1] [K.sup.-1]

[Re.sub.W] - Wall Reynolds number, d [[phi].sub.W] / [[mu].sub.i] or d [[phi].sub.W] / [[micro].sub.o], -

[Re.sub.Z] - Axial Reynolds number, [d.sub.i][[phi].sub.i]/[[micro].sub.i] or [d.sub.o][[phi].sub.o] / [[micro].sub.o], -

S - Soot loading, kg [m.sup.-3]

t - Time, s

[T.sub.i], [T.sub.o] - Mixing cup temperature of the gas in the inlet and outlet channels, K

[T.sub.ln] - Gas inlet temperature, K

[T.sub.W] - Temperature of the PF wall and soot cake (and of gas passing through these), K

u - Velocity of gas in the channels as a function of position, m [s.sup.-1]

V - Mean velocity of gas in inlet channel, i.e. <[u.sub.i]>, m [s.sup.-1]

[V.sub.o] - Mean velocity of gas in outlet channel, i.e. <[u.sub.o]> , m [s.sup.-1]

[V.sub.W] - Superficial gas velocity through the soot cake and PF wall referenced to d, m [s.sup.-1]

[w.sub.S], [w.sub.W] - Thickness of the soot cake and the PF wall, m

y - Distance or coordinate in the through-wall direction, m

z - Axial distance or coordinate, m

[alpha] - Momentum flux correction factor, <[u.sup.2]> / [<u>.sup.2], -

[[alpha].sub.e] - Kinetic energy flux correction factor, <[u.sup.3]> / [<u>.sup.3], -

[[alpha].sub.j] - Momentum flux correction factor at location j in Fig. 6, -

[[beta].sub.S], [[beta].sub.W] - Forchheimer constant for the soot cake and PF wall, [m.sup.-1]

[delta]z - An element of length in the axial direction, m

[epsilon] - Open cross-sectional area of the soot-free PF in the region away from the plugs divided by the cross-sectional area of the PF, -

[[epsilon].sub.i] - Open area fraction of the inlet face of the PF, i.e. open area of front face divided by total face area, -

[[epsilon].sub.o] - Open area fraction of the outlet face of the PF, i.e. open area of rear face divided by total face area, -

[[zeta].sub.i] - Contraction loss coefficient, -

[[zeta].sub.o] - Expansion loss coefficient, -

[[lambda].sub.i], [[lambda].sub.o] - Thermal conductivity of the gas in the inlet and outlet channels, W [m.sup.-1] [K.sup.-1]

[[lambda].sub.S], [[lambda].sub.W] - Thermal conductivity of the soot cake and the PF wall, W [m.sup.-1] [K.sup.-1]

[[micro].sub.i], [[micro].sub.o] - Viscosity of the gas in the inlet and outlet channels, Pa s

[[micro].sub.W] - Viscosity of the gas passing through the soot cake and the PF wall, Pa s

[rho] - Density, kg [m.sup.-3]

[[rho].sub.Cell] - Cell density, i.e. number of channels per unit cross-sectional area of the whole PF, [m.sup.-2]

[[rho].sub.g,W] - Density of the gas passing through the soot cake and the PF wall, kg [m.sup.-3]

[[rho].sub.i], [[rho].sub.0] - Density of the gas in the inlet and outlet channels, kg [m.sup.-3]

[[rho].sub.s] - Soot packing density, kg [m.sup.-3]

[[rho].sub.W] - Density of the PF wall, kg [m.sup.-3]

[[alpha].sub.i] - Contraction ratio, i.e. ratio of open area of the PF front face divided by the cross-sectional area for flow before the gas enters the PF, -

[[sigma].sub.0] - Expansion ratio, i.e. ratio of open area of the PF rear face divided by the cross-sectional area for flow after the gas leaves the PF, -

[[phi].sub.i] - Mean mass flux of gas in inlet channel, [[rho].sub.i][V.sub.i], kg [m.sup.-2] [s.sup.-1]

[[phi].sub.o] - Mean mass flux of gas in outlet channel, [[rho].sub.o][V.sub.o], kg [m.sup.-2] [s.sup.-1]

[[phi].sub.W] - Through-wall mass flux referenced to d, [[rho].sub.g, W] [V.sub.W], kg [m.sup.-2] [s.sup.-1]

< > - Average of the bracketed quantity over a cross section

APPENDIX

APPENDIX A: THE VISCOUS LOSS TERM IN THE MOMENTUM BALANCE EQUATIONS

The friction factor for fluid flow in a conduit is defined by [66]:

Force = (Friction factor)(Wetted surface area)(1/2[rho][V.sup.2]) (A.1)

Note that V is the average velocity across the cross section of the conduit. The friction factor for laminar flow has the form Constant/Re [66], where [Re.sub.z] = d[rho]V/[micro]. The wetted surface area of an element of square conduit of length [delta]z is 4d [delta]z, where d is the width of the conduit. Substituting in these values gives:

Force = Constant[mu] / d[rho]V (4d[delta]z)(1/2[rho][V.sup.2]) = 2 (Constant) [micro]V[delta]z = F[micro]V[delta]z (A.2)

The right hand term above is the viscous loss term from the momentum balance equation (Eq. (18) and (19)). Note that F is twice the constant from the expression for the friction factor or 2(Friction factor)[Re.sub.z].

APPENDIX B: PHYSICAL PROPERTIES

Table B.1. Physical properties used for the simulations. Some values have been converted into SI units from the values given in the original reference. Where a function has been fitted to measured data, this is noted; otherwise the functions listed were taken directly from the reference. T is the temperature in K. Quantity Value [C.sub.pg,i], [C.sub.pg,0] 934.5 + 0.1970 T + 1.133x[10.sup.6] [T.sup.-2] J [kg.sup.-1] [K.sup.-1] (a,b) [C.sub.pW] 070 + 0.156 T - 3.44 x [10.sup.7] [T.sup.-2] J [kg.sup.-1] [K.sup.-1] (c) F 28.454d M 0.028 kg [mol.sup.-1] Nu 2.98 (d,e) [alpha] 1.377 (d) [[alpha].sub.2], [[alpha].sub.5] 1.377 [[alpha].sub.6] 1 (f) [[alpha].sub.e1] V 2.717 x [10.sup.-4] [T.sup.0.8023] W [m.sup.-1][K.sup.-1] (a,b) [lambda]w 0.8 W [m.sup.-1] [K.sup.-1] (c) [micro] 3.531 x [10.sup.-7] [T.sup.0.6968] Pa s (a,b) pw 925 kg [m.sup.-3] (c,g) Quantity Reference [C.sub.pg,i], [C.sub.pg,0] [67] [C.sub.pW] [15] F [21,48,68] M - Nu [68] [alpha] [28] [[alpha].sub.2], [[alpha].sub.5] [28] [[alpha].sub.6] - [[alpha].sub.e1] - [60] [lambda]w [15] [micro] [59] pw [69] (a) Function fitted by to data from the reference (Fig. B. 1, B.2). (b) Value for air. (c) Value for cordierite. (d) In section 6, Rew-dependant values will be used instead of these. (e) Value for laminar flow in a square conduit with constant wall temperature when axial fluid heat conduction is negligible. (f) Fully developed turbulent flow has been assumed, so this is the value for a flat velocity profile. (g) Based on an intrinsic density for cordierite of 2500 kg [m.sup.-3] [69] and an assumed porosity of 63%.

Timothy C. Watling and Maya R. Ravenscroft

Johnson Matthey Technology Centre

Jason P.E. Cleeton, Ian D. Rees, and David A.R. Wilkins

Johnson Matthey ECT

(1) From a standard identity [33], which is easily proved by multiplying out the terms:

[nabla][rho]uu = [rho]u [nabla]u + u[nabla] * [rho]u

Under steady-state conditions the mass balance is [nabla] * [rho]u = 0 [32] and hence [nabla] * [rho]<uu = [rho]u * [nabla]u.

(2) The full version of this equation is as below [34], but the second term is zero for an ideal gas

dH = ([partial derivative]H / [[partial derivative]T).sub.P]dT + [([partial derivative]H / [partial derivative]P).sub.T] dP = [C.sub.p]dT + 1 / [rho] - T[([partial derivative][[rho].sup.-1] / [partial derivative]T]).sub.P]dP

(3) Compared to the modern form of Darcy's law, Darcy's original equation [37] lacks the viscosity term.

(4) In Forchheimer's original paper [38], pressure drop per unit length of column for water flowing through a column of sand was modelled by two equations, viz. aV + b[V.sup.2] and m[V.sup.n], where a, b, m and n are constants. The former equation would appear to be the origin of the modern Forchheimer extension to Darcy's law.

(5) The Forchheimer term is also non-linear in velocity, but not listed here as its contribution is negligible. The pressure drop of all the nonlinear terms increase with the square of the velocity

doi:10.4271/2017-01-0974

Table 1. Comparison of literature values for [[zeta].sub.i] and [[zeta].sub.0] with those from Eq. (45) and (50). "Method" is the way in which [[zeta].sub.i] and [[zeta].sub.0] values were determined. For [[zeta].sub.i] and [[zeta].sub.0] defined by equations, calculated values over the range relevant for PF are given (0<[[epsilon].sub.i]<0.4, 0.2<[[epsilon].sub.0]<0.4). Where a single value is given for [[zeta].sub.i] and [[zeta].sub.0] this is a combined value for [[zeta].sub.i] + [[zeta].sub.0]. Values in parenthesis after the values of [[zeta].sub.i] or [[zeta].sub.0] are the corresponding value for [[epsilon].sub.i] or [[epsilon].sub.0], as appropriate. Reference Method [[zeta].sub.i] [[zeta].sub.0] This work Equation 1.9 (0.4)-2.2 (0) 0.5 (0.2)-0.8 (0.4) [17] Equation (a) 0.9 (0.4)-1.1 (0) 0.4 (0.4)-0.6 (0.2) [19] CFD 3.0 (0.16-0.32) 0.2 (0.32) [19] Empirical 2.5 (0.32) (b) [36] Empirical 2.9 (0.32) - 7.2 0.33) (b) [35] CFD 2.2 (0.37) - 8.2 (0.37) (b) [48] Literature (c) 0.2 - 0.8 0.2-0.8 (a) Standard equation from literature for heat exchangers [54]. (b) Combined value for [[zeta].sub.i] + [[zeta].sub.0]. (c) Values obtained from standard equation in literature. Table 2. PF wall permeability and Forchheimer constant obtained by fitting the data in Fig. 13. All simulations were run with these values. Parameter Value kw 5.95x [10.sup.-13] [m.sup.2] [beta]w 4.48 x [10.sup.7] [m.sup.-1]