Oscillatory Convection Onset in a Porous Rectangle with Non-analytical Corners

Abstract

The analytical theory on Darcy–Bénard convection is dominated by normal-mode approaches, which essentially reduce the spatial order from four to two. This paper goes beyond the normal-mode paradigm of convection onset in a porous rectangle. A handpicked case where all four corners of the rectangle are non-analytical is therefore investigated. The marginal state is oscillatory with one-way horizontal wave propagation. The time-periodic convection pattern has no spatial periodicity and requires heavy numerical computation by the finite element method. The critical Rayleigh number at convection onset is computed, with its associated frequency of oscillation. Snapshots of the 2D eigenfunctions for the flow field and temperature field are plotted. Detailed local gradient analyses near two corners indicate that they hide logarithmic singularities, where the displayed eigenfunctions may represent outer solutions in matched asymptotic expansions. The results are validated with respect to the asymptotic limit of Nield (Water Resour Res 11:553–560, 1968).

1 Introduction

Normal modes are in common use as analytical tools for solving eigenvalue problems of the fourth order and higher. These convenient tools have the disadvantage of overruling the nature of a higher-order eigenvalue problem, since normal modes are eigenfunctions of the second-order Helmholtz equation. Normal modes can properly be applied only to fourth-order eigenvalue problems with some degeneracy that reduces their effective spatial order from fourth down to second order only.

The Darcy–Bénard onset problem for convection in a porous medium is the simplest thermo-mechanical eigenvalue problem. Horton and Rogers (1945) solved this onset problem for a horizontal layer of unlimited horizontal extent. Their work was confirmed by Lapwood (1948), and the onset problem for a rectangular geometry is called the HRL (Horton–Rogers–Lapwood) problem.

Normal modes in a rectangular geometry can be understood by the concept of convection cell structures. It follows that the normal-mode compatible boundary conditions are not true external boundary conditions but rather internal conditions of continuity across a cell wall. Accordingly, it is important to distinguish between a flow cell and a thermal cell, as two parallel but interwoven thermo-mechanical structures. In the horizontal direction, there is a quarter of a wavelength phase shift between a flow cell and its associated thermal cell. No such phase shift exists in the vertical direction.

The set of 2D normal-mode compatible boundary conditions along a horizontal boundary is either two Dirichlet conditions for the temperature and streamfunction, or two Neumann conditions for these eigenfunctions. The standard set of 2D normal-mode compatible conditions along a vertical boundary is a Neumann temperature condition combined with a Dirichlet condition for the streamfunction, expressing the quarter of a wavelength phase shift between these eigenfunctions. A Dirichlet temperature condition combines with a Neumann condition for the streamfunction to another valid set of normal-mode compatible conditions along a vertical wall.

A corner in a HRL problem with rectangular geometry is analytical for the eigenfunctions only as far as at least one of the sides meeting in the corner has normal-mode compatible boundary conditions. The present work is the first theoretical work on a HRL problem where all four corners of the rectangle are non-analytical. In essence, no eigenfunctions will have physical meaning if they are analytically extended on the other side of a boundary. Moreover, it means that no valid Taylor expansions exist around a corner in the rectangle. Thereby high complexity is induced in the eigenfunctions, as the present paper will show.

The present paper terminates a journey of departure away from the normal-mode paradigm, leading to the highest possible complexity for a HRL problem in 2D. Nield (1968) wrote the pioneering paper that terminated the tacit monopoly of normal-mode solutions in the HRL problem. Nield studied all combinations of Dirichlet and Neumann thermo-mechanical conditions at the top and bottom of a horizontal layer. The majority of Nield’s tabulated cases are not compatible with normal modes. His three normal-mode solutions give critical Rayleigh numbers of $4 \pi ^2$, $\pi ^2$ and zero. The first non-normal-mode solutions for the horizontal direction were presented by Nilsen and Storesletten (1990). They studied the HRL problem in 2D with conducting and impermeable sidewalls. Tyvand (2002) tabulated the HRL onset problem in 2D with all different combinations of Dirichlet and Neumann conditions for the horizontal sidewalls. The single example with no closed-form solution lets one sidewall be impermeable and conducting, while the other wall is open and thermally insulating. This RT problem was solved by Rees and Tyvand (2004a). It is of particular mathematical interest since it has no steady-state solution, and the thermo-mechanical eigenfunction is oscillatory. It represents a wave that originates at one sidewall and is advected through the layer to be transmitted out through the other end-wall.

If the boundary conditions for a HRL problem are compatible with normal modes in all spatial directions except for one, the eigenvalue problem can in most cases be solved analytically. The eigenvalue problem is much more complicated when boundary conditions are incompatible with normal modes in two perpendicular directions. Such eigenvalue problems do not separate in the spatial variables, and it is impossible to derive analytical eigenfunctions. Tyvand and Nøland (2019) designed and solved numerically such a 2D problem in a rectangle, with non-normal mode eigenfunctions where the horizontal and vertical dependencies cannot be separated from one another. There are no known analytical solution methods, as there is at least one corner with conflicting thermo-mechanical conditions along with the boundaries meeting in that corner.

The present contribution is based on the RT problem, but we add extra complications by introducing boundary conditions at the upper plane that are incompatible with normal modes. Thereby we remove the possibilities of developing locally valid analytical solutions, as we probably introduce weak singularities in all four corners of the porous rectangle. There is no concept of a wave number, and the internal cell walls for the sets of eigenfunctions will neither be exactly vertical nor exactly horizontal. This situation contrasts the RT model, where all cell walls are vertical or horizontal. The rectangular shape of all cell walls in the RT paper is a trivial fact, following from the separation of variables in the horizontal and vertical directions, being dictated by the normal-mode conditions at the lower and upper boundary.

In this paper, the highest possible complexity of a HRL eigenvalue problem in 2D will be exposed, with the simple combinations of either Dirichlet or Neumann conditions at all four walls of a porous rectangle with horizontal and vertical sides.

2 Mathematical Formulation

A porous medium is considered, with 2D flow in the x, z plane. Its physical counterpart is a 3D porous medium with impermeable and insulating walls $y=0$ and $y=\Delta y$, where the distance $\Delta y$ is sufficiently small compared with the vertical length scale, implying that the preferred mode of free convection will be 2D in the x, z plane. The z axis is directed vertically upwards. The porous medium is homogeneous and isotropic. The assumption of 2D flow in the x, z plane is valid also for a 3D porous medium with thickness $\Delta y$ in the y direction, if two constraints are met: (i) The thermo-mechanical wall conditions at $y=0$ and $y=\Delta y$ are those of impermeable and thermally insulating walls. (ii) The cell width in the x direction must be sufficiently large, compared with $\Delta y$. From now on, we assume these conditions for 2D flow in the x, z plane to be met without specifying them further. Storesletten and Tveitereid (1991) investigated a problem where the range of validity of this assumption of 2D flow has been scrutinized.

The velocity vector ${\mathbf {v}}$ has Cartesian components (u, w). The temperature field is T(x, z, t), with t denoting time. The undisturbed motionless state has a uniform vertical temperature gradient. In the undisturbed state, the lower plane $z=0$ has the temperature $T=T_0$, and the upper plane $z=h$ has a temperature $T=T_0-\Delta T$, where $\Delta T$ is the temperature difference across the porous layer. The gravitational acceleration g is written in vector form as ${\mathbf {g}}$.

The standard Darcy–Boussinesq equations for convection in a homogeneous and isotropic porous medium can be written

$$\begin{aligned}&\nabla p +\frac{\mu }{K} {\mathbf {v}} + \rho _0 \beta ~ (T - T_0) {\mathbf {g}}=0,\end{aligned}$$

(1)

$$\begin{aligned}&\nabla \cdot {\mathbf {v}}=0,\end{aligned}$$

(2)

$$\begin{aligned}&(\rho c_p)_{\mathrm{m}}\frac{\partial T}{\partial t}+(\rho c_p)_{\mathrm{f}}~ {\mathbf {v}}\cdot \nabla T = \lambda _{\mathrm{m}} \nabla ^2 T. \end{aligned}$$

(3)

In these equations, p is the dynamic pressure, $\beta$ is the coefficient of thermal expansion, $\rho =\rho _0$ is the fluid density at the reference temperature $T_0$, $\mu$ is the dynamic viscosity of the saturating fluid, K is the permeability, $c_p$ is the specific heat at constant pressure, and $\lambda _{\mathrm{m}}$ is the thermal conductivity of the saturated porous medium. The subscript ‘m’ refers to an average over the solid/fluid mixture, while the subscript ‘f’ refers to the saturating fluid alone.

We consider a 2D porous medium with a rectangular shape, bounded by the four sides

$$\begin{aligned} x=0,\quad x= l,\quad z=0,\quad z=h, \end{aligned}$$

(4)

which means that the 2D rectangular porous box has height h and width l.

2.1 Dimensionless Equations

From now on we work with dimensionless variables. We reformulate the mathematical problem in dimensionless form by means of the transformations

$$\begin{aligned}&\frac{1}{h} (x, z) \rightarrow (x, z),\quad \frac{h}{\kappa _{\mathrm{m}} }(u, w)\rightarrow (u, w),\quad h \nabla \rightarrow \nabla ,\nonumber \\&\frac{1}{\Delta T} (T - T_0) \rightarrow T,\quad \frac{K }{\mu \kappa _{\mathrm{m}}} (p - p_0) \rightarrow p,\quad \frac{(\rho c_p)_{\mathrm{f}} \kappa _{\mathrm{m}} }{(\rho c_p)_{\mathrm{m}} h^2} t \rightarrow t, \end{aligned}$$

(5)

where $\kappa _{\mathrm{m}} = \lambda _{\mathrm{m}}/(\rho _0 c_p)_{\mathrm{f}}$ is the thermal diffusivity of the saturated porous medium. We denote the vertical unit vector by ${\mathbf {k}}$, directed upwards.

The dimensionless governing equations can then be written

$$\begin{aligned}&{\mathbf {v}}+\nabla p - R~ T {\mathbf {k}}=0.\end{aligned}$$

(6)

$$\begin{aligned}&\nabla \cdot {\mathbf {v}}=0\end{aligned}$$

(7)

$$\begin{aligned}&\frac{\partial T}{\partial t} +{\mathbf {v}} \cdot \nabla T = \nabla ^2 T. \end{aligned}$$

(8)

The Rayleigh number R has been introduced. It is defined as

$$\begin{aligned} R =\frac{\rho _0 g \beta K \Delta T h}{\mu \kappa _{\mathrm{m}} }. \end{aligned}$$

(9)

2.2 Basic Solution

The stationary basic solution of Eqs. (6)–(8) is given subscript “b”

$$\begin{aligned} {\mathbf {v}}_{\mathrm{b}} = 0,\quad T_{\mathrm{b}} = - z,\quad p_{\mathrm{b}} = - \frac{1}{2} R~ z^2. \end{aligned}$$

(10)

2.3 Linearized Perturbation Equations

In our stability analysis, we disturb the static state (10) with perturbed fields

$$\begin{aligned} {\mathbf {v}} = {\mathbf {v}}_{\mathrm{b}} + {\mathbf {v}},\quad T = T_{\mathrm{b}}(z) + \varTheta ,\quad p = p_{\mathrm{b}}(z) + p'. \end{aligned}$$

(11)

with perturbations ${\mathbf {v}}, \varTheta , p'$ that are functions of x, z and t. Linearizing Eqs. (6)–(8) with respect to perturbations and eliminating the pressure gives

$$\begin{aligned}&\nabla ^2 w = R \frac{\partial ^2 \varTheta }{\partial x^2},\end{aligned}$$

(12)

$$\begin{aligned}&\frac{\partial \varTheta }{\partial t} - w = \nabla ^2 \varTheta . \end{aligned}$$

(13)

We introduce the streamfunction $\varPsi$ by the relationships

$$\begin{aligned} (u, w) = \left( \frac{\partial \varPsi }{\partial z}, - \frac{\partial \varPsi }{\partial x}\right) , \end{aligned}$$

(14)

satisfying the continuity equation (7).

2.4 The Selection of Boundary Conditions

We will design the present problem for achieving the highest possible complexity of a HRL problem. This will be done by starting from the classical normal-mode eigenfunctions and gradually introduce boundary conditions that are incompatible with normal modes.

2.4.1 Normal-Mode Boundary Conditions for a Vertical Rectangle

Normal modes for the thermo-mechanical problem of convection onset has qualitative differences between the horizontal and vertical directions. Let us look at the preferred steady onset mode of the classical HRL problem

$$\begin{aligned} \varPsi = A \sin (\pi x) \sin (\pi z),\quad \varTheta = B \cos (\pi x) \sin (\pi z), \end{aligned}$$

(15)

written in a form valid for a square cavity filled with porous medium. A and B are amplitudes. This set of eigenfunctions satisfies the following boundary conditions at the lower and upper boundaries

$$\begin{aligned}&\varPsi = \varTheta = 0,\quad z=0,\end{aligned}$$

(16)

$$\begin{aligned}&\varPsi = \varTheta = 0,\quad z=1. \end{aligned}$$

(17)

The boundary conditions at the vertical sidewalls are

$$\begin{aligned}&\varPsi = \frac{\partial \varTheta }{\partial x} = 0,\quad x=0,\end{aligned}$$

(18)

$$\begin{aligned}&\varPsi = \frac{\partial \varTheta }{\partial x} = 0,\quad x=L. \end{aligned}$$

(19)

Here, L denotes the aspect ratio of the rectangle $L = l/h$. The normal modes (15) are valid for the square ($L=1$), but can be generalized to the rectangle by the transformation $x \rightarrow m x/L$. Here, m is a positive integer denoting the mth normal mode in the horizontal direction.

Normal modes in a rectangular geometry take care of continuity across internal cell walls. A normal mode formula for the x direction represents a possible physical reality on the other side of a boundary $x=x_0$. This means that a Taylor expansion in x is valid around $x=x_0$.

The four sides of the square represent cell walls for a flow cell. There is a horizontal phase shift of a quarter of a period between a flow cell and a thermal cell. Therefore, the eigenfunctions (15) represent thermal cells that are split into halves by the vertical walls.

There is another legal set of boundary conditions that produce normal modes in a vertical porous rectangle. It is comprised of the following boundary conditions at the lower and upper boundaries

$$\begin{aligned}&\frac{\partial \varPsi }{\partial z} = \frac{\partial \varTheta }{\partial z}= 0,\quad z=0,\end{aligned}$$

(20)

$$\begin{aligned}&\frac{\partial \varPsi }{\partial z} = \frac{\partial \varTheta }{\partial z}= 0,\quad z=1. \end{aligned}$$

(21)

The alternative normal-mode compatible boundary conditions at the vertical sidewalls are

$$\begin{aligned}&\frac{\partial \varPsi }{\partial x} = \varTheta = 0,\quad x=0,\end{aligned}$$

(22)

$$\begin{aligned}&\frac{\partial \varPsi }{\partial x} = \varTheta = 0,\quad x=L. \end{aligned}$$

(23)

These alternative boundary conditions compatible with normal modes allow the following set of eigenfunctions

$$\begin{aligned} \varTheta = A \sin (\pi x) \cos (\pi z),\quad \varPsi = B \cos (\pi x) \cos (\pi z), \end{aligned}$$

(24)

Now, the rectangular boundaries represent cell walls for the thermal cell, while the flow cells are split into halves by the vertical walls.

2.4.2 Non-normal-Mode Boundary Conditions for a Vertical Rectangle

We have now exposed the classes of Dirichlet- and Neumann-type boundary conditions compatible with normal modes in a rectangle. In the vertical direction, the two conditions for the eigenfunctions $\varPsi$ and $\varTheta$ had to be either two Dirichlet conditions or two Neumann conditions. In the horizontal direction, the two normal-mode compatible conditions for the eigenfunctions $\varPsi$ and $\varTheta$ are of opposite types, one of them a Dirichlet condition and the other one a Neumann condition.

Now, we have the background for constructing boundary conditions that are incompatible with normal modes. In the horizontal direction, such conditions for $\varPsi$ and $\varTheta$ may be both of the Dirichlet type or both of the Neumann type. In the vertical direction, we need a Dirichlet condition for $\varPsi$ and a Neumann condition for $\varTheta$, or vice versa.

2.4.3 The Choice of Boundary Conditions Implying Non-analytical Corners

We have now outlined boundary conditions that give non-normal modes in the horizontal direction. The most complicated eigenfunctions arise with Dirichlet conditions for both $\varPsi$ and $\varTheta$ at the left-hand wall $x=0$

$$\begin{aligned} \varPsi = \varTheta = 0,\quad x=0, \end{aligned}$$

(25)

in combination with Neumann conditions for both $\varPsi$ and $\varTheta$ at the right-hand wall $x=L$

$$\begin{aligned} \frac{\partial \varPsi }{\partial x} = \frac{\partial \varTheta }{\partial x} = 0,\quad x=L. \end{aligned}$$

(26)

This open-wall condition for $\varPsi$ corresponds to a hydrostatic reservoir surrounding the porous medium on the right-hand side, combined with a thermal condition of zero heat flux. This combination of kinematic and thermal condition may seem artificial, but it appears as a consistent limit case of the general thermo-mechanical Robin conditions for vertical walls, derived by Nygård and Tyvand (2010). Barletta and Storesletten (2015) applied this kinematic open-wall condition to the onset problem for a circular cylinder, where the wall was a heat conductor with zero perturbation temperature. The present set of end-wall conditions are the same as those applied in the RT paper, as we intend to formulate the most complicated non-normal-mode counterpart of that paper.

The chosen boundary conditions in the horizontal direction are mutually asymmetric to cause oscillatory thermo-mechanical eigenfunctions, known from the RT paper, where the problem was solved with normal modes in the vertical direction. In the recent TN paper, the upper boundary condition was chosen incompatible with normal modes. This gave a complicated time-dependent eigenfunction, yet the lower corners of the rectangle had a valid Taylor expansion for the vertical direction. This was because the lower boundary $z=0$ was chosen to be compatible with normal modes in the TN paper. Maximum complexity in the vertical direction is achieved when the boundary conditions are not compatible with normal modes, neither at the lower nor the upper boundary. We keep the same non-normal mode condition as in the TN paper at the upper boundary

$$\begin{aligned} \varPsi = \frac{\partial \varTheta }{\partial x} = 0,\quad z=1. \end{aligned}$$

(27)

The important novelty in the present paper is our choice of the second type of non-normal mode condition at the lower boundary

$$\begin{aligned} \frac{\partial \varPsi }{\partial x} = \varTheta = 0,\quad z=0. \end{aligned}$$

(28)

Note that we avoid the case of two equal conditions of the non-normal type. This is because such combinations imply symmetry or antisymmetry that will reduce the complexity of the eigenfunctions. It was shown by Nilsen and Storesletten (1990) that two equal Dirichlet conditions would give combinations of symmetric and antisymmetric eigenfunctions, which is indeed plausible because of the symmetry in the conditions. Such a right-left symmetry of the boundary conditions makes oscillations impossible, and it may allow internal cell walls that are parallel to the outer walls of the rectangle.

We sum up this formulation stage by concluding that the set of boundary conditions (25)–(28) gives the maximal complexity in the eigenfunctions for a vertical rectangle. The local solutions for each of the rectangle corners O, A, B and C are non-analytical with no valid Taylor expansions around these four corners. Any solution can only be valid within the local angular sector of $\pi /2$ around each corner, because of conflicting boundary conditions along the meeting sides of the rectangle.

3 The Thermo-Mechanical Eigenvalue Problem and Its Solution

The porous rectangle with the chosen boundary condition has no analytical corners, which complicates the solution procedure for several reasons: (i) We lack a good benchmarking for finite values of the aspect ratio L. There exists an analytical solution in the limit $L \rightarrow \infty$, given below. (ii) There may be weak integrable singularities in each of the four corners, but this is a challenging topic that we leave for future work. (iii) The numerical solution by the finite element method needs a very fine mesh near all the boundaries and around the four corners. An easy visual check of the convergence is that all streamlines or isotherms must be either tangential to a boundary (with a Dirichlet condition there) or perpendicular to a boundary (with a Neumann condition there).

With no available local results valid near the four corners of the rectangle, there is only one set of analytical results that we can use for direct comparison with our numerical results. It is the asymptotic onset criterion for the limit $L \rightarrow \infty$

$$\begin{aligned} R_{\mathrm{c}} = 17.65,\quad k_{\mathrm{c}} = 1.75, \end{aligned}$$

(29)

with zero frequency of oscillation. This onset criterion for a horizontal layer of infinite horizontal extent with the present choice of boundary conditions for the vertical direction was derived by Nield (1968).

In Fig. 1, we present a definition sketch in terms of dimensionless variables. Here, the four corners of the rectangle are denoted by O, A, B, C. Their respective coordinates are (0, 0), (L, 0), (L, 1) and (0, 1). The local eigenfunction fields at the onset of convection behave similarly near each of the four corners, and no analytical solutions are known for the local fields near these corners. We will attempt to use the concept of a local gradient angle (LGA) developed in our previous TN paper, to extract some analytical information from the numerical solutions around the four non-analytical corners.

Fig. 1

Illustration of the governing equations and corresponding boundary conditions of the coupled eigenvalue problem for a porous rectangle with four non-analytical corners

3.1 The Eigenfunctions at Marginal Stability

Following TN, the eigenfunctions $\varPsi (x,y,t)$ and $\varTheta (x,y.t)$ can be rewritten as

$$\begin{aligned} (\varPsi , \varTheta ) = (\psi (x,z), \theta (x,z)) e^{\sigma t}, \end{aligned}$$

(30)

with complex growth rate $\sigma = \sigma _r + i \sigma _i$ and i denoting the imaginary unit. The coupled governing equations can now be written

$$\begin{aligned}&\nabla ^2 \theta - \sigma \theta = \frac{\partial \psi }{\partial x},\end{aligned}$$

(31)

$$\begin{aligned}&\nabla ^2 \psi + R \frac{\partial \theta }{\partial x} = 0, \end{aligned}$$

(32)

with the boundary conditions

$$\begin{aligned}&\frac{\partial \psi }{\partial z} = \theta = 0,\quad z=0,\end{aligned}$$

(33)

$$\begin{aligned}&\psi = \frac{\partial \theta }{\partial z} = 0,\quad z=1,\end{aligned}$$

(34)

$$\begin{aligned}&\psi = \theta = 0,\quad x = 0,\end{aligned}$$

(35)

$$\begin{aligned}&\frac{\partial \psi }{\partial x} = \frac{\partial \theta }{\partial x} = 0,\quad x = L. \end{aligned}$$

(36)

Marginal stability is defined by $\sigma _r = 0$, defining $\omega = \sigma _i$. $\omega$ arise as a new eigenvalue in the problem, where real and imaginary parts of the eigenfunctions are defined as

$$\begin{aligned} \psi (x,z) = \psi _r + i \psi _i,\quad \theta (x,z) = \theta _r + i \theta _i. \end{aligned}$$

(37)

The governing equations can now be split in real and imaginary parts as follows

$$\begin{aligned}&\nabla ^2 \theta _r + \omega \theta _i = \frac{\partial \psi _r}{\partial x},\end{aligned}$$

(38)

$$\begin{aligned}&\nabla ^2 \theta _i - \omega \theta _r = \frac{\partial \psi _i}{\partial x},\end{aligned}$$

(39)

$$\begin{aligned}&\nabla ^2 \psi _r + R \frac{\partial \theta _r}{\partial x} = 0,\end{aligned}$$

(40)

$$\begin{aligned}&\nabla ^2 \psi _i + R \frac{\partial \theta _i}{\partial x} = 0. \end{aligned}$$

(41)

3.2 Numerical Solutions

The full solution for the marginal state of convection must be found numerically, since no analytical solutions are known for this fourth-order eigenvalue problem without degeneracies in 2D. We use the commercial finite element (FE) package Comsol Multiphysics.

The numerical solution may represent the outer solution in the sense of matched asymptotics, with respect to all the four corners of the rectangle. The two eigenvalues R and $\omega$ have to be solved simultaneously. The solution is found by iterations with respect to the growth rate $\sigma _r$, in order to achieve a valid result that fulfills $\sigma _r = 0$ with sufficient accuracy. The imaginary part that arises numerically for the eigenvalue $\omega$, is equal to $-\sigma _r$. The numerical search is based on an optimization problem to be solved, where the solver finally picks the correct $\omega$ to reach an imaginary part ($-\sigma _r$) less than $10^{-5}$.

A single convection cell as we define it here is either a flow cell or a thermal cell. A flow cell is a closed structure of corotating fluid with the same direction of vorticity. A thermal cell is a closed structure of positive (hot fluid) or negative perturbation temperature (cold fluid), The classical normal modes of the HRL problem has sinusoidal variation with x, where one may consider a full spatial period as a cell width, enclosing two neighboring structures with opposite sign. Here, we operate outside the paradigm of normal modes, and the appropriate way of looking at cell structures is to consider one single corotating structure as an individual flow cell. We therefore consider a hot structure with only positive perturbation temperature as an individual thermal cell.

3.2.1 The Porous Square $(L=1)$

The natural start is the square domain, which is displayed in Fig. 2. The four eigenfunctions with the lowest Rayleigh numbers are shown. In each case, five snapshots of the eigenfunctions are displayed, representing half an oscillation period. The isotherms are represented by solid black lines, supported by color coding. The streamlines are conveniently represented by white graphs on the colored background. There are constant intervals in the displayed values of the perturbation temperature and streamfunction.

Fig. 2

Streamlines (white) and isotherms (black) at equally spaced time intervals over one-half of a period for convection in a square (aspect ratio $L=1$) for the four first eigenfunctions. Time increases from the left to the right

The preferred mode for convection has a critical Rayleigh number $R = 49.48$, about 25% above the classical value of $4 \pi ^2$ (Horton and Rogers 1945). We observe the curved cell walls, both for flow cells and for thermal cells. The concept of individual cells is harder to define than in the previous models of oscillatory convection (RT and TN), as fragmented cell structures abound. A dilemma is that a Neumann condition does not expose a kinematic cell wall (a dividing streamline) among all other isolines that are also normal to the boundary. Dividing streamlines that express cell walls appear only at the upper boundary. The dividing isotherms (of zero perturbation temperature) can be seen only along the lower boundary.

The lateral walls of the rectangle are easier to interpret, as far as the cell structure is concerned. The left-hand boundary $x=0$ represents an external wall for both the streamlines and the isotherms, so it may have the role of a cell wall, even though it is in principle different from an internal cell wall. The right-hand boundary seems to take the role of a vertical cut through an individual cell, since both the streamlines and the isotherms are horizontal at the right-hand boundary $x>L$, directed perpendicularly to it.

There are no possible lateral extrapolations of the mathematical eigenfunctions to be valid solutions for $x<0$ or $x>L$. The present problem is unique in the sense that no valid extrapolations exist beyond the four sides of the rectangle, neither from a mathematical nor a physical point of view. Therefore, all the four corners are non-analytical, which means that no Taylor expansions are valid in the infinitesimal neighborhood of these corners. The numerical solutions that we find are at best valid outer solutions in the sense of matched asymptotics, and in a later section, we will try to extract additional information about possible inner solutions close to the corners. There is no spatial periodicity within the square geometry of Fig. 2, yet we can estimate roughly the number of half-cells in each direction.

The first eigenmode shown in Fig. 2 has a one-cell structure horizontally, and the snapshot in the middle shows only cold fluid, which means that it will have only hot fluid half an oscillation period later. In the vertical direction, there is almost a full cell structure.

The second eigenmode in Fig. 2 has a one-cell structure vertically. The horizontal direction covers more than two cell widths. The third eigenmode in Fig. 2 increases the cell number vertically from one to two. The fourth node in Fig. 2 preserves the two-cell structure horizontally and increases the number of cells horizontally to roughly two and a half cell.

3.2.2 A Narrow Porous Domain $(L=0.5)$

Figure 3 shows the four eigenfunctions with the lowest Rayleigh numbers for a narrow rectangle with an aspect ratio $L=0.5$. We observe how the approximate number of cells vertically starts with one, increasing to two and three and then back to one again for the fourth eigenmode. The number of horizontal cells stays slightly below two for the first three eigenmodes, and it jumps to three horizontal cells for the fourth eigenmode. A peculiarity of the fourth eigenmode is that there are stages where much of the flow goes directly through the porous cavity, entering at the bottom and leaving at the open right-hand wall. When this happens, the thermal cells are mostly closed structures. In other stages during one period of oscillation, the flow is recirculating while the thermal field does not make closed cells.

Fig. 3

Streamlines (white) and isotherms (black) at equally spaced time intervals over one-half of a period for convection in a narrow rectangle (aspect ratio $L=0.5$) for the 4 first eigenfunctions. Time increases from the left to the right

3.2.3 A Moderately Wide Porous Domain $(L=2)$

Figure 4 displays the two eigenfunctions with the lowest Rayleigh numbers for a narrow rectangle with the aspect ratio $L=2$. It is followed up by Fig. 5, which shows the third and fourth eigenfunctions for $L=2$. As in the previous figures, we see how the number of cells in the horizontal and vertical directions tend to increase with increasing Rayleigh numbers.

The first mode in Fig. 4 shows isolated blobs of closed isolines near the lower and upper isolines, surrounded by isolines that go perpendicularly through the boundary. These blobs are of two types: (i) The pointed oval blobs surrounded by isolines that penetrate a boundary where the blob meets that boundary in a single apex point. Such thermal blobs are always attached to the upper boundary, while similar streamline blobs are attached to the lower boundary. (ii) The almost triangular blobs that are confined by two right-angle isolines meeting a boundary. Such thermal blobs are attached to the lower boundary, while similar streamline blobs are attached to the upper boundary.

Fig. 4

Streamlines (white) and isotherms (black) at equally spaced time intervals over one-half of a period for convection in a wide rectangle (aspect ratio $L=2$) for the 2 first eigenfunctions. Time increases from the top to the bottom

Fig. 5

Streamlines (white) and isotherms (black) at equally spaced time intervals over one-half of a period for convection in a wide rectangle (aspect ratio $L=2$) for the third and fourth eigenfunction. Time increases from the top to the bottom

3.2.4 The Local Gradient Angle for $L=2$

Our recent TN paper introduced the concept of a local gradient angle (LGA). It is relevant near a Dirichlet corner, where both the boundary conditions for the two sides meeting in the corner are of the Dirichlet type. The present LGA concept comes with only one of the two coupled eigenfunctions having Dirichlet conditions at a corner: The streamfunction has its Dirichlet corner at the origin O, while the temperature perturbation has C as its only Dirichlet corner. These corners are non-analytical, with no valid Taylor series around them, different from our previous TN model with valid Taylor series in the z direction around the corner where we computed the LGA.

Figure 6 shows an intuitive illustration of the LGA concept, based on how a ray of light emitted from the origin is reflected back (with the isoline acting as a mirror surface) along the same direction as it came from, while all other rays directed differently are scattered from this surface. The LGA for the temperature perturbation is computed at the origin, while the streamfunction LGA is computed at the corner C. As in our previous TN paper, we pick the aspect ratio $L=2$ for the LGA computations. This chosen moderate value of L aims to prevent too much disturbance from the singularities of other corners in this local LGA analysis.

Let us recall what an LGA analysis for a normal-mode solution would be like. A full normal-mode solution in 2D would always give LGA equal to $\pi /4$, which is seen from the classical HRL solution. However, the solution from Nilsen and Storesletten (1990) shows a Dirichlet corner for the streamfunction where the LGA could suddenly become zero as one internal cell wall happens to coincide with the outer cell wall, as explained further by Rees and Tyvand (2004b). The degeneracy of the rectangle problem with a conducting outer wall causes exceptional point values of zero for the LGA to occur, but its regular value will remain fixed as $\pi /4$.

Fig. 6

Local gradient angle (LGA) at different corner radii, to be computed as a function of the phase of the oscillation period for the first eigenfunction with aspect ration $L=2$. This concept works only at corners where the relevant boundary conditions along the two meeting sides are of the Dirichlet type. a Corner O for the thermal eigenfunction. b Corner C for the stream eigenfunction

These insights from earlier work are kept in mind for analyzing the present LGA computations. The variation of the LGA in two corners during an oscillation period for the case $L=2$ are shown in Fig. 7. The upper plot (a) depicts the LGA for the temperature field at the corner O. The lower plot (b) depicts the LGA for the streamfunction at the corner C. First, we show computations for the radius 0.1, which is considered too large to give convergence for the local solution near each corner. Next, we show results with the same radius 0.05 for calculating LGA as in our previous TN paper, noting that this is almost exactly the normal-mode type of solution that comes out from the model of Nilsen and Storesletten (1990), if the internal cell wall moves steadily. The only difference is that our LGA needs a finite phase interval to shift from $\pi /4$ to zero, and then back to $\pi /4$ again, but this is a necessity due to the chosen finite value for the radius of computation around the corner.

Fig. 7

Local gradient angle (LGA) at different corner radii as a function of the phase of the oscillation period for the first eigenfunction with aspect ration $L=2$. a Corner O for the thermal eigenfunction. b Corner C for the stream eigenfunction

These results tell us that abrupt changes in the LGA happen when an internal cell wall merges with the outer physical wall, which also clarifies the similar behavior observed in the TN paper.

We go further in investigating the local convergence of the LGA at the corners O and C, but what we find is that the results do not converge. When we reduce the radius of computation from 0.05 down to 0.025, and further down to 0.01, the LGA curves drift away from the normal-mode-like behavior that we observe at a radius of computation 0.05. These systematic deviations are almost the same for the two corners O and A. As a result, they can hardly be due to numerical errors, so there should be mathematical reasons for them to appear. A plausible reason is that logarithmic singularities exist in each corner, which can in principle be studied by matched asymptotic expansions.

Logarithmic singularities can be integrable and should not affect the numerical convergence for the onset criterion governed by the outer flow in terms of matched asymptotic expansions. To solve the outer problem analytically by matched asymptotic expansions is no option because of its complexities. The LGA computations indicate that our numerical solution represents the leading-order outer solution of a matched asymptotic expansion.

In the “Appendix,” we give a broader discussion on the non-analytical corners and their possible singularities. An analytical corner would give a converged, constant value for its local gradient angle (LGA).

3.2.5 A Very Wide Porous Domain $(L=6)$

There is no wave number concept emerging from our numerical results so far. The four non-analytical corners of the rectangle forbid regular cell structures to occur in their neighborhood. Nevertheless, we may identify a horizontal cell width along a horizontal Dirichlet boundary, if we evaluate it sufficiently far from the corners. Such calculations require a relatively large value of L. We have therefore chosen the value $L=6$ for plotting the preferred set of eigenfunctions in Fig. 8, where the values for the central cell width extracted from this figure are tabulated in Table 1. The above asymptotic result (29) derived by Nield (1968) predicts a cell width of $\pi /1.75 = 1.795$ length units, as $L \rightarrow \infty$. This result gives a reasonable benchmarking for the thermal cell widths in Table 1. Nield’s result does not fit so well with the tabulated widths of flow cells. This is not surprising since asymmetries in the mutual couplings of thermal and stream eigenfunctions prevent their asymptotic behavior from being identical.

Fig. 8

Streamlines (white) and isotherms (black) at equally spaced time intervals over one-half of a period for convection in a very wide rectangle (aspect ratio $L=6$) for the first eigenfunction. Time increases from the top to the bottom

The higher eigenfunctions for $L=6$ fail to give converged cell widths, seen from the snapshots of the second eigenfunction shown in Fig. 9.

Fig. 9

Streamlines (white) and isotherms (black) at equally spaced time intervals over one-half of a period for convection in a very wide rectangle (aspect ratio $L=6$) for the second eigenfunction. Time increases from the top to the bottom

Table 1 Central cell width at a Dirchlet-type boundary for the first thermal and stream eigenfunctions with a wide aspect ratio ($L=6$) at different phase instants of the oscillation period

3.3 Critical Rayleigh Number and Oscillation Frequency

The Rayleigh number at marginal stability is calculated as a function of the aspect ratio and is plotted in Fig. 10, and the associated angular frequency is plotted in Fig. 11. At each value of L, the five successive modes with the lowest Rayleigh numbers are being selected for these plots. In Fig. 10, we are able to see some slope discontinuities at points where one onset mode is replaced by a different mode. These switches from one onset mode to another are much more dramatic in Fig. 11, where great jumps in the oscillation frequencies are observed. The ordering of unstable modes is based on their Rayleigh numbers only, not on their oscillation frequencies, so it should not be surprising that such jumps occur. In Fig. 10, we have included a dashed line marking the analytical limit for the Rayleigh number derived by Nield (1968), given by Eq. (29) above.

Fig. 10

The critical Rayleigh number for the onset of convection is plotted as a function of the aspect ratio. At each value of L, the five modes with the lowest values of R are displayed. The lowest portion of the figure is magnified, with an inclusion of Nield’s asymptotic value $R=17.65$, valid as $L \rightarrow \infty$ for the lowest mode

Fig. 11

The angular frequency of oscillation $\omega$ at marginal stability as a function of the aspect ratio L, plotted with logarithmic scale. The five onset modes with the lowest Rayleigh numbers are plotted. These are not ordered according to the magnitude of the frequency, which is why a number of discontinuities appear in this figure. The curves are plotted with straight lines connecting the computed points $(L, \omega )$ obtained by successive increments 0.05 in the aspect ratio L

In the figures, we have not included comparisons with earlier work on oscillating convection. The critical Rayleigh number and the associated angular frequency of oscillation are not far away from the results for the benchmark RT problem (Rees and Tyvand 2004a).

4 Discussion and Outlook

The Darcy–Bénard onset problem of convection in a porous medium with rectangular geometry enjoys the status of being a highly mature problem, where little is left for further analysis. This is true only as far as the conventional normal-mode paradigm is concerned (Horton and Rogers 1945), allowing spatially separable extrapolations from normal modes. Little is still known on non-separable eigenfunctions of non-normal mode type (Tyvand and Nøland 2019; Tyvand et al. 2019). The shortcomings of normal modes can no longer be overlooked, as these second-order tools are not adapted to the full fourth-order physics of a Darcy–Bénard problem without degeneracies.

Qualitative comparisons can be made between the fourth-order Darcy–Bénard problem and the eigenvalue problem of vibrating elastic plates (Reddy 1999). The standard boundary condition of a clamped elastic plate falls outside the normal-mode paradigm, posing standing mathematical challenges even for rectangular plates. With exterior forcing, a classical non-normal mode solution for wedge-shaped elastic plates was found by Moffatt (1964). He derived this singular solution for the mathematically equivalent problem of a driven viscous cavity with Stokes flow.

Thermo-mechanical normal modes for 2D Darcy–Bénard convection in a vertical porous rectangle has a peculiarity compared with 2D elastic plates. The role of buoyancy in free convection creates asymmetries between the temperature field and the flow variables. Asymmetries between a normal-mode streamfunction and its coupled temperature perturbation are manifested as phase shifts. In the horizontal direction, these coexisting cells have a mutual phase shift of $\pi /2$, while they are synchronized with a zero phase shift in the vertical direction. These asymmetries give the basis for designing the present eigenvalue problem.

Our motivation was to depart as far as possible from the normal-mode paradigm for a 2D Darcy–Bénard problem. We have designed the richest eigenvalue problem available for a vertical porous rectangle with Dirichlet and/or Neumann boundary conditions. The problem that we established and solved numerically has four corners that are all non-analytical. No Taylor expansions are valid extrapolations out from the corners, and no extrapolations beyond the boundaries of the rectangle can have physical meaning.

Singularities should exist at non-analytical corners, since a local Taylor series around a corner would have a finite radius of convergence in the absence of singularities. The mildest singularities are of a logarithmic type, and we cannot rule them out.

In this paper, it is inevitably more difficult to classify the individual convection cells in comparison with earlier contributions reporting oscillatory convection. The perturbation isotherms penetrate the upper boundary, while the streamlines penetrate the lower boundary. The non-analytical corners tend to forbid regular cell structures. The chosen boundary conditions allow no spatial periodicity inside the rectangle. The closed-cell concept is fragmented for a narrow rectangle.

The concept of a local gradient angle (LGA) was introduced in our previous work (TN), and it is further extended in the present paper. The numerical LGA sensitivity analysis suggests that there exist logarithmic singularities in each corner. The wave number concept does not apply for most of the studied cases. A sufficiently wide rectangle (with an aspect ratio of 6) allows an approximate cell width to emerge from careful evaluation sufficiently far from the corners. Thereby the large aspect ratio limit of the preferred onset mode is assessed against the analytical predictions of Nield (1968).

We have designed an eigenvalue problem with boundary conditions incompatible with normal modes both horizontally and vertically, resulting in a complexity in time-dependent eigenfunctions beyond the known HRL solutions from the literature. The presented numerical solutions for the coupled eigenfunctions are considered as representing outer solutions of matched asymptotic expansions. A difficult challenge is to derive valid inner solutions for the present fourth-order thermo-mechanical problem, where a major obstacle is that the singularities in the outer solution are still unknown.

Appendix: On Corner Singularities due to Conflicting Boundary Conditions

Eigenvalue problems of fourth order in rectangular domains may possess corner singularities for several reasons: (i) The fourth-order biharmonic operator does not separate in orthogonal coordinate systems. (ii) A double set of boundary conditions are needed along each wall, enhancing possibilities of conflict at any corner where walls with different conditions meet.

A typical eigenvalue problem of second-order is based on the Laplace operator (Helmholtz equation), and it does not possess corner singularities. At any corner in a rectangular geometry, a local Taylor expansion will be valid for an eigenvalue problem with the Helmholtz equation.

We consider a 2D rectangular domain with a corner in origin O where the two sides $z=0$ and $x=0$ meet. An analytical corner in the origin O for an eigenfunction $\theta$ can be of three different types: (i) Two sides with Dirichlet conditions meet, and the leading Taylor expansion out from the corner is $\theta = A x^m z^n$, where m and n are positive integers. (ii) Two sides with Neumann conditions meet, and the leading Taylor expansion is $\theta = (A + B x^{2m})(1 + C z^{2n})$, where again m and n are positive integers. (iii) The x axis has a Dirichlet condition while the z axis has a Neumann condition (or vice versa), which gives a leading Taylor expansion $\theta = (A + B x^{2m})z^n$, where as before m and n are positive integers.

The normal-mode eigenfunctions for the Helmholtz equation provide a basis for series expansions and integral transforms. Analytical methods intended for higher-order problems may be based on normal modes. Such tools are not versatile, as they may restrict physical modeling when mathematical degeneracies are induced. Normal modes cannot be fully representative of fourth-order problems.

The present eigenvalue problem does not allow normal modes to be used since all four corners are non-analytical. No local Taylor expansion is valid around any of these corners since there are conflicting boundary conditions along with any pair of sides meeting in the four corners.

Before a further discussion of these non-analytical corners, we remark that conflicts between boundary conditions may arise even for the Laplace equation. An instructive example is an impulsive wavemaker introduced by Peregrine (1972). A piston forces the liquid into impulsive horizontal motion, while the free surface allows only vertical motion. At the waterline where the piston and the free surface meets, an inevitable singularity arises for the surface particle which has an infinite ratio between its free vertical velocity and its forced horizontal velocity. This is a logarithmic singularity for the velocity potential, but it is an integrable singularity representing a finite dimensionless virtual mass 0.5428. This is the effective liquid mass impulsively forced into motion by the piston, measured by the unit depth and the liquid density. Needham et al. (2007) resolved the logarithmic singularity of the impulsive wavemaker by developing an inner solution in terms of matched asymptotic expansions. Corner singularities for second-order problems may thus occur by forcing on one side of a corner which is in conflict with a homogeneous boundary condition on the other side of the corner.

There are several ways of obtaining singularities for fourth-order problems involving the biharmonic operator: (i) A forcing problem in a rectangle, where a forcing condition on one side of a corner is in conflict with a homogeneous boundary condition on the other side of the corner. (ii) Forcing problems where the non-separability of the biharmonic operator causes singularities. The sequence of Moffatt vortices Moffatt (1964) is a known example of this type of corner singularity. (iii) Eigenvalue problems with corner singularities. These may be caused by the non-separability of the biharmonic operator in combination with a mutual conflict between homogeneous boundary conditions along with the two sides meeting in the corner.

The set of coupled thermo-mechanical eigenfunctions $\psi$ and $\theta$ cannot be separated when the corners are non-analytical. Even though $\psi$ and $\theta$ can be uncoupled to form two separate fourth-order equations, it is impossible to formulate a full set of boundary conditions for $\psi$ alone or for $\theta$ alone. Thus, the eigenvalue problem has to consist of two coupled second-order differential equations, with one Dirichlet or Neumann condition for $\psi$ and one such condition for $\theta$ along each side of the rectangle.

Non-analytical corners have no locally valid Taylor expansions. They arise when one of the thermo-mechanical eigenfunctions switches between Dirichlet type and Neumann type at a corner while the other eigenfunction maintains its condition (either Dirichlet or Neumann type) through the same corner. At all four corners of the present rectangle, there are similar mutual conflicts in the set of boundary conditions: At the corners O and B the streamfunction conditions switch, while the temperature perturbation maintains its Dirichlet condition (at O) or Neumann condition (at B). At the two other corners, A and C the temperature conditions switch, while the streamfunction maintains its Neumann condition (at A) or Dirichlet condition (at C) through these corners. In all these cases, the local streamlines are perpendicular to the local isotherms along one of the sides meeting in a corner, while these streamlines are parallel to the isotherms along with the other side meeting in the same corner. This is a conflict that leads to great mathematical difficulties. Corner singularities seem unavoidable, in light of the previous work for the simpler uncoupled problems. The oscillatory time dependence of the eigenfunctions adds to the difficulties of this coupled fourth-order problem.

The fact that we find satisfactory numerical convergence for the essential global parameter (the Rayleigh number eigenvalue) indicates that the singularities are weak, and they are probably integrable. Again, we recall the integrable logarithmic singularity from Needham et al. (2007), giving a finite virtual mass at the end of a semi-infinite fluid column put impulsively into motion.