Adapting a Fourier pseudospectral method to Dirichlet boundary conditions for Rayleigh-Benard convection

 

I. C. Ramos,¹ C. B. Briozzo¹ *

* E-mail: briozzo@famaf.unc.edu.ar

¹ Facultad de Matemática, Astronomía y Física, Universidad Nacional de Cárdoba, X5000HUA Cárdoba, Argentina.

---

We present the adaptation to non-free boundary conditions of a pseudospectral method based on the (complex) Fourier transform. The method is applied to the numerical integration of the Oberbeck-Boussinesq equations in a Rayleigh-Benard cell with no-slip boundary conditions for velocity and Dirichlet boundary conditions for temperature. We show the first results of a 2D numerical simulation of dry air convection at high Rayleigh number (R ~ 10⁹). These results are the basis for the later study, by the same method, of wet convection in a solar still.

Keywords Rayleigh-Benard convection; Pseudospectral method;

---

 

I. Introduction

Experimental observations 1 show that the onset of a turbulent convective flux can significantly enhance the efficiency of a basin-type solar still, but until now a theoretical explanation is lacking. Any adequate hydrodynamical simulation must incorporate the effects of moisture and condensation. Recent works 2, 3 show that this can be achieved through a Boussinesq-like approximation, which simplifies considerably the problem. However, realistic simulations are still demanding, given the need to resolve fine flux details and cope with Rayleigh numbers up to ~ 10⁹ 4.

Spectral methods 5 are well suited for this kind of tasks, and have many attractive features: they are simple to implement, show much better resolution and accuracy properties than finite difference or finite volume methods, and are highly efficient in large-scale simulations 6. Fourier-based pseudospectral methods are the simplest and fastest, since the discretized spatial differential operators are local, nonlinear terms can be computed through Fast Fourier Transform (FFT) convolutions, and solving the Poisson equation originating from the incompressibility (divergence-free) condition is almost trivial. Nevertheless, they usually work only for free (in fact, periodic) boundary conditions (BCs).

The presence of non-free (e.g., Dirichlet or Neumann) BCs introduces additional complications. For example, two-dimensional Rayleigh-Benard convection with laterally periodic BCs can be treated by using a spectral Galerkin-Fourier technique in the horizontal coordinate and a collocation-Chebyshev method in the vertical one 7, but vertical derivatives must then be computed by matrix multiplication. On a grid with N horizontal and M vertical points, this needs the solution of a linear system of dimension M for each of the N horizontal Fourier modes, at each timeintegration step.

Another complication with non-free BCs arises from the need to fulfill numerically the divergence-free condition, leading mainly to two different groups of methods (see Ref. 6 and references therein). In a first group the velocity field is written in terms of scalar potentials such that the divergence-free condition is satisfied by construction, e.g., in the 2D streamfunction-vorticity formulation or the 3D decomposition into toroidal and poloidal velocity potentials. In these methods, pressure is not present in the equations, but they lead to systems of higher-order partial differential equations with coupled BCs. In a second group, a primitive variable formulation of the equations is adopted and projection methods 8 are used to decouple velocity and pressure. These methods use a specific splitting of the equation system based on the chosen time-integration scheme, and determine pressure by projecting an appropriate velocity field onto a divergence-free space, leading to predictor-corrector algorithms. Besides the problem of correctly specifying the pressure BCs 9, 10, these methods require solving a Poisson equation for the pressure at each time-integration step. On a N = N x M grid, the best Fourier-based Fast Poisson Solvers (FPS) have operation counts O(Nlog₂ N) for the lowest (second) order discretization (and significantly worse for higher orders) 11,12, and those using GMRES are KO(N) with K > 100 13-15.

In this work, we will show how a Fourier-based pseudospectral method can be adapted to simple non-free (but periodic) BCs without losing its more appealing features. This is a first step towards building a pseudospectral simulation of wet air convection inside a basin-type solar still, and must be considered just as a proof of concept.

II. System

In these equations, the dynamical variables are the density p, the velocity u, the temperature T, and the pressure P. The parameters are the shear (or dynamic) viscosity n, the constant pressure specific heat capacity c_p, the thermal conductivity K, and the gravitatioinal acceleration .

The Boussinesq approximation 17 consists in discarding the dependence of n, c_p, and K on temperature and density, keeping

in the buoyancy term (-pgz in Eq. 1) but otherwise setting p = p elsewhere. Here, T = ¹ (T_h + T_c) is a reference temperature, p is a reference density (that of air at normal temperature and pressure), and a is the thermal expansion coefficient. Dropping the bars signifying reference quantities, absorbing some constants into the pressure gradient term, and defining the viscous diffusivity (or kinematic viscosity) v = n/p and the thermal diffusivity k = K/(pc_p), we obtain the (dimensional) Oberbeck-Boussinesq equations 17

where P is a reference (e.g. normal atmospheric) pressure, and z^' = z - H/2.

We now scale lengths with the cell height H, times with the characteristic vertical thermal diffusion time t_c = H²/k, and temperatures with the temperature difference AT. The nondimensional lengths, times, and velocities are then

is the Rayleigh number and a = v/k is Prandtl's number (~ 0.7 for dry air). The BCs we adopt are periodic in the horizontal direction, and homogeneous Dirichlet for both velocity u (no-slip BCs) and temperature 9 (perfect thermal contact) on the lower and upper plates.

Note that Eq. (13) is not a differential equation but a constitutive relationship, expressing the in-compressibility of the flux. In fact, the pressure term in Eq. (11) is computed by enforcing Eq.(13), which gives

where », j = x, z. In primitive variable integration schemes, this Poisson equation must be solved with adequate BCs at each time-step 6,9, to insure V . (dtu) = 0.

III. Helmholtz decomposition

Given a vector field f, twice continuously differen-tiable, Helmholtz's Theorem 18 states that it can be decomposed as

 

However, on a finite domain, the surface terms cannot be ignored, since they are essential for f_|_ to have the correct BCs, which by Eq. (22) are the same as those for u in Eq. (11). Using Eq. (24), the field in Eq. (23) can be seen to be a particular solution (the free-space solution) of the Poisson equation

where the last equation is needed to insure the transversality of w and hence of fj. This requirement can be automatically fulfilled by writing w explicitly as

which shows explicitly the transversality of fj_. The determination of the value of c, and the treatment of possible divergences in v at k = 0, are closely related and will be dealt with in the next section. The BCs for w at the lower and upper plates can be obtained from those for fj_, and are w = -v; the BCs on the horizontal direction are periodic but otherwise free, and will be automatically fulfilled by the constructive procedure for w given in the next section.

IV. Ultra-fast Laplace solver

We start by solving Eq. (30) with c = 0, that is Laplace's equation, on the rectangular domain 0 < x < L, 0 < z < H, which is an elementary problem in harmonic analysis. Over an unbounded domain, the solutions have the form e^iAxe^Az, where A is an arbitrary separation constant; the particular solutions for the case A = 0 are 1, x, z, and xz. Periodicity in x on 0, L imposes A = 2np/L with p G Z; the general solution is then

where a_p and b_p are constants, and z' = z - H/2. Here, we have used that cosh(0) = 1 to absorb all constant terms in ao, and used that sinh(0) = 0 to absorb the linear term in z^' as the particular term of the sum for p = 0, leaving Eq. (34) explicitly in the form of a (complex) Fourier series in x. The hyperbolic and linear functions in z^' can, in turn, be rather trivially expanded as complex Fourier series in z on the interval 0, H, leaving w in the form of the double Fourier series

Equations (44) and (31) show that the transverse field fj^_pq on the wavenumber grid can be obtained directly in terms of the non-transverse field ^f_pq without needing, at any point, to return to the coordinate grid. Moreover, Eq. (45) shows that c_pq and s_pq are given matrices that can be computed just once at the start of the simulation, as is the denominator k² _p + k² in Eq. (31). The only tasks to be performed at each time-step are then: first, computing the scalars v_pq from ^f_pq, requiring three multiplications per grid point; second, computing the sums in Eq. (44), which requires one multiplication per grid point; third, multiplying them by c_pq and s_pq, costing two multiplications per grid point; and fourth, multiplying by k_{z , q} and k_{x p} at each grid point. The total cost of obtaining fj_^_q is then 8N multiplications on a N = Nx M grid, thus outperforming even the best FPS by a significant factor on large grids.

It must be noted that the method introduced here has some similarity to the streamfunction-vorticity formulation 6, in the sense that the scalar fields v and W play a role similar to these potentials. However, in our method they are not taken as dynamical variables, and the evolution equations are not formulated in terms of them but of primitive variables. The method presented here shows also a strong similarity with projection methods 8, 9, but differently from them, pressure is not computed along the time evolution and, in fact, does no longer appear in the evolution equations.

V. Algorithm outline

We present now an outline of the numeric algorithm as we implemented it in the simulations. Initialization:

. Take zero velocity and Gaussian white noise for the temperature (with amplitude of the order of thermal noise), discretized on the coordinate grid, and take their FFT to get w^_pq, w^_pq, 9_pq on the wavenumber grid.

. Compute convolutions of U FFT with 2/3 rule.

Here, we have denoted by (f * g)_pq the convolution product of fields f_pq and g_pq discretized on the wavenumber grid, which is performed by FFT.

It must be noted that the spatial discretization of the evolution equations has been performed in a closed form, independent of the time-stepping algorithm to be employed to solve the resulting set of ordinary differential equations, outlined above. We must also note that this system does not involve multiplication by matrices with dimension N, the only nonlocal parts being the convolution products (handled through FFT) and the elimination of the longitudinal component of the velocity.

The time-stepping can then be performed by any algorithm designed to solve systems of ordinary differential equations. In our case, we opted for an adaptive-stepsize fifth order Runge-Kutta-Cash-Karp algorithm 19, which in our previous experience we have found efficient, stable, and flexible.

VI. Test runs

. Pre-compute (just once) the matrices c_pq and Even for a simple system like the one we are study- present only a brief outline. All results are given in terms of the laterally-infinite cell crytical Rayleigh number R_c — 1701, and the characteristic vertical diffusion time t_c which for our cell and medium is — 11797s. Coordinates and fields are in the dimen-sionless variables of Eqs. (11)-(13).

 

Figure 1: Temperature and velocity fields for R = 5R_c at t = 2t_c on a 32x 16 grid.

Figure 2: Temperature and velocity fields for R = 50R_c at t = t_c on a 64x32 grid.

 

At low R, a stationary regime state, consisting in two counter-rotating rolls, is reached in times - t_c or less. This time falls rapidly with increasing R, to - 0.1t_c at R - 1000R_c (see Figs. 1, 2 and 3).

Around R - 5000R_c, these rolls develop lateral oscillations, and the first "secondary structures" (small whirlpools) appear near the base of the ascending and descending plumes (see Fig. 4).

Above R - 5000R_c, the regime state becomes disordered and aperiodic, consisting of intermittent plumes and whirlpools in a wide size range (see Fig. 5).

At R - 5x 10⁵R_c, the temperature difference is - 65K; the smallest whirlpools are - 1cm wide, and the typical wind speeds are - 1m/s, in agreement

Figure 3: Temperature and velocity fields for R = 500R_c at t = 0.25t_c on a 128x64 grid (velocity decimated to a 64x 32 grid).

Figure 4: Temperature and velocity fields for R = 5000R_c at t = 0.05t_c on a 192x96 grid (velocity decimated to a 64x 32 grid).

with experimental observations 1 (see Fig. 6). The time to reach this regime state is rather short, ~0.001tc.

Note that in all figures, except for Figs. 1 and 2, the velocity grid has been decimated to enhance clarity.

VII. Code performance

Over the full range of R tested here (more than five orders of magnitude), the code maintained the typical velocity divergence and the field values at the bottom and top boundaries at essentially machine-precision zero, showing that the implemented method is sound.

Also over this range of R, the grid spacing needed to achieve "smooth" fields (i.e., to capture all the physical detail down to the smallest present scales) is consistent with the width of the (thermal) boundary layer. However, for coarser grids the code still gives qualitatively sound results; typically a checkerboard-like instability develops, but the algorithm keeps it quenched, showing very good stability even in presence of a severe accuracy loss.

Figure 5: Temperature and velocity fields for R = 5 x 10⁴R_c at t = 0.01t_c on a 384x192 grid (velocity decimated to a 64x 32 grid).

The algorithm is also fast: the simulation for R - 10⁹ (see Fig. 6), on a 512x256 grid, took less than one day per simulated minute on a single core of the 3GHz PentiumD CPU on which all our runs were performed, with no code optimization.

It is difficult to find in the literature a directly comparable simulation for more accurate benchmarking. However, from Ref. 6 we can see that, for example, the relaxation to a (steady) regime state in Rayleigh-Benard convection at R < 30000 - 20R_c on a - 20000-node grid takes 46 minutes on a similar processor at similar speed (3.2GHz Pentium 4), by the method implemented there. In the case of our code, at R = 20R_c on a 200 x 100 grid, and with the sole optimization of grouping the real FFTs in pairs (see the 2FFT algorithm in Ref. 19), the equivalent relaxation took 8 minutes, which is shorter by a factor of - 6. But it must be taken into account that with our initial conditions (zero velocities and thermal noise in temperatures) the convection onset is slow, and is followed by a transient stage with strong and disordered convec-tive patterns that decay to the regime state very slowly.

Figure 6: Temperature and velocity fields for R = 5x 10⁵R_c at t = 0.002t_c on a 512x256 grid (velocity decimated to a 64x 32 grid).

VIII. Conclusions and outlook

We have been able to show that a Fourier-based pseudospectral method can be adapted to a (admittedly simple) non-free BC setting, at the cost of moderate analytical work on the solutions of Pois-son's and Laplace's equations. The method is formulated in primitive variables, but the pressure is not explicitly computed nor referenced, like in a streamfunction-vorticity formulation. It also shares some properties with projection methods, but it decouples the implementation of the incompressibility condition from the time-stepping scheme, allowing great flexibility in the selection of the last. The resulting code is fast and stable, and implements the BCs and the incompressibility condition essentially to machine-precision.

Work on the extension of this scheme to a fully closed Rayleigh-Benard cell (i.e., with non-free BCs also on the lateral walls) is currently under course.

1 I De Paul, Evidence of chaotic heat enhancement in a solar still, Appl. Thermal Eng. 29, 1840 (2009).         [ Links ]

2 O Pauluis, Thermodynamic consistency of the anelastic approximation for a moist atmosphere, J. Atmos. Sci. 65, 2719 (2008).

3 O Pauluis, J Schumacher, Idealized moist Rayleigh-Benard convection with piecewise linear equation of state, Commun. Math. Sci. 8, 295 (2010).

4 T Weidauer, J Schumacher, Moist turbulent Rayleigh-Benard convection with Neumann and Dirichlet boundary conditions, Phys. Fluids 24, 076604 (2012)

5 B Fornberg, A practical guide to pseudospectral methods, Cambridge monographs on applied and computational mathematics, Cambridge University Press, Cambridge, UK (1999).

6 I Mercader, O Batiste, A Alonso, An efficient spectral code for incompressible flows in cylindrical geometries, Computers & Fluids 39, 215 (2010).

7 I Mercader, O Batiste, A Alonso, Continuation of travelling-wave solutions of the Navier-Stokes equations, Int. J. Numer. Meth. Fluids 52, 707 (2006).

8 A J Chorin, Numerical solution of the Navier-Stokes equations, Math. Comput. 22, 745 (1968).

9 P M Gresho, R L Sani, On the pressure boundary conditions for the incompressible Navier-Stokes equations, Int. J. Numer. Meth. Fluids 7, 1111 (1987).

10 R L Sani, J Shen, O Pironneau, P M Gresho, Pressure boundary condition for the time-dependent incompressible Navier-Stokes equations, Int. J. Numer. Meth. Fluids 50, 673 (2006).

11 V Fuka, PoisFFT - A free parallel fast Poisson solver, Appl. Math. Comput., 267, 356 (2015).

12 E Braverman, M Israeli, A Averbuch, L Vo-zovoiy, A fast 3D Poisson Solver of arbitrary order accuracy, J. Comput. Phys. 144, 109 (1988).

13 P M Gresho, D F Griffiths, D J Silvester, Adaptive time-stepping for incompressible flow. Part I: Scalar advection-diffusion, SIAM J. Sci. Comput. 30, 2018 (2008).

14 P M Gresho, D F Griffiths, D J Silvester, Adaptive time-stepping for incompressible flow. Part II: Navier-Stokes equations, SIAM J. Sci. Comput. 32, 111 (2010).

15 Y Saad, M H Schultz, GMRES: A generalized minimal residual algorithm for solving non-symmetric linear systems, SIAM J. Sci. Com-put. 7, 856 (1986).

16 L D Landau, E M Lifshitz, Fluid mechanics, Pergamon Press, Oxford (1959).

17 M C Cross, P C Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. 65, 851 (1993).

18 H von Helmholtz, Uber integrale der hydrodynamischen gleichungen, welche den wirbelbewegungen entsprechen, Celles J. 55, 25 (1858).

19 WH Press, S A Teukolsky, W T Vetterling, B P Flannery, Numerical recipes, Cambridge University Press, Cambridge, UK (1996).