Skip to main content

Modelling, asymptotic analysis and simulation of the gas dynamics in a chimney

Abstract

This paper presents a one dimensional model for the gas dynamics in a chimney. This is a prototype example of a small Mach number flow with strong heat sources. Due to the small Mach number of the gas flow an asymptotic model is derived from a fully compressible model, which then is compared to the original model and to the often used Boussinesq approximation. The Boussinesq approximation is shown to be inappropriate for this application, whereas the small Mach number asymptotic model we propose is shown to be a very good approximation. In particular it allows very fast numerical simulations. Finally, all this is underlined by numerical simulations where we validate the various models.

MSC:35C20, 35Q35, 76M45, 76N15, 80A20.

1 Introduction

Chimneys are typical representatives of the beginning of the industrial era in the 18th century. However there is a much longer story behind. Starting from the 13th century chimneys were build and used to channel away smoke and hot air. Nowadays chimneys are common presences in our everyday life, from industrial applications to the heating systems in our houses. Just to give some examples, we can think to modern heating systems, whit different kinds of fuel, and also to Controlled Ventilation Systems, where the basic principle of chimney is used, often with the help of a fan, to keep an healthy level of oxygen in modern buildings [1].

Here by chimney we mean a vertical structure for venting into the outside atmosphere hot flue gases produced e.g. in a combustion process like in a boiler, in a stove, in a furnace or in a fireplace. The underlying main physical principles are well understood. Therefore it is considered only a technical issue to give precise and quantitative predictions of the gas dynamics in a chimney. On one hand there is huge technological know how (for example see [2] and [3]) with empirical formulas for very specific situations based on restrictive (often unrealistic) assumptions. On the other hand one can find sophisticated CFD gasdynamics simulations based on time expensive numerical simulations. There seems to be a peculiar lack of models in between these two extreme areas, i.e. models which are physically reasonably realistic and numerically cheap at the same time. This paper goes exactly in this direction.

When describing the gas dynamics in a chimney a few crucial facts can be identified.

  • The flow is driven by density gradients (buoyancy force).

  • The flow velocities are slow compared to the speed of sound.

  • The temperature differences are big and decisive for the flow behaviour.

Small flow velocities indicate an incompressible model. In combination with the buoyancy as driving force a model with variable density is necessary. In such a context a widely used approach is the Boussinesq approximation going back originally to Boussinesq [4] with later justifications by Bois, Zeytounian etc. (see e.g. [57] and citations therein). However we will show that here a Boussinesq type approach is not appropriate. As an alternative we use a variable density low Mach number model (see [812] for the non isentropic or the non isothermal cases) adapted to this approximation which results to be realistic and numerically cheap at the same time.

The aim of the paper is to study the flow dynamics both qualitatively and quantitatively. A common measure for the quality of the full process in a chimney is the draft in the chimney. And the draft in the chimney depends e.g. on the heat source, on the chimney geometry, on the chimney friction and heat loss parameters, on the outside meteorological situation etc.

In Section 2 we describe the underlying equations. In Section 3 we will scale the equation and identify various asymptotic regions. In Section 4 we discuss the non validity of a Boussinesq type model and propose a variable density low Mach number model. Finally in Section 5 we will provide several numerical examples, first using simulations on a full system of equations as a reference to underline the validity of our approach, and then analysing in more details the dependence on the parameters.

2 The model

We describe the motion of smoke and flue gases produced by combustion trough a chimney. An ideal example that we can have in mind is that of a home fireplace (Figure 1), but the model applies to a great variety of other cases (modern heating and ventilation systems, industrial chimneys etc.).

Figure 1
figure 1

Basic setting of the chimney and the fireplace.

The fresh cold air enters at the base of the structure and, after flowing through a small section of the chimney, it reaches the fireplace. There a combustion process takes place and the air gains some heat and mass (in the combustion process the air mixture typically looses oxygen and gains carbon dioxide and water vapor, the total variation resulting in an net increase of the air mass). As a result the density decreases and the velocity increases. The hot air smoke mixture enters the main part of the chimney, where the so called chimney effect acts as the main driving force. Clearly, in that part there are also other effects like friction and heat losses. Finally the hot air smoke mixture exits at the top of the structure.

We assume that the main dynamics acts in vertical direction, as we do not expect significant flow phenomena in other directions. Then almost all gas particles experience the same when passing through the structure. Therefore we describe the gas dynamics by a one dimensional model, where the single dimension x starts at the base of the structure and follows the air flow towards the top of the chimney. All quantities under consideration have to be interpreted as mean values over the cross sectional area. However, this approach still allows to consider variable cross sections e.g. a larger cross section in the fire place area.

Keeping in mind the above considerations and simplifications we start with a one dimensional fully compressible model. Let ρ ˜ , u ˜ , p ˜ and T ˜ be respectively the density, the velocity, the pressure and the temperature of the air flow, all averaged on the cross section A ˜ and functions of space x ˜ and time t ˜ . The space variable x ˜ is zero at the bottom of the structure and its positive direction points toward the chimney top. We consider an unscaled version of the Euler equations of gas dynamics with a balance equation for mass

( A ˜ ρ ˜ ) t ˜ + ( A ˜ ρ ˜ u ˜ ) x ˜ = A ˜ m ˜ ,
(1)

a balance equation for momentum

( A ˜ ρ ˜ u ˜ ) t ˜ + ( A ˜ ( ρ ˜ u ˜ 2 + p ˜ ) ) x ˜ = A ˜ x ˜ p ˜ A ˜ g ρ ˜ λ 2 D A ˜ ρ ˜ u ˜ | u ˜ |+ Δ ˜ p f a n A ˜ L ,
(2)

and a balance equation for energy

( A ˜ ( c v p ˜ R + ρ ˜ u ˜ 2 2 + g ρ ˜ h ˜ ) ) t ˜ + ( A ˜ u ˜ ( c v p ˜ R + ρ ˜ u ˜ 2 2 + g ρ ˜ h ˜ ) + A ˜ p ˜ u ˜ ) x ˜ = A ˜ q ˜ k 2 ( T ˜ T ˜ w ) A ˜ σ ( T ˜ T ˜ w ) 4 .
(3)

This system is completed by the ideal gas law

p ˜ =R ρ ˜ T ˜ .

In (1) the term on the right hand side m ˜ is due to the additionally produced gas in the combustion process as a consequence of the oxygen consumption. We will propose a model for this mass variation in the next section.

The term A ˜ g ρ ˜ on the right hand side of equation (2) describes the gravitational force which is volume dependent and so linearly dependent on the cross section. The friction on the surface of the chimney walls is described by λ 2 D A ˜ ρ ˜ u ˜ | u ˜ | where λ is the friction factor and

D=2 A ˜ / π

is the so called hydraulic diameter. Realistic values for λ can be found in literature, see, for example [13].

In case of an additional fan (for example in low temperature chimneys) a pressure gain term Δ ˜ p f a n /L is introduced, where L is the reference value for the length of the chimney.

In the energy equation (3) h=h(x) is the height of the structure and c v is the specific heat of dry air at constant volume. On the right hand side of (3) we have the heat source term due to the combustion process A ˜ q ˜ , the heat loss trough the chimney walls k 2 ( T ˜ T ˜ w ) and the energy loss due to thermal radiation A ˜ σ ( T ˜ T ˜ w ) 4 . Here T ˜ w is the wall temperature, k depends on the cross section and on the heat loss coefficient α according the relation

k=α2 A ˜ π .

A typical description for the heat loss coefficient is

α=St( c v +R) ρ ˜ u ˜ ,

where St=λ/8 is the Stanton number [14]. Using average density and velocity values we can compute α. Note that in Section 5 we will use slightly smaller α values to take into account the presence of insulating materials typically used in chimney walls [2, 3, 15, 16].

Finally σ=5.670400× 10 8  J s 1  m 2  K 4 is the thermal radiation coefficient (Stefan’s constant).

2.1 Modelling the mass loss and the heat source

We assume both terms m ˜ and q ˜ depending on the air flow ρ ˜ u ˜ . Assuming that in the fireplace there is always enough fuel, the reaction rate depends on the amount of oxygen available for the combustion process. The amount of oxygen (being a fixed percentage in the air) is linearly dependent on the mass flow. Thus we set

q ˜ = c ˜ q ρ ˜ | u ˜ |

and

m ˜ = c ˜ m ρ ˜ | u ˜ |,

where c ˜ m and c ˜ q are a mass variation and a heat production coefficient. To determine reasonable values to these coefficients we recall that burning wood produces 106 J/kg of thermal energy ([17, 18] and [19]), while chimneys used for heating in common handholds have usually a power of 102 kW [15, 16]. Considering this and taking into account the reference values for density and velocity, we will use a heat production coefficient c ˜ q 10 5 10 6  kg m/s 2 .

As for the mass variation in the air, there are different possible approaches [20]. We can consider the burning reaction and compare the wood mass loss with the mass gain in the air, and then use the results from [18] and [19]. Otherwise we can use the experimental results on the oxygen variation in air after combustion. We will use c ˜ m 10 2  1/m.

3 Scaling

Now we scale the problem. Every quantity z ˜ is associated by z ˜ = z r z to a dimensionless quantity z, where we use the sub-index r to indicate the reference values. Similar, every function f ˜ is related by f r f(x,t)= f ˜ ( x ˜ , t ˜ ) to a dimensionless function f. Realistic reference values are listed in Table 1.

Table 1 Scaling table with typical reference values

In this way we obtain the scaled model

ρ t + ( ρ u ) x = A x A ρu+ c m ρ|u|,
(4)
ρ u t +ρu u x + 1 γ M 2 p x = c m ρ|u|u 1 F r 2 ρ ξ A ρu|u|+Δ p f a n ,
(5)
( p + γ M 2 ( γ 1 ) ρ u 2 2 + γ M 2 ( γ 1 ) F r 2 ρ h ) t + 1 A [ A u ( p + γ M 2 ( γ 1 ) ρ u 2 2 + γ M 2 ( γ 1 ) F r 2 ρ h ) + ( γ 1 ) A u p ] x = c q ρ | u | a ( T T w ) A r T 4 ,
(6)

while the state equation becomes

p=ρT.
(7)

(Here ρ r := p r /R T r , see Table 1.)

We introduced the dimensionless constants

γ= c v + R c v , 1 γ M 2 = p r ρ r u r 2 , 1 F r 2 = L u r 2 g,

where M is the Mach-number and Fr is the Froude-number. Moreover we define

ξ= λ π L 4 A r ,Δ p f a n = 1 ρ r u r 2 Δ ˜ p f a n ,r= t r T r 4 σ p r (γ1)

and

c m = c ˜ m L, c q = c ˜ q ( γ 1 ) ρ r L p r ,a= α T r π t r A r p r (γ1).

In Table 2 we remember the values of the constants used in the model.

Table 2 Constants

We note that in our setting M 10 3 and therefore 1 M 2 10 6 . In addition we have Fr 10 1 and therefore 1 F r 2 10 2 . These facts suggest to perform an asymptotic analysis with the intention to further simplify the model.

3.1 Initial and boundary conditions

Equations (4)-(7) need to be endowed with the appropriate initial and boundary conditions. We assign initial data for the velocity u, the density ρ and the pressure p:

u(x,0)= u 0 (x),ρ(x,0)= ρ 0 (x),p(x,0)= p 0 (x)
(8)

and boundary condition for p:

p(0,t)= p l ,p(1,t)= p r ,
(9)

where p l and p r are, respectively, pressure at the bottom and pressure at the top of the structure.

Finally, as the density ρ satisfies a transport equation, we prescribe inflow conditions:

ρ ( 0 , t ) = ρ l if  u ( 0 , t ) > 0 , ρ ( 1 , t ) = ρ r if  u ( 1 , t ) < 0 .
(10)

We still need to assign reasonable values to p l , p r , ρ l and ρ r . Thus we study the gasdynamics outside the structure. Assuming a stationary state without sources and chimney effects we have from the equations (4)-(6)

( ρ h u ) x = 0 , ρ h u u x + 1 ε ( p h ) x = 1 F r 2 ρ h , [ u ( p h + ε ( γ 1 ) ρ h u 2 2 + ε ( γ 1 ) F r 2 ρ h h ) + ( γ 1 ) u p h ] x = 0 ,
(11)

where ρ h and p h are the outside density and pressure. If in addition we assume no flow (i.e. u=0) we obtain the scaled adiabatic formulas for density, pressure and temperature [13, 21]

ρ h = ρ h 0 ( 1 ( γ 1 ) M 2 F r 2 h ( x ) ) 1 γ 1 ,
(12)
p h = p h 0 ( 1 ( γ 1 ) M 2 F r 2 h ( x ) ) γ γ 1 ,
(13)
T h = T h 0 ( 1 ( γ 1 ) M 2 F r 2 h ( x ) ) .
(14)

Remember that in a vertical setting the height is h(x)=x, x=0 at the bottom, x axis oriented toward the top of the chimney. In our scaling the values at the bottom are given by ρ h 0 =1 and p h 0 =1.

Now we can complete (9) and (10) imposing:

p l = p h | x = 0 =1, p r = p h | x = 1 = ( 1 ( γ 1 ) M 2 F r 2 ) γ γ 1 ,
(15)
ρ l = ρ h | x = 0 =1, ρ r = ρ h | x = 1 = ( 1 ( γ 1 ) M 2 F r 2 ) 1 γ 1 .
(16)

Thus the problem to solve consists of the equations (4)-(7), the initial conditions (8), the boundary condition (9) and (10) specified in (15) and (16).

4 Asymptotics

In the previous section we encountered two small parameters, the Mach number M 10 3 and the Froude number Fr 10 1 . In our application they seem to be of a different order of magnitude. Therefore we do not consider the so called quasistatic case

MFr1.

However, the cases

F r 2 M1,
(17)

and

Fr fixed,M1
(18)

should not be ruled out. The case (17) leads to the so called Boussinesq approximation, the case (18) is known as the small Mach number limit. Both approximations are widely used, see [47] and citations therein for the Boussinesq approximation. For the small Mach number approximation in general (non isentropic and non isothermal) see [812, 22], for numerical issues see [2326] and references therein. For a variety of one dimensional applications similar to the one in this paper see [2730]. In the literature (for example in [31]) the Boussinesq number is defined as B= M 2 /F r 2 . Both cases (17) and (18) refer to the limit B0, whereas the quasistatic case would correspond to a limit with B B 0 0.

In the following for sake of simplicity we neglect the terms involving the fan and the radiation energy. Also, we deal with a vertical chimney such that h(x)=x holds. We start with the classical Boussinesq approximation.

4.1 Boussinesq approximation

In this section we assume MF r 2 and set c= M F r 2 =O(1). The following derivation mimics the asymptotic derivation of the Boussinesq approximation presented by [7] extended to the case with mass losses, heat sources etc. (see also Remark 2 at the end of this section). Note that from the state equation we get T= p ρ so that we consider the equations (4)-(6) for the unknowns ρ, p and u. For z=ρ,p,u, we start with the following asymptotic expansion in the small Mach number M

z(x,t)= z ( 0 ) (x,t)+ z ( 1 ) (x,t)M+ z ( 2 ) (x,t) M 2 +O ( M 3 ) ,

to obtain in (4)-(6) a hierarchy of equations.

From (4) in O( M 0 ) and O( M 1 ) we obtain

ρ t ( 0 ) + ( ρ ( 0 ) u ( 0 ) ) x = A x A ρ 0 u 0 + c m ρ ( 0 ) | u ( 0 ) |,
(19)
ρ t ( 1 ) + ( ρ ( 1 ) u ( 0 ) + ρ ( 0 ) u ( 1 ) ) x = ( c m A x A ) ( ρ ( 1 ) | u ( 0 ) | + ρ ( 0 ) | u ( 1 ) | ) .
(20)

From (5) in O( M 2 ), O( M 1 ), O( M 0 ) we get

p x ( 0 ) =0,
(21)
p x ( 1 ) =γc ρ ( 0 ) ,
(22)
ρ ( 0 ) u t ( 0 ) + ρ ( 0 ) u ( 0 ) u x ( 0 ) + 1 γ p x ( 2 ) = c m ρ ( 0 ) | u ( 0 ) | u ( 0 ) c ρ ( 1 ) ξ A ρ ( 0 ) u ( 0 ) | u ( 0 ) |.
(23)

From (6) in O( M 0 ) and O( M 1 ) we obtain

( p ( 0 ) ) t + ( γ u ( 0 ) p ( 0 ) + c m ρ ( 0 ) u ( 0 ) T ( 0 ) ) x = A x A γ ρ ( 0 ) p ( 0 ) c q ρ ( 0 ) | u ( 0 ) |a ( p ( 0 ) ρ ( 0 ) T w ) A ,
(24)
( p ( 1 ) + c ( γ 1 ) γ ρ ( 0 ) x ) t + 1 A ( A γ ( u ( 0 ) p ( 1 ) + u ( 1 ) p ( 0 ) ) ) x = c q ( ρ ( 0 ) | u ( 1 ) | + ρ ( 1 ) | u ( 0 ) | ) a p ( 1 ) ρ ( 0 ) ρ ( 1 ) ( ρ ( 0 ) ) 2 A .
(25)

Expanding the boundary conditions (15) we find

p ( 0 , t ) = 1 , p ( 1 , t ) = ( 1 ( γ 1 ) c M ) γ γ 1 = 1 γ c M + γ c 2 M 2 ,

and we deduce

p ( 0 ) ( 0 , t ) = 1 , p ( 0 ) ( 1 , t ) = 1 , p ( 1 ) ( 0 , t ) = 0 , p ( 1 ) ( 1 , t ) = c γ , p ( 2 ) ( 0 , t ) = 0 , p ( 2 ) ( 1 , t ) = + γ c 2 .

In a similar way, the inflow condition for ρ ( 0 ) and ρ ( 1 ) are given by

ρ ( 0 ) ( 0 , t ) = 1 if  u ( 0 , t ) > 0 , ρ ( 0 ) ( 1 , t ) = 1 if  u ( 1 , t ) < 0 , ρ ( 1 ) ( 0 , t ) = 0 if  u ( 0 , t ) > 0 , ρ ( 1 ) ( 1 , t ) = c if  u ( 1 , t ) < 0 .

The main assumption in the classical Boussinesq approximation is a constant leading order density ρ ( 0 ) =1, where the value is fixed to 1 due to the (scaled) boundary data. From (21), (22) and the boundary and inflow conditions we conclude

p ( 0 ) (x,t)1, p ( 1 ) = p ( 1 ) ( ρ ( 0 ) ) =γc ρ ( 0 ) x.

This leads in (19) and (24) to the contradiction

u x ( 0 ) = c m ρ ( 0 ) | u ( 0 ) | , u x ( 0 ) = c q ρ ( 0 ) | u ( 0 ) | a ( p ( 0 ) ρ ( 0 ) T w ) A .

This means that a mass loss or a heat source of order one is not compatible with the Boussinesq assumption of a constant leading order density. This is reasonable since order one mass loss or heat source terms are supposed to induce significant order one changes in the density (and therefore in temperature). We conclude that the Boussinesq approximation is NOT applicable in our case of significant mass losses or heat source. Therefore in Section 4.2 we pass to the small Mach number asypmtotics.

Remark 1 Only for small mass losses and heat sources, i.e. in the case

c m =O(M), c q =O(M),a=O(M)

the Boussinesq approximation would be applicable and we obtain the well known ‘incompressibility’ condition u x ( 0 ) =0 and a u ( 0 ) = u ( 0 ) (t). Then (20) and (25) become the two equations for the perturbations u ( 1 ) and ρ ( 1 ) . Finally ρ ( 1 ) enters in the buoyancy term on the right hand side of (23).

For simplicity let’s assume the very simple case c m = c q =a=0 corresponding to the outside case with no sources. Then equation (25) reduces to

u x ( 1 ) =c u ( 0 )

which gives in (20)

ρ t ( 1 ) + u ( 0 ) ( ρ ( 1 ) ) x =c u ( 0 ) .

Assuming positive velocities u ( 0 ) and the boundary condition ρ ( 1 ) (0,t)=0 gives the long time behaviour

ρ ( 1 ) (x,t)=cx.

Thus

ρ(x,t)=1Mcx+O ( M 2 ) =1 M 2 F r 2 x+O ( M 2 )

which are exactly the first two terms in the expansion of the adiabatic formula (12) for the density. In Section 5.2 we come back to these asymptotics when doing numerical comparisons.

Remark 2 We said that this derivation of the Boussinesq approximation is in the line of [7] which is different from the more classical asymptotic derivations going back to [5, 6]. In fact in [5, 6] a multiple scale analysis is used to obtain the final equations. However, the non applicability of the Boussinesq approximation in our application is independent of the chosen approach. For details on the differences between the approaches see [7].

4.2 Small Mach number approximation

Here we keep the Froude number Fr fixed and consider only the Mach number to be small. This limit is well known (see citations at the beginning of this section). Note again that from the state equation we get T=p/ρ so that we consider only the equations (4)-(6) for the unknowns ρ, p and u.

Since the equations (4)-(6) only contain M 2 (and not M itself) for z=ρ,p,u, we make an Ansatz of an asymptotic expansion in the small parameter M 2

z(x,t)= z ( 0 ) (x,t)+ z ( 2 ) (x,t) M 2 +O ( M 4 ) ,

which leads again to a hierarchy of equations. From (4) in O( M 0 ) we have

ρ t ( 0 ) + ( ρ ( 0 ) u ( 0 ) ) x = A x A ρ ( 0 ) u ( 0 ) + c m ρ ( 0 ) | u ( 0 ) |.
(26)

From (5) in O( M 2 ) and O( M 0 ) we obtain

p x ( 0 ) =0,
(27)
ρ ( 0 ) u t ( 0 ) + ρ ( 0 ) u ( 0 ) u x ( 0 ) + 1 γ p x ( 2 ) = c m ρ ( 0 ) | u ( 0 ) | u ( 0 ) 1 F r 2 ρ ( 0 ) ξ A ρ ( 0 ) u ( 0 ) | u ( 0 ) |.
(28)

Finally (6) gives in O( M 0 )

( p 0 ) t + 1 A ( γ A u ( 0 ) p ( 0 ) ) x = c q ρ ( 0 ) | u ( 0 ) |a ( p ( 0 ) ρ ( 0 ) T w ) A ,

where we used T ( 0 ) = p ( 0 ) / ρ ( 0 ) .

The boundary conditions for p can be expanded

p ( 0 , t ) = 1 , p ( 1 , t ) = ( 1 ( γ 1 ) M 2 F r 2 ) γ γ 1 1 γ M 2 F r 2 ,

from which we obtain

p ( 0 ) (0,t)=1, p ( 0 ) (1,t)=1,
(29)
p ( 2 ) (0,t)=0, p ( 2 ) (1,t)= γ F r 2 .
(30)

In a similar way we get the inflow conditions for ρ ( 0 ) :

ρ ( 0 ) (0,t)=1if u(0,t)>0, ρ ( 0 ) (1,t)=1if u(1,t)<0.
(31)

From the boundary condition (29) on p ( 0 ) and from the O( M 2 ) relation (27) obtained from the momentum equation, we conclude p ( 0 ) (x,t)1. Therefore we are left with three unknowns, ρ ( 0 ) , u ( 0 ) , p ( 2 ) satisfying a set of three equations

ρ t ( 0 ) + ( ρ ( 0 ) u ( 0 ) ) x = A x A ρ ( 0 ) u ( 0 ) + c m ρ ( 0 ) | u ( 0 ) |,
(32)
ρ ( 0 ) u t ( 0 ) + ρ u ( 0 ) u x ( 0 ) + 1 γ ( p ( 2 ) ) x = c m ρ ( 0 ) | u ( 0 ) | u ( 0 ) ρ ( 0 ) F r 2 ξ A ρ ( 0 ) u ( 0 ) | u ( 0 ) | ,
(33)
γ u x ( 0 ) =γ A x A u ( 0 ) + c q ρ ( 0 ) | u ( 0 ) | a A ( 1 ρ ( 0 ) T w ) .
(34)

These equations have to be completed with initial data for ρ ( 0 ) , u ( 0 ) and boundary conditions (30) for p ( 2 ) and (31) for ρ ( 0 ) .

Note that ρ ( 0 ) is not constant and describes eventually big density changes which imply big temperature changes. The same holds for u ( 0 ) . Contrary to the Boussinesq approximation we encounter no contradiction with order one mass loss or heat source terms. Therefore we proceed with this set of equations in the following sections.

As we can see in (32)-(34), the model contains 4 parameters which have to be chosen appropriately, i.e. c m , c q , ξ and a. The former two parameters depend on the chemistry of the fire, the latter two (friction and heat loss) only on the construction of the chimney. The quality of quantitative results obtained by this model depends on the precision of knowledge of the 4 mentioned parameters. By the way, this problem is in common with all approaches, in particular also with all the sophisticated CFD ones. The model, being 1 dimensional and not fully compressible, allows - as we will see also in the examples - fast simulations. Fast direct simulations are the basis of any optimisation strategy. In the planning stage of a chimney the diameter of the chimney, the height of the chimney, the size of the allowed heat loss etc. are typical parameters to be optimised.

5 Numerical simulations

5.1 A reformulation of the small Mach number equations

We are now interested in solving numerically equations (32)-(34) with initial data (8) for ρ ( 0 ) and u ( 0 ) , boundary conditions (30) for p ( 2 ) and inflow conditions (31) for ρ ( 0 ) .

First we notice that it is more convenient to consider the air pressure in relation to the atmospheric pressure. So, if we call p a c = p a c ( 0 ) + M 2 p a c ( 2 ) +O( M 4 ) the altitude corrected pressure, given by

p a c =p p h = ( p ( 0 ) 1 ) + M 2 ( p ( 2 ) + γ F r 2 x ) +O ( M 4 )

we obtain

ρ t + ( ρ u ) x = A x A ρu+ c m ρ|u|,
(35)
ρ u t +ρu u x + 1 γ ( p a c ( 2 ) ) x = c m ρ|u|u ( ρ 1 ) F r 2 ξ A ρu|u|+Δ p f a n ,
(36)
γ u x =γ A x A u+ c q ρ|u| a A ( 1 ρ T w ) r ( 1 ρ ) 4 ,
(37)

where we dropped the index ( 0 ) for ρ ( 0 ) and u ( 0 ) . The boundary conditions for the altitude corrected pressure are given by p a c (0,t)= p a c (1,t)=0. From now on we drop the superscript ac.

Also, here we included again the additional fan and the radiation term, that were previously neglected for sake of simplicity.

We reformulate (35)-(37) in order to simplify and in particular to obtain a system of evolution equations. First we integrate (37) in space on (0,x), and define a new function v(t):=A(0)u(0,t); then we rewrite (35) using (37); finally we integrate (36) in space on (0,1).

In this way we get

u ( x , t ) v ( t ) A ( x ) = 1 γ A ( x ) 0 x A ( y ) c q ρ ( y , t ) | u ( y , t ) | d y u ( x , t ) v ( t ) A ( x ) = 1 γ A ( x ) 0 x { a A ( y ) ( 1 ρ ( y , t ) T w ) + r ( 1 ρ ( y , t ) ) 4 } d y ,
(38)
ρ t +u ρ x = 1 γ ρ c q ρ|u|+ c m ρ|u|+ a γ A ρ ( 1 ρ T w ) + r γ ρ ( 1 ρ ) 4 ,
(39)
v t = ( 0 1 ρ A 1 ) 1 { 0 1 ρ u u x d x + 0 1 c m ρ | u | u d x + 1 F r 2 0 1 ( ρ 1 ) d x v t = + 0 1 ξ A ρ u | u | d x + Δ p ( 2 ) γ 0 1 Δ p f a n d x + 0 1 ρ γ A ( 0 x A c q ( ρ | u | ) t d y v t = + 0 x a A ( 1 ρ 2 ρ t ) d y + 4 0 x r A ( 1 ρ ) 3 ( 1 ρ 2 ρ t ) d y ) d x } ,
(40)

where Δ p ( 2 ) =( p ( 2 ) (1,t) p ( 2 ) (0,t))=0.

Thus we have to solve a PDE for the density ρ(x,t), an ODE for v(t), both evolution equations, and an integral equation for u(x,t).

Finally we have to discuss how the initial, boundary and inflow conditions translate in terms of these new equations.

The initial data (8) can be transformed in initial data for equations (39)-(40), while the boundary conditions for the pressure are included in equation (40). The inflow conditions (31) for the density stay unaltered.

To solve equations (38), (39) and (40) we use

  • an explicit upwind for the PDE for ρ in (39),

  • an explicit time step for the ODE for v in (40),

  • an explicit time step for the integral equation for u in (38).

In what follows i will be the index for the space discretization and n the one for the time discretization; moreover, for the sake of brevity, we summarize equations (38), (39) and (40) as

u = v A + F u [ ρ , u ] , ρ t = u ρ x + F ρ [ ρ , u ] , v t = F v [ ρ , ρ t , u , u x , u t ] ,

where F u , F ρ , F v are the right hand sides of (38), (39), (40), respectively.

With this notation, the explicit upwind for the PDE for ρ in (39) is given by:

ρ i n + 1 ρ i n Δ t + 1 2 ( u i n | u i n | ) ρ i + 1 n ρ i n Δ x + 1 2 ( u i n + | u i n | ) ρ i n ρ i 1 n Δ x = F ρ [ ρ i n , u i n ] .

To solve the explicit time step for the ODE for v in (40), we use:

v n + 1 v n Δ t = F [ ρ i n , 1 2 ( u i n | u i n | ) ρ i + 1 n ρ i n Δ x 1 2 ( u i n + | u i n | ) ρ i n ρ i 1 n Δ x + F ρ [ ρ i n , u i n ] , u i n , 1 2 ( u i n | u i n | ) u i + 1 n u i n Δ x 1 2 ( u i n + | u i n | ) u i n u i 1 n Δ x , u i n u i n 1 Δ t ] .

Finally, we compute

u i n + 1 = v n + 1 A i + F u [ ρ i n , u i n ] .

The integrals involved are solved by summation. To use ( u i n u i n 1 ) Δ t at the starting time (level n=0), we prescribe (the same) initial values for u also at a level n=1.

To assure stability of the scheme we only need to check CFL condition for the upwind scheme. This method is the simplest reasonable choice (from a numerical point of view) and promises to be fast. We will see that this simple choice already guarantees numerical robustness. In the future we plan to study more sophisticated numerical approaches to simulate our set of equations.

To be more precise we expect to be (much) faster than the corresponding simulations on the full set of compressible equations (4)-(6). This is due to the fact that the CFL condition in the compressible equations (4)-(6) is very restrictive due to big sound speed (due to the small Mach number). Contrary to that the CFL condition in (40)-(39) is determined by the flow velocity (which is of order one) and allows much bigger time steps. In addition every time step in the asymptotic model consists of an ODE and a single PDE timestep, whereas in the fully compressible model we have a time step in a coupled system of three PDE’s. In a similar context a precise comparison of the numerical effort was presented [29].

In the following examples we will use as geometrical data for our structure those given in Table 3; we will study the dependence on these data in Section 5.4.

Table 3 Geometrical data

We will run our simulations with final time T=20 s since within this time range the stationary state is reached in all our examples.

5.2 Example 1: Basic description - comparisons with other models

At first we would like to compare the results obtained from our model, with those obtained from the Boussinesq approximation (19)-(25) and from the full system (4)-(7), which will be used as our reference.

For the Boussinesq approximation, we consider equations (19)-(25) and assume

c m =M c ˜ m , c q =M c ˜ q ,a=M a ˜ ,

with c ˜ m , q ˜ m , a ˜ of order 1. If we also assume ρ ( 0 ) =1, trough computation similar to those used for the small Mach number case in Section 5.1, we get:

u t ( 0 ) = 1 0 1 d x A { Δ p ( 2 ) γ 0 1 ( c ρ ( 1 ) + ξ A u ( 0 ) | u ( 0 ) | A x A ( u ( 0 ) ) 2 ) d x } ,
(41)
ρ t ( 1 ) + u ( 0 ) ρ x ( 1 ) = ( c ˜ m c ˜ q γ ) | u ( 0 ) |γc u ( 0 ) + a ˜ γ A (1 T w ).
(42)

Here Δ p ( 2 ) =( p ( 2 ) (1,t) p ( 2 ) (0,t))=0, where we used the altitude corrected pressure; the boundary conditions are those obtained in Section 4.1. To solve these equations we use an explicit upwind scheme for ρ ( 1 ) and an explicit time step for u ( 0 ) .

Solving (41)-(42) we will be able to compute the unknowns in the zero order approximation, i.e.

ρ ρ ( 0 ) ,p p ( 0 ) ,T T ( 0 ) = p ( 0 ) ρ ( 0 ) ,u= u ( 0 ) ,
(43)

where only the velocity is not constant. In the first order approximation we have

ρ ρ ( 0 ) +M ρ ( 1 ) ,p p ( 0 ) +M p ( 1 ) ,T= T ( 0 ) + T ( 1 ) M= p ( 0 ) + p ( 1 ) M ρ ( 0 ) + ρ ( 1 ) M .
(44)

In order to solve the full system (4)-(7), we use the standard technique of splitting the problem into two sub-problems, solving first the hyperbolic conservation law by an upwind scheme, and then solving the ODE for the source term with an explicit Euler scheme. In this case, we choose to use a first order scheme in order to have a reasonable comparison with the result obtained from the other two models.

For simplicity, we consider a basic example with no heat loss through the chimney walls and constant cross section. We consider the values for our parameters as given in Table 4.

Table 4 Data for basic examples

As for the initial data, we consider u 0 (x)0, ρ 0 (x)= ρ h (x) and, in the case of the full system, p 0 (x)= p h (x).

After a simulation time of few seconds, all three models reach stationary values, and we can compare the result for velocity, density and temperature of air through the chimney.

For both the small Mach number and the full model we see that the temperature reaches a value of 250C (with ΔT 10 2 ) and a velocity of 3 m/s (Figures 2-4). We have very good agreement although the simulation time of the small Mach approach is about three orders of magnitude lower (and faster than real time).

Figure 2
figure 2

Temperature. Simulations for data as in Table 4. While the results obtained with the full system and the system of equations reduced via small Mach number approximation, coincide, the results obtained with the Boussinesq approximation are unreasonable.

Meanwhile, the Boussinesq approximation gives unrealistic results, involving negative density (Figure 3).

Figure 3
figure 3

Density. Simulations for data as in Table 4. While the results obtained with the full system and the system of equations reduced via small Mach number approximation, coincide, the results obtained with the Boussinesq approximation are unreasonable.

Figure 4
figure 4

Velocity. Simulations for data as in Table 4. While the results obtained with the full system and the system of equations reduced via small Mach number approximation, coincide, the results obtained with the Boussinesq approximation are unreasonable.

In Figure 2 after a first period where the density is still positive (and positive temperature), also the temperature becomes negative due to (44). This is not a numerical problem, but a problem of the wrong asymptotics. The reason is that for strong heat sources the Boussinesq approximation is a bad approximation since the higher order terms in the asymptotics are not of higher order but of the same order as the leading order terms.

On the other hand, lowering the heat sources, the Boussinesq approximation becomes reasonable, in the sense that the density and the temperature are no more negative.

When we finally compare the three models for no heat source, the solutions given from Boussinesq approximation, that computes the result up to order one, overlap those given by the full system (Figure 5).

Figure 5
figure 5

Temperature and density for no heat source. The results obtained with the (first order) Boussinesq approximation coincide with those obtained with the full system, while the small Mach number approximation coincides with the zeroth order Boussinesq approach.

5.3 Example 1.2: adding heat loss

Now we concentrate on the results obtained with our model.

We consider a variable cross section, and add the effect of the heat loss along the last 4 meters of the chimney walls setting the heat loss coefficient α=20W/( m 2 K) (see Figure 6).

Figure 6
figure 6

Adding heat loss. Data as in Table 4, but with α=20W/( m 2 K).

This choice is motivated by the fact that many chimneys are not isolated over the last meters (at least in older buildings) since there the chimney passes under the roof (which may be a not heated area) or even outside of the roof. In the heated part of the building chimneys are typically well isolated.

As we can see from Figure 6, when the air flow passes the fireplace, temperature increases (as density decreases). The wider cross section in the fireplace reduces the velocity. The stationary state is reached very rapidly.

The maximal velocity and maximal temperature are, respectively, 2 m/s and 240C, while at the top of the chimney, the exit velocity is about 1.8 m/s and the exit temperature is about 155C. All these are reasonable values. In modern chimneys, thanks to the use of special materials, the friction factor and the heat loss coefficient are kept considerably smaller, such that for small heat sources there is still a good draft.

We finally note that adding the radiation term in the energy equation does not affect the previous results since the Stefan constant is too small. Or in other words, the reached temperature is too low to induce significant losses due to radiation.

5.4 Example 2: effects of the geometry

Here we show how the draft velocity depends - for fixed cross sectional area at the fireplace - on the geometry of the structure.

Figure 7 shows the dependence on the radius of the chimney (for three different values of the heat loss coefficient α). For a small diameter the surface effects dominate. Therefore for a small diameter the friction losses are dominant. For bigger diameters the surface effects become less important. Then increasing the diameter the draft is decreasing since the total buoyancy force is distributed over a bigger cross section. As a final result there is an optimal diameter for which a maximal value for the velocity is reached.

Figure 7
figure 7

Exit velocity as a function of the chimney radius for three different heat loss coefficients.

On the other side we can vary the effective length of the chimney i.e. the length of the chimney above the fire place. In Figure 8 the result for two different values of heat loss coefficient is given. We assume heat loss only in the last part of the chimney (1/4 of the total effective length), all the other data are as in the previous example. In a short chimney the buoyancy force is too low, in a very high chimney the friction losses dominate. Again, there is an optimal effective length for which the velocity reaches a maximum.

Figure 8
figure 8

Exit velocity as a function of the effective length for two different heat loss coefficients.

5.5 Example 3: dependence on the friction factor

In Examples 1 and 2 we noticed the dependence of the draft on the heat loss coefficient α. Here we discuss the dependence on another surface effect, the friction factor λ. It is not surprising that increasing the friction factor the draft velocity decreases. In Figure 9 we choose the logarithm of the exit velocity, to underline the sensitivity on small variations in the friction factor value.

Figure 9
figure 9

Exit velocity as a function of the friction factor (on a logarithmic scale).

5.6 Example 4: external influences

Finally we model the presence of strong winds on the chimney top that may interfere with the hot gas flowing out of the structure.

We assume that the pressure variation at the top of the chimney due to wind refers mainly to horizontal wind directions. In case of a vertical wind component there would be an influence on the direction of the vertical part of the velocity but, thanks to rotational symmetry, no influence on the direction of the horizontal part is expected.

So we simulate the outside wind by increasing the external pressure.

The results are shown in Figure 10 where we draw the gas temperature as a function of the fire strength (represented by the heat production coefficient) over the added external pressure. In the dark blue area in the lower right part the outside pressure is too high for the given fire strength and no updraft in the chimney is possible. In the upper left part the chimney is active. There, a higher fire strength induces a higher gas temperature. A higher gas temperature induces more buoyancy. Therefore the chimney works even for stronger winds (i.e. higher outside pressures).

Figure 10
figure 10

Exit temperature as a function of the external pressure on top of the chimney. In the dark blue area in the lower right part no updraft in the chimney is possible. In the upper left part the chimney is active.

Another outside parameter is the external temperature. As we are working with the assumption of constant pressure, when modifying the external temperature we modify the reference values T r and ρ r (see Table 1). Consequently we modify the heat production coefficient c ˜ q (see definition of c ˜ q in Section 3). Thus, increasing the external temperature decreases the air velocity (Figure 11). This result agrees with empirical experience.

Figure 11
figure 11

Exit velocity as a function of the external temperature for three different heat loss parameters.

6 Conclusions

This paper proposes a small Mach number model to describe the gas dynamics in a chimney. It is shown that in this small Mach number application with strong heat sources the Boussinesq approach is not adequate. However, other small Mach number asymptotic models are shown to be a very good and fast alternative to a fully compressible model. A numerical scheme is proposed and completed with numerical simulations. Being a compromise between highly sophisticated CFD simulations and basic empirical formulas this approach is appropriate to perform parameter studies or to apply parameter studies (e.g. diameter, height, heat loss etc.) or to apply control procedures.

References

  1. McDermott H: Handbook of Ventilation for Contaminant Control. 1985.

    Google Scholar 

  2. [http://www.eca-europe.org].

  3. [http://www.nace.org.uk].

  4. Boussinesq J 2. In Théorie analytique de la chaleur: mise en harmonie avec la thermodynamique et avec la théorie mécanique de la lumière. Gauthier-Villars, Paris; 1903.

    Google Scholar 

  5. Bois PA: Asymptotic aspects of the Boussinesq approximation for gases and liquids. Geophys Astrophys Fluid Dyn 1991,58(1–4):45–55. 10.1080/03091929108227330

    Article  Google Scholar 

  6. Zeytounian RK: Joseph Boussinesq and his approximation: a contemporary view. C R, Méc 2003,331(8):575–586.

    Article  MATH  Google Scholar 

  7. Principe J, Codina R: Mathematical models for thermally coupled low speed flows. Adv Theor Appl Mech 2009, 2: 93–112.

    MATH  Google Scholar 

  8. Rehm RG, Baum HR: The equations of motion for thermally driven, buoyant flows. J Res Natl Bur Stand 1978,83(3):297–308. 10.6028/jres.083.019

    Article  MATH  Google Scholar 

  9. Paolucci S: Filtering of sound from the Navier-Stokes equations. NASA STI/Recon Technical Report N 1982, 83:26036.

  10. Majda A, Sethian J: The derivation and numerical solution of the equations for zero Mach number combustion. Combust Sci Technol 1985,42(3–4):185–205. 10.1080/00102208508960376

    Article  Google Scholar 

  11. Fedorchenko AT: A model of unsteady subsonic flow with acoustics excluded. J Fluid Mech 1997,334(1):135–155.

    Article  MATH  MathSciNet  Google Scholar 

  12. Métivier G, Schochet S: The incompressible limit of the non-isentropic Euler equations. Arch Ration Mech Anal 2001, 158: 61–90. 10.1007/PL00004241

    Article  MATH  MathSciNet  Google Scholar 

  13. Kuchling H: Taschenbuch der Physik. Fachbuchverl, Leipzig; 1994.

    Google Scholar 

  14. Majda A 53. In Compressible Fluid Flow and Systems of Conservation Laws in Several Space Variables. Springer, Berlin; 1984.

    Chapter  Google Scholar 

  15. [http://www.klb-klimaleichtblock.de].

  16. [http://www.schiedel.de].

  17. Heskestad G: Fire plumes, flame height and air entrainment. In Sfpe Handbook of Fire Protection Engineering. National Fire Protection Association, Quincy; 2002.

    Google Scholar 

  18. Hu LH, Li YZ, Wang HB, Huo R: An empirical equation to predict the growth coefficient of burning rate of wood cribs in a linear growth model. J Fire Sci 2006.,24(2): Article ID 153 Article ID 153

  19. Mizuno T, Kawagoe K: Burning behaviour of upholstered chair: Part1. Fire Sci Technol 1984,4(1):37–45. 10.3210/fst.4.37

    Article  Google Scholar 

  20. Styahar A: Modelling and simulation of chimney using one-dimensional solar updraft tower model. Master thesis. Department of Mathematics, University of Hamburg, 2010.

  21. Lions PL Oxford Lecture Series in Mathematics and Its Applications 10. Mathematical Topics in Fluid Mechanics, Volume 2, Compressible Models 1998. reprinted 2007 reprinted 2007

    Google Scholar 

  22. Schochet S: The compressible Euler equations in a bounded domain: existence of solutions and the incompressible limit. Commun Math Phys 1986, 104: 61–90.

    Article  MathSciNet  Google Scholar 

  23. Birken P, Meister A: Stability of preconditioned finite volume schemes at low Mach numbers. BIT Numer Math 2005,45(3):463–480. 10.1007/s10543-005-0009-0

    Article  MATH  MathSciNet  Google Scholar 

  24. Klein R: Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: one-dimesional flow. J Comput Phys 1995,121(2):213–237. 10.1016/S0021-9991(95)90034-9

    Article  MATH  MathSciNet  Google Scholar 

  25. Meister A: Asymptotic single and multiple scale expansions in the low Mach number limit. SIAM J Appl Math 1999,60(1):256–271. 10.1137/S0036139998343198

    Article  MATH  MathSciNet  Google Scholar 

  26. Müller B: Low mach number asymptotics of the Navier Stokes equations. J Eng Math 1998, 34: 97–109. 10.1023/A:1004349817404

    Article  MATH  Google Scholar 

  27. Gasser I: Modelling and simulation of a solar updraft tower. Kinet Relat Models 2009,2(1):191–204.

    Article  MATH  MathSciNet  Google Scholar 

  28. Bauer M, Gasser I: Modeling, asymptotic analysis, and simulation of an energy tower. SIAM J Appl Math 2012,72(1):362–381. 10.1137/100797977

    Article  MATH  MathSciNet  Google Scholar 

  29. Gasser I, Rybicki M: Modeling and simulation of gas dynamics in an exhaust pipe. Appl Math Model 2013,73(5):2747–2764.

    Article  MathSciNet  Google Scholar 

  30. Gasser I, Struckmeier J: An asymptotic-induced one-dimensional model to describe fires in tunnels. Math Methods Appl Sci 2002,25(14):1231–1249. 10.1002/mma.337

    Article  MATH  MathSciNet  Google Scholar 

  31. Batchelor GK: The conditions for dynamical similarity of motions of a frictionless perfect-gas atmosphere. Q J R Meteorol Soc 1953,79(340):224–235. 10.1002/qj.49707934004

    Article  Google Scholar 

Download references

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Ingenuin Gasser.

Authors’ original submitted files for images

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Felaco, E., Gasser, I. Modelling, asymptotic analysis and simulation of the gas dynamics in a chimney. J.Math.Industry 3, 3 (2013). https://doi.org/10.1186/2190-5983-3-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/2190-5983-3-3

Keywords