On the relevance of oscillatory convection regimes in porous media — review and numerical experiments

Computers & Fluids 30 (2001) 189±209

www.elsevier.com/locate/compfluid

On the relevance of oscillatory convection regimes in porous media Ð review and numerical experiments Ekkehard Holzbecher Department of Ecohydrology, IGB-Institute of Freshwater Ecology and Inland Fisheries, Rudower Chaussee 6A, Building 21.2, 12484 Berlin, Germany Received 12 July 1999; received in revised form 3 April 2000; accepted 3 April 2000

Abstract Using spectral and bifurcation methods, several studies have found the critical margin for the onset oscillatory convection in porous medium at a Rayleigh number 1390. The threshold is derived for ®rst mode 2D rolls, but can hardly be observed in experimental setups. This paper shows that aside from experimental diculties the reason for this phenomenon is that for supercritical Rayleigh numbers in the given range (100±1000) the occurrence of ®rst mode square rolls is unlikely. Steady convection in second and third mode is preferred to oscillatory convection in the ®rst mode. Using a branch tracing technique, it is shown that for higher modes the transition to oscillatory regime occurs at higher Rayleigh numbers (1510 for the second mode, 1970 for the third mode). For Rayleigh numbers 11000 oscillatory convection in the third mode is preferred to stable convective regimes. Results of the paper are obtained by numerical experiments. Nusselt numbers representing mean heat transfer through the system are taken as indicators for oscillatory convection. Model design, set up using the FAST-C(2D) code, is described. FAST-C(2D) code, developed by the author, is introduced and its derivation from basic principles is brie¯y outlined. Results are discussed with respect to the relevance or irrelevance of oscillatory convection. 7 2000 Elsevier Science Ltd. All rights reserved.

1. Introduction Convective motions have been the subject of scienti®c research since the experiments of BeÂnard at the start of the century. BeÂnard [1] heated a box ®lled with ¯uid from below and E-mail address: [email protected] (E. Holzbecher). 0045-7930/01/$ - see front matter 7 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 0 0 ) 0 0 0 0 8 - 6

190

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

Nomenclature Coux , Couz Courant numbers D thermal di€usivity g gravity vector with length equal to acceleration due to gravity pointing to the earth center H system and model height eddy height Heddy k tensor of permeabilities k scalar permeability permeabilities in principal directions (x: horizontal; z: vertical) kx , kz L system and model length eddy length (horizontal extension) Leddy m mode number Neu Neumann number Nu Nusselt number Pe grid PeÂclet number Ra Rayleigh number critical Rayleigh number Racrit critical Rayleigh number for transition from conduction to steady convection Racrit1 critical Rayleigh number for transition from steady to oscillatory convection Racrit2 t time T temperature v velocity vector components of velocity vector in direction of both coordinate axes vx , vz x, z coordinate axes (x: horizontal; z: vertical) w automatic timestepping parameter Dt timestep Dx, Dz grid spacing in both coordinate axes e accuracy required by iterative solution algorithm g ratio of speci®c heat for ¯uid and porous medium j porosity m dynamical viscosity r density W convection cell aspect ratio cold water reference density r0 Dr density di€erence between hot and cold water C streamfunction y normalized temperature or salinity o vorticity

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

191

discovered what we today call convection cells. Moreover, BeÂnard found that for low heater temperatures the system remained static, but ¯ow appeared when temperatures at the heated surface were increased. It is now well understood that convective motions emerge when a certain critical value of a parameter is exceeded. It is common to use a dimensionless combination of parameters, the so called Rayleigh-number Ra. Lord Rayleigh [2] was the ®rst who provided a rigorous analytical treatment of BeÂnard convection by studying the exchange of stability in a pure ¯uid experiment. The critical Rayleigh number Racrit1 marks the transition from the conductive no-¯ow state to steady convective ¯ow patterns. For ¯ow in porous media the same is true: there is a critical margin where conductive and convective regimes exchange stability. Despite the fact that Lord Rayleigh treated the pure ¯uid case only, the dimensionless combination of parameters that determines the ¯ow pattern in scienti®c literature is referred to as a Rayleigh number. The de®nition of the `®ltration Rayleigh number' in porous media is given below in Eq. (9). Lapwood [3] was the ®rst to derive the critical margin Racrit1 ˆ 4p 2 for the onset of convection in porous media. Convective motions are relevant in various situations. The ¯ow pattern is most important in studies on heat transfer because the transition from conductive regime to convective regime is accompanied by a signi®cant increase of heat transfer through the system. The Nusselt number Nu is the dimensionless parameter which characterises mean heat transfer through the system from the heated to the cooled surface. For an understanding of the heat transfer the dependency of the Nusselt number on the Rayleigh number is most important and is usually given in a Ra±Nu diagram. The Ra±Nu curve is not unique. Nusselt number and heat transfer depend additionally on the convection mode. The term mode stems from the spectral approach used in the classical treatment of the problem for which functions of the type exp…imx=p†exp…iz=p†

…1†

are studied in order to ®nd solutions for the system of di€erential equations Ð technique applied already by Lord Rayleigh [2] and by Lapwood [3]. For porous media it is found that the amplitude responsible for the onset of convection at minimum Rayleigh number Racrit1 ˆ 4p 2 is given for m ˆ 1 (1st mode). The next possible solution in the unit square is given for m ˆ 2 (2nd mode). This notation is adopted from the classical papers and used in the remainder of this paper. The mode is characterised by the convection cell aspect ratio W ˆ Leddy =Heddy , which should not be mistaken with system or model aspect ratio L/H. For a given system aspect ratio, only certain cell aspect ratios are possible. In a developed convection pattern length L and height H are multiples of cell length Leddy and cell height Heddy : The aspect ratio for the ®rst mode is W ˆ 1; for the second mode W ˆ 1=2; and for the third mode W ˆ 1=3: With a decreasing aspect ratio eddies become more slender. Di€erent modes represent di€erent branches of solutions for a nonlinear set of di€erential equations, which describes density driven ¯ow in porous media. When a branch is traced in the direction of increasing Rayleigh numbers, a transition to oscillatory convection may be observed. There is a second critical Rayleigh number Racrit2 , at which the steady convective regime becomes unstable and changes to a ¯uctuating convective regime. The phenomenon has been studied using various numerical techniques and a short review is provided in the following.

192

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

Most studies on oscillatory convection combine analytical and numerical methods. Caltagirone [4] studied the relation between Ra, Nu and the cell aspect ratio W using numerical simulations. The Galerkin technique led him to a critical Rayleigh number of Racrit2 ˆ 38425 for the square convection cell. Straus and Schubert [5] used the spectral method to treat 2D and 3D con®gurations and obtained a range for the change to unsteady motions for Ra between 300 and 320. Later the same authors corrected that result using an improved numerical technique on advanced hardware with an increased number of terms in the approximating functions. Straus and Schubert [6] studied Rayleigh numbers up to 650 in a square 2D-cavity and found several margins with transitions between di€erent distinct unsteady states. According to the improved results, the onset of oscillatory convection in the unit cell is given for Ra between 380 and 400. At a (third) critical value Racrit3 between 480 and 500, there is a switch to a two frequency oscillatory quasi-periodic mode. If Ra is increased further (500±520), a single frequency becomes dominant again Ð having a di€erent frequency than all frequencies observed before. Caltagirone and Fabrie [7] observed the second single frequency mode for Rayleigh numbers up to 600. In contrast, calculations by Kimura et al. [8] showed a switch to a second single frequency periodic regime at Racrit4 between 560 and 570. Moreover, Kimura et al. [8] found the transition to a chaotic regime for Rayleigh numbers between 850 and 1000 and wrote: ``routes to chaos in porous-medium convection are fundamentally di€erent than those in ordinary ¯uids''. They point to the fact that the step from a 2-frequency regime back to a single frequency regime (with increasing Ra ) jumps from higher disorder to lower disorder, while in pure ¯uid convection there is a straight line with increasing disorder from conduction to chaos. Caltagirone and Fabrie [7] traced the dynamic regimes of unicellular convection in a square system up to Ra ˆ 1500 and basically con®rmed the previous ®ndings. They identi®ed the emergence of a second, quasi-periodic pattern for Rayleigh numbers above 1000 …Racrit5 † on the route to chaos. The references cited so far used the spectral approach in combination with the Galerkin technique. The set of amplitude equations Ð a set of ordinary di€erential equations Ð is treated as an initial value problem. There are some di€erences in the numerical treatment: Kimura et al. [8] and Caltagirone and Fabrie [7] applied a modi®ed technique, the `pseudospectral' method, with higher computational eciency than spectral methods used in former studies. Other authors used techniques from bifurcation theory to treat the set of ordinary di€erential equations. Stable and oscillatory convection are di€erent types of solutions, which can be represented by a point and a circle in the phase space. The Hopf bifurcation marks the transition from a situation with an equilibrium point to a qualitatively di€erent situation with a limit cycle. Bifurcations can be characterized by the eigenvalues of the linear system: it is an indicator for the Hopf bifurcation when a pair of imaginary eigenvalues crosses the imaginary axis. Aidun and Steen [9] used an eigenvalue expansion and a branch-tracing technique to treat the bifurcation problem. They obtained Racrit2 ˆ 390:7 for 2D problems. Aidun [10] analysed the 3D case and was able to con®rm the previous result under the condition that the extension in the third dimension is small enough. Using bifurcation techniques, Riley and Winters [11] obtained a value of 390.72 for the second critical Rayleigh number.

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

193

Oscillatory motions have hardly been observed in laboratory experiments with a porous medium. Caltagirone [4] reported experiments in a horizontal layer with rolls large in horizontal size. The results of these experiments do not coincide with numerical ®ndings. In a 2D experiment, local instabilities were found to disappear in the main vortex. Schneider [12], Masuoka [13], Buretta and Berman [14], Jonsson and Cattan [15] performed experiments with Rayleigh numbers in the given range without mentioning oscillatory convection. All these experimental studies observed stable convection for Rayleigh numbers above 390. In a cylindrical setup, Bau and Torrance [16] found the critical margin for oscillatory convection for Ra=Racrit1 15:5:

2. Theory and numerical implementation 2.1. Principles and di€erential equations Convective motions result from an interaction of ¯ow and transport processes. Fluid ¯ow in¯uences transport in any non-trivial case. There is an in¯uence from a transported component on the ¯ow when ¯uid characteristics depend on the component. Components with an in¯uence on density and viscosity of the ¯uid are heat and salt. In this paper the focus is on the thermal case only, but results are valid for the saline case as well in which convective motions are induced by salt concentration gradients. Double di€usive regimes with a combined in¯uence of salt concentration and heat show a di€erent behaviour and are not touched by the results of this study. In the following, di€erential equations are derived from basic principles such as mass and energy conservation. Well-established linear relations are taken into account: Darcy's law for ¯uid ¯ow in porous media and Fick's law for thermal di€usion. The set of di€erential equations describes the main processes in porous media convection under certain simplifying conditions. In the absence of sources or sinks, ¯uid ¯ow is governed by the principle of mass conservation and Darcy's law: @ …jr † ˆ ÿr
…rv † @t k v ˆ ÿ …rp ÿ rg† m

…2a,b†

The Oberbeck±Boussinesq assumption leads to a substantial simpli®cation of the di€erential equations. The most simple formulation of that assumption is that changes in density can be neglected except from the buoyancy term rg in Darcy's law (2b). Then Eq. (2a,b) reduce to the system: @j ˆr
v @t

k v ˆ ÿ …rp ÿ rg† m

…3a,b†

194

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

The assumption that porosity does not change in time is mostly valid and thus holds: r
v ˆ 0: Here the streamfunction is used as main ¯ow variable in 2D cross-sections. Under the condition of vanishing divergence r
v ˆ 0 the streamfunction C in a (x, z )-coordinate system is given by the two de®ning equations: @C ˆ vz @x

@C ˆ ÿvx @z

…4a,b†

For the 2D-case using Eq. (4) and Darcy's law (2b), the following di€erential equation evolves for the streamfunction:         @ m @C @ m @C @ m @ m @r vz ÿ vx ˆ ÿg ‡ ˆ …5† @x kz @x @z kx @z @x kz @z kx @x The mixed derivative terms of pressure vanish and a simple equation for the streamfunction remains, where the density on the right side stems from the buoyancy term. Eqs. (4) and (5) are a set of three di€erential equations describing the ¯ow of water with variable properties. It is completed by the explicit formulae for ¯uid viscosity m and ¯uid density r: Note that the entire set is derived with no assumptions about the speci®c parameter dependencies of r, m: In addition to mass conservation, the principle of energy conservation provides another equation which is necessary to describe convective ¯ow. The energy equation is understood as a di€erential equation with temperature T as the dependent variable. It is a transport equation describing advective and dispersive ¯uxes of heat. For the following, it is assumed that dispersion can be described by Fick's law with a scalar di€usion coecient D. A more detailed derivation is provided by Combarnous and Bories [17] or Holzbecher [18]. It is useful to move to dimensionless variables. The lengths are transformed using a typical length unit (here H ). Time is transformed with respect to di€usivity D: x4x=H z4z=H t4t
D=H 2 v4v
g
H=D C4C
g=D

…6†

In the following, the notation for the transformed variables will be the same as for the original variables, as it is mostly clear which of both alternatives is meant. A new notation is introduced only for the transport variable: y denotes normalized temperature: y ˆ …T ÿ Tmin †=…Tmax ÿ Tmin †: With the given transformations, ¯ow and transport equations come into the form:

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

    @ 1 @C @ 1 @C ggH @r @y ‡ ˆÿ @x kz @x @z kx @z Dm @y @x @y ˆ ÿv
ry ‡ r 2 y @t

195

…7a,b†

Some further simpli®cations and assumptions are generally made in order to derive the classical formulation with a characterising dimensionless number. In homogeneous and isotropic porous media, Eq. (7a) can be multiplied by the permeability k. Dynamic viscosity m is taken as constant in order to enable comparison with other publications in which the same assumption is made. Note that the assumption is not ful®lled when wide temperature ranges have to be considered. A linear state equation is assumed to hold: r ˆ r0 ÿ yDr

…8†

The di€erential equation for the streamfunction then becomes: @ 2C @ 2C @y ‡ ˆ ÿRa 2 2 @x @z @x

with Ra ˆ

kgDrgH Dm

…9†

The dimensionless constant Ra is the so called Rayleigh number. In his classical publication from 1916, Lord Rayleigh treated variable density ¯ow of pure ¯uid systems and derived a dimensionless number as a characteristic of the free ¯uid system, which as a parameter combination is not identical to the de®nition given in Eq. (9). In a similar manner, a set of equations can be derived for free ¯uid ¯ow. When instead of Darcy's law the Navier±Stokes equations describe ¯uid motions, another di€erential equation for vorticity needs to be taken into account (see Elder [19]). Note that for porous media ¯ow, the right-hand side of Eq. (9) is equal to vorticity (see Elder [20]). o ˆ Ra

@y @x

…10†

De Josselin de Jong [21] and Wooding [22] were probably the ®rst who used an equivalent formulation for transient convection in porous media as given above.

2.2. Software (FAST-C(2D)) FAST-C(2D) is designed for modelling density-driven ¯ow in porous media. Holzbecher [18] presents a detailed description and several application cases: heating from the side, geothermal ¯ow, saltwater intrusion, saltwater upconing, ¯ow above a saltdome. Not all features of the software are relevant for convection modelling. Only some basic features which are important for the oscillatory convection case are mentioned in the following. FAST-C(2D) works with the streamfunction C as the dependent ¯ow variable and temperature T as the transport variable. The code FAST-C(2D) is based on a block centered ®nite di€erence discretisation of Eqs. (4) and (7). Boundary conditions are to be given in terms of streamfunction and temperature. It is possible to specify Dirichlet (1st order) and Neumann

196

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

(2nd order) boundary conditions at single block edges at the outer boundaries of the model region. For temporal discretisation, the Crank±Nicolson procedure is recommended. Helpful for the modeller is the automatic timestepping option. With this option selected, timesteps are calculated during the simulation run to ful®l all generalised Neumann- and Courant-criteria with user-de®ned timestepping parameter w in the entire model area:   vx Dt vz Dt 1 1 DtRw …11† Rw and Couz ˆ Rw and Neu ˆ D ‡ Coux ˆ Dx Dz Dx 2 Dz 2 For wR1 the well-known stability criteria Coux R1, Couz R1, NeuR1 are ful®lled everywhere in the model. This may be quite restrictive when one of the criteria is violated only at a few locations in the model area. Thus, it is sometimes reasonable to use values w > 1: FAST-C(2D) is equipped with a variety of output options, numerical and graphical. For convective motions, the Nusselt number output is relevant. There are three approximations of Nusselt numbers which are calculated according to the following formula: Nu…Z † ˆ

1 X @y…Z † @z L

Z 2 f1, 2, 3g

…12†

with: … †

@y 1 2 …y1 ÿ y † ˆ Dz @z  …3 †    … † @y 2 1 1 @y 1 15 5 3      …y1 ÿ y † ÿ …y2 ÿ y † ‡ …y3 ÿ y † 3…y1 ÿ y † ÿ …y2 ÿ y † ˆ ˆ Dz 3 Dz 4 6 20 @z @z Here y denotes the pre-set value of y at the boundary and y1 , y2 , y3 are numerically derived values in the three blocks next to the boundary. In convection studies, the sum in Eq. (12) is taken over all blocks on a horizontal boundary. As the calculations are made for the upper and lower edges separately, it is possible to check the quality of numerically-computed steady states. The quadratic approximation is always greater than the linear, and the cubic lies between both lower order approximations. Very useful for the simulation runs of this study are the RESTART options implemented in FAST-C(2D). At the end of the simulation, the ®nal state can be stored in a `restart'-®le from which it can be taken to start another run with the same or modi®ed parameters. There are two RESTART options: one to read the restart-®le and one to write the restart-®le. Both options can be set in the same run, which is extremely useful in a branch-tracing simulation series (see below). It is possible with FAST-C(2D) software to tackle cases with nonlinear density (see Holzbecher [23]) and nonlinear viscosity (see Holzbecher [24]). But these options are not used in this study in order to obtain results which are comparable with former ®ndings obtained using the classical concept of a constant Rayleigh number. Applications of FAST-C(2D), among others, are published by Baumann and Moser [25], Holzbecher and Baumann [26], Holzbecher and Heinl [27], and Holzbecher and Yusa [28].

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

197

2.3. Model and experiment design Using FAST-C(2D) (see Holzbecher [18]) a model for the 2D rectangular system of unit height H is set up. The length L is changed depending on the mode to be traced (see below). At the top and bottom of the model Dirichlet (1st order) boundary conditions are prescribed: constant high temperature at the bottom …y ˆ 1† and constant low temperature at the top …y ˆ 0). Adiabatic conditions are assumed to be valid at the sidewalls that are achieved by the Neumann (2nd order) boundary condition @y=@x ˆ 0: For the streamfunction, the Dirichlet boundary condition C ˆ 0 is required at all boundaries of the modelled area, i.e., the boundary of the rectangular model area is a streamline. The discretization type for 1st order derivatives in the transport Eq. (7b) is CIS (central in space). The alternative option BIS (backward in space) is unsuitable to model oscillatory convection. Holzbecher [18] demonstrated that the ¯uctuations of the oscillatory regime are totally smoothed out by numerical dispersion which accompanies the use of the BIS scheme. The simulation using BIS did not reveal the oscillatory regime, while ¯uctuations were clearly visible in the same model run with the CIS option. Automatic timestepping is selected in the model simulations with parameter w ˆ 1: Simulation runs are stopped at dimensionless time tmax ˆ 100 or at a numerical steady state, i.e., when according to the accuracy requirement of the iterative algorithm the system does not change during a timestep simulation. Two di€erent modes of applications are used: branch tracing, and preferred mode simulation. In branch tracing, a simulation series is studied according to the following recipe: starting from a solution with a certain number of eddies simulations with a slightly changed Rayleigh number are initiated in which the ®nal state of the former run is taken as the initial state of the following run. In FAST-C(2D) this can easily be performed using the RESTART options as explained above. The steady state of a model is a good initial state for another simulation with a slightly changed parameter value when it is intended to keep the mode. In preferred mode simulations, the same initial state is used for a series of simulation runs. Usually it is not obvious to which convection mode the initial state leads. Initial states can be random distributions of any variable or peak disturbances in any variable. For the study, an equidistant grid spacing with 20  20 blocks is used in branch tracing simulations. Preferred mode simulations are run on a 40  40 grid. The following part gives a review of former ®ndings and new results for the branch tracing simulation series with particular interest on the 2nd mode convection cell. The last part of the paper tackles preferred mode simulations with particular focus on heat transfer. 3. Results and discussion 3.1. Oscillatory convection A preliminary study on oscillatory convection was undertaken within the software testing phase for the FAST-C(2D) code. Using a coarse grid Holzbecher [18] already showed that the critical value Racrit2 can be reproduced with the model quite accurately. The model was set up

198

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

Plate 1.

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

199

Fig. 1. Snapshot from oscillatory convection; streamfunction surface for eddies with opposite rotation.

on a 20  20 grid in a square box with spatial extension H ˆ L ˆ 1: The initial state of model runs for Ra ˆ 360, 380 and 400 was the steady state one-cell solution calculated for Ra ˆ 100 (1st mode) which had been saved in a previous run using the RESTART option. For Ra ˆ 360 and 400, Holzbecher [18] studied temporal changes of the maximum of velocity. The di€erence between the two simulations is striking: while for the lower Rayleigh number a constant value is approached, the regime for the higher Rayleigh number is clearly oscillatory. For Ra ˆ 400 the numerical simulation reveals an oscillatory ¯ow regime. For the lower Rayleigh number of 360 there are ¯uctuations as well, but these are much less pronounced than for the other Rayleigh numbers and the amplitude decreases. Apparently, the numerical model reproduces nicely the qualitative change of the solution from the steady to the oscillating regime. Table 1 Maximum and minimum of Nusselt numbers in the oscillation cycle for Ra ˆ 400, depending on boundary and approximation order (superscript)

Maximum Minimum

Nu (1)-top

Nu (1)-bottom

Nu (2)-top

Nu (2)-bottom

Nu (3)-top

Nu (3)-bottom

5.6269 5.2921

5.6033 5.2733

6.0014 5.6205

5.9553 5.5808

5.9241 5.5471

5.8685 5.4992

200

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

As the computations of Holzbecher [18] show, the ¯ow switches between two types of unicellular regimes during an oscillation cycle. In one, the ¯ow in the interior away from the boundaries is almost circular. In the other extreme, an elliptical shape of ¯uid ¯ow can be observed. The ellipse is lying diagonally in the square system. The change of the ¯ow pattern can best be viewed in the twelve snapshots of one cycle in the colour plate. The development of the ¯ow patterns coincides very well with motion patterns presented by Kimura et al. [8] and Caltagirone and Fabrie [7]. In the FAST-C(2D) model for Ra ˆ 400, the timeperiod of the oscillations is approximately 0.0128 in dimensionless time-units, which corresponds to a frequency around 491. This is close to the value of 471 derived by Riley and Winters [11] for Ra ˆ 400:2: Note that the ¯ow direction is not unique. If there is a convection cell with upward ¯ow on the right side and downward ¯ow on the left side, another cell with upward ¯ow on the left side and downward ¯ow on the right side is a valid solution as well. Streamfunction has positive values in one situation and negative in the other (see Fig. 1). In the mentioned study of Holzbecher [18] maximum velocity was taken as the indicator for the transition from steady to oscillatory convection. Other characteristics of the ¯ow pattern can be chosen to reveal the qualitative change in ¯ow regime, such as streamfunction extremes. Extreme values of the streamfunction can be found in the center of the convection eddies, while C approaches zero near all boundaries (see Fig. 1). Thus, extremes of C are good indicators for the inner cell ¯ow pattern. The maximum velocity is found near the edges of the modelled system and is an appropriate indicator for the behaviour of the solutions near the boundaries. Indicators which can be used alternatively are heat transfer or dimensionless Nu. The Nusselt number is selected in this study as the indicator of the hydraulic regime because it represents increased heat transfer additionally. A second reason is that Nusselt numbers of di€erent approximation order are better comparable on models with di€erent grid spacing than maximum velocity. In the simulations of Holzbecher [18] Nusselt numbers were calculated after each timestep for upper and lower boundaries using the three di€erent approximation schemes (1st order, 2nd order and 3rd order) mentioned above. Minima and maxima in one cycle as presented in Table 1 show di€erences between upper and lower boundary and changes with approximation order. The absolute value of heat transfer is delivered by the models with some uncertainty. Computed Nusselt numbers for upper and lower boundary are never equal Ð even for steady convection regimes. Nevertheless, for oscillatory regimes minima and maxima of Nu are observed at the same simulated time. Absolute values of Nu depend on the approximation order (see Eq. (12)). The uncertainty increases with increasing Rayleigh numbers, but could be reduced at least to some extent by using higher accuracy and ®ner grids. As the focus of this study is on the onset of oscillatory convection this point is not investigated further. In order to check the dependency of results on spatial discretization, a coarse grid …20  20† and a ®ne grid …40  40† were used. An oscillating regime is obtained in both ®ne and coarse grid simulations. Nevertheless, oscillations in the re®ned model are not as regular as in runs with the coarse grid. Nusselt numbers are overestimated in the coarse grid simulations. The di€erences are in the range between 0.3 and 0.6. Oscillations of the streamfunction extreme are smaller in the ®ne grid. Nevertheless, the mean value of the streamfunction peak is almost the same. The maximum velocity cannot be taken as an indicator of eddy strength near the

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

201

boundaries, when models with di€erent grid spacing are to be compared. vmax is smaller in the coarse grid, because velocities at grid points nearer to the model boundary are taken into account. The grid-PeÂclet-number provides a criterion if the grid is suciently ®ne. With maxima of velocity components vx, max and vz, max the maximum value for the grid-PeÂclet-number is computed as   vx, max Dx vz, max Dz …13† , Pemax ˆ Max D D For dimensionless variables, the expressions in the brackets become simply vx, max Dx and vz, max , because the time scale is normalized with respect to the di€usivity D (see transformation (6)). When velocity takes maximum values in the range of 100 the grid spacing needs to come down to 1/50, to ful®l the grid-PeÂclet-criterion in the entire system. In the FAST-model with Ra ˆ 400 velocities in that range are obtained. For higher Rayleigh numbers even higher maxima can be observed. Thus, even the ®ne grid models presented here usually have the grid PeÂclet-criterion not ful®lled in the entire model area. The grid-PeÂclet-criterion is violated in some of the blocks near the boundary where velocities are the highest. Note that the maximum velocities in a re®ned model are higher because grid points in closer vicinity to the boundaries are computed. That means that in contrast with common transport models, the grid-PeÂcletnumber is not necessarily lower on a ®ner grid. For this study further experiments are made for 2nd mode patterns with eddy aspect ratio

Fig. 2. Transient development of Nusselt numbers (2nd order app.) at bottom boundary for 2nd mode cells, modelled in a cavity with aspect ratio L=H ˆ 1=2 (FAST-C(2D), 20  20 grid, e ˆ 10 ÿ5 ); indication for a critical value Racrit; 2 ‰500, 520†:

202

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

  W ˆ Leddy Heddy ˆ 1 2: In order to start on a branch with a 2nd mode eddy, a single eddy initial state was chosen in a half cell …L ˆ 0:5, H ˆ 1† and a two-eddy initial state in a unit cell …L ˆ H ˆ 1). The initial situations are analytically equivalent. Nevertheless, there are di€erences from the numerical aspect, because both systems are modelled using a 20  20 grid. Thus the halfcell system is ®ner discretised in the horizontal direction than the unit cell system. Another principal di€erence is that transitions between the 2nd mode to 1st mode can only be studied in the unit cell. Results in the halfcell system indicate that there is a critical margin for the Rayleigh number which is between 500 and 520 (see Fig. 2). Oscillations can be observed in all simulations, but for the lower Rayleigh numbers ¯uctuations decay gradually while they remain at the same level for higher Ra between 510 and 520. Riley and Winters [11] found that `as the cavity becomes ¯atter (larger aspect ratio), the ¯ow becomes less stable'. The cited authors refer to systems with Leddy rHeddy in which, starting from aspect ratio 1.0, they observed the critical margin with increasing W to move to lower Rayleigh numbers. For aspect ratio W ˆ 3=2 Riley and Winters [11] obtained the critical number Racrit1 ˆ 310: As the preceding simulations with FAST-C(2D) show, the cited statement holds for slender convection cells, i.e., for Leddy Heddy RHeddy as well. Vice versa the formulation is: as the cavity becomes more slender (smaller aspect ratio), the ¯ow becomes more stable. A critical margin of Racrit2 ˆ 510 for W ˆ 1=2 ®ts well when the curve of Riley and Winters [11] is extrapolated for slender rolls. Experiments with the unit square system reveal that there is no transition from 2nd to 1st mode! In contrast, transitions in the opposite direction can be observed in some numerical experiments with FAST-C(2D). This again indicates that slender rolls are more stable. As the following shows, convective patterns with slender eddies are preferred for higher Rayleigh numbers. 3.2. Heat transfer in preferred modes In a dimensionless formulation, mean heat transfer through the system from the heated to the cooled boundary is given by the Nusselt number Nu. The Ra±Nu diagram illustrates how heat transfer changes with Rayleigh number. It is often forgotten to note that the curve is not unique. Nusselt number and heat transfer depend on one additional parameter Ð that is, the convection mode. If there is only one curve in the Ra±Nu diagram it should be mentioned which mode is drawn, or, if the mode is not constant, under which conditions the plot has been obtained. In the following, it is questioned which mode is likely or preferred for a given Rayleigh number. It is clear that in a deterministic approach, as it is presented here, the ®nal convection state is determined by the initial state of the system. Initial states which are by ®rst view near to certain modes are arti®cial and unlikely in natural systems. Convective patterns presented by Holzbecher [18] were studied on unicellular motions in square cells (1st mode). Starting from the one cell situation in the square cross-section, there remained one cell in all simulations. These type of cells had been chosen in order to compare with ®ndings published by other researchers. Using di€erent initial patterns, it turns out that the 1st mode convection is not `natural' under all circumstances. One example: when the simulation starts with a layered temperature

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

203

distribution which is disturbed at a single location, for Ra ˆ 500 a three-eddy steady pattern is approached as steady state. In fact it turns out that the stability of one-eddy ¯ow regimes for high Rayleigh numbers strongly depends on the unicellular initial condition. Numerical experiments with FAST-C(2D) show that for Ra  Racrit1 unicellular motions can only be reached when the initial state is near the 1-eddy branch. Following the numerical experiments 1st mode convection is not even likely to exist when the Rayleigh number is signi®cantly above the critical value and slender rolls are preferred. This view is supported by the series of model runs presented in the following. The background for the presented series of modeling is that inhomogeneities cause deviations from the homogeneous situation simulated in all numerical experiments which are reported so far. This is obvious under natural ®eld conditions, but even in laboratory experiments the observed states will show irregularities and disturbances at least at a minor scale. When this is true for pure ¯uid systems it is even more important in porous media where inhomogeneities can hardly be excluded. Even in controlled experiments inhomogeneities cannot be avoided m not even in the most carefully set up laboratory. Thus, if an initial state is to resemble a natural situation, it should contain some disturbances. This becomes most obvious when a numerical simulator is started from the pure conductive state without any disturbances: even for supercritical parameter combinations, the code does not start to show convective ¯ow regimes. The conductive state is a solution of the nonlinear Table 2 Nusselt numbers for steady convection; comparison of results from di€erent studies Mode

Rayleigh number Caltagirone [4] Schubert and Straus [29] Steen and Aidun [30] FAST-C(2D) (this study)

First

100 150 200 250 300

2.651

Second 100 150 200 250 300 400 500

2.136

Third

4.662b

a b

300 400 500 600 800

3.813 4.199 4.523

2.651 3.322 3.808 4.514

4.026

4.022

5.016

5.005 5.677

6.287a

7.034b 8.747b

4.947 6.108 7.489 8.312

Linear interpolated from Nu numbers for nearby Ra numbers. Aspect ratio y ˆ 0:3 (instead of 1/3).

2.65 3.32 3.81 4.19 4.52

2.620 3.363a 3.809 4.229a 4.603

2.14 3.24 4.03 4.58 5.02 5.69 6.20

2.136 3.261a 4.109 4.614 5.019 5.789 6.324 5.030 6.192 7.075 7.674 8.596

204

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

Fig. 3. Preferred mode simulation; streamfunction contours at dimensionless times t ˆ 1:e-4, 5.e-3, 0.01, 0.02, 0.04; left column: Ra ˆ 250, right column: Ra ˆ 480:

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

205

Fig. 4. Ra±Nu diagram for preferred mode simulations; 1st and 2nd mode results for comparison.

Fig. 5. Onset of oscillatory convection for 3rd mode cell in the unit square; Nu values calculated with FAST-C(2D).

206

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

set of di€erential equations and there is no change to other values. The numerical solution under these conditions does not re¯ect the physical instability. In that situation the modeler applies the simple trick to add a single disturbance in a state variable. If temperature, for example, has a single disturbance peak in the linear conductive state, the simulation moves towards a convective state. 1st and 2nd order oscillatory convection modes, which have been studied in the previous part, emerge out of a very arti®cial situation in which a parameter is changed gradually over a range of approximately one order of magnitude. It is the aim of the experiment series to ®nd the preferred mode for various Rayleigh numbers between 100 and 1000. The FAST-C(2D) code with automatic timestepping …w ˆ 1† is applied on a grid with 40  40 equidistant blocks. The accuracy criterion within the iterations is e ˆ 10 ÿ4 : Initial state to start the simulation is a random distribution of (normalized) temperature with mean value ymean ˆ 0:5 and variance 0.1. The same random distribution is used for all simulation runs. Intercomparison of results is possible for some Rayleigh numbers and some modes with three studies on steady convection. In Table 2, all possible combinations are listed which were presented by Caltagirone [4], Schubert and Straus [29] and Steen and Aidun [30]. From FASTC(2D) output the 1st order approximation for Nu on the top of the model are extracted throughout for this comparison. Taking into account that very di€erent numerical approaches are applied, the correspondence of results is very reasonable, for single values remarkable. As expected with increasing Rayleigh number, the di€erences between Nusselt numbers increase. For illustration purposes the development of the ¯uid pattern is portrayed for two simulation runs. Fig. 3 depicts the development of convective ¯ow patterns from the initial random distribution for two Rayleigh numbers: Ra ˆ 250 and 480. In one simulation the system ®nally runs into a two-eddy pattern while a three-roll pattern is approached in the other case. As an overview, the preferred modes are presented in a phase diagram. Fig. 4 portrays Nusselt number increase in dependence of Rayleigh number. For Rayleigh numbers between 100 and 300 the second mode of stable convection obviously dominates, giving the regime to the third mode for values above that interval. For comparison, ®rst mode and second mode branches are depicted additionally. The bifurcation between second and third mode obviously occurs where the Nusselt numbers are equal Ð so there is no discontinuity in the curve representing the preferred modes. In contrast, a similar switch between ®rst and second mode could not be observed. Moreover the numerical experiments show that inhomogeneities and local disturbances of Table 3 Maximum and minimum of Nusselt numbers in the oscillation cycle, depending on boundary and approximation order (for Ra ˆ 980)

Maximum Minimum

Nu (1)-top

Nu (1)-bottom

Nu (2)-top

Nu (2)-bottom

Nu (3)-top

Nu (3)-bottom

9.4003 9.1825

9.3420 9.2045

9.8611 9.6125

9.7910 9.6384

9.6960 9.4549

9.6268 9.4800

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

207

the ideal ¯ow pattern lead to a switch into another stable convective mode rather than to oscillatory convection. Thus the second critical Rayleigh number of the 1st mode Racrit2 1390, which marks the transition to oscillatory regimes, is practically irrelevant. At approximately Ra1970, the third mode switches to an oscillatory regime. Second order Nusselt numbers calculated from numerical output for the lower boundary of the model show the transition most clearly. They are presented in Fig. 5. As a decline of the ¯uctuations is obvious for Ra ˆ 960, the oscillating regime for Ra ˆ 980 is clearly visible. As Table 3 shows, computed Nu numbers depend on the approximation order. Apparently, ¯uctuations have a higher amplitude at the top boundary than at the bottom boundary. The values are given here in order to present the approximate magnitude and uncertainty of heat transfer to be expected. Di€erences in Nusselt numbers (for di€erent approximation orders and di€erent boundaries) can be taken as an indicator for the quality of the numerical results. Using this measure, the results obtained by a numerical simulation are accurate for low Rayleigh numbers and become less accurate with increasing Ra. The reason is that grid-PeÂclet-criterion PeR2 becomes less or not ful®lled with increasing Rayleigh number when the grid remains ®xed. Note that in contrast to mere transport calculations, the velocity is not known before, but is computed during the simulation run by a ®nite di€erence approximation of Eq. (4). Thus, the grid cannot be constructed on the basis of the grid-PeÂclet-criterion prior to the ®rst model run. It can be constructed a posteriori, i.e., after the ®rst model run. But even then it has to be realized that the grid-PeÂclet-number depends heavily on the grid re®nement at the boundaries, as just mentioned. For oscillatory convection in the 2nd mode Schubert and Straus [29] reported to have found the critical margin between Rayleigh numbers 650 and 700. Caltagirone [4] as well as Steen and Aidun [30] found a stable situation even for Ra ˆ 800: Both results are quite di€erent. Moreover, they di€er from the critical margin of 510 obtained in this study. There is better agreement for the 3rd mode. Caltagirone found that for Ra ˆ 1000, the cell with aspect ratio W ˆ 0:3 is stable, while the cell with W ˆ 0:4 is oscillatory. This can be taken as an indication that the 3rd mode (with W ˆ 1=3† is near the margin between stable and ¯uctuating regimes for Rayleigh numbers near 1000. This corresponds very well with the critical value of Ra1970; obtained in this study.

4. Conclusion A numerical model based on a nonlinear system describing density driven ¯ow in porous media is set up. For the veri®cation of the model some well-known characteristics of oscillatory convection, derived by spectral or bifurcation methods, can be reproduced. The critical margin of 1390 for oscillatory convection of the 1st mode pattern in the unit square is con®rmed by numerical experiments. For the 2nd mode a value of 1510 is obtained, for the 3rd mode it is found at 1970. Obviously critical Rayleigh numbers are increased for slender eddies; i.e., convection patterns with slender eddies are more stable. A series of experiments with FAST-C(2D) in the unit square with a random initial state are carried out in order to ®nd preferred modes. It turns out that the ®rst mode is relevant for low

208

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

supercritical Rayleigh numbers only. The second mode is relevant up to Rayleigh numbers of 300, when the third mode takes over the regime. For Rayleigh numbers up to 1970 stable convection patterns can be observed, before the 3rd mode oscillatory convection gets preference. It is unlikely that oscillatory convection can be observed in controlled experiments in porous media where inhomogeneities favour the transition to other stable modes which coexist with oscillatory convection in the studied Rayleigh number range. As Rayleigh numbers in natural geological formations are hardly supercritical in relation to the lowest critical margin Racrit1 (see [31]), oscillations can hardly play a role in ®eld situations either. Nevertheless, oscillatory convection is relevant from the mathematical point of view. It helps understanding the transition from steady to oscillatory and from oscillatory to chaotic regimes which is fundamentally di€erent in porous media than in free ¯uid convection. From the viewpoint of numerical modelling, there is another reason to study oscillatory regimes as it is done in this paper: well established results which are accepted in the scienti®c community can be used in veri®cation exercises for numerical codes.

References [1] BeÂnard, H. Les tourbillons cellulaires dans nappe liqiude transportant de la chaleur par convections en regime permanent. Rev. Gen. Sci. Pures Appl. Bull. Assoc. 11, 1900. 1261±1271, 1309±1328. [2] Rayleigh Lord. On convection currents in a horizontal layer of ¯uid when the higher temperature is on the under side. Phil Mag 1916;XXXII:529±46. [3] Lapwood ER. Convection of a ¯uid in a porous medium. Proc Camb Phil Soc A 1948;225:508±21. [4] Caltagirone J-P. Thermoconvective instabilities in a horizontal porous layer. J Fluid Mech 1975;72(2):269±87. [5] Straus JM, schubert G. Three-dimensional convection in cubic box of ¯uid-saturated porous material. J Fluid Mech 1979;121:301±13. [6] Schubert G, Straus JM. Transitions in time-dependent thermal convection in ¯uid-saturated porous media. J Fluid Mech 1982;121:301±13. [7] Caltagirone J-P, Fabrie P. Natural convection in a porous medium at high Rayleigh numbers. Eur. J. Mech. B/ Fluids 1989;(3):207±227. [8] Kimura S, Schubert G, Straus JM. Route to chaos in porous medium thermal convection. J Fluid Mech 1986;166:305±24. [9] Aidun CK, Steen PH. Transition to oscillatory convection in a ¯uid-saturated porous medium. J Thermophys 1987;1(3):268±73. [10] Aidun CK. Stability of convection rolls in porous media. In: Bifurcation phenomena in thermal processes and convection, HTD-vol. 94/AMD-vol. 89. ASME, 1987. p. 31±6. [11] Riley DS, Winters KH. Time-periodic convection in porous media: the evolution of Hopf-bifurcations with aspect ratio. J Fluid Mech 1991;223:457±74. [12] Schneider KI. Investigations of the in¯uence of free convection on heat transfer through granular material. In: 11th Int. Congr. on Refrigeration, Proc., 1963. p. 247±54. [13] Masuoka T. Heat transfer by free convection in a porous layer heated from below. Heat Transfer - Jap Res 1972;1(1):39±45. [14] Buretta RJ, Berman AS. Convective heat transfer in a liquid saturated porous layer. ASME J of Appl Mech 1976;43:249±53. [15] Jonsson T, Cattan I. Prandtl number dependence of natural convection in porous media. ASME J of Heat Transfer 1987;109:371±7.

E. Holzbecher / Computers & Fluids 30 (2001) 189±209

209

[16] Bau HH, Torrance KE. Low Rayleigh number thermal convection in a vertical cylinder ®lled with porous materials and heated from below. ASME J Heat Transfer 1982;104:166±72. [17] Combarnous MA, Bories SA. Hydrothermal convection in saturated porous media. Adv Hydrosci 1975;10:231± 307. [18] Holzbecher E. Modeling density-driven ¯ow in porous media. Berlin: Springer, 1998. [19] Elder JW. The unstable thermal interface. J Fluid Mech 1968;32(1):69±96. [20] Elder JW. Transient convection in a porous medium. J Fluid Mech 1967;27(3):609±23. [21] de Josselin de Jong G. Singularity distributions for the analysis of multiple-¯uid ¯ow through porous media. J Geophys Res 1960;65(11):3739±58. [22] Wooding RA. Rayleigh instability of a thermal boundary layer in ¯ow through a porous medium. J Fluid Mech 1960;9:183±92. [23] Holzbecher E. Numerical studies on thermal convection in cold groundwater. Int J Heat Mass Transfer 1997;40(3):605±12. [24] Holzbecher E. The in¯uence of variable viscosity on thermal convection in porous media. In: Int. Conf. Heat Transfer 98, Proc. 1998. p. 115±24. [25] Baumann R, Moser H. Modellierung der Meerwasserinvasion im Delta arider und semiarider Gebiete am Beispiel des Nildeltas. Z dt Geol Ges 1992;143:316±24. [26] Holzbecher E, Baumann R. Numerical simulations of saltwater intrusion into the Nile Delta Aquifer. In: Peters A, Wittum G, Herrling B, Meissner U, Brebbia CA, Gray WG, Pinder GF, editors. Comp Meth in Water Res X, Proc, vol. 2. Dordrecht: Kluwer Academic Publishers, 1994. p. 1011±8. [27] Holzbecher E, Heinl M. Anisotropy and dispersivity e€ects on saltwater upconing. Comp. Methods and Water Ressources, Beirut 1995:117±126. [28] Holzbecher E, Yusa Y. Numerical experiments on free and forced convection in porous media. Int J Heat Mass Transfer 1995;38(11):2109±15. [29] Schubert G, Straus JM. Three-dimensional and multicellular steady and unsteady convection in ¯uid-saturated porous media at high Rayleigh-numbers. J Fluid Mech 1979;94:25±38. [30] Steen PH, Aidun CK. Time periodic convection in porous media: transition mechanism. J Fluid Mech 1988;196:263±90. [31] Bjùrlykke K, Mo A, Palm E. Modelling of thermal convection in sedimentary basins and its relevance to diagenetic reactions. Marine and Petroleum Geology 1988;5:338±51.