1. Introduction
Pattern formation in large systems that are driven away from equilibrium is ubiquitous and is observed in diverse contexts. Important fluid examples include turbulence (Pope Reference Pope2000; Duraisamy et al. Reference Duraisamy, Iaccarino and Xiao2019), fluid convection (Cross & Hohenberg Reference Cross and Hohenberg1993; Bodenschatz et al. Reference Bodenschatz, Pesch and Ahlers2000), cellular flame fronts (Sivashinsky Reference Sivashinsky1983), chemical reaction fronts (Lee et al. Reference Lee, McCormick, Ouyang and Swinney1993; Rongy et al. Reference Rongy, Trevelyan and De Wit2008; Mukherjee & Paul Reference Mukherjee and Paul2019, Reference Mukherjee and Paul2020, Reference Mukherjee and Paul2022; Gao et al. Reference Gao, Paul, Chen, Hong, Chamorro, Sottos and Geubelle2023) and droplet freezing (Kant et al. Reference Kant, Müller-Groeling and Lohse2020). Pattern formation is also important in many biological systems such as microbial colony organisation (Levine & Reynolds Reference Levine and Reynolds1991; Tai et al. Reference Tai, Mukherjee, Nero, Olson, Tithof, Nadell and Yan2022; Chuang et al. Reference Chuang, Rao and Leibler2025), leaf and flower arrangement in plants (Pennybacker & Newell Reference Pennybacker and Newell2013), wound healing (Sherratt & Murray Reference Sherratt and Murray1990; Maini et al. Reference Maini, McElwain and Leavesley2004) and neuronal depolarisation waves in the brain (Somjen Reference Somjen2004; Dahlem et al. Reference Dahlem, Graf, Strong, Dreier, Dahlem, Sieber, Hanke, Podoll and Schöll2010; Mukherjee et al. Reference Mukherjee, Mirzaee and Tithof2023).
In many of these systems, the patterns emerge behind the propagation of a front separating the patterned and unpatterned regions in space (Ben-Jacob et al. Reference Ben-Jacob, Brand, Dee, Kramer and Langer1985; Cross & Hohenberg Reference Cross and Hohenberg1993; van Saarloos Reference van Saarloos2003). Examples are proliferating band formation of differing cell densities in bacterial colonies subjected to oxygen and nutrient gradients (Chuang et al. Reference Chuang, Rao and Leibler2025), radially propagating chemical fronts coupled to advection (Maharana et al. Reference Maharana, Negrojević, Comolli and De Wit2026), cellular clump formation in the wake of ‘active’ cell swarms (Ford et al. Reference Ford, Celora, Westbrook, Dalwadi, Walker, Baumann, Weijer, Pearce and Chubb2025), sporulation behind expanding biofilm matrix fronts (Srinivasan et al. Reference Srinivasan, Vladescu, Koehler, Wang, Mani and Rubinstein2018) and oscillating structures in the wake of bimolecular chemical reactions in fluids influenced by buoyancy and Marangoni instabilities (Budroni et al. Reference Budroni, Upadhyay and Rongy2019; Bigaj et al. Reference Bigaj, Budroni, Escala and Rongy2023).
The pattern formation that occurs behind a front has been carefully studied in the laboratory using well-controlled fluid experiments. This includes propagating vortex fronts in Taylor–Couette flow (Ahlers & Cannell Reference Ahlers and Cannell1983; Ben-Jacob et al. Reference Ben-Jacob, Brand, Dee, Kramer and Langer1985), vortex fronts in the wake of a cylinder (Yang & Zebib Reference Yang and Zebib1989), patterns behind a front in a layer unstable to the Rayleigh–Taylor instability (Fermigier et al. Reference Fermigier, Limat, Wesfreid, Boudinet and Quilliet1992) and propagating fronts of vortices (Ben-Jacob et al. Reference Ben-Jacob, Brand, Dee, Kramer and Langer1985; Fineberg & Steinberg Reference Fineberg and Steinberg1987). The velocity of the front, and the spatial scales of the fluid structures left in its wake, can depend upon how far the system is driven from equilibrium, the details of the underlying dynamics, the boundary conditions, the mechanism of front initiation and the size of the domain for finite systems.
Considerable effort has been devoted to building a physical understanding of pattern-forming fronts (Dee & Langer Reference Dee and Langer1983; Dominguez-Lerma et al. Reference Dominguez-Lerma, Ahlers and Cannell1984; Ben-Jacob et al. Reference Ben-Jacob, Brand, Dee, Kramer and Langer1985; Fineberg & Steinberg Reference Fineberg and Steinberg1987). Near critical, an analysis using the amplitude equation yields a prediction for the velocity of the front (Dee & Langer Reference Dee and Langer1983). A useful prediction of the selected pattern is often the one which maximises the growth rate of perturbations in the linear regime (Dominguez-Lerma et al. Reference Dominguez-Lerma, Ahlers and Cannell1984). These theoretical predictions have been tested numerically (Ben-Jacob et al. Reference Ben-Jacob, Brand, Dee, Kramer and Langer1985; Lücke et al. Reference Lücke, Mihelcic and Kowalski1987; Kockelkoren et al. Reference Kockelkoren, Storm and van Saarloos2003) and experimentally (Fineberg & Steinberg Reference Fineberg and Steinberg1987). While the predicted front velocity is generally in good agreement with experimental and numerical results, discrepancies remain in the wavenumber selected by the convection rolls behind the propagating front (Cross & Hohenberg Reference Cross and Hohenberg1993; Kockelkoren et al. Reference Kockelkoren, Storm and van Saarloos2003).
Rayleigh–Bénard convection (RBC) is a canonical system for studying pattern formation in a controlled setting that is accessible to experiments (Bodenschatz et al. Reference Bodenschatz, Pesch and Ahlers2000) and numerical simulations (cf. Paul et al. Reference Paul, Chiam, Cross, Fischer and Greenside2003). RBC is typically studied as the bulk instability that results in the buoyancy-driven fluid motion of a shallow fluid layer of depth
$d$
when heated uniformly from below in a gravitational field. When the temperature difference,
$\Delta T = T_h - T_c$
, between the hot bottom surface
$T_h$
and the cold top surface
$T_c$
equals a critical value
$\Delta T_c$
, yielding the critical Rayleigh number
${Ra}_c$
, convective motion arises. The fluid motion occurs over the entire fluid layer simultaneously in the form of convection rolls. The convection rolls form with a critical wavenumber
$q_c$
(and critical wavelength
$\lambda _c = 2 \pi /q_c$
). For an infinite layer of fluid, with no-slip top and bottom surfaces,
${Ra}_c = 1707.6$
and
$q_c = 3.117$
(
$\lambda _c = 2.016$
) where
$q_c$
and
$\lambda _c$
have been non-dimensionalised using
$d$
(Chandrasekhar Reference Chandrasekhar1961).
In this study, we do not explore the bulk instability where convective fluid motion occurs everywhere in the domain simultaneously. Instead, we consider an initially quiescent layer of fluid which has been locally perturbed to initiate the formation of a convection roll which yields a front propagating into still fluid leaving convection rolls in its wake. The details of our process for accomplishing this are the following.
For time
$t \lt 0$
, we set
$\Delta T = 0$
and the entire fluid layer is motionless. At
$t = 0$
we prescribe that the temperature of the fluid layer varies linearly from the hot bottom surface to the cold top surface as required by heat conduction in the absence of fluid motion. For
$t \ge 0$
we set
$\Delta T \ge \Delta T_c$
and the temperature difference then remains at this constant value for all time. A local perturbation is introduced at
$t = 0$
to initiate the formation of a convection roll which results in a propagating front of emerging convection rolls. It is useful to note that the timescale for the spontaneous nucleation of convection rolls from the bulk instability scales inversely with the reduced Rayleigh number, near onset (Ahlers et al. Reference Ahlers, Cross, Hohenberg and Safran1981; Fineberg & Steinberg Reference Fineberg and Steinberg1987). However, the front dynamics we study here occurs in the window after front initiation and before this bulk instability develops. Once the bulk instability occurs, the entire fluid layer undergoes convective motion which annihilates the propagating front. Our analysis is of the formation and propagation of the pattern-forming front prior to its destruction by the bulk instability.
The front velocity
$v_{\kern-1pt f}$
and the wavenumber
$q$
of the convection rolls that form behind the front depend upon the reduced Rayleigh number,
$\epsilon$
, which is a measure of how far the fluid layer is away from threshold where
$\epsilon = ({Ra} - {Ra}_c) / {Ra}_c$
. As the fluid layer approaches critical from above,
$\epsilon \rightarrow 0$
, the front velocity vanishes and the wavenumber of the forming convection rolls approaches the finite value
$q_0$
.
Linear stability theory predicts that the front velocity will grow as
$v_{\kern-1pt f} \sim \epsilon ^{1/2}$
near threshold (Dee & Langer Reference Dee and Langer1983). The front-selected wavenumber at
$\epsilon = 0$
,
$q_0$
, must be contrasted with the critical wavenumber,
$q_c$
, which describes the wavenumber of convection rolls that form in an infinite fluid layer due to the bulk fluid instability at
$\epsilon =0$
. Linear stability theory predicts
$q_0 = q_c$
for a front of convection rolls in an infinite layer of fluid and for
$\epsilon \ll 1$
the wavenumber increases as
$q \sim \epsilon$
(Dominguez-Lerma et al. Reference Dominguez-Lerma, Ahlers and Cannell1984). As we will discuss in detail,
$q_0 \ne q_c$
for all of the conditions we explore.
In a series of pioneering experiments by Fineberg & Steinberg (Reference Fineberg and Steinberg1987) propagating fronts of convection rolls were investigated over a wide range of conditions. These experiments yielded
$v_{\kern-1pt f} \sim \epsilon ^{1/2}$
in agreement with predictions using the amplitude equation. However, the measured values of wavenumber were significantly different from the theoretical predictions. The wavelength of the front-selected convection rolls in experiment is described by
$\lambda = \lambda _0 (1-b\sqrt {\epsilon })$
, where
$\lambda _0 = 2.29$
(
$q_0 = 2.744$
) and
$b = 0.18$
. It is important to emphasise that
$\lambda _0$
is the wavelength of the front-selected convection rolls at
${Ra}_c$
.
The experimental measurements of the convection rolls disagree with theory in several important respects. The experiments yield
$\lambda _0 \gt \lambda _c$
, where the wavelength of the front-selected convection rolls measured in experiment for
$\epsilon \rightarrow 0$
is 13.6 % larger than
$\lambda _c$
. In addition, the variation of the wavelength with
$\epsilon$
in experiment is different from the theoretical prediction. Expressing the experimental result in terms of wavenumber, expanding for small
$\epsilon$
and keeping only the leading term yields
$q/q_0 = 1 + 0.18 \epsilon ^{1/2}$
. In summary, theory predicts
$q_0 = q_c$
and
$q/q_0 \sim \epsilon$
, while experimental measurement yields
$q_0 \lt q_c$
and
$q/q_0 \sim \epsilon ^{1/2}$
.
Previous numerical simulations by Lücke et al. (Reference Lücke, Mihelcic and Kowalski1987) reproduced the theoretically predicted linear scaling
$q/q_0 \sim \epsilon$
with
$q_0 = q_c$
. The difference between
$q_0$
and
$q_c$
in experiment was attributed to the details of the apparatus used (Lücke et al. Reference Lücke, Mihelcic and Kowalski1987; Cross & Hohenberg Reference Cross and Hohenberg1993) and to the slow relaxation of the fronts to their asymptotic values which makes measurements very difficult (Kockelkoren et al. Reference Kockelkoren, Storm and van Saarloos2003).
In this work, we numerically investigate the propagation of convection rolls over a range of experimentally relevant conditions, including domain size, domain geometry, method of front initiation and boundary conditions. We explore how these factors influence the front velocity and the wavenumber of the convection rolls. We use a highly flexible spectral element approach (Nek5000 2026) for a range of conditions, including the addition of specific geometrical features of the experiments of Fineberg & Steinberg (Reference Fineberg and Steinberg1987), to quantitatively explore fundamental features of pattern-forming fronts.
The remainder of the paper is organised as follows. In § 2, we discuss the general approach, including the governing equations and the details of the computational domains, boundary conditions and methods of front initiation that are used. In § 3, the numerical results are discussed. We first discuss the propagation of fronts which leave behind straight parallel convection rolls. We begin with a two-dimensional (2-D) domain and then compare with three-dimensional (3-D) box domains which include a range of features and boundary conditions. We then discuss fronts of concentric rolls that occur in a large cylindrical domain. Lastly, our conclusions are presented in § 4.
2. Approach
2.1. Convection domains
We use several computational domains to explore the propagation of pattern-forming fronts. Our intention with this section is to clearly present the details of the different simulations we conduct prior to discussing the numerical results. We begin our study with a 2-D domain which we then extend to become a 3-D box domain. With the box domain, we explore a range of experimentally motivated features including different front initiation approaches and the influence of thin fin structures attached to the sidewalls to inhibit fluid motion near the walls. Finally, we explore fronts of concentric rolls formed in a large cylindrical domain. These different cases are all discussed in turn below.
The 2-D domain, shown in figure 1, has a length
$L$
in the
$x$
direction where
$d$
is the depth of the fluid layer in the
$z$
direction. The bottom surface is hot (red) with temperature
$T_h$
, the top surface is cold (blue) with temperature
$T_c$
and gravity opposes the
$z$
direction. All material boundaries are no-slip surfaces. For most of our results, we have used a domain with aspect ratio
$\varGamma = L/d = 30$
. However, we found it necessary to extend the domain to
$\varGamma = 60$
to study fronts near onset,
$\epsilon \ll 1$
, which required a longer time and larger domain for the fronts to approach their asymptotic state. For all of the 2-D simulations, the fronts are initiated using a constant-temperature hot sidewall at
$x = 0$
(red) and the fronts propagate from left to right.
The 2-D domain used to study the propagation of fronts of straight parallel rolls. A fluid layer of length
$L$
and depth
$d$
with a hot bottom wall at
$T_h$
(red) and cold top wall at
$T_c$
(blue) where
$T_h \gt T_c$
. The front is initiated using a hot wall located at
$x = 0$
(red) and propagates from left to right. The aspect ratios
$\varGamma = L/d$
used in the simulations are
$\varGamma = 30$
and 60. Gravity opposes the
$z$
direction.
We also explore the propagation of fronts of straight parallel rolls in 3-D box domains (see figure 2). The bottom surface is hot at
$T_h$
, the top surface is cold
$T_c$
and gravity opposes the
$z$
direction. The first box domain explored, not shown explicitly in figure 2, has finite lengths in the
$x$
and
$y$
directions of
$L_x$
and
$L_y$
, respectively, to yield aspect ratios
$\varGamma _x = L_x/d$
and
$\varGamma _y = L_y/d$
. We have used
$\varGamma _x = 27.3$
and
$\varGamma _y = 6.54$
to align with the apparatus used in the experiments of Fineberg & Steinberg (Reference Fineberg and Steinberg1987). For this domain,
$\varGamma _x \gg \varGamma _y$
, and as a result the convection rolls orient their axes along the direction of the shorter side of the domain. As a result, straight parallel rolls form with their axes aligned in the
$y$
direction and the formation of new rolls propagates in the positive
$x$
direction as indicated by the arrow.
A box domain with fins attached to the sidewalls. (a) A fluid layer of depth
$d$
with a rectangular cross-section of lengths
$L_x$
and
$L_y$
in the
$x$
and
$y$
directions, respectively to yield
$\varGamma _x = L_x/d$
and
$\varGamma _y = L_y/d$
. The bottom wall is hot (red) and the top wall is cold (blue). Thin fins are attached to the sidewalls (shown in grey). The front is initiated at the wall located at
$x = 0$
(red) and propagates in the positive
$x$
direction. (b) A
$y$
–
$z$
cross-section highlighting the fin geometry. The fins are no-slip surfaces of finite thickness in the
$y$
and
$z$
directions and are perfect thermal conductors. We also use a box domain without fins (not shown) which is this domain with the grey fins removed. Gravity opposes the
$z$
direction and the schematics are not drawn to scale.
We also investigate this box domain with the addition of thin fins attached to the sidewalls (see figure 2), as used in the experiments of Fineberg & Steinberg (Reference Fineberg and Steinberg1987). The fins are attached to the sidewalls, at midheight (
$z = 1/2$
), and extend from the sidewalls into the domain. This is shown in figure 2(a) where the fins are the grey surfaces. The wall where the front is initiated, at
$x = 0$
, is shown in red and the front propagates in the positive
$x$
direction.
A
$y$
–
$z$
cross-section of the domain is shown in figure 2(b) highlighting the fin geometry. The fins have a thickness, in the
$z$
direction, of 0.1 and they extend from the sidewall into the domain by a distance of 1 where distance has been non-dimensionalised using the layer depth
$d$
. All material surfaces, including the fins, are no-slip surfaces and the fins are perfect thermal conductors.
In experiment, it is often desired to include a ‘soft’ boundary (Kramer et al. Reference Kramer, Ben-Jacob, Brand and Cross1982) on the sidewalls that is less restrictive than a no-slip surface. Examples include ramped and finned boundaries (Daviaud & Pocheau Reference Daviaud and Pocheau1989; de Bruyn et al. Reference de Bruyn, Bodenschatz, Morris, Cannell and Ahlers1996; Bajaj et al. Reference Bajaj, Mukolobwiez, Currier and Ahlers1999; Paul et al. Reference Paul, Cross and Fischer2002, Reference Paul, Chiam, Cross, Fischer and Greenside2003). The intention behind using a finned boundary is often to provide a sidewall boundary condition with reduced thermal forcing and reduced viscous shear. The fins are no-slip surfaces extending into the domain that quench convection in the region above and below the fins. The quenching of the fluid motion can be traced to the cubic depth dependence of the Rayleigh number. However, the influence of the finned boundaries on the convection rolls, and on the propagating front, is quite complex (Paul et al. Reference Paul, Chiam, Cross, Fischer and Greenside2003). A thorough theoretical understanding of how the fins affect the fluid dynamics is not currently available and we use experimentally accurate numerical simulations to explore this further.
Figure 3 shows a cylindrical domain of aspect ratio
$\varGamma = r_0/d$
, where
$r_0$
is the radius of the domain. The front is initiated at the centre (
$x = y = 0$
) by a locally imposed thermal perturbation causing the formation of a concentric convection roll. A concentric convection roll is interesting in light of our study because it represents an experimentally accessible 3-D case of a convection roll without a sidewall boundary due to its axisymmetric structure. The concentric roll at the centre causes the formation of additional rolls resulting in a front propagating radially outward leaving a field of concentric rolls (or target pattern) in its wake.
Cylindrical domain used to study the propagation of concentric convection rolls. A fluid layer of depth
$d$
, radius
$r_0$
and aspect ratio
$\varGamma = r_0/d$
with a hot bottom wall (red) at
$T_h$
and cold top wall (blue) at
$T_c$
. All material boundaries are no-slip surfaces. The front is initiated at the origin by introducing a thermal perturbation at
$t = 0$
. The front on concentric convection rolls propagates radially outward from the origin towards the sidewalls. Gravity opposes the
$z$
direction and in our study
$\varGamma = 40$
.
2.2. Computational approach
The fluid motion due to RBC is described by the non-dimensional Boussinesq equations
\begin{eqnarray} \textit{Pr}^{-1} \left ( \frac {\partial \boldsymbol u}{\partial t} + \boldsymbol u \boldsymbol{\cdot } \boldsymbol{\nabla }\boldsymbol u \right ) &=& - \boldsymbol{\nabla }p + {\nabla} ^2 \boldsymbol u + {Ra} T \hat z, \\[-28pt] \nonumber \end{eqnarray}
\begin{eqnarray} \frac {\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{\nabla }T &=&{\nabla} ^2 T , \\[-28pt] \nonumber \end{eqnarray}
\begin{eqnarray} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} &=& 0 , \\[0pt] \nonumber \end{eqnarray}
which represent the conservation of momentum, energy and mass, respectively where
$\boldsymbol{u}$
is the fluid velocity vector,
$T$
is the temperature,
$p$
is the pressure,
$\hat {z}$
is a unit vector in the
$z$
direction and
$ \textit{Pr}$
is the Prandtl number. The non-dimensionalisation is done in the typical manner using
$d$
as the length scale, the vertical diffusion of heat
$d^2/\kappa$
as the timescale, where
$\kappa$
is the thermal diffusivity of the fluid, and
$\Delta T$
as the temperature scale. In the following, we assume that all variables are in non-dimensional form.
We integrate (2.1)–(2.3) using a parallel spectral-element approach (Deville et al. Reference Deville, Fischer and Mund2002; Nek5000 2026) that has been used extensively to study open questions regarding fluid convection (Paul et al. Reference Paul, Cross, Fischer and Greenside2001; Paul et al. Reference Paul, Cross and Fischer2002; Scheel & Cross Reference Scheel and Cross2006; Karimi & Paul Reference Karimi and Paul2012; Scheel et al. Reference Scheel, Emran and Schumacher2013; Mehrvarzi & Paul Reference Mehrvarzi and Paul2014; Xu & Paul Reference Xu and Paul2016; Mukherjee & Paul Reference Mukherjee and Paul2019). All material surfaces are no slip,
$\boldsymbol{u} = 0$
, and the temperatures of the bottom and the top surfaces are
$T(z = 0) = 1$
and
$T(z = 1) = 0$
, respectively. All sidewalls, other than the sidewall where front initiation occurs in the rectangular domains, and the fins, when included, are perfect thermal conductors. As a result, these sidewalls and the fins are held at a constant temperature for all time given by the thermal conduction profile
$T(z) = 1 - z$
.
In the 2-D domain (figure 1) the front is initiated using a hot sidewall
$T(x = 0) = T_0$
where we use
$T_0 = 1$
. In the box domains, we initiate the fronts using a constant-temperature hot sidewall, as well as a constant-heat-flux
$q''$
sidewall such that
$q''(x = 0) = q_0''$
, where
$q_0''$
is a constant. The heat flux has been non-dimensionalised using
$k_{\kern-1pt f} \Delta T/d$
, where
$k_{\kern-1pt f}$
is the thermal conductivity of the fluid. Using Fourier’s law of heat conduction this boundary condition can also be expressed as
$\partial T / \partial x|_{x=0} = - q_0''$
. We have used
$q_0'' = 0.014$
, which aligns with the wall heating used in the experiments of Fineberg & Steinberg (Reference Fineberg and Steinberg1987). We have found that
$v_{\kern-1pt f}$
and
$q$
do not vary significantly for
$1\times 10^{-3} \lesssim q_0'' \lesssim 1$
.
In the cylindrical domain (figure 3) the front is initiated at the centre of the domain by imposing a localised Gaussian perturbation to the temperature field which can be expressed as
\begin{align} T(x,y,z,t = 0) = {\rm e}^{- \tfrac{(x^2 + y^2)}{\eta ^2}} ,\end{align}
where
$\eta ^2 = 1/2$
. This yields a rapid, localised drop-off of the initial disturbance.
2.3. Simulation cases
The simulation cases we use are summarised in table 1. Simulation 1 uses the 2-D domain (figure 1) with a constant-temperature wall boundary condition to initiate the front. We use
$\varGamma = 30$
except for simulations near the onset of convection,
$\epsilon \lesssim 0.1$
, where we use
$\varGamma = 60$
to allow for the front to settle towards its asymptotic state. We use both
$ \textit{Pr} = 1$
and
$ \textit{Pr} = 5.373$
for
${Ra}_c \lesssim {Ra} \leq 10^4$
. For
${Ra}{\gtrsim }10^4$
, the bulk instability of the fluid layer occurs too rapidly to study propagating fronts using our approach. Simulations 2 and 3 use box domains without and with fins, respectively. For the box domains we use both constant-temperature and constant-heat-flux boundary conditions to initiate the fronts. Simulation 4 uses the cylindrical domain (figure 3) with
$ \textit{Pr} = 1$
and
$3 \times 10^3 \leq {Ra} \leq 10^4$
where the front is initiated by a thermal perturbation at the centre of the domain.
Summary of the simulations including the aspect ratio, domain geometry, Rayleigh number, Prandtl number and the method to initiate the front. Here
$T_0$
indicates a hot wall with a constant temperature,
$q_0''$
is a wall with a constant heat flux and ‘Gaussian’ represents a localised Gaussian disturbance of the temperature given by (2.4).
Table 1. Long description
The table presents a summary of four simulations, detailing their aspect ratios, domain geometries, Rayleigh numbers, Prandtl numbers, and initiation methods. It has four rows and six columns. The columns are labeled Simulation, Aspect ratio, Domain, Ra, Pr, and Initiation. Row 1: Simulation 1, Aspect ratio Γ = 30, 60, Domain 2-D, figure 1, Ra 1708 ≤ Ra ≤ 10^4, Pr 1, 5.373, Initiation T0. Row 2: Simulation 2, Aspect ratio Γx = 27.3, Γy = 6.54, Domain Box without fins, Ra 1708 ≤ Ra ≤ 1800, Pr 1, 5.373, Initiation T0, q’’0. Row 3: Simulation 3, Aspect ratio Γx = 27.3, Γy = 6.54, Domain Box with fins, figure 2, Ra 1708 ≤ Ra ≤ 1800, Pr 5.373, Initiation T0, q’0. Row 4: Simulation 4, Aspect ratio Γ = 40, Domain Cylindrical, figure 3, Ra 3 × 10^3 ≤ Ra ≤ 10^4, Pr 1, Initiation Gaussian.
3. Results and discussion
3.1. Fronts of straight parallel convection rolls
We first discuss fronts which leave in their wake a field of straight and parallel convection rolls using the 2-D domain (simulation 1). Figure 4(a) shows colour contours of
$T(x,z)$
for
$0 \leq x \leq 20$
in a domain with aspect ratio
$\varGamma = 60$
where
${Ra} = 1745$
. The front is initiated using
$T(x = 0) = 1$
and is travelling from left to right. Red and blue contours represent hot and cold fluid, respectively. The temperature field is shown at four instances of time with time increasing from top to bottom.
The temperature field can be used to identify locations with, and without, fluid motion. This is illustrated in figures 4 and 5 where we show temperature contours, temperature profiles and contours of the stream function. In the absence of fluid motion,
$T(x,z)$
varies linearly in the
$z$
direction due to heat conduction to yield
$T = 1 - z$
. This is shown on the far right-hand side of the top panel of figure 4(a) where the fluid is motionless since the front has not yet arrived at this location and, as a result,
$T(x,z)$
varies uniformly from red to blue. However, near the left wall, where convective motion is present,
$T(x,z)$
deviates from its linear variation indicating the presence of hot rising fluid and cool falling fluid due to the convection rolls.
The temperature variation of a fluid layer containing a front in a 2-D domain forming a chain of counter-rotating convection rolls. Simulation parameters: simulation 1,
$\varGamma = 60$
,
${Ra} = 1745$
,
$ \textit{Pr} = 1$
. (a) Colour contours of
$T(x,z)$
at four times, the front propagating from left to right. The region
$0 \leq x \leq 20$
is shown for clarity. Red is hot rising fluid and blue is cool falling fluid. Time increases from top to bottom. (b) Temperature profiles
$T(x,z = 1/2)$
as a function of time
$t$
. Time increases from bottom to top. The bottom curve is at
$t = 0$
, the interval between the remaining profiles being 5.68 time units.
Figure 4. Long description
Panel A: This panel shows a series of color contour plots representing temperature variations in a fluid layer over time. The front propagates from left to right. The region displayed is highlighted for clarity. Red indicates hot rising fluid, while blue indicates cool falling fluid. Time increases from top to bottom. Panel B: This panel presents temperature profiles as a function of time. The profiles are plotted with time increasing from bottom to top. The bottom curve corresponds to a specific time, with subsequent profiles spaced at intervals of 5.68 time units.
We use the temperature field at the midplane,
$T(x,z = 1/2)$
, referred to as the temperature profile, to quantify the position of the front
$x_{\kern-1pt f}$
and the wavelength
$\lambda$
of the convection rolls. The variation of the temperature profile at several times is shown in figure 4(b). The horizontal line at the bottom is the initial state of no fluid motion, with a linear conduction temperature variation that yields
$T(x,z = 1/2) = 0.5$
over the entire domain. Each curve is the temperature profile at a different instant of time where time increases from the bottom to the top. The different temperature profiles are translated vertically in order to show them clearly on one plot.
The front location
$x_{\kern-1pt f}$
is identified as the position where the temperature profile first deviates from a value of 1/2 when examining the temperature profile from right to left (for example, see figure 4
b). Numerically, we quantify
$x_{\kern-1pt f}$
as the largest value of
$x$
such that
$|T(x,z = 1/2) - 1/2| \gt \delta$
for
$\delta = 1 \times 10^{-3}$
. We have verified that our results do not vary significantly upon the value of the threshold used over the range
$1 \times 10^{-4} \leq \delta \leq 1 \times 10^{-2}$
. The front location separates quiescent fluid (
$x \gt x_{\kern-1pt f}$
) from convective fluid motion (
$x \le x_{\kern-1pt f}$
). The front velocity is the time rate of change of the front location,
$v_{\kern-1pt f} = \dot {x}_{\kern-1pt f}$
.
Peaks, or local maxima, in the temperature profiles shown in figure 4(b) indicate rising hot fluid and troughs, or local minima, indicate descending cold fluid. The
$x$
location of the centre of a roll pair can be identified as a local maximum in the temperature profile. The local maximum indicates the upflow which occurs at the shared roll boundary at the centre of the roll pair. This roll pair is then bounded on its left, and right, by adjacent local minima indicating the region of downflow for each roll of the pair. A similar argument can be made using a local minimum in the temperature profile to locate the centre of a roll pair where the adjacent local maxima locate the outer edges of the roll pair. The centres of individual convection rolls occur at locations where
$T(x,z = 1/2) = 0.5$
in regions where convective motion is occurring. We quantify
$\lambda$
using the location of the roll centres.
The connection between
$T(x,z)$
and
$\boldsymbol{u}(x,z)$
is shown in figure 5 where each row shows
$T(x,z$
) contours on the left and contours of the stream function
$\psi (x,z)$
on the right at the time when
$x_{\kern-1pt f} = 24$
. The stream function is defined in the usual manner as
$u = \partial \psi /\partial z$
and
$w = - \partial \psi /\partial x$
, where
$u(x,z)$
and
$w(x,z)$
are the
$x$
and
$z$
components of the fluid velocity vector, respectively. Red indicates counterclockwise fluid motion, blue indicates clockwise fluid motion and green represents negligible fluid motion. The contours of
$\psi (x,z)$
yield the structure of the counter-rotating convection rolls. Only the spatial region
$13 \leq x \leq 26$
is shown to highlight the region near the front.
Colour contours of temperature (a–d) and normalised stream function (e–h) for fronts in a 2-D domain. Each row is for a different value of
${Ra}$
and shows contours at the time when
$x_{\kern-1pt f} = 24$
. The spatial region
$13 \leq x \leq 26$
is shown for clarity. (a,e)
${Ra} = 2000$
; (b,f)
${Ra} = 3300$
; (c,g)
${Ra} = 5000$
; (d,h)
${Ra} = 10^4$
. Simulation parameters: simulation 1,
$\varGamma = 30$
,
$ \textit{Pr} = 1$
.
Figures 5(a) and 5(e) show
$T$
and
$\psi$
contours, respectively, for
${Ra} = 2000$
. As expected for this low value of the Rayleigh number, the temperature exhibits a decaying sinusoidal structure as the front is approached from the left, and the variation of
$\psi$
indicates the presence of counter-rotating convection rolls in the wake of the front.
As
${Ra}$
increases, the spatial structure of
$T$
and
$\psi$
deviate from this description. Most striking is the transition of
$T(x,z)$
towards a plume-like structure. In figure 5(b–d), the variation of the sinusoidal temperature contours towards mushroom-shaped structures is evident. In addition, the spacing between successive peaks in the temperature profile decreases with increasing
${Ra}$
indicating a decrease in the wavelength of the convection rolls. As
${Ra}$
is increased, the flow field begins to exhibit asymmetries in the convection rolls near the leading edge of the front as shown on the far right of figures 5(d) and 5(h).
The variation of
$x_{\kern-1pt f}$
with scaled time
$t^*$
is shown in figure 6 where
$t^* = \epsilon t/\tau _0$
with the characteristic time
$\tau _0$
given by
$\tau _0^{-1} = 19.65 {Pr}/({Pr} + 0.5117)$
for RBC with rigid boundaries (Cross Reference Cross1980). The convergence of
$v_{\kern-1pt f}(t)$
to its steady asymptotic value
$\bar {v}_{\kern-1pt f}$
is algebraically slow (Ebert & van Saarloos Reference Ebert and van Saarloos2000). It has been shown that the convergence is approximately achieved when
$t^* \gtrsim 5$
(Kockelkoren et al. Reference Kockelkoren, Storm and van Saarloos2003). The upper curve (blue) is for
${Ra} = 1745$
with
$\varGamma = 60$
and the lower curve (red) is for
${Ra} = 2000$
with
$\varGamma = 30$
. The slightly jagged nature of the variation of
$x_{\kern-1pt f}(t^*)$
is an artefact of how we numerically determine its value from the temperature profiles.
The Ginzburg–Landau, or amplitude, equation serves as a model system that generalises properties of complex non-equilibrium systems like RBC near threshold,
$\epsilon \to 0$
. The amplitude equation is a reduction of the the full Boussinesq equations through an asymptotic expansion in the small parameter
$\epsilon$
(Newell & Whitehead Reference Newell and Whitehead1969; Segel Reference Segel1969), and has been compared with experiments and simulations, with very good agreement (Cross Reference Cross1980; Fineberg & Steinberg Reference Fineberg and Steinberg1987; Lücke et al. Reference Lücke, Mihelcic and Kowalski1987). The equation describing the growth of the finite amplitude,
$A(x,t)$
, of a slowly varying and spatially periodic state near threshold is
\begin{align} \tau _0 \frac {\partial A}{\partial t} = \epsilon A + \xi _0^2 \frac {\partial ^2 A}{\partial x^2} - g |A|^2 A. \end{align}
A theoretical prediction of the variation of the asymptotic front velocity with
$\epsilon$
can be determined by solving the amplitude equation to yield
\begin{align} \bar {v}_{f,0} = 2\xi _0 \tau _0^{-1} \epsilon ^{1/2} ,\end{align}
where
$\xi _0$
is the correlation length (van Saarloos Reference van Saarloos2003). For RBC with rigid boundaries the correlation length is
$\xi _0^2 = 0.148$
(Cross Reference Cross1980). The front velocity given by (3.2) is the pulled reaction–diffusion front velocity. Depending on their initiation and the front velocity, propagating fronts can be classified as either pushed or pulled. Pulled fronts propagate with a velocity that is determined by the linearised dynamics at the leading edge of the front. In contrast, pushed fronts propagate with a front velocity that is greater than this value and is governed by the nonlinearities behind the leading edge of the front (van Saarloos Reference van Saarloos2003).
It is interesting to note the similarity between autocatalytic reaction–diffusion fronts and the amplitude equation (3.1) (van Saarloos Reference van Saarloos2003). The term
$\xi _0^2 \partial _x^2 A$
is analogous to the diffusion term in an autocatalytic reaction–diffusion equation. The term
$\epsilon A - g |A|^2 A$
is analogous to the autocatalytic production term, which, when linearised about the base state of
$A=0$
, yields
$\epsilon A$
. As in a pulled reaction–diffusion front, the propagation speed is set by the linear spreading rate of the leading edge, which yields the particular form of the front velocity given by (3.2). For linearised dynamics to prevail during the front evolution, the spatial variation of the initial condition must be steeper than
${\rm e}^{-\sqrt {\epsilon }/\xi _0}$
. We have ensured that this condition is satisfied in our simulations.
The variation of
$x_{\kern-1pt f}$
with
$t^*$
for fronts of straight parallel rolls propagating in a 2-D domain with
$ \textit{Pr} = 1$
where
$t^* = \epsilon t/\tau _0$
. Upper two curves:
${Ra} = 1745$
,
$\varGamma = 60$
, numerical results (blue) and the prediction (black) given by (3.4). Lower two curves:
${Ra} = 2000$
,
$\varGamma = 30$
, numerical results (red) and the prediction (black) given by (3.4).
The convergence of the front velocity towards its asymptotic value can be expressed as (Kockelkoren et al. Reference Kockelkoren, Storm and van Saarloos2003)
\begin{align} \frac {v_{\kern-1pt f}(t^*)}{\bar {v}_{f,0}} = 1 - \frac {3}{4} t^{*-1} + \frac {3 \sqrt {\pi }}{4} t^{*-3/2} + \mathcal{O} \big(t^{*-2} \big). \end{align}
It will be convenient to integrate (3.3) to yield an expression for the front position:
\begin{align} x_{\kern-1pt f}(t^*) = \bar {v}_{f,0} \left ( t^* - \frac {3}{4} \ln t^* - \frac {3 \sqrt {\pi }}{ 2} t^{*-1/2} \right ) + \mathcal{O} \big(t^{*-1} \big), \end{align}
where we have set the integration constant to zero assuming the front is initially at the origin.
The solid lines (black) in figure 6 are the predictions given by (3.4). For
${Ra} = 1745$
(upper two curves) the front position asymptotically converges to the theoretical value for
$t^* \gtrsim 5$
. For
${Ra} = 2000$
(lower two curves) the asymptotic state is approximately reached for
$t^* \gtrsim 1.5$
. As expected, the time for convergence increases as
$\epsilon \rightarrow 0$
.
For all of our 2-D results we have ensured that the fronts have significantly approached their asymptotic state prior to quantifying the asymptotic front velocity
$\bar {v}_{\kern-1pt f}$
. It is important to highlight that in experiment it is often very difficult to reach the asymptotic state due to the large-aspect-ratio domain that would be required. In this respect, the numerical simulations provide direct access to the long-time asymptotic dynamics. In the experiments of Fineberg & Steinberg (Reference Fineberg and Steinberg1987) the largest times accessible were
$3 \lesssim t^* \lesssim 4$
.
The variation of
$\bar {v}_{\kern-1pt f}$
with
$\epsilon$
is shown in figure 7 for several cases and for an
$\epsilon$
variation of over four orders of magnitude. The front velocity for rolls in a 2-D domain with
$ \textit{Pr} = 1$
are shown using red squares. The solid lines are the predictions of the front velocity using (3.2). The lower line is the theoretically predicted front velocity with
$ \textit{Pr} = 1$
, which yields
$\bar {v}_{f,0} = 10 \epsilon ^{1/2}$
. The front velocity is well described by the theoretical prediction over the entire range explored. The front velocity for
$ \textit{Pr} = 5.373$
is shown using green diamonds where the upper line is the predicted front velocity for
$ \textit{Pr} = 5.373$
, which yields
$\bar {v}_{f,0} = 13.8 \epsilon ^{1/2}$
. For
$\epsilon \lesssim 1$
the front velocity follows the
$\epsilon ^{1/2}$
trend.
Variation of
$\bar {v}_{\kern-1pt f}$
with
$\epsilon$
, plotted on logarithmic axes. Red squares: 2-D domain,
$ \textit{Pr} = 1$
; green diamonds: 2-D domain,
$ \textit{Pr} = 5.373$
; grey triangles: box domain with fins,
$ \textit{Pr} = 5.373$
; blue circles: cylindrical domain,
$ \textit{Pr} = 1$
. Solid lines are the predicted front velocity using (3.2):
$ \textit{Pr} = 1$
(lower),
$ \textit{Pr} = 5.373$
(upper).
We observe a deviation from the
$\epsilon ^{1/2}$
trend for
$\epsilon \gtrsim 1$
when Pr = 5.373, where the measured front velocities (green diamonds) exceed the theoretical prediction. The fronts are faster than predicted for a pulled front under these conditions. This deviation could indicate a transition from pulled- to pushed-front behaviour, where the nonlinearities behind the leading edge increase the front velocity (van Saarloos Reference van Saarloos2003). It has been found that adding a symmetry-breaking term in the Swift–Hohenberg model leads to a transition from pulled to pushed fronts (van Saarloos Reference van Saarloos1988, Reference van Saarloos1989, Reference van Saarloos2003). In the context of RBC, the symmetry breaking could be the transition to the plume-like structure of the flow field, away from the sinusoidal structure at higher Rayleigh numbers, as shown in figure 5(a–d). This transition of the convective rolls occurs at a lower value of the Rayleigh number as the Prandtl number is increased, which could explain why the Pr = 1 results follow the theoretical expectations. We, however, do not explore this transition in detail further here.
Our results show that
$\bar {v}_{\kern-1pt f}$
increases with increasing
$ \textit{Pr}$
. The
$\bar {v}_{\kern-1pt f} \propto \epsilon ^{1/2}$
trend in the front velocity is derived from the amplitude equation in the limit of small
$\epsilon$
and it is not expected to be valid for larger values
$\epsilon \gtrsim 1$
. It is interesting to note that it continues to describe the fronts well for
$\epsilon$
as large as
$\epsilon \approx 10$
for our results with
$ \textit{Pr} = 1$
.
The variation of the asymptotic wavenumber with
$\epsilon$
is shown in figure 8. We quantify the time variation of the average wavenumber of the rolls in the following manner. At each time, we determine the wavelengths of the convection rolls that are present as the distance between three adjacent roll centres using the temperature profile. We average the measured wavelengths to obtain an average value of the wavelength for the entire fluid layer at time
$t$
which we refer to as
$\lambda (t)$
. The average wavenumber is
$q(t) = 2 \pi /\lambda (t)$
and we estimate
$\bar {q}$
using the fit
$q(t) = \bar {q} + b/t$
, where
$b$
is a constant. Both the asymptotic values
$\bar {q}$
and
$\bar {v}_{\kern-1pt f}$
are obtained from fits to the data which return their values within a 95 % confidence interval. We define the uncertainty as the half-width of this interval, which for the values reported here is within 0.01 % of the fitted values.
Variation of
$\hat {q}$
with
$\epsilon$
, plotted on logarithmic axes. Red squares: 2-D domain,
$ \textit{Pr} = 1$
,
$q_0 = 3.114$
; green diamonds: 2-D domain,
$ \textit{Pr} = 5.373$
,
$q_0 = 3.114$
; grey triangles: box domain with fins,
$ \textit{Pr} = 5.373$
,
$q_0 = 3.026$
; blue circles: cylindrical domain,
$ \textit{Pr} = 1$
,
$q_0 = 3.117$
. Solid line (lower): curve fit using
$\hat {q} = 0.28 \epsilon$
for
$\epsilon \lesssim 0.1$
. Solid line (upper): curve fit using
$\hat {q} = 0.12\epsilon ^{1/4}$
for
$\epsilon \gtrsim 1$
. The dashed line represents the experimental measurements of Fineberg & Steinberg (Reference Fineberg and Steinberg1987) where
$\hat {q} \propto \epsilon ^{1/2}$
with
$q_0 = 2.74$
.
In order to compare the wavenumber trends for multiple cases on one plot, we use the normalised wavenumber
$\hat {q} = (\bar {q}-q_0)/q_0$
. We use
$q_0$
rather than
$q_c$
in this normalisation since
$q_0 \ne q_c$
for all of the cases we study. The value of
$q_0$
is determined by computing
$\bar {q}(\epsilon )$
for decreasing
$\epsilon$
and using a linear curve fit to obtain a wavenumber value at
$\epsilon = 0$
. For our 2-D simulations this resulted in
$q_0 = 3.114$
(essentially yielding
$q_0 \approx q_c$
) where
$q_0$
is independent of
$ \textit{Pr}$
.
The variation of
$\hat {q}$
for the 2-D domain, with
$ \textit{Pr} = 1$
, is shown by the red squares. Near threshold,
$\epsilon \lesssim 0.1$
,
$\hat {q}$
varies linearly as indicated by the lower solid line. Further away from threshold,
$\epsilon \gtrsim 1$
, the variation transitions to a
$\epsilon ^{1/4}$
dependence as indicated by the fit given by the solid line passing through the results for larger
$\epsilon$
. The wavenumber variation in the 2-D domain with
$ \textit{Pr}=5.373$
is shown by the green diamonds. Overall, the trends are similar with an
$\epsilon$
scaling transitioning to
$\epsilon ^{1/4}$
further from threshold. The wavenumbers of the rolls are larger for larger
$ \textit{Pr}$
.
The linear scaling of
$\hat {q}$
near threshold agrees with the wavenumber
$q_m$
which maximises the growth rate of perturbations in the linear regime (Dominguez-Lerma et al. Reference Dominguez-Lerma, Ahlers and Cannell1984). This can be expressed as
$\hat {q}_m = \alpha \epsilon$
where
$q_0 = q_c$
. The constant
$\alpha$
accounts for the
$ \textit{Pr}$
dependence and is given by
$\alpha = 0.0494 + 0.295 {Pr} ({Pr} + 0.509)^{-1}$
(Dominguez-Lerma et al. Reference Dominguez-Lerma, Ahlers and Cannell1984). For
$ \textit{Pr}=1$
this yields
$\alpha = 0.245$
; for
$ \textit{Pr} = 5.373$
this yields
$\alpha = 0.319$
. It is interesting to note that the linear scaling that fits our data,
$\hat {q} = 0.28 \epsilon$
, is the mean of these two values of
$\alpha$
. A linear scaling of the wavenumber was also reported in a 2-D numerical study for
$0.01 \le \epsilon \le 0.2$
by Lücke et al. (Reference Lücke, Mihelcic and Kowalski1987).
It is insightful to compare these results with the wavenumbers generated by pattern-forming fronts using the Swift–Hohenberg equation. In this case the wavenumber selected behind the front
$q_{\textit{SH}}$
is given by (van Saarloos Reference van Saarloos1989)
\begin{align} \frac {q_{\textit{SH}}}{q_c} = \frac {3 (3+ \sqrt {1 + 6 \epsilon })^{3/2}}{8 (2+\sqrt {1 + 6 \epsilon })}. \end{align}
In the limit
$\epsilon \ll 1$
this becomes
$\hat {q}_{\textit{SH}} = \epsilon /8$
where
$q_0 = q_c$
which again recovers the linear variation with
$\epsilon$
. Additionally, in the limit of
$\epsilon \gg 1$
, equation (3.5) yields
$\hat {q}_{\textit{SH}} \approx 0.587 \epsilon ^{1/4}$
in agreement with the large-
$\epsilon$
trends of the 2-D numerical results. This is expected since the Swift–Hohenberg equation is a model equation that shares a qualitative structure with the full Boussinesq equations (Swift & Hohenberg Reference Swift and Hohenberg1977).
The experimental results of Fineberg & Steinberg (Reference Fineberg and Steinberg1987) are shown in figure 8 by the dashed line. The experimentally measured wavenumbers are significantly different from the wavenumbers found in the 2-D domain. Furthermore, a striking difference between experiment and the 2-D results is the variation of
$q_0$
. In experiment (Fineberg & Steinberg Reference Fineberg and Steinberg1987)
$q_0 = 2.74$
and in our 2-D numerics
$q_0 = 3.114$
. The wavelength of the convection rolls at critical in experiment is over 13 % larger than what is found in the 2-D numerics. We re-emphasise that linear theory of bulk convective instability in an infinite fluid layer predicts
$q_c=3.117$
, whereas the front-selected wavenumbers at critical, in both experiments by Fineberg & Steinberg (Reference Fineberg and Steinberg1987) and our numerical results, satisfy
$q_0 \lt q_c$
. Our 2-D numerical results yield
$q_0=3.114$
which is very close to
$q_c$
in magnitude, approximately 0.1 % smaller, while
$q_0$
from Fineberg & Steinberg (Reference Fineberg and Steinberg1987) is approximately 12 % smaller than
$q_c$
. In simulation 3, where we reproduce the experimental domain with finned sidewalls and a constant flux at the left wall, the front-selected wavenumber
$q_0$
is still about 2 % less than
$q_c$
.
We illustrate this more clearly in figure 9 where the variation of the asymptotic wavelength with
$\epsilon$
is shown near critical. The experimental results of Fineberg & Steinberg (Reference Fineberg and Steinberg1987) for
$ \textit{Pr} = 5.373$
are shown by the dashed line. Our 2-D numerical results using
$ \textit{Pr} = 5.373$
are the green diamonds. The significant difference in
$\lambda _0$
is clearly evident. For reference, we also include the 2-D results with
$ \textit{Pr} = 1$
as the red squares illustrating an insensitivity of
$\lambda _0$
on
$ \textit{Pr}$
.
Variation of the wavelength
$\bar {\lambda }$
with
$\epsilon$
. Red squares: 2-D domain,
$ \textit{Pr} = 1$
; green diamonds: 2-D domain,
$ \textit{Pr} = 5.373$
; pink right triangles: box domain without fins and a left wall with constant temperature
$T_0$
;
$ \textit{Pr} = 5.373$
for this and all the subsequent cases; orange right triangles: box domain without fins and a left wall with constant heat flux
$q_0''$
; white left triangles: 3-D box domain with fins and a left wall held at constant temperature
$T_0$
; grey left triangles: box domain with fins and a left wall with constant heat flux
$q_0''$
. The solid line through the 2-D results is
$\bar {\lambda }(\epsilon ) = \lambda _0(1+\alpha \epsilon )^{-1}$
, where
$\lambda _0 = 2.018$
and
$\alpha = 0.28$
. The dashed line denotes the experimental measurements (Fineberg & Steinberg Reference Fineberg and Steinberg1987) given by
$\lambda = \lambda _0 (1-b\sqrt {\epsilon })$
with
$\lambda _0 = 2.29$
and
$b = 0.18$
. In all cases shown,
$T_0 = 1$
and
$q_0'' = 0.014$
. Inset: a close-up view to clearly show the differences in the numerical results for the different 3-D cases shown.
In comparison with the 2-D results, the asymptotic wavelengths
$\bar {\lambda }$
obtained from the simulations using box domains (simulations 2 and 3) are larger. We explore box domains with and without fins, and initiate convection with either a constant temperature
$T_0$
or a constant flux
$q_0''$
at the left wall. For ease of comparison, the inset of figure 9 shows a close-up view of the results using the box domains.
The pink right triangles represent results for the box domain without fins, where we have used a constant temperature
$T_0$
at the left wall to initiate the front. This box domain can be thought of as a lateral extension of the 2-D domain in the
$y$
direction (red squares and green diamonds). The convection rolls, in this case, select larger
$\bar {\lambda }$
in comparison with the 2-D case. When the same configuration is driven by
$q_0''$
on the left wall instead of constant
$T_0$
, the wavelength of the convection rolls increases further, as shown by the orange right triangles.
Including fins on the sidewalls of the box domain (simulation 4) further modifies the asymptotic wavelength. The white left triangles represent results from the box domain with fins using a constant
$T_0$
to initiate the rolls at the left wall. The presence of fins increases the wavelength of the convection rolls when compared with the box domain without fins (pink right triangles). Finally, a box domain with fins using a constant flux initiation
$q_0''$
yields the largest values of
$\bar {\lambda }$
in our study as shown by the grey left triangles.
The influence of adding fins to the box domain is directly reflected in the resulting fluid dynamics. Flow-field images are shown in figure 10(a) using colour contours of the fluid velocity in the
$z$
direction at the horizontal midplane,
$w(x,y,z = 1/2)$
. The front is travelling from left to right where red is rising fluid, blue is falling fluid and green indicates negligible fluid motion. It is clear from these flow fields that the propagating front of convection rolls is a 3-D phenomenon. For example, the growth of the new convection rolls initially occurs in the centre of the domain (
$y \approx \varGamma _y/2$
) away from the no-slip sidewalls located at
$y = 0$
and
$y = \varGamma _y$
.
Propagating fronts in box domains (
${Ra} = 1710$
,
$ \textit{Pr} = 5.373$
). Contours of the fluid velocity in the
$z$
direction at the horizontal midplane,
$w(x,y,z = 1/2)$
. Red is rising fluid, blue is falling fluid and green is negligible fluid motion. (a) Box domain (simulation 3). (b) Box domain with fins (simulation 4). In both panels the front is initiated using a constant heat flux at
$x = 0$
of
$q_0''=0.014$
.
The variations of
$x_{\kern-1pt f}$
and
$q$
with the scaled time
$t^*$
are shown by the solid curves in figure 11 for
${Ra} = 1750$
and
$ \textit{Pr} = 5.373$
. Only the time window
$2 \lesssim t^* \lesssim 5$
is shown. For
$t^* \lesssim 2$
, the front position
$x_{\kern-1pt f}$
deviates significantly from the asymptotic state. For
$t^* \gtrsim 5$
, the convection rolls in the box domains interact with the far right boundary of the domain.
The front position is consistent with the prediction given by (3.4) as shown by the black dashed lines in figure 11(a). The colour conventions used here are consistent with those of figure 9. The red, pink and black solid curves correspond to the 2-D domain, box domain and box domain with fins, respectively, with each initiated by using a constant temperature
$T_0$
at the left wall. The orange and grey curves correspond to the box domain and box domain with fins, respectively, which have been initiated with a constant heat flux
$q_0''$
.
It is interesting to note that the front position, when initiated with the constant-temperature sidewall, is always larger than the front position when the initiation is done with the constant-heat-flux sidewall despite each having the same asymptotic front speed. The difference results from a delay in the initiation of the fronts when initiated with a constant flux. The initial heat flux into the fluid for the constant-temperature-sidewall case is much larger than the value of
$q_0''$
that we use. As a result, more time is required for the convection roll initiation in the constant-heat-flux case. Despite this initial offset, the asymptotic front velocity selected by all the cases remains identical, as indicated by the slopes of the two black dashed lines. We anticipate that this delay in initiation is a function of
$q_0''$
and we have not explored this further here.
The variation of the wavenumber in figure 11(b) is quite interesting for several reasons. The numerical results are shown by the coloured curves (using the conventions of figure 11
a) and the experimental measurement is the dash-dotted line. These results indicate that the wavenumber in the 2-D domain, with initiation using
$T_0$
, is the largest (red). A box domain yields a smaller wavenumber (pink) when using
$T_0$
to initiate the front. The wavenumber reduces further for a box domain which uses a constant-heat-flux sidewall for front initiation (orange). The inclusion of fins to the box domain reduces the wavenumber further, the black curves using
$T_0$
and the grey curve using
$q_0''$
for initiation. Overall, we find that the smallest wavenumber, on average, occurs for the configuration that most aligns with the experiments of Fineberg & Steinberg (Reference Fineberg and Steinberg1987). However, the wavenumbers in all of the numerical simulations are significantly larger than those measured experimentally as indicated by the dash-dotted line.
The variation of
$x_{\kern-1pt f}$
and
$q$
with scaled time for fronts of straight parallel rolls in 2-D and box domains (
${Ra} = 1750$
,
$ \textit{Pr} = 5.373$
). The variation of (a)
$x_{\kern-1pt f}(t^*)$
and (b)
$q(t^*)$
. Red: 2-D domain with constant-temperature sidewall; pink: box domain with constant-temperature sidewall; orange: box domain with constant-flux sidewall; black: finned box domain with constant-temperature sidewall; grey: finned box domain using a constant-heat-flux sidewall (
$q_0'' = 0.014$
). All cases with a hot sidewall use
$T_0 = 1$
. The black dashed lines in (a) are obtained from (3.4) and the dash-dotted line in (b) is
$\bar {q}$
from experiment (Fineberg & Steinberg Reference Fineberg and Steinberg1987).
It is useful to note that for values of
${Ra}$
near onset, the accessible time window for the study of front propagation is limited to
$t^* \lesssim 3$
due to the size of the box domains. Therefore, it is expected that the fronts will not have reached their asymptotic states during this time. For these cases, we fit the available data for
$q(t)$
and
$x_{\kern-1pt f}(t)$
to determine the asymptotic values of wavenumber and front speed.
The front velocity in the box domain with fins is similar to what is found using a 2-D domain. This is illustrated by comparing the green diamonds and the grey triangles in figure 7. Similarly, the variation of the scaled wavenumber of the rolls in the box domain with fins is similar to the results from the 2-D domain as illustrated in figure 8. However, the wavelength of the rolls at critical is larger for the box domain with fins as indicated in figure 9. Although the wavelength has increased, it still remains smaller than the experimental values indicated by the dash-dotted line in figure 11(b).
3.2. Fronts of concentric rolls
We next discuss the propagation of a front which forms concentric convection rolls in a cylindrical domain (simulation 4). Figure 12 shows colour contours of the temperature field at the horizontal midplane,
$T(x,y,z = 1/2)$
, for a front at several instances of time. The front propagates radially outward towards the boundary as seen in figure 12(a–c). At later times, shown in figure 12(d–f), the bulk convective instability leads to the emergence of the spiral defect chaos state (Morris et al. Reference Morris, Bodenschatz, Cannell and Ahlers1993; Bodenschatz et al. Reference Bodenschatz, Pesch and Ahlers2000; Vitral et al. Reference Vitral, Mukherjee, Leo, Viñals, Paul and Huang2020) in the region beyond the front. For very long times (not shown) the concentric rolls are annihilated and the entire fluid layer exhibits spiral defect chaos.
The propagation of concentric convection rolls in a large cylindrical domain:
${Ra} = 3000$
(
$\epsilon = 0.76$
),
$ \textit{Pr} = 1$
,
$\varGamma = 40$
, simulation 4. The front is initiated at the centre of the domain and propagates radially outward towards the sidewalls. Colour contours of
$T(x,y,z = 1/2)$
are shown at the times: (a)
$t = 0.28$
, (b)
$t = 1.2$
, (c)
$t = 2.3$
, (d)
$t = 2.73$
, (e)
$t = 4.62$
, (f)
$t = 11$
. Red is hot rising fluid and blue is cool falling fluid. For these conditions, the fluid layer is unstable to the spiral defect chaos state which appears at later times (see d–f).
Figure 12. Long description
Panel A: A circular diagram showing a small red core at the center, surrounded by a green background. Panel B: A circular diagram with a larger red spiral pattern at the center, surrounded by a green background. Panel C: A circular diagram with a more complex red spiral pattern at the center, surrounded by a green background with small disturbances. Panel D: A circular diagram with a highly disturbed red spiral pattern at the center, surrounded by a green background with more pronounced disturbances. Panel E: A circular diagram with a chaotic red spiral pattern at the center, surrounded by a background with significant disturbances. Panel F: A circular diagram with an even more chaotic red spiral pattern at the center, surrounded by a background with extensive disturbances.
The velocity of the front leaving concentric rolls in its wake is shown in figure 7 by the blue circles. The front velocity again follows the
$\epsilon ^{1/2}$
scaling as predicted by (3.2). The scaled wavenumbers of the concentric rolls are shown in figure 8 by the blue circles where we find
$\hat {q} = 0.12 \epsilon ^{1/4}$
in agreement with the
$\epsilon ^{1/4}$
scaling predicted by the Swift–Hohenberg equation away from threshold. It is interesting to point out that the wavelength of the propagating concentric rolls is smaller than the wavelength of the disordered rolls present in the spiral defect region (see figure 12
e,f). It is useful to note that the concentric convection roll front is initiated with the perturbation given by (2.4) which has been chosen to ensure a rapid, localised drop-off of the initial disturbance. We have not evaluated the dependence of the selected wavenumbers against the details of this perturbation; this would be an interesting direction of future work.
It is insightful to compare our results for the wavenumber of concentric rolls formed behind a front with the wavenumber selected by concentric rolls in the long-time limit due to the bulk instability. The asymptotic wavenumber selected by concentric rolls which completely fill a cylindrical domain has been studied experimentally and theoretically in detail (Cross & Hohenberg Reference Cross and Hohenberg1993; Bodenschatz et al. Reference Bodenschatz, Pesch and Ahlers2000). It has been shown that the asymptotic wavenumber selected by concentric rolls due to the bulk instability decreases with increasing
$\epsilon$
for
$ \textit{Pr} \gtrsim 1$
(Buell & Catton Reference Buell and Catton1986). This must be contrasted with the findings for the concentric rolls that form behind a front which yield an increase in the selected wavenumber with increasing
$\epsilon$
as indicated in figure 8 and the
$\hat {q} \sim \epsilon ^{1/4}$
scaling that is shown. The values of the wavenumbers are also significantly different, the concentric rolls from the bulk instability yielding much smaller wavenumbers than what is found for the front-selected wavenumbers. For example, the wavenumber selected by the concentric convection rolls shown in figure 12 is
$\bar {q} = 3.47$
. However, for these conditions (
${Ra} = 3000$
,
$ \textit{Pr} = 1$
) the bulk instability selects a wavenumber of
$\bar {q} \approx 2.82$
(Buell & Catton Reference Buell and Catton1986). The average selected wavenumber for the spiral defect chaos state is even smaller with a value of
$\bar {q} \approx 2.5$
(Chiam et al. Reference Chiam, Paul, Cross and Greenside2003).
In a recent study, we explored the formation of convection rolls behind a reaction front that added heat to the fluid while also changing its density (Mukherjee & Paul Reference Mukherjee and Paul2022). The propagating reaction front temporarily annihilated the spiral defect chaos state as it travelled and in its wake a front forming concentric convection rolls emerged (see figures 15 and 16 in Mukherjee & Paul (Reference Mukherjee and Paul2022)). The wavenumber of the forming concentric rolls was larger than the average wavenumber of convection rolls undergoing spiral defect chaos in the bulk. The average wavenumber selected by the forming convection rolls in the wake of the reaction front for
${Ra} = 6 \times 10^3$
was
$ \bar {q} \approx 3.9$
. This value of the wavenumber is comparable with the wavenumbers we find here for concentric rolls that remain behind a propagating front in an initially quiescent fluid domain.
4. Conclusion
We have explored propagating fronts that form convection rolls in RBC across a wide range of conditions. When the Rayleigh number is just above the convective instability, a heated sidewall is used to initiate the propagation of a chain of convection rolls. The pattern-forming front selects a unique wavenumber and front velocity that depend on the Rayleigh number, on the Prandtl number and on the details of the convection domain.
We have quantified the front velocity by tracking the leading edge of the temperature profile. We find that
$v_{\kern-1pt f} \sim \epsilon ^{1/2}$
for most of the cases we have explored. The scaling agrees with the theoretical expression for the velocity of pulled fronts. Deviations from the square-root scaling are observed only for
$ \textit{Pr} = 5.373$
when
$\epsilon \gtrsim 1$
which we anticipate is due to the increasing role of nonlinearities as
$\epsilon$
increases.
We have explored the slow asymptotic convergence of the front velocity with time. Near onset, the fronts must evolve for a sufficient time such that
$t^* \gtrsim 5$
. The aspect ratio of the domain must be tailored with this slow convergence in mind. If a domain is not large enough, the front will make contact with the far sidewall prior to reaching its asymptotic state.
Near the convective threshold, the asymptotic wavenumber scales linearly with the reduced Rayleigh number. This is in agreement with the theoretical predictions based on the wavenumber of the maximum growth in the linear regime. For larger
$\epsilon$
the wavenumber transitions to a
$\epsilon ^{1/4}$
scaling which is predicted by an analysis using the Swift–Hohenberg equation.
We have compared our numerical results with the experimental measurements of Fineberg & Steinberg (Reference Fineberg and Steinberg1987) which reported a
$\epsilon ^{1/2}$
dependence of the wavenumber near onset. An important comparison is with our simulation 3 which uses a convection domain constructed to match the geometric details and boundary conditions of the experiment. We find that the wavenumbers selected in the simulations exhibit a linear scaling with
$\epsilon$
near threshold.
We find that the wavenumber selection is influenced by the geometry of the convection domain, boundary conditions and method of front initiation. In particular, box domains produce rolls with larger wavelengths than what is found using a 2-D domain. Similarly, initiating convection using a constant-heat-flux boundary condition results in larger wavelength rolls than when using a constant-temperature boundary to initiate the front. Additionally, incorporating fins on the sidewalls of the box domains results in convection rolls with larger wavelength.
Simulations designed to incorporate details of the experimental set-up of Fineberg & Steinberg (Reference Fineberg and Steinberg1987) produced rolls with the largest wavelengths among all of the cases we explored, yet the wavelengths remained
$ \sim 10\,\%$
smaller than the experimentally measured values. The front velocity
$v_{\kern-1pt f}$
and the scaling behaviour of the normalised wavenumber
$\hat {q}$
with
$\epsilon$
do not depend significantly upon the domain geometry, boundary conditions or method of front initiation that are used. We, however, note that this universality concerns the scaling exponents of the normalised wavenumbers. The prefactors in the relations and the absolute value of the selected wavenumbers vary appreciably between cases, as shown in figure 9.
The physical origin of the large-wavelength convection rolls found in the experiments of Fineberg & Steinberg (Reference Fineberg and Steinberg1987) for
$\epsilon \ll 1$
remains an open question. Our results suggest this is not due to the use of a specific geometry, the use of finned sidewalls or the constant-heat-flux sidewall that was used to initiate the fronts. A factor that could contribute to the discrepancy is the finite thermal conductivity of the sidewalls and fins that are used in the experiment, which we include as perfect thermal conductors. Quantifying the role of bounding surfaces with finite thermal conductivity on the wavenumber selection would be an interesting direction in the future.
Acknowledgements
We acknowledge many fruitful interactions with P. Fischer and the Nek5000 user group. S.M. acknowledges B. Kumar for help with mesh generation. Portions of the numerical computations were conducted with generous support from the Advanced Research Computing Center at Virginia Tech and the High Performance Computing Center at Iowa State University.
Declaration of interests
The authors report no conflict of interest.
L
d
Th
Tc
Th>Tc
x=0
Γ=L/d
Γ=30
z
d
Lx
Ly
x
y
Γx=Lx/d
Γy=Ly/d
x=0
x
y
z
y
z
z
d
r0
Γ=r0/d
Th
Tc
t=0
z
Γ=40
T0
q0″
Γ=60
Ra=1745
Pr=1
T(x,z)
0≤x≤20
T(x,z=1/2)
t
t=0
Ra
xf=24
13≤x≤26
Ra=2000
Ra=3300
Ra=5000
Ra=104
Γ=30
Pr=1
xf
t∗
Pr=1
t∗=ϵt/τ0
Ra=1745
Γ=60
Ra=2000
Γ=30
v¯f
ϵ
Pr=1
Pr=5.373
Pr=5.373
Pr=1
Pr=1
Pr=5.373
q^
ϵ
Pr=1
q0=3.114
Pr=5.373
q0=3.114
Pr=5.373
q0=3.026
Pr=1
q0=3.117
q^=0.28ϵ
ϵ≲0.1
q^=0.12ϵ1/4
ϵ≳1
q^∝ϵ1/2
q0=2.74
λ¯
ϵ
Pr=1
Pr=5.373
T0
Pr=5.373
q0″
T0
q0″
λ¯(ϵ)=λ0(1+αϵ)−1
λ0=2.018
α=0.28
λ=λ0(1−bϵ)
λ0=2.29
b=0.18
T0=1
q0″=0.014
Ra=1710
Pr=5.373
z
w(x,y,z=1/2)
x=0
q0″=0.014
xf
q
Ra=1750
Pr=5.373
xf(t∗)
q(t∗)
q0″=0.014
T0=1
q¯
Ra=3000
ϵ=0.76
Pr=1
Γ=40
T(x,y,z=1/2)
t=0.28
t=1.2
t=2.3
t=2.73
t=4.62
t=11
Open access
