On the use and misuse of the Oberbeck-Boussinesq approximation

The Oberbeck-Boussinesq approximation is the most widely employed theoretical scheme for the study of natural or mixed convection flows. However, the misunderstanding of this approximated framework is a possibility that may cause the emergence of paradoxes or, at least, incorrect conclusions. In this note, the basic features of the Oberbeck-Boussinesq approximation are briefly recalled and three simple examples where this theoretical scheme may be misused are provided. Such misuses of the approximation lead to erroneous conclusions that, in the examples presented in this note, entail violations of the principle of mass conservation. A discussion about the Oberbeck-Boussinesq approximation as an asymptotic theory obtained by letting the product of the thermal expansion coefficient and the reference temperature difference tend to zero is also presented.


Introduction
The study of natural or mixed convection flows either in fluids or in fluid-saturated porous media is, with a limited number of exceptions, modelled theoretically by claiming the validity of the Oberbeck-Boussinesq approximation.There are several thorough and comprehensive analyses of how this approximation can be established starting from a general formulation of the local balance equations of mass, momentum and energy for a fluid.Beyond the many textbooks of fluid dynamics and convection heat transfer, we mention the analyses of this topic presented in chapter 8 of Zeytounian [1], in Rajagopal et al. [2] and in Zeytounian [3].Such discussions on the origin and on the range of applicability of the Oberbeck-Boussinesq approximation stem from pioneering papers such as Spiegel and Veronis [4], Gray and Giorgini [5] and Hills and Roberts [6].A recent interesting review on the topic has been presented by Mayeli and Sheard [7].
The aim of this short paper is to highlight an aspect of the Oberbeck-Boussinesq approximation that may be the source of pitfalls, i.e., the duality of variable fluid density and constant fluid density.If it is recognised that the approximation scheme predicates a variable density which is pressure independent and linearly varying with the temperature, in some instances one may forget that such variable density serves only to define the buoyancy force within the local momentum balance equation.Utilising the variable density outside this very specific context may lead to unphysical predictions and, hence, to incorrect conclusions.Such conclusions are incorrect as they usually conflict with the principle of mass conservation.A final discussion about the interpretation of the Oberbeck-Boussinesq approximation as a limiting case of the general local balance equations of mass, momentum and energy is presented.

A minimalistic survey of the Oberbeck-Boussinesq model
There are several detailed and thorough descriptions of the Oberbeck-Boussinesq model for natural and mixed convection flows.The basic system of partial differential equations, expressing the local mass, momentum and energy balances, is given by where u is the velocity, P is the pressure, T is the temperature, t is time, µ is the dynamic viscosity and α is the thermal diffusivity.There are two densities in Eqs.(1).One is the reference density ρ 0 , i.e., the fluid density evaluated at the reference temperature T 0 .On the other hand, ρ denotes the fluid density evaluated at the local temperature T through the linear equation of state, where β is the thermal expansion coefficient.We mention that the linear equation of state (2) is to be replaced by a quadratic equation of state in special situations such as pure water close to 4 • C. It must be stressed that, in Eqs. ( 1) and ( 2), (ρ 0 , µ, α, β) are constant, pressure-independent, fluid properties evaluated at the reference temperature T 0 .In order to maximise the reliability of the approximation, a judicious choice of the reference temperature is an average temperature over the flow spatial domain and over the time interval of the flow process.We also implicitly assume, for the local energy balance (1c), that the viscous heating effect is negligible and that no internal heat source is present.
When the Oberbeck-Boussinesq scheme is to be applied to the seepage flow in a porous medium, then the local momentum balance equation may be modelled through Darcy's law, namely where u now denotes the seepage, or Darcy's, velocity and K is the permeability of the medium.Then, Eq. (3) supersedes Eq. (1b), in this case.Also Eq. (1c) is to be reformulated when we study the seepage flow in a porous medium.In fact, the local energy balance reads where σ is the ratio between the average heat capacity of the saturated porous medium and the heat capacity of the fluid, while α m is the average thermal diffusivity of the saturated porous medium.

Hydrostatic pressure, buoyancy force
The gravitational force term, ρg in either Eq.(1b) or in Eq. ( 3), is the one and only quantity where the use of Eq. ( 2) is allowed.Thus, we can write where is the hydrostatic pressure, with x = (x, y, z) denoting the position vector, and is the buoyancy force.Both P and P h depend on the spatial position.It is a common practice to use a specific notation for their local difference, We call p hydrodynamic pressure, for the sake of brevity, to mark the distinction from the pressure and the hydrostatic pressure.In Eq. ( 6), P e is an arbitrary constant value which, physically, can be assigned in order to fix a possible overall external pressurisation of the fluid.In most usual cases, P e is assumed to coincide with the atmospheric pressure.
Out of the present context, where the hydrostatic pressure and the buoyancy force are defined on the basis of a chosen reference temperature, T 0 , the fluid density is considered to be a constant equal to ρ 0 .If one violates this simple rule, then erroneous conclusions could be drawn which will typically lead to violations of the mass conservation principle.

A rectangular cavity with side heating
Let us consider the classical problem of two-dimensional natural convection in a rectangular cavity with rigid and impermeable side boundaries kept at uniform, but different, temperatures, T 1 and T 2 , while the upper and lower sides are kept adiabatic.It is not restrictive to assume T 1 > T 2 .A sketch of the system is provided in Fig. 1.
It is well-known that, with sufficiently small differences T 1 − T 2 , a steady-state natural convection flow occurs in the cavity.In non-dimensional terms, this restriction is equivalent to a sufficiently small Rayleigh number.The steady-state flow is cellular in character, with one or more convective cells.For a given fluid, the number of stationary convective cells depends on the aspect as an obvious consequence of Eq. (1a).In fact, one can just integrate Eq. (1a) over the upper half-domain {x ∈ [−L, L], y ∈ [0, H]} and then employ Gauss' theorem by recalling that S is the only permeable boundary of such a domain.Here, v is the y component of u.Then, we could wonder how we could evaluate the mass flow rate, ṁ, per unit depth (in the z direction) across S.
There is a correct way and an incorrect way.The correct way is by applying the principle that the fluid density is to be considered constant and equal to ρ 0 , so that Eq. ( 9) yields The incorrect way is by employing the variable density expressed by Eq. ( 2), so that Eq. ( 9) The reason why Eq. ( 11) yields ṁ < 0 is that, at least with a sufficiently small temperature difference T 1 − T 2 , both v and T − T 0 are odd functions of x with x ∈ [−L, L], positive for x ∈ (−L, 0) and negative for x ∈ (0, L) (see, for instance, de Vahl Davis [8]).Thus, their product is an even and positive function of x throughout the domain of integration x ∈ [−L, L], so that Eq. ( 11) leads to the conclusion that ṁ < 0. The numerical solution discussed by de Vahl Davis [8] is complemented by a simple analytical solution which predicts the same symmetry for v and T − T 0 and which holds for a very tall cavity, H L. Such an asymptotic solution is briefly outlined in Appendix A. We mention that Eq. ( 9) holds also for any other y = constant plane S, as a consequence of Eq. (1a).We said that Eq. ( 11) expresses the incorrect way to evaluate ṁ since the conclusion ṁ = 0 is an evident violation of the mass conservation within the upper half- In fact, both for the upper and the lower half-domains, S would be the only permeable boundary and it would be crossed by a net mass flow rate.Such a situation, in a stationary regime, yields a violation of the principle of mass conservation.

Mixed convection duct flow
Let us consider the internal mixed convection in a duct with an impermeable wall having an increasing temperature along the streamwise direction.Such a behaviour is observed, for instance, when the duct wall is subject to an incoming uniform heat flux.We consider the case of stationary flow.
As sketched in Fig. 2, we consider the region V delimited by the cross-sections S 1 and S 2 .Since, the fluid is heated in the streamwise direction, we have an average temperature at the cross-section S 1 , denoted by T 1 , smaller than the average temperature T 2 evaluated at the cross-section S 2 .Then, a judicious choice of the reference temperature T 0 for the Oberbeck-Boussinesq approximation in the domain V is the volume-averaged temperature, Such a volume-averaged temperature value, T 0 , is larger than T 1 and smaller than T 2 .Let us evaluate the average velocities across S 1 and S 2 , where n1 and n2 are the unit vectors of the surfaces S 1 and S 2 , oriented in the streamwise direction as shown in Fig. 2. Equation (1a), after an integration over V and use of Gauss' theorem yields the equality S 1 u m1 = S 2 u m2 .One can wonder how one can evaluate the mass flow rates across S 1 and S 2 , i.e., ṁ1 and ṁ2 , respectively.The right way is by assuming that the fluid density is to be considered constant and equal to ρ 0 all over V , so that Since S 1 u m1 = S 2 u m2 , Eqs. ( 13) and ( 14) allow one to conclude that ṁ1 = ṁ2 , which is in perfect agreement with the principle of mass conservation.
The incorrect way to evaluate ṁ1 and ṁ2 is by employing the variable density ρ, Eq. ( 2).In this case, we have From Eq. ( 15) and from the equality S 1 u m1 = S 2 u m2 , one can write By employing again the equality S 1 u m1 = S 2 u m2 and by recalling that the average temperature over S 1 is smaller than the average temperature over S 2 , one can immediately conclude that ṁ1 − ṁ2 > 0, i.e., the mass rate flowing through S 1 is larger than the mass rate flowing through S 2 which, in a stationary regime, means a violation of the principle of mass conservation.

A vertical porous slab separating two fluid reservoirs
Let us consider an infinitely wide wall separating two fluid reservoirs kept at different temperatures, T 1 and T 2 .As shown in Fig. 3, the wall has a porous insertion bounded by two planes S 1 and S 2 .
The fluid in the left-hand reservoir is the same as that in the right-hand reservoir and they are both in a rest state.The same fluid saturates the porous slab.The Oberbeck-Boussinesq approximation can be applied assuming the reference temperature T 0 = (T 1 + T 2 )/2.One can evaluate the pressure distribution on the boundaries S 1 and S 2 as the hydrostatic pressure.By employing Eq. ( 6), the pressure distribution on the plane S 1 is given by while, on S 2 , we have Equations ( 17) and ( 18) have been used to formulate the pressure conditions at the porous slab boundaries by Barletta [9,10], Barletta and Celli [11] and by Barletta and Rees [12].The pressurisation constants, P e1 and P e2 , can be equal or not depending on the conditions externally imposed on the two reservoirs.For instance, one may have a situation where the reservoir at temperature T 1 is compressed on its top, while the reservoir at temperature T 2 is open to the atmospheric pressure.
In such a situation, one has P e1 > P e2 .Equations ( 17) and (18) reveal that the pressure distribution on S 1 and S 2 is y dependent, but the pressure difference across the porous slab is a constant, If none of the two reservoirs is pressurised, then ∆P = 0.If ∆P = 0, one may have a seepage throughflow across the porous slab [10].
As sketched in Fig. 3 the porous insertion is bounded above and below by an impermeable material, so that application of Eq. (1a) yields where u is the x component of the seepage velocity, u, across the porous material.By recalling that the fluid density is to be intended as uniform with the value ρ 0 , Eq. ( 20) implies that the mass flow rates across S 1 and S 2 are equal as required by mass conservation, This is the correct and consistent application of the Oberbeck-Boussinesq scheme to this sample case.
There is always the possibility of introducing errors as in the previous examples.One may enforce the validity of Eq. ( 2) in the evaluation of the hydrostatic pressures in the two reservoirs so that Eqs. ( 17) and ( 18) are replaced by and Equations ( 22) and ( 23) yield an y dependent pressure difference across the porous layer Even in the absence of pressurisation in one of two reservoirs (P e1 = P e2 ), Eq. ( 24) entails a pressure difference across the porous layer which, in turn, leads to the prediction of a horizontal throughflow across the porous layer [13].Such a stationary throughflow is unphysical as it would lead to a violation of mass conservation.In fact, consistently with the assumption of a variable density, as given by Eq. ( 2), the mass flow rate across S 1 is expressed as while the mass flow rate across S 2 is given by The violation of mass conservation is quantified by a relative error where Eqs. ( 20), ( 25) and ( 26) have been employed.Interestingly enough, should δ be considered as negligible, δ 1, then the pressure distributions given by Eqs. ( 22) and ( 23), and erroneously employed in the study carried out by Vynnycky and Mitchell [13], would turn out to match perfectly the correct pressure distributions expressed by Eqs. ( 17) and (18).More precisely, the corrective terms introduced in Eqs. ( 22) and (23), i.e. ∓ δ/2, can be taken into account consistently only if the solenoidal constraint for the velocity, ∇ • u = 0, is replaced by the variable-density local mass balance equation, ∇ • (ρ u) = 0.The intermediate approximation used by Vynnycky and Mitchell [13], where Eqs. ( 22) and (23) are employed in combination with ∇ • u = 0, is flawed as it leads to a violation of the principle of mass conservation.

The Froude number and the Rayleigh number
Let us reconsider Eqs. ( 5)- (7).The pressure gradient and gravitational force contributions to the local momentum balance are given by the force per unit volume where we are assuming that the y axis is vertical and orientated upward.We can rescale F in order to obtain a dimensionless formulation of the local momentum balance equation.The dimensionless F denoted with an asterisk is given by where the constant ∆P r is a reference pressure difference and the constant L is a reference length.One can also define the dimensionless coordinates pressure and temperature as where ∆T r is a reference temperature difference.On account of Eqs. ( 29) and (30), Eq. ( 28) yields where the Froude number is defined as and δ is given by It is to be mentioned that the usual definition is Fr = U 2 0 /(g L) (see, for instance, Mayeli and Sheard [7]), where U 0 is a reference velocity.In fact, Eq. (32) matches perfectly the usual definition provided that one chooses U 0 = ∆P r /ρ 0 .It is also to be mentioned that, according to other authors (e.g.Zeytounian [1]), Fr should be defined as a velocity ratio, so that it is given by the square root of the quantity identified by Mayeli and Sheard [7] as the Froude number.In our discussion, we will rely on the definition given by Eq. (32).
In the case of natural convection, where buoyancy alone causes the flow, a typical choice of ∆P r is Thus, the ratio δ/Fr coincides with the Rayleigh number, In the discussions available in the literature about the rigorous derivation of the Oberbeck-Boussinesq approximation, the approximate governing equations (1) are considered as a limiting case of the local balance equations for a fully-compressible (variable density) flow when This double limit is stated, although in slightly different terms, by Hills and Roberts [6] and reported by Zeytounian [3,14].It is also implicitly employed in the analysis carried out by Rajagopal et al. [2].Also the recent study by Vynnycky and Mitchell [13] relies on this limiting scheme as a basis for the Oberbeck-Boussinesq approximation.However, this approach leads to a singular behaviour of F * , as one can easily infer from Eq. ( 31), since the term −(1/Fr) êy blows up in the limit given by Eq. ( 36).The singular nature of the Oberbeck-Boussinesq limit of F * is unavoidable, even if it could have been concealed should one have employed, on account of Eqs. ( 6) and ( 8), the gradient of the dynamic pressure, ∇p, instead of the pressure gradient, ∇P , on writing the expression of F in Eq. ( 28).However, this trick would have simply swept the dust under the carpet without actually solving the problem.In fact, there is just one way to avoid a singular limiting bahaviour of F * in the limit defined by Eq. (36).One should constrain ∇P to be equal to −ρ 0 g êy whenever the Oberbeck-Boussinesq approximation is used.Unfortunately, such a constraint is unphysical except for a quite limited number of special cases.Interestingly enough, the derivation of the Oberbeck-Boussinesq approximation presented by Rajagopal et al. [2] leads to the same unphysical conclusion drawn above, namely that the pressure gradient must always coincide with the hydrostatic pressure gradient.We mention that the treatment presented by Rajagopal et al. [2] is based on a dimensionless scaling of the governing balance equations where the reference length is O −1 and the reference velocity is O( ), where is the perturbation parameter.The Oberbeck-Boussinesq approximation is defined as the asymptotic case where → 0. Thus, the dimensionless scaling is singular in this limit.Incidentally, is proportional to δ 1/3 , with δ given by Eq. (33).The paper by Rajagopal et al. [2] discusses in detail the serious drawbacks of the previous theoretical studies that define the Oberbeck-Boussinesq approximation as a limiting case obtained by letting one or more perturbation parameters to zero.Examples are the papers by Spiegel and Veronis [4], by Gray and Giorgini [5] and by Hills and Roberts [6].As a consequence, one can say that an asymptotic theory, based on a suitable perturbation scheme, which is aimed at a rigorous deduction of the Oberbeck-Boussinesq set of governing equations ( 1) is still lacking.It is also possible that the existing approaches that define the Oberbeck-Boussinesq approximation as a limiting case of the fully compressible set of local balance equations are intrinsically biased.It is the authors' opinion that, until a rigorous and nonsingular theoretical scheme will be set up to justify the approximation as an asymptotic regime, its validity relies entirely on the widely-documented experimental validations available in the literature for a very broad range of flow regimes.

Conclusions
The Oberbeck-Boussinesq approximation for the local balance equations of mass, momentum and energy of buoyancy-induced fluid flows has been briefly outlined.It has been stressed that the linear temperature-dependent density expression is to be used, within the local momentum balance equation, only to transform the combined pressure gradient force and gravity force into a combination of terms involving the dynamic pressure gradient and the buoyancy force.Out of this very specific use, the variable density expression cannot be employed consistently with the approximated Oberbeck-Boussinesq scheme.Three examples have been discussed, relative to either natural or mixed convection flows, where it has been shown that an improper use of the variable density unavoidably leads to violations of the principle of mass conservation.The third of these examples points to a misuse of the Oberbeck-Boussinesq approximation presented in a recent paper [13].A final section has been devoted to the idea of the Oberbeck-Boussinesq approximation as a limiting case of the fully-compressible local balance equations where the dimensionless parameter δ = β ∆T r tends to zero, with β the thermal expansion coefficient and ∆T r the reference temperature difference.It has been pointed out that such an idea leads either to a singular behaviour for the gradient of the dynamic pressure, or to the unphysical constraint that such gradient be identically zero.We conclude that, so far, the Oberbeck-Boussinesq approximation is to be classified as a phenomenologically based model.

Figure 1 :
Figure 1: Sketch of the rectangular cavity with side heating.

Figure 3 :
Figure 3: Sketch of of the porous slab separating two fluid reservoirs.